Commit 5d96e144 authored by David Frank's avatar David Frank Committed by Tobias Lasser

#64 Move geometry to core, use Strong typing in Geometry interface

- Move geometry to core module
- Use strong types for Geometry (Old constructors are kept for now, but are just wrappers, most of them use constexpr)
- Added some convience overloads for RotationAngles3D and GeometryData
- Fixed tests for removing old constructors (CircluarTrajectoryGenerator, Projectors)
parent 4d703b69
Pipeline #277596 passed with stages
in 59 minutes and 23 seconds
......@@ -15,6 +15,8 @@
#include <cstdlib>
using namespace elsa;
using namespace elsa::geometry;
static const index_t dimension = 2;
void iterate2D(const Geometry& geo)
......@@ -48,18 +50,20 @@ TEST_CASE("Ray generation for 2D")
{
IndexVector_t volCoeff(2);
volCoeff << 5, 5;
VolumeDescriptor ddVol(volCoeff);
IndexVector_t detCoeff(1);
detCoeff << 5;
VolumeDescriptor ddDet(detCoeff);
IndexVector_t detCoeff(2);
detCoeff << 5, 1;
real_t s2c = 10;
real_t c2d = 4;
auto volData = VolumeData2D{Size2D{volCoeff}};
auto sinoData = SinogramData2D{Size2D{detCoeff}};
GIVEN("Geometry without offset and rotation")
{
Geometry g(s2c, c2d, 0, ddVol, ddDet);
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{0},
std::move(volData), std::move(sinoData));
// test outer + central pixels
iterate2D(g);
......@@ -67,8 +71,8 @@ TEST_CASE("Ray generation for 2D")
GIVEN("Geometry with offset but no rotation")
{
real_t offset = 2;
Geometry g(s2c, c2d, 0, ddVol, ddDet, offset);
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{0},
std::move(volData), std::move(sinoData), PrincipalPointOffset{2});
// test outer + central pixels
iterate2D(g);
......@@ -76,8 +80,8 @@ TEST_CASE("Ray generation for 2D")
GIVEN("Geometry at 90°, but no offset")
{
real_t angle = pi_t / 2; // 90 degrees
Geometry g(s2c, c2d, angle, ddVol, ddDet);
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{90},
std::move(volData), std::move(sinoData));
// test outer + central pixels
iterate2D(g);
......@@ -85,10 +89,9 @@ TEST_CASE("Ray generation for 2D")
GIVEN("Geometry at 45° with offset")
{
real_t angle = pi_t / 4; // 45 degrees
real_t cx = -1;
real_t cy = 2;
Geometry g(s2c, c2d, angle, ddVol, ddDet, 0, cx, cy);
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{45},
std::move(volData), std::move(sinoData), PrincipalPointOffset{0},
RotationOffset2D{-1, 2});
// test outer + central pixels
iterate2D(g);
......@@ -101,16 +104,20 @@ TEST_CASE("Ray generation for 3D")
volCoeff << 5, 5, 5;
VolumeDescriptor ddVol(volCoeff);
IndexVector_t detCoeff(2);
detCoeff << 5, 5;
IndexVector_t detCoeff(3);
detCoeff << 5, 5, 1;
VolumeDescriptor ddDet(detCoeff);
real_t s2c = 10;
real_t c2d = 4;
auto volData = VolumeData3D{Size3D{volCoeff}};
auto sinoData = SinogramData3D{Size3D{detCoeff}};
GIVEN("Geometry without offset and rotation")
{
Geometry g(s2c, c2d, ddVol, ddDet, 0);
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
std::move(volData), std::move(sinoData), RotationAngles3D{Gamma(0)});
// test outer + central pixels
iterate3D(g);
......@@ -118,9 +125,9 @@ TEST_CASE("Ray generation for 3D")
GIVEN("Geometry with offset but no rotation")
{
real_t px = -1;
real_t py = 3;
Geometry g(s2c, c2d, ddVol, ddDet, 0, 0, 0, px, py);
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
std::move(volData), std::move(sinoData), RotationAngles3D{Gamma{0}},
PrincipalPointOffset2D{0, 0}, RotationOffset3D{-1, 3, 0});
// test outer + central pixels
iterate3D(g);
......@@ -128,8 +135,8 @@ TEST_CASE("Ray generation for 3D")
GIVEN("Geometry at 90°, but no offset")
{
real_t angle = pi_t / 2;
Geometry g(s2c, c2d, ddVol, ddDet, angle);
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
std::move(volData), std::move(sinoData), RotationAngles3D{Gamma{pi_t / 2}});
// test outer + central pixels
iterate3D(g);
......@@ -139,11 +146,11 @@ TEST_CASE("Ray generation for 3D")
{
real_t angle1 = pi_t / 4;
real_t angle2 = pi_t / 2;
RealVector_t offset(3);
offset << 1, -2, -1;
Geometry g(s2c, c2d, ddVol, ddDet, angle1, angle2, 0, 0, 0, offset[0], offset[1],
offset[2]);
Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
std::move(volData), std::move(sinoData),
RotationAngles3D{Gamma{angle1}, Beta{angle2}}, PrincipalPointOffset2D{0, 0},
RotationOffset3D{1, -2, -1});
// test outer + central pixels
iterate3D(g);
}
......@@ -173,4 +180,4 @@ TEST_CASE("Ray generation for 3D")
// test outer + central pixels
iterate3D(g);
}
}
\ No newline at end of file
}
......@@ -24,7 +24,9 @@ set(MODULE_HEADERS
Handlers/DataHandlerMapCPU.h
LinearOperator.h
Expression.h
ExpressionPredicates.h)
ExpressionPredicates.h
Geometry.h
StrongTypes.h)
# list all the code files of the module
set(MODULE_SOURCES
......@@ -37,7 +39,9 @@ set(MODULE_SOURCES
DataContainer.cpp
Handlers/DataHandlerCPU.cpp
Handlers/DataHandlerMapCPU.cpp
LinearOperator.cpp)
LinearOperator.cpp
Geometry.cpp
StrongTypes.cpp)
if(${ELSA_CUDA_VECTOR})
set(MODULE_HEADERS
......
......@@ -7,10 +7,12 @@
namespace elsa
{
Geometry::Geometry(real_t sourceToCenterOfRotation, real_t centerOfRotationToDetector,
real_t angle, const DataDescriptor& volumeDescriptor,
const DataDescriptor& sinoDescriptor, real_t offset,
real_t centerOfRotationOffsetX, real_t centerOfRotationOffsetY)
using namespace geometry;
Geometry::Geometry(SourceToCenterOfRotation sourceToCenterOfRotation,
CenterOfRotationToDetector centerOfRotationToDetector, Radian angle,
VolumeData2D&& volData, SinogramData2D&& sinoData,
PrincipalPointOffset offset, RotationOffset2D centerOfRotOffset)
: _objectDimension{2},
_P{RealMatrix_t::Identity(2, 2 + 1)},
_Pinv{RealMatrix_t::Identity(2 + 1, 2)},
......@@ -20,35 +22,36 @@ namespace elsa
_S{RealMatrix_t::Identity(2 + 1, 2 + 1)},
_C{RealVector_t::Zero(2)}
{
auto [volSpacing, volOrigin] = std::move(volData);
auto [sinoSpacing, sinoOrigin] = std::move(sinoData);
// setup rotation matrix _R
real_t c = std::cos(angle);
real_t s = std::sin(angle);
_R << c, -s, s, c;
// setup scaling matrix _S
_S << volumeDescriptor.getSpacingPerDimension()[0], 0, 0, 0,
volumeDescriptor.getSpacingPerDimension()[1], 0, 0, 0, 1;
_S << volSpacing[0], 0, 0, 0, volSpacing[1], 0, 0, 0, 1;
// set the translation _t
RealVector_t centerOfRotationOffset(_objectDimension);
centerOfRotationOffset << centerOfRotationOffsetX, centerOfRotationOffsetY;
// auto centerOfRotationOffset = static_cast<RealVector_t>(centerOfRotOffset);
// centerOfRotationOffset << centerOfRotationOffsetX, centerOfRotationOffsetY;
_t = _R * (-centerOfRotationOffset - volumeDescriptor.getLocationOfOrigin());
_t = _R * (-centerOfRotOffset.get() - volOrigin);
_t[_objectDimension - 1] += sourceToCenterOfRotation;
// set the intrinsic parameters _K
real_t alpha = sinoDescriptor.getSpacingPerDimension()[0];
real_t alpha = sinoSpacing[0];
_K << (sourceToCenterOfRotation + centerOfRotationToDetector) / alpha,
(sinoDescriptor.getLocationOfOrigin()[0] / alpha + offset), 0, 1;
(sinoOrigin[0] / alpha + offset), 0, 1;
buildMatrices();
}
Geometry::Geometry(real_t sourceToCenterOfRotation, real_t centerOfRotationToDetector,
const DataDescriptor& volumeDescriptor, const DataDescriptor& sinoDescriptor,
real_t gamma, real_t beta, real_t alpha, real_t px, real_t py,
real_t centerOfRotationOffsetX, real_t centerOfRotationOffsetY,
real_t centerOfRotationOffsetZ)
Geometry::Geometry(SourceToCenterOfRotation sourceToCenterOfRotation,
CenterOfRotationToDetector centerOfRotationToDetector,
VolumeData3D&& volData, SinogramData3D&& sinoData, RotationAngles3D angles,
PrincipalPointOffset2D offset, RotationOffset3D centerOfRotOffset)
: _objectDimension{3},
_P{RealMatrix_t::Identity(3, 3 + 1)},
_Pinv{RealMatrix_t::Identity(3 + 1, 3)},
......@@ -58,6 +61,13 @@ namespace elsa
_S{RealMatrix_t::Identity(3 + 1, 3 + 1)},
_C{RealVector_t::Zero(3)}
{
auto [volSpacing, volOrigin] = std::move(volData);
auto [sinoSpacing, sinoOrigin] = std::move(sinoData);
real_t alpha = angles.alpha();
real_t beta = angles.beta();
real_t gamma = angles.gamma();
// setup rotation matrix
real_t ca = std::cos(alpha);
real_t sa = std::sin(alpha);
......@@ -72,26 +82,24 @@ namespace elsa
* (RealMatrix_t(3, 3) << ca, 0, sa, 0, 1, 0, -sa, 0, ca).finished();
// setup scaling matrix _S
_S << volumeDescriptor.getSpacingPerDimension()[0], 0, 0, 0, 0,
volumeDescriptor.getSpacingPerDimension()[1], 0, 0, 0, 0,
volumeDescriptor.getSpacingPerDimension()[2], 0, 0, 0, 0, 1;
_S << volSpacing[0], 0, 0, 0, 0, volSpacing[1], 0, 0, 0, 0, volSpacing[2], 0, 0, 0, 0, 1;
// setup the translation _t
RealVector_t centerOfRotationOffset(_objectDimension);
centerOfRotationOffset << centerOfRotationOffsetX, centerOfRotationOffsetY,
centerOfRotationOffsetZ;
// RealVector_t centerOfRotationOffset(_objectDimension);
// centerOfRotationOffset << centerOfRotationOffsetX, centerOfRotationOffsetY,
// centerOfRotationOffsetZ;
_t = _R * (-centerOfRotationOffset - volumeDescriptor.getLocationOfOrigin());
_t = _R * (-centerOfRotOffset.get() - volOrigin);
_t(_objectDimension - 1) += sourceToCenterOfRotation;
// setup the intrinsic parameters _K
real_t alpha1 = sinoDescriptor.getSpacingPerDimension()[0];
real_t alpha2 = sinoDescriptor.getSpacingPerDimension()[1];
real_t alpha1 = sinoSpacing[0];
real_t alpha2 = sinoSpacing[1];
_K << (sourceToCenterOfRotation + centerOfRotationToDetector) / alpha1, 0,
sinoDescriptor.getLocationOfOrigin()[0] / alpha1 + px, 0,
sinoOrigin[0] / alpha1 + offset[0], 0,
(sourceToCenterOfRotation + centerOfRotationToDetector) / alpha2,
sinoDescriptor.getLocationOfOrigin()[1] / alpha2 + py, 0, 0, 1;
sinoOrigin[1] / alpha2 + offset[1], 0, 0, 1;
buildMatrices();
}
......
......@@ -2,19 +2,21 @@
#include "elsaDefines.h"
#include "DataDescriptor.h"
#include "StrongTypes.h"
#include <utility>
#include <cassert>
namespace elsa
{
/**
* \brief Class representing 2d/3d projective camera geometry for use in CT projectors.
* @brief Class representing 2d/3d projective camera geometry for use in CT projectors.
*
* \author Matthias Wieczorek - initial code
* \author Maximilian Hornung - modularization, redesign
* \author David Frank - bugfixes
* \author Nikola Dinev - refactoring
* \author Tobias Lasser - refactoring, modernization
* @author Matthias Wieczorek - initial code
* @author Maximilian Hornung - modularization, redesign
* @author David Frank - bugfixes, strong typing
* @author Nikola Dinev - refactoring
* @author Tobias Lasser - refactoring, modernization
*
* The location of X-ray source, volume (typically containing the center of rotation), and X-ray
* detector are encoded using projection matrices (see A. Zissermann, "Multiple View Geometry in
......@@ -24,67 +26,80 @@ namespace elsa
{
public:
/**
* \brief Constructor for 2D projective geometry
* @brief Constructor for 2D projective geometry
*
* \param[in] sourceToCenterOfRotation distance from source to the center of rotation (along
* y axis) \param[in] centerOfRotationToDetector distance from center of rotation to
* detector (along y axis) \param[in] angle rotation angle (in radians) \param[in]
* volumeDescriptor descriptor for the 2d volume \param[in] sinoDescriptor descriptor for
* the sinogram \param[in] offset offset of the principal point [default 0] \param[in]
* centerOfRotationOffsetX offset of the center of rotation in x direction [default 0]
* \param[in] centerOfRotationOffsetY offset of the center of rotation in y direction
* [default 0]
* @param[in] sourceToCenterOfRotation distance from source to the center of rotation (along
* y axis)
* @param[in] centerOfRotationToDetector distance from center of rotation to detector (along
* y axis)
* @param[in] angle rotation angle (in radians)
* @param[in] volData descriptor for the 2d volume
* @param[in] sinoData descriptor for the sinogram
* @param[in] offset offset of the principal point [default 0]
* @param[in] centerOfRotOffset offset of the center of rotation
*
*
* VolumeData2D and SinogramData2D are taken as r-value references, as it's cheaper to move,
* them in, and they are only intended as temporary objects. o construct a Geometry with
* VolumeData2D{...}/SinogramData2D{...} as temporary or move the object in, but be aware of
* reusing it
*/
Geometry(real_t sourceToCenterOfRotation, real_t centerOfRotationToDetector, real_t angle,
const DataDescriptor& volumeDescriptor, const DataDescriptor& sinoDescriptor,
real_t offset = static_cast<real_t>(0.0),
real_t centerOfRotationOffsetX = static_cast<real_t>(0.0),
real_t centerOfRotationOffsetY = static_cast<real_t>(0.0));
Geometry(geometry::SourceToCenterOfRotation sourceToCenterOfRotation,
geometry::CenterOfRotationToDetector centerOfRotationToDetector,
geometry::Radian angle, geometry::VolumeData2D&& volData,
geometry::SinogramData2D&& sinoData,
geometry::PrincipalPointOffset offset = geometry::PrincipalPointOffset{0},
geometry::RotationOffset2D centerOfRotOffset = geometry::RotationOffset2D{0, 0});
/**
* \brief Constructor for 3D projective geometry using Euler angles
* @brief Constructor for 3D projective geometry using Euler angles
*
* \param[in] sourceToCenterOfRotation distance from source to the center of rotation (along
* z axis) \param[in] centerOfRotationToDetector distance from center of rotation to
* detector (along z axis) \param[in] volumeDescriptor descriptor for the 3d volume
* \param[in] sinoDescriptor descriptor for the sinogram
* \param[in] gamma third rotation angle around y''-axis (in radians)
* \param[in] beta second rotation angle around z' axis (such that y->y'', in radians)
* [default 0] \param[in] alpha first rotation angle around y axis (such that z->z', in
* radians) [default 0] \param[in] px offset of the principal point in x direction [default
* 0] \param[in] py offset of the principal point in y direction [default 0] \param[in]
* centerOfRotationOffsetX offset of the center of rotation in x direction [default 0]
* \param[in] centerOfRotationOffsetY offset of the center of rotation in y direction
* [default 0] \param[in] centerOfRotationOffsetZ offset of the center of rotation in z
* direction [default 0]
* @param[in] sourceToCenterOfRotation distance from source to the center of rotation (along
* z axis)
* @param[in] centerOfRotationToDetector distance from center of rotation to detector (along
* z axis)
* @param[in] volData descriptor for the 3d volume
* @param[in] sinoData descriptor for the sinogram
* @param[in] angles (gamma -> around y''-axis, beta -> around z' axis, alpha -> around y
* axis) in radians
* @param[in] offset offset of the principal point
* param[in] centerOfRotOffset offset of the center of rotation
*
* Alpha, beta, gamma are Euler rotation angles using the YZY convention. They are specified
* in radians. In standard circular trajectory CT settings, we would have alpha = beta = 0,
* while gamma is the angle of rotation)
*
* VolumeData3D and SinogramData3D are taken as r-value references, as it's cheaper to move,
* them in, and they are only intended as temporary objects. So construct a Geometry with
* VolumeData3D{...}/SinogramData3D{...} as temporary or move the object in, but be aware of
* reusing it
*/
Geometry(real_t sourceToCenterOfRotation, real_t centerOfRotationToDetector,
const DataDescriptor& volumeDescriptor, const DataDescriptor& sinoDescriptor,
real_t gamma, real_t beta = static_cast<real_t>(0.0),
real_t alpha = static_cast<real_t>(0.0), real_t px = static_cast<real_t>(0.0),
real_t py = static_cast<real_t>(0.0),
real_t centerOfRotationOffsetX = static_cast<real_t>(0.0),
real_t centerOfRotationOffsetY = static_cast<real_t>(0.0),
real_t centerOfRotationOffsetZ = static_cast<real_t>(0.0));
Geometry(geometry::SourceToCenterOfRotation sourceToCenterOfRotation,
geometry::CenterOfRotationToDetector centerOfRotationToDetector,
geometry::VolumeData3D&& volData, geometry::SinogramData3D&& sinoData,
geometry::RotationAngles3D angles,
geometry::PrincipalPointOffset2D offset = geometry::PrincipalPointOffset2D{0, 0},
geometry::RotationOffset3D centerOfRotOffset = geometry::RotationOffset3D{0, 0,
0});
/**
* \brief Constructor for 3D projective geometry using a 3x3 rotation matrix
* @brief Constructor for 3D projective geometry using a 3x3 rotation matrix
*
* \param[in] sourceToCenterOfRotation distance from source to the center of rotation (along
* z axis) \param[in] centerOfRotationToDetector distance from center of rotation to
* detector (along z axis) \param[in] volumeDescriptor descriptor for the 3d volume
* \param[in] sinoDescriptor descriptor for the sinogram
* \param[in] R a 3x3 rotation matrix
* \param[in] px offset of the principal point in x-direction [default 0]
* \param[in] py offset of the principal point in y-direction [default 0]
* \param[in] centerOfRotationOffsetX offset of the center of rotation in x direction
* [default 0] \param[in] centerOfRotationOffsetY offset of the center of rotation in y
* direction [default 0] \param[in] centerOfRotationOffsetZ offset of the center of rotation
* in z direction [default 0]
* @param[in] sourceToCenterOfRotation distance from source to the center of rotation (along
* z axis)
* @param[in] centerOfRotationToDetector distance from center of rotation to
* detector (along z axis)
* @param[in] volumeDescriptor descriptor for the 3d volume
* @param[in] sinoDescriptor descriptor for the sinogram
* @param[in] R a 3x3 rotation matrix
* @param[in] px offset of the principal point in x-direction [default 0]
* @param[in] py offset of the principal point in y-direction [default 0]
* @param[in] centerOfRotationOffsetX offset of the center of rotation in x direction
* [default 0]
* @param[in] centerOfRotationOffsetY offset of the center of rotation in y direction
* [default 0]
* @param[in] centerOfRotationOffsetZ offset of the center of rotation in z direction
* [default 0]
*/
Geometry(real_t sourceToCenterOfRotation, real_t centerOfRotationToDetector,
const DataDescriptor& volumeDescriptor, const DataDescriptor& sinoDescriptor,
......@@ -95,42 +110,42 @@ namespace elsa
real_t centerOfRotationOffsetZ = static_cast<real_t>(0.0));
/**
* \brief Compute a ray (ray origin and ray direction) that hits the specified point p on
* @brief Compute a ray (ray origin and ray direction) that hits the specified point p on
* the detector
*
* \param[in] p point p on detector
* @param[in] p point p on detector
*
* \returns pair of <ro, rd>, where ro = ray origin and rd = ray direction
* @returns pair of <ro, rd>, where ro = ray origin and rd = ray direction
*
* Computation are done using the projection matrix.
*/
std::pair<RealVector_t, RealVector_t> computeRayTo(const RealVector_t& p) const;
/**
* \brief Return the projection matrix
* @brief Return the projection matrix
*
* \returns projection matrix
* @returns projection matrix
*/
const RealMatrix_t& getProjectionMatrix() const;
/**
* \brief Return the inverse of the projection matrix
* @brief Return the inverse of the projection matrix
*
* \returns the inverse of the projection matrix
* @returns the inverse of the projection matrix
*/
const RealMatrix_t& getInverseProjectionMatrix() const;
/**
* \brief Return the camera center corresponding to the projection matrix
* @brief Return the camera center corresponding to the projection matrix
*
* \returns the camera center (as a coordinate vector)
* @returns the camera center (as a coordinate vector)
*/
const RealVector_t& getCameraCenter() const;
/**
* \brief Return the rotation matrix corresponding to the projection matrix
* @brief Return the rotation matrix corresponding to the projection matrix
*
* \returns the rotation matrix
* @returns the rotation matrix
*/
const RealMatrix_t& getRotationMatrix() const;
......
#include "StrongTypes.h"
namespace elsa::geometry::detail
{
// Explicit template instantiation
template class StaticVectorTemplate<1, RealVector_t>;
template class StaticVectorTemplate<2, RealVector_t>;
template class StaticVectorTemplate<3, RealVector_t>;
template class StaticVectorTemplate<4, RealVector_t>;
template class StaticVectorTemplate<1, IndexVector_t>;
template class StaticVectorTemplate<2, IndexVector_t>;
template class StaticVectorTemplate<3, IndexVector_t>;
template class StaticVectorTemplate<4, IndexVector_t>;
template class GeometryData<1>;
template class GeometryData<2>;
template class GeometryData<3>;
template class GeometryData<4>;
} // namespace elsa::geometry::detail
\ No newline at end of file
#pragma once
#include "elsaDefines.h"
#include "DataDescriptor.h"
#include <utility>
#include <cassert>
#include <array>
namespace elsa
{
namespace geometry
{
class Degree;
/**
* @brief Class describing angles in Radians. Mostly used for parameter passing.
* Has an implicit constructor from Degree, and converts to real_t implicitly as well
*/
class Radian
{
public:
/// Default constructor (sets 0 radians)
constexpr Radian() : _radian(0) {}
/// Constructor from real_t
constexpr explicit Radian(real_t radian) : _radian(radian) {}
/// Constructor (implicit) from Degree
constexpr Radian(Degree d);
/// Conversion to degree (as real_t)
constexpr real_t to_degree() const { return _radian * 180 / pi_t; }
/// Conversion (implicit) to real_t
constexpr operator real_t() const { return _radian; }
private:
real_t _radian;
};
/**
* @brief Class describing angles in Degree. Mostly used for parameter passing.
* Has an implicit constructor from Radians, and converts to real_t implicitly as well
*/
class Degree
{
public:
/// Default constructor (sets 0 degrees)
constexpr Degree() : _degree(0) {}
/// Constructor from real_t
constexpr explicit Degree(real_t degree) : _degree(degree) {}
/// Constructor (implicit) from Radian
constexpr Degree(Radian r);
/// Conversion to radian (as real_t)
constexpr real_t to_radian() const { return _degree * pi_t / 180; }
/// Conversion (implicit) to real_t
constexpr operator real_t() const { return _degree; }
private:
real_t _degree;
};
// Defined out of class, to access the to_... functions
constexpr Radian::Radian(Degree d) : _radian(d.to_radian()) {}
constexpr Degree::Degree(Radian r) : _degree(r.to_degree()) {}
/// Strong type class to a specific angle alpha
class Alpha : Radian
{
public:
using Base = Radian;
using Base::Base;
using Base::to_degree;
using Base::operator real_t;
};
/// Strong type class to a specific angle beta
class Beta : Radian
{
public:
using Base = Radian;
using Base::Base;
using Base::to_degree;
using Base::operator real_t;
};
/// Strong type class to a specific angle gamma
class Gamma : Radian
{
public:
using Base = Radian;
using Base::Base;
using Base::to_degree;
using Base::operator real_t;
};
namespace detail
{
/**
* @brief Class to store a fixed length array of Radians. Ensures, that all entries are
* convertible to radians
*/
template <index_t Size>
class RotationAngles
{
public:
// Alias for Radian
using Type = Radian;
private:
/// Alias for, enable if, all (the conjunction of) Ts... are convertible to Radian
template <typename... Ts>
using AllRadian = typename std::enable_if_t<
std::conjunction<std::is_convertible<Ts, Type>...>::value>;
public:
RotationAngles(const RotationAngles&) = default;
RotationAngles& operator=(const RotationAngles&) = default;
RotationAngles(RotationAngles&&) = default;
RotationAngles& operator=(RotationAngles&&) = default;
/// Construct array (expect it to be all convertible to Radians and number of
/// arguments, equal to size)
template <typename... Ts, typename = AllRadian<Ts...>,
typename = std::enable_if_t<(Size > 0) && (Size == (sizeof...(Ts)))>>
constexpr RotationAngles(Ts&&... ts)
: RotationAngles(std::index_sequence_for<Ts...>{}, std::forward<Ts>(ts)...)
{
}
/// Access operator
constexpr Type operator[](index_t i) const { return _angles[i]; }
/// get function (non-const reference) to enable structured bindings
template <size_t I>
Type& get() &
{
return (std::get<I>(_angles));
}
/// get function (const reference) to enable structured bindings
template <size_t I>
const Type& get() const&
{
return (std::get<I>(_angles));
}
/// get function (r-value reference) to enable structured bindings
template <size_t I>
Type&& get() &&
{
return std::move(std::get<I>(_angles));
}
private:
/// Helper constructor to init array
template <std::size_t... Is, typename... Ts, typename = AllRadian<Ts...>>
constexpr RotationAngles(std::index_sequence<Is...>, Ts&&... vals)
{
// Fold expression expands to _angles[0] = vals0, _angles[1] = vals1, ...,
// _angles[k] = valsk;
(static_cast<void>(_angles[Is] = vals), ...);
}
std::array<Type, Size> _angles;
};
} // namespace detail
/**
* @brief Strong type class for 3D geometry, which expects 3 angles (alpha, beta, gamma).
* Class is constructible from the previously defined Alpha, Beta, Gamma and only subsets of
* it (i.e. only with Beta and Alpha)
*
*/
class RotationAngles3D : detail::RotationAngles<3>
{
public:
using Base = detail::RotationAngles<3>;
private:
using Base::Base;
public:
using Base::get;
using Base::operator[];