JosephsMethod.h 5.47 KB
Newer Older
Nikola Dinev's avatar
Nikola Dinev committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#pragma once

#include "LinearOperator.h"
#include "Geometry.h"
#include "BoundingBox.h"

#include <vector>
#include <utility>

#include <Eigen/Geometry>

namespace elsa
{
    /**
     * \brief Operator representing the discretized X-ray transform in 2d/3d using Joseph's method.
     *
     * \author Christoph Hahn - initial implementation
     * \author Maximilian Hornung - modularization
     * \author Nikola Dinev - fixes
Jens Petit's avatar
Jens Petit committed
20
     *
Nikola Dinev's avatar
Nikola Dinev committed
21
     * \tparam data_t data type for the domain and range of the operator, defaulting to real_t
Jens Petit's avatar
Jens Petit committed
22
     *
Nikola Dinev's avatar
Nikola Dinev committed
23
24
25
26
     * 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.
Jens Petit's avatar
Jens Petit committed
27
28
29
30
     *
     * The geometry is represented as a list of projection matrices (see class Geometry), one for
     * each acquisition pose.
     *
Nikola Dinev's avatar
Nikola Dinev committed
31
32
33
     * 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)
Jens Petit's avatar
Jens Petit committed
34
35
     * or the nearest 4 voxels (in 3D).
     *
Nikola Dinev's avatar
Nikola Dinev committed
36
37
38
39
     * Forward projection is accomplished using apply(), backward projection using applyAdjoint().
     * This projector is matched.
     */
    template <typename data_t = real_t>
Jens Petit's avatar
Jens Petit committed
40
41
    class JosephsMethod : public LinearOperator<data_t>
    {
Nikola Dinev's avatar
Nikola Dinev committed
42
43
44
45
46
47
48
49
50
51
    public:
        /// Available interpolation modes
        enum class Interpolation { NN, LINEAR };

        /**
         * \brief Constructor for Joseph's traversal method.
         *
         * \param[in] domainDescriptor describing the domain of the operator (the volume)
         * \param[in] rangeDescriptor describing the range of the operator (the sinogram)
         * \param[in] geometryList vector containing the geometries for the acquisition poses
52
         * \param[in] interpolation enum specifying the interpolation mode
Jens Petit's avatar
Jens Petit committed
53
         *
Nikola Dinev's avatar
Nikola Dinev committed
54
55
56
57
         * 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).
         */
        JosephsMethod(const DataDescriptor& domainDescriptor, const DataDescriptor& rangeDescriptor,
Jens Petit's avatar
Jens Petit committed
58
59
                      const std::vector<Geometry>& geometryList,
                      Interpolation interpolation = Interpolation::LINEAR);
Nikola Dinev's avatar
Nikola Dinev committed
60
61
62
63
64

        /// default destructor
        ~JosephsMethod() = default;

    protected:
65
66
67
        /// default copy constructor, hidden from non-derived classes to prevent potential slicing
        JosephsMethod(const JosephsMethod<data_t>&) = default;

Nikola Dinev's avatar
Nikola Dinev committed
68
        /// apply Joseph's method (i.e. forward projection)
69
        void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
Nikola Dinev's avatar
Nikola Dinev committed
70

Nikola Dinev's avatar
Nikola Dinev committed
71
        /// apply the adjoint of Joseph's method (i.e. backward projection)
72
73
        void applyAdjointImpl(const DataContainer<data_t>& y,
                              DataContainer<data_t>& Aty) const override;
Nikola Dinev's avatar
Nikola Dinev committed
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

        /// implement the polymorphic clone operation
        JosephsMethod<data_t>* cloneImpl() const override;

        /// implement the polymorphic comparison operation
        bool isEqual(const LinearOperator<data_t>& other) const override;

    private:
        /// the bounding box of the volume
        BoundingBox _boundingBox;

        /// the geometry list
        std::vector<Geometry> _geometryList;

        /// the interpolation mode
        Interpolation _interpolation;

        /// the traversal routine (for both apply/applyAdjoint)
        template <bool adjoint>
Jens Petit's avatar
Jens Petit committed
93
94
        void traverseVolume(const DataContainer<data_t>& vector,
                            DataContainer<data_t>& result) const;
Nikola Dinev's avatar
Nikola Dinev committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

        /// convenience typedef for ray
        using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;

        /**
         * \brief computes the ray to the middle of the detector element
         *
         * \param[in] detectorIndex the index of the detector element
         * \param[in] dimension the dimension of the detector (1 or 2)
         *
         * \returns the ray
         */
        Ray computeRayToDetector(index_t detectorIndex, index_t dimension) const;

        /**
         * \brief  Linear interpolation, works in any dimension
         *
112
113
114
115
116
117
         * \param[in] vector the input DataContainer
         * \param[out] result DataContainer for results
         * \param[in] fractionals the fractional numbers used in the interpolation
         * \param[in] adjoint true for backward projection, false for forward
         * \param[in] domainDim number of dimensions
         * \param[in] currentVoxel coordinates of voxel for interpolation
Jens Petit's avatar
Jens Petit committed
118
119
120
         * \param[in] intersection weighting for the interpolated values depending on the incidence
         * angle \param[in] from index of the current vector position \param[in] to index of the
         * current result position \param[in] mainDirection specifies the main direction of the ray
Nikola Dinev's avatar
Nikola Dinev committed
121
         */
122
        void linear(const DataContainer<data_t>& vector, DataContainer<data_t>& result,
Jens Petit's avatar
Jens Petit committed
123
124
125
                    const RealVector_t& fractionals, bool adjoint, int domainDim,
                    const IndexVector_t& currentVoxel, float intersection, index_t from, index_t to,
                    int mainDirection) const;
Nikola Dinev's avatar
Nikola Dinev committed
126
127
128
129
130
131

        /// lift from base class
        using LinearOperator<data_t>::_domainDescriptor;
        using LinearOperator<data_t>::_rangeDescriptor;
    };
} // namespace elsa