elsa projectors

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

bool operator==(const Geometry &other) const

comparison operator

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>
class elsa::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>
void traverseVolume(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>
class elsa::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>
void traverseVolume(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>
class elsa::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

enum Interpolation

Available interpolation modes.

Values:

enumerator NN
enumerator LINEAR

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>
void traverseVolume(const DataContainer<data_t> &vector, DataContainer<data_t> &result) const

the traversal routine (for both apply/applyAdjoint)

template<bool adjoint>
void linear(const DataContainer<data_t> &vector, DataContainer<data_t> &result, const RealVector_t &fractionals, index_t domainDim, const IndexVector_t &currentVoxel, 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

Public Members

index_t _dim

the number of dimensions (2 or 3)

RealVector_t _min = {_dim}

the front corner of the box

RealVector_t _max = {_dim}

the back corner of the box

IndexVector_t _voxelCoordToIndexVector = {_dim}

helper to convert coordinates to indices

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

RealVector_t _entryPoint = {_aabb._dim}

the entry point parameters of the ray in the aabb

RealVector_t _stepDirection = {_aabb._dim}

the step direction of the traverser

RealVector_t _currentPos = {_aabb._dim}

the current position of the traverser in the aabb

RealVector_t _tMax = {_aabb._dim}

the current maximum step parameter along the ray

RealVector_t _tDelta = {_aabb._dim}

the step sizes for the step parameter along the ray

bool _isInAABB = {false}

flag if traverser still in bounding box

real_t _tExit = {0.0}

the current step parameter exiting the current voxel

const RealVector_t _EPS{RealVector_t(_aabb._dim).setConstant(std::numeric_limits<real_t>::epsilon())}

constant vector containing epsilon

const RealVector_t _MAX{RealVector_t(_aabb._dim).setConstant(std::numeric_limits<real_t>::max())}

constant vector containing the maximum number

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
using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>

convenience type def for ray

Private Functions

void initFractionals(const RealVector_t &currentPosition)

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 &currentPosition)

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 _entryPoint = {_aabb._dim}

the entry point parameters of the ray in the aabb

RealVector_t _exitPoint = {_aabb._dim}

the exit point parameters of the ray in the aabb

RealVector_t _currentPos = {_aabb._dim}

the current position of the traverser in the aabb

RealVector_t _stepDirection = {_aabb._dim}

the step direction of the traverser

RealVector_t _nextStep = {_aabb._dim}

the step sizes for the next step along the ray

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

STAGE _stage = {FIRST}

stage at which the traversal currently is