elsa projectors¶
Table of Contents
Geometry¶
-
class
elsa
::
Geometry
¶ Class representing 2d/3d projective camera geometry for use in CT projectors.
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
Computer Vision”). Detectors are assumed to be flat.
- 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
Public Functions
-
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})¶ Constructor for 2D projective geometry.
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- Parameters
[in] sourceToCenterOfRotation
: distance from source to the center of rotation (along y axis)[in] centerOfRotationToDetector
: distance from center of rotation to detector (along y axis)[in] angle
: rotation angle (in radians)[in] volData
: descriptor for the 2d volume[in] sinoData
: descriptor for the sinogram[in] offset
: offset of the principal point [default 0][in] centerOfRotOffset
: offset of the center of rotation
-
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})¶ Constructor for 3D projective geometry using Euler angles.
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)
- Parameters
[in] sourceToCenterOfRotation
: distance from source to the center of rotation (along z axis)[in] centerOfRotationToDetector
: distance from center of rotation to detector (along z axis)[in] volData
: descriptor for the 3d volume[in] sinoData
: descriptor for the sinogram[in] angles
: (gamma -> around y’’-axis, beta -> around z’ axis, alpha -> around y axis) in radians[in] offset
: offset of the principal point param[in] centerOfRotOffset offset of the center 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, const RealMatrix_t &R, 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))¶ Constructor for 3D projective geometry using a 3x3 rotation matrix.
- Parameters
[in] sourceToCenterOfRotation
: distance from source to the center of rotation (along z axis)[in] centerOfRotationToDetector
: distance from center of rotation to detector (along z axis)[in] volumeDescriptor
: descriptor for the 3d volume[in] sinoDescriptor
: descriptor for the sinogram[in] R
: a 3x3 rotation matrix[in] px
: offset of the principal point in x-direction [default 0][in] py
: offset of the principal point in y-direction [default 0][in] centerOfRotationOffsetX
: offset of the center of rotation in x direction [default 0][in] centerOfRotationOffsetY
: offset of the center of rotation in y direction [default 0][in] centerOfRotationOffsetZ
: offset of the center of rotation in z direction [default 0]
-
const RealMatrix_t &
getProjectionMatrix
() const¶ Return the projection matrix.
- Return
projection matrix
-
const RealMatrix_t &
getInverseProjectionMatrix
() const¶ Return the inverse of the projection matrix.
- Return
the inverse of the projection matrix
-
const RealVector_t &
getCameraCenter
() const¶ Return the camera center corresponding to the projection matrix.
- Return
the camera center (as a coordinate vector)
-
const RealMatrix_t &
getRotationMatrix
() const¶ Return the rotation matrix corresponding to the projection matrix.
- Return
the rotation matrix
Private Functions
-
void
buildMatrices
()¶ build the projection matrix, its inverse and the camera center
Private Members
-
index_t
_objectDimension
¶ the dimension of the object space / volume (either 2 or 3)
-
RealMatrix_t
_P
¶ the projection matrix (= [_K|0] * [_R|_t] * _S)
-
RealMatrix_t
_Pinv
¶ the inverse of the projection matrix
-
RealMatrix_t
_K
¶ the intrinsic parameters _K
-
RealMatrix_t
_R
¶ the rotation matrix
-
RealVector_t
_t
¶ the translation in object space
-
RealMatrix_t
_S
¶ the scaling in object space
-
RealVector_t
_C
¶ the camera center _C
BinaryMethod¶
-
template<typename
data_t
= real_t>
classelsa
::
BinaryMethod
: public elsa::LinearOperator<real_t>¶ Operator representing the discretized X-ray transform in 2d/3d using a simplistic binary hit/miss method.
The volume is traversed along the rays as specified by the
Geometry. Each ray is traversed in a continguous fashion (i.e. along long voxel borders, not diagonally) and each traversed voxel is counted as a hit with weight 1.- Author
Tobias Lasser - initial code, modernization
- Author
David Frank - rewrite and fixes
- Author
Maximilian Hornung - modularization
- Author
Nikola Dinev - fixes
- Template Parameters
data_t
: data type for the domain and range of the operator, defaulting to real_t
The geometry is represented as a list of projection matrices (see class Geometry), one for each acquisition pose.
Forward projection is accomplished using apply(), backward projection using applyAdjoint(). This projector is matched.
Warning: This method is not particularly accurate!
Public Functions
-
BinaryMethod
(const VolumeDescriptor &domainDescriptor, const DetectorDescriptor &rangeDescriptor)¶ Constructor for the binary voxel traversal method.
The domain is expected to be 2 or 3 dimensional (volSizeX, volSizeY, [volSizeZ]), the range is expected to be matching the domain (detSizeX, [detSizeY], acqPoses).
- Parameters
[in] domainDescriptor
: describing the domain of the operator (the volume)[in] rangeDescriptor
: describing the range of the operator (the sinogram)
-
~BinaryMethod
() override = default¶ default destructor
Protected Functions
-
BinaryMethod
(const BinaryMethod<data_t>&) = default¶ default copy constructor, hidden from non-derived classes to prevent potential slicing
-
void
applyImpl
(const DataContainer<data_t> &x, DataContainer<data_t> &Ax) const override¶ apply the binary method (i.e. forward projection)
-
void
applyAdjointImpl
(const DataContainer<data_t> &y, DataContainer<data_t> &Aty) const override¶ apply the adjoint of the binary method (i.e. backward projection)
-
BinaryMethod<data_t> *
cloneImpl
() const override¶ implement the polymorphic clone operation
-
bool
isEqual
(const LinearOperator<data_t> &other) const override¶ implement the polymorphic comparison operation
Private Types
-
using
Ray
= Eigen::ParametrizedLine<real_t, Eigen::Dynamic>¶ convenience typedef for ray
Private Functions
-
template<bool
adjoint
>
voidtraverseVolume
(const DataContainer<data_t> &vector, DataContainer<data_t> &result) const¶ the traversal routine (for both apply/applyAdjoint)
Private Members
-
BoundingBox
_boundingBox
¶ the bounding box of the volume
-
DetectorDescriptor &
_detectorDescriptor
¶ Reference to DetectorDescriptor stored in LinearOperator.
-
VolumeDescriptor &
_volumeDescriptor
¶ Reference to VolumeDescriptor stored in LinearOperator.
SiddonsMethod¶
-
template<typename
data_t
= real_t>
classelsa
::
SiddonsMethod
: public elsa::LinearOperator<real_t>¶ Operator representing the discretized X-ray transform in 2d/3d using Siddon’s method.
The volume is traversed along the rays as specified by the
Geometry. Each ray is traversed in a continguous fashion (i.e. along long voxel borders, not diagonally) and each traversed voxel is counted as a hit with weight according to the length of the path of the ray through the voxel.- Author
David Frank - initial code
- Author
Nikola Dinev - modularization, fixes
- Template Parameters
data_t
: data type for the domain and range of the operator, defaulting to real_t
The geometry is represented as a list of projection matrices (see class Geometry), one for each acquisition pose.
Forward projection is accomplished using apply(), backward projection using applyAdjoint(). This projector is matched.
Public Functions
-
SiddonsMethod
(const VolumeDescriptor &domainDescriptor, const DetectorDescriptor &rangeDescriptor)¶ Constructor for Siddon’s method traversal.
The domain is expected to be 2 or 3 dimensional (volSizeX, volSizeY, [volSizeZ]), the range is expected to be matching the domain (detSizeX, [detSizeY], acqPoses).
- Parameters
[in] domainDescriptor
: describing the domain of the operator (the volume)[in] rangeDescriptor
: describing the range of the operator (the sinogram)
-
~SiddonsMethod
() override = default¶ default destructor
Protected Functions
-
SiddonsMethod
(const SiddonsMethod<data_t>&) = default¶ default copy constructor, hidden from non-derived classes to prevent potential slicing
-
void
applyImpl
(const DataContainer<data_t> &x, DataContainer<data_t> &Ax) const override¶ apply Siddon’s method (i.e. forward projection)
-
void
applyAdjointImpl
(const DataContainer<data_t> &y, DataContainer<data_t> &Aty) const override¶ apply the adjoint of Siddon’s method (i.e. backward projection)
-
SiddonsMethod<data_t> *
cloneImpl
() const override¶ implement the polymorphic clone operation
-
bool
isEqual
(const LinearOperator<data_t> &other) const override¶ implement the polymorphic comparison operation
Private Types
-
using
Ray
= Eigen::ParametrizedLine<real_t, Eigen::Dynamic>¶ convenience typedef for ray
Private Functions
-
template<bool
adjoint
>
voidtraverseVolume
(const DataContainer<data_t> &vector, DataContainer<data_t> &result) const¶ the traversal routine (for both apply/applyAdjoint)
Private Members
-
BoundingBox
_boundingBox
¶ the bounding box of the volume
-
DetectorDescriptor &
_detectorDescriptor
¶ Reference to DetectorDescriptor stored in LinearOperator.
-
VolumeDescriptor &
_volumeDescriptor
¶ Reference to VolumeDescriptor stored in LinearOperator.
JosephsMethod¶
-
template<typename
data_t
= real_t>
classelsa
::
JosephsMethod
: public elsa::LinearOperator<real_t>¶ Operator representing the discretized X-ray transform in 2d/3d using Joseph’s method.
The volume is traversed along the rays as specified by the
Geometry. For interior voxels the sampling point is located in the middle of the two planes orthogonal to the main direction of the ray. For boundary voxels the sampling point is located at the center of the ray intersection with the voxel.- Author
Christoph Hahn - initial implementation
- Author
Maximilian Hornung - modularization
- Author
Nikola Dinev - fixes
- Template Parameters
data_t
: data type for the domain and range of the operator, defaulting to real_t
The geometry is represented as a list of projection matrices (see class Geometry), one for each acquisition pose.
Two modes of interpolation are available: NN (NearestNeighbours) takes the value of the pixel/voxel containing the point LINEAR performs linear interpolation for the nearest 2 pixels (in 2D) or the nearest 4 voxels (in 3D).
Forward projection is accomplished using apply(), backward projection using applyAdjoint(). This projector is matched.
Public Types
Public Functions
-
JosephsMethod
(const VolumeDescriptor &domainDescriptor, const DetectorDescriptor &rangeDescriptor, Interpolation interpolation = Interpolation::LINEAR)¶ Constructor for Joseph’s traversal method.
The domain is expected to be 2 or 3 dimensional (volSizeX, volSizeY, [volSizeZ]), the range is expected to be matching the domain (detSizeX, [detSizeY], acqPoses).
- Parameters
[in] domainDescriptor
: describing the domain of the operator (the volume)[in] rangeDescriptor
: describing the range of the operator (the sinogram)[in] interpolation
: enum specifying the interpolation mode
-
~JosephsMethod
() = default¶ default destructor
Protected Functions
-
JosephsMethod
(const JosephsMethod<data_t>&) = default¶ default copy constructor, hidden from non-derived classes to prevent potential slicing
-
void
applyImpl
(const DataContainer<data_t> &x, DataContainer<data_t> &Ax) const override¶ apply Joseph’s method (i.e. forward projection)
-
void
applyAdjointImpl
(const DataContainer<data_t> &y, DataContainer<data_t> &Aty) const override¶ apply the adjoint of Joseph’s method (i.e. backward projection)
-
JosephsMethod<data_t> *
cloneImpl
() const override¶ implement the polymorphic clone operation
-
bool
isEqual
(const LinearOperator<data_t> &other) const override¶ implement the polymorphic comparison operation
Private Types
-
using
Ray
= Eigen::ParametrizedLine<real_t, Eigen::Dynamic>¶ convenience typedef for ray
Private Functions
-
template<bool
adjoint
>
voidtraverseVolume
(const DataContainer<data_t> &vector, DataContainer<data_t> &result) const¶ the traversal routine (for both apply/applyAdjoint)
-
template<bool
adjoint
>
voidlinear
(const DataContainer<data_t> &vector, DataContainer<data_t> &result, const RealVector_t &fractionals, index_t domainDim, const IndexVector_t ¤tVoxel, float intersection, index_t from, index_t to, int mainDirection) const¶ Linear interpolation, works in any dimension.
- Template Parameters
adjoint
: true for backward projection, false for forward
- Parameters
[in] vector
: the input DataContainer[out] result
: DataContainer for results[in] fractionals
: the fractional numbers used in the interpolation[in] domainDim
: number of dimensions[in] currentVoxel
: coordinates of voxel for interpolation[in] intersection
: weighting for the interpolated values depending on the incidence angle[in] from
: index of the current vector position[in] to
: index of the current result position[in] mainDirection
: specifies the main direction of the ray
Private Members
-
BoundingBox
_boundingBox
¶ the bounding box of the volume
-
DetectorDescriptor &
_detectorDescriptor
¶ Reference to DetectorDescriptor stored in LinearOperator.
-
VolumeDescriptor &
_volumeDescriptor
¶ Reference to VolumeDescriptor stored in LinearOperator.
-
Interpolation
_interpolation
¶ the interpolation mode
Implementation Details¶
BoundingBox¶
-
struct
elsa
::
BoundingBox
¶ Helper struct defining a 2d or 3d axis-aligned bounding box (or in short: AABB).
- Author
David Frank - initial code
- Author
Tobias Lasser - minor changes
Public Functions
-
BoundingBox
(const IndexVector_t &volumeDimensions)¶ Construct AABB of particular size.
- Parameters
[in] volumeDimensions
: the number of coefficients per volume dimension
Intersection¶
-
class
elsa
::
Intersection
¶ The intersection class computes intersections between rays and axis-aligned bounding boxes (AABBs)
- Author
Tobias Lasser - initial code, modernization
- Author
David Frank - various fixes
- Author
Maximilian Hornung - modularization
- Author
Nikola Dinev - various fixes
Public Static Functions
-
std::optional<IntersectionResult>
withRay
(const BoundingBox &aabb, const Ray &r)¶ Compute entry and exit point of ray in a volume (given as an AABB)
If the ray is running along a border of the bounding box, the lower bound will be counted as in the bounding and the upper bound will be counted as outside.
- Return
nullopt if the volume is not hit, otherwise IntersectionResult with entry/exit parameters tmin/tmax
- Parameters
[in] aabb
: the volume specified through an axis-aligned bounding box[in] r
: the ray which we test for intersection with aabb
Method adapted from https://tavianator.com/fast-branchless-raybounding-box-intersections-part-2-nans/
Private Types
-
using
Ray
= Eigen::ParametrizedLine<real_t, Eigen::Dynamic>¶ the type for a ray using Eigen
TraverseAABB¶
-
class
elsa
::
TraverseAABB
¶ Class implementing a voxel traversal of a volume (AABB) along a ray.
This traversal always proceeds along “long” voxel edges, it will never “jump diagonally”. The method is based on J. Amantides, A. Woo: A Fast Voxel Traversal Algorithm for Ray Tracing.
- Author
Tobias Lasser - initial code
- Author
David Frank - major rewrite
- Author
Maximilian Hornung - modularization
- Author
Nikola Dinev - fixes
Public Functions
-
TraverseAABB
(const BoundingBox &aabb, const Ray &r)¶ Constructor for traversal, accepting bounding box and ray.
- Parameters
[in] aabb
: axis-aligned boundary box describing the volume[in] r
: the ray to be traversed
-
void
updateTraverse
()¶ Update the traverser status by taking the next traversal step.
-
real_t
updateTraverseAndGetDistance
()¶ Update the traverser status by taking the next traversal step, return the distance of step.
- Return
the distance of the step taken
-
bool
isInBoundingBox
() const¶ Return whether the traversal is still in the bounding box.
- Return
true if still in bounding box, false otherwise
-
IndexVector_t
getCurrentVoxel
() const¶ Return the current voxel in the bounding box.
- Return
coordinate vector
Private Types
-
using
Ray
= Eigen::ParametrizedLine<real_t, Eigen::Dynamic>¶ convenience type def for ray
Private Functions
-
void
calculateAABBIntersections
(const Ray &r)¶ compute the entry and exit points of ray r with the volume (aabb)
-
void
initStepDirection
(const RealVector_t &rd)¶ setup the step directions (which is basically the sign of the ray direction rd)
-
void
selectClosestVoxel
()¶ select the closest voxel to the entry point
-
void
initDelta
(const RealVector_t &rd)¶ setup the step sizes considering the ray direction rd
-
void
initMax
(const RealVector_t &rd)¶ setup the maximum step parameters considering the ray direction rd
-
bool
isCurrentPositionInAABB
(index_t index) const¶ check if the current index is still in the bounding box
Private Members
-
BoundingBox
_aabb
¶ the volume / axis-aligned bounding box
-
bool
_isInAABB
= {false}¶ flag if traverser still in bounding box
-
real_t
_tExit
= {0.0}¶ the current step parameter exiting the current voxel
TraverseAABBJosephsMethod¶
-
class
elsa
::
TraverseAABBJosephsMethod
¶ Special traversal wrapper for Joseph’s method.
Public Functions
-
TraverseAABBJosephsMethod
(const BoundingBox &aabb, const Ray &r)¶ Constructor for traversal, accepting bounding box and ray.
- Parameters
[in] aabb
: axis-aligned boundary box describing the volume[in] r
: the ray to be traversed
-
void
updateTraverse
()¶ Update the traverser status by taking the next traversal step.
-
bool
isInBoundingBox
() const¶ Return whether the traversal is still in the bounding box.
- Return
true if still in bounding box, false otherwise
-
const RealVector_t &
getFractionals
() const¶ Returns the fractional part of the current position wrt. the center of the pixel/voxel.
For example: if the current position is (15.4, 21.7, 23.0), then the current voxel is (15,21,23) and the fractional parts are (15.4, 21.7, 23.0) - (15.5, 21.5, 23.5) = (-0.1, 0.2, -0.5)
- Return
RealVector_t fractionals
-
int
getIgnoreDirection
() const¶ Get direction to be ignored during interpolation.
The sampling points are chosen so that one of the directions can be ignored during interpolation, effectively reducing the number of operations needed for the calculation of the interpolated value by half. For interior points, the direction to be ignored is always the main direction. For the entry and exit points it’s the entry resp. exit direction.
- Return
int index of the direction to be ignored
-
real_t
getIntersectionLength
() const¶ Get the intersection length for the current step.
- Return
real_t intersection length
-
IndexVector_t
getCurrentVoxel
() const¶ Get the coordinates of the voxel which contains the current sampling point.
- Return
IndexVector_t coordinates of the voxel which contains the current sampling point
Private Types
-
enum
STAGE
¶ Possible stages of the traversal.
FIRST: currently handling the first sampling INTERIOR: currently handling an interior point LAST: currently handling the last sampling point
If the ray only intersects a single pixel, the sole sampling point is considered last.
Values:
-
enumerator
FIRST
¶
-
enumerator
INTERIOR
¶
-
enumerator
LAST
¶
-
enumerator
-
using
Ray
= Eigen::ParametrizedLine<real_t, Eigen::Dynamic>¶ convenience type def for ray
Private Functions
-
void
initFractionals
(const RealVector_t ¤tPosition)¶ initialize fractionals from a coordinate vector
-
void
initStepDirection
(const RealVector_t &rd)¶ setup the step directions (which is basically the sign of the ray direction rd)
-
real_t
calculateAABBIntersections
(const Ray &ray)¶ compute the entry and exit points of ray r with the volume (aabb), returns the length of the intersection
-
void
selectClosestVoxel
(const RealVector_t ¤tPosition)¶ select the closest voxel to the entry point (considering the current position)
-
bool
isCurrentPositionInAABB
(index_t index) const¶ check if the current index is still in the bounding box
-
void
moveToFirstSamplingPoint
(const RealVector_t &rd, real_t intersectionLength)¶ advances the traversal algorithm to the first sampling point
Private Members
-
BoundingBox
_aabb
¶ the volume / axis-aligned bounding box
-
RealVector_t
_fractionals
= {_aabb._dim}¶ the fractional parts of the current position coordinates (actually frac(_currentPos) - 0.5)
-
bool
_isInAABB
= {false}¶ flag if traverser still in bounding box
-
real_t
_intersectionLength
= {0}¶ length of ray segment currently being handled
-
int
_ignoreDirection
= {-1}¶ index of direction for which no interpolation needs to be performed
-
int
_exitDirection
= {-1}¶ index of direction in which the ray exits the volume
-