ConeBeamGeometry

class odl.applications.tomo.geometry.conebeam.ConeBeamGeometry(apart, dpart, src_radius, det_radius, det_curvature_radius=None, pitch=0, axis=(0, 0, 1), src_shift_func=None, det_shift_func=None, **kwargs)[source]

Bases: DivergentBeamGeometry, AxisOrientedGeometry

Cone beam geometry with circular/helical source curve.

The source moves along a spiral oriented along a fixed axis, with radius src_radius in the azimuthal plane and a given pitch. The detector reference point is opposite to the source, i.e. in the point at distance src_rad + det_rad on the line in the azimuthal plane through the source point and axis.

The motion parameter is the 1d rotation angle parameterizing source and detector positions simultaneously.

In the standard configuration, the rotation axis is (0, 0, 1), the initial source-to-detector vector is (0, 1, 0), and the initial detector axes are [(1, 0, 0), (0, 0, 1)].

For details, check the online docs.

__init__(apart, dpart, src_radius, det_radius, det_curvature_radius=None, pitch=0, axis=(0, 0, 1), src_shift_func=None, det_shift_func=None, **kwargs)[source]

Initialize a new instance.

Parameters

apart1-dim. RectPartition: Partition of the angle interval.
dpart2-dim. RectPartition: Partition of the detector parameter rectangle.
src_radiusnonnegative float: Radius of the source circle.
det_radiusnonnegative float: Radius of the detector circle. Must be nonzero if src_radius is zero.
det_curvature_radius2-tuple of nonnegative floats, optional: Radius or radii of the detector curvature. If None, a flat detector is used. If (r, None) or (r, float('inf')), a cylindrical detector is used. If (r1, r2), a spherical detector is used.
pitchfloat, optional: Constant distance along axis that a point on the helix traverses when increasing the angle parameter by 2 * pi. The default case pitch=0 results in a circular cone beam geometry.
axisarray-like, shape (3,), optional: Vector defining the fixed rotation axis of this geometry.
src_shift_funccallable, optional: Function with signature src_shift_func(angle) -> shift returning a source shift for a given angle. Each shift is interpreted as a vector [shift_d, shift_t, shift_r], where “d”, “t” and “r” denote shifts along the following directions: detector-to-source, tangent to the rotation (projected on plane perpendicular to rotation axis), rotation axis.
det_shift_funccallable, optional: Function with signature det_shift_func(angle) -> shift returning a detector shift for a given angle. Each shift is interpreted as a vector [shift_d, shift_t, shift_r], where “d”, “t” and “r” denote shifts along the following directions: source-to-detector, tangent to the rotation (projected on plane perpendicular to rotation axis), rotation axis.

Other Parameters

offset_along_axisfloat, optional: Scalar offset along the axis at angle=0, i.e., the translation along the axis at angle 0 is offset_along_axis * axis. Default: 0.
src_to_det_initarray-like, shape (3,), optional: Initial state of the vector pointing from source to detector reference point. The zero vector is not allowed. The default depends on axis, see Notes.
det_axes_initlength-2-sequence of array-like’s, optional: Initial axes defining the detector orientation, provided as arrays with shape (3,). Default: depends on axis, see Notes.
translationarray-like, shape (3,), optional: Global translation of the geometry. This is added last in any method that computes an absolute vector, e.g., det_refpoint, and also shifts the axis of rotation. Default: (0, 0, 0)
check_boundsbool, optional: If True, methods computing vectors check input arguments. Checks are vectorized and add only a small overhead. Default: True

Notes

In the default configuration, the rotation axis is (0, 0, 1), the initial source-to-detector direction is (0, 1, 0), and the default detector axes are [(1, 0, 0), (0, 0, 1)]. If a different axis is provided, the new default initial position and the new default axes are the computed by rotating the original ones by a matrix that transforms (0, 0, 1) to the new (normalized) axis. This matrix is calculated with the rotation_matrix_from_to function. Expressed in code, we have

init_rot = rotation_matrix_from_to((0, 0, 1), axis)
src_to_det_init = init_rot.dot((0, 1, 0))
det_axes_init[0] = init_rot.dot((1, 0, 0))
det_axes_init[1] = init_rot.dot((0, 0, 1))

Examples

Initialization with default parameters and some (arbitrary) choices for pitch and radii:

>>> apart = odl.uniform_partition(0, 4 * np.pi, 10)
>>> dpart = odl.uniform_partition([-1, -1], [1, 1], (20, 20))
>>> geom = ConeBeamGeometry(
...     apart, dpart, src_radius=5, det_radius=10, pitch=2)
>>> geom.src_position(0)
array([ 0., -5.,  0.])
>>> geom.det_refpoint(0)
array([ 0., 10.,  0.])
>>> np.allclose(geom.src_position(2 * np.pi),
...             geom.src_position(0) + (0, 0, 2))  # z shift by pitch
True

Checking the default orientation:

>>> geom.axis
array([ 0.,  0.,  1.])
>>> geom.src_to_det_init
array([ 0.,  1.,  0.])
>>> geom.det_axes_init
array([[ 1.,  0.,  0.],
       [ 0.,  0.,  1.]])

Specifying curvature of the cylindrical detector:

>>> apart = odl.uniform_partition(0, 4 * np.pi, 10)
>>> dpart = odl.uniform_partition(
...     [-np.pi / 2, -1], [np.pi / 2, 1], (20, 20))
>>> geom = ConeBeamGeometry(apart, dpart,
...     src_radius=5, det_radius=10, det_curvature_radius=(10, None))
>>> # (10*sin(pi/2), 10*cos(pi/2), 1)
>>> np.round(geom.det_point_position(0, [ np.pi / 2, 1] ), 2)
array([ 10., 0., 1.])

Specifying curvature of the spherical detector:

>>> apart = odl.uniform_partition(0, 4 * np.pi, 10)
>>> dpart = odl.uniform_partition([-np.pi / 2, -np.pi / 4],
...                               [ np.pi / 2,  np.pi / 4], (20, 20))
>>> geom = ConeBeamGeometry(apart, dpart,
...     src_radius=5, det_radius=10, det_curvature_radius=(10, 10))
>>> # 10*( cos(pi/4), 0, sin(pi/4))
>>> np.round(geom.det_point_position(0, [ np.pi / 2, np.pi / 4] ), 2)
array([ 7.07, 0.  , 7.07])

Specifying an axis by default rotates the standard configuration to this position:

>>> e_x, e_y, e_z = np.eye(3)  # standard unit vectors
>>> geom = ConeBeamGeometry(
...     apart, dpart, src_radius=5, det_radius=10, pitch=2,
...     axis=(0, 1, 0))
>>> np.allclose(geom.axis, e_y)
True
>>> np.allclose(geom.src_to_det_init, -e_z)
True
>>> np.allclose(geom.det_axes_init, (e_x, e_y))
True
>>> geom = ConeBeamGeometry(
...     apart, dpart, src_radius=5, det_radius=10, pitch=2,
...     axis=(1, 0, 0))
>>> np.allclose(geom.axis, e_x)
True
>>> np.allclose(geom.src_to_det_init, e_y)
True
>>> np.allclose(geom.det_axes_init, (-e_z, e_x))
True

The initial source-to-detector vector and the detector axes can also be set explicitly:

>>> geom = ConeBeamGeometry(
...     apart, dpart, src_radius=5, det_radius=10, pitch=2,
...     src_to_det_init=(-1, 0, 0),
...     det_axes_init=((0, 1, 0), (0, 0, 1)))
>>> np.allclose(geom.axis, e_z)
True
>>> np.allclose(geom.src_to_det_init, -e_x)
True
>>> np.allclose(geom.det_axes_init, (e_y, e_z))
True

Specifying a flying focal spot and detector offset:

>>> apart = odl.uniform_partition(0, 2 * np.pi, 4)
>>> geom = ConeBeamGeometry(
...     apart, dpart,
...     src_radius=1, det_radius=5,
...     src_shift_func=lambda angle: odl.applications.tomo.flying_focal_spot(
...         angle, apart=apart, shifts=[(0, 0.1, 0), (0, 0, 0.1)]),
...     det_shift_func=lambda angle: [0.0, 0.05, 0.03])
>>> geom.src_shift_func(geom.angles)
array([[ 0. , 0.1, 0. ],
       [ 0. , 0. , 0.1],
       [ 0. , 0.1, 0. ],
       [ 0. , 0. , 0.1]])
>>> geom.det_shift_func(geom.angles)
[0.0, 0.05, 0.03]