SiddonsMethodCUDA.h 4.78 KB
Newer Older
Nikola Dinev's avatar
Nikola Dinev committed
1
2
#pragma once

3
#include <cuda_runtime.h>
Nikola Dinev's avatar
Nikola Dinev committed
4

5
#include "elsaDefines.h"
Nikola Dinev's avatar
Nikola Dinev committed
6
7
8
9
#include "LinearOperator.h"
#include "Geometry.h"
#include "BoundingBox.h"

10
#include "TraverseSiddonsCUDA.cuh"
Nikola Dinev's avatar
Nikola Dinev committed
11
12
13
14

namespace elsa
{
    /**
Jens Petit's avatar
Jens Petit committed
15
16
     * \brief GPU-operator representing the discretized X-ray transform in 2d/3d using Siddon's
     * method.
17
18
19
20
21
     *
     * \author Nikola Dinev
     *
     * \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
23
24
25
     * 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.
26
     *
Jens Petit's avatar
Jens Petit committed
27
28
     * The geometry is represented as a list of projection matrices (see class Geometry), one for
     * each acquisition pose.
29
30
31
     *
     * Forward projection is accomplished using apply(), backward projection using applyAdjoint().
     * This projector is matched.
Jens Petit's avatar
Jens Petit committed
32
33
34
     *
     * Currently only utilizes a single GPU. Volume and images should both fit in device memory at
     * the same time.
Nikola Dinev's avatar
Nikola Dinev committed
35
     */
36
    template <typename data_t = real_t>
Nikola Dinev's avatar
Nikola Dinev committed
37
38
39
40
    class SiddonsMethodCUDA : public LinearOperator<data_t>
    {
    public:
        /**
Jens Petit's avatar
Jens Petit committed
41
42
43
44
45
46
47
48
49
50
51
52
         * \brief Constructor for Siddon's method traversal.
         *
         * \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
         *
         * 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).
         */
        SiddonsMethodCUDA(const DataDescriptor& domainDescriptor,
                          const DataDescriptor& rangeDescriptor,
                          const std::vector<Geometry>& geometryList);
Nikola Dinev's avatar
Nikola Dinev committed
53

54
        /// destructor
Nikola Dinev's avatar
Nikola Dinev committed
55
56
57
        ~SiddonsMethodCUDA();

    protected:
58
59
60
        /// default copy constructor, hidden from non-derived classes to prevent potential slicing
        SiddonsMethodCUDA(const SiddonsMethodCUDA<data_t>&) = default;

61
        /// apply Siddon's method (i.e. forward projection)
62
        void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
Nikola Dinev's avatar
Nikola Dinev committed
63

64
        /// apply the adjoint of Siddon's method (i.e. backward projection)
65
66
        void applyAdjointImpl(const DataContainer<data_t>& y,
                              DataContainer<data_t>& Aty) const override;
Nikola Dinev's avatar
Nikola Dinev committed
67
68
69
70
71
72
73
74

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

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

    private:
75
76
77
78
79
80
81
        /// the bounding box of the volume
        BoundingBox _boundingBox;

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

        /// threads per block used in the kernel execution configuration
82
        const int _threadsPerBlock = TraverseSiddonsCUDA<data_t>::MAX_THREADS_PER_BLOCK;
Nikola Dinev's avatar
Nikola Dinev committed
83

84
        /// inverse of of projection matrices; stored column-wise on GPU
Nikola Dinev's avatar
Nikola Dinev committed
85
        cudaPitchedPtr _projInvMatrices;
86
87

        /// ray origins for each acquisition angle
Nikola Dinev's avatar
Nikola Dinev committed
88
89
        cudaPitchedPtr _rayOrigins;

90
        /// sets up and starts the kernel traversal routine (for both apply/applyAdjoint)
Jens Petit's avatar
Jens Petit committed
91
        template <bool adjoint>
Nikola Dinev's avatar
Nikola Dinev committed
92
93
        void traverseVolume(void* volumePtr, void* sinoPtr) const;

94
95
        /**
         * \brief Copies contents of a 3D data container between GPU and host memory
Jens Petit's avatar
Jens Petit committed
96
         *
97
98
         * \tparam direction specifies the direction of the copy operation
         * \tparam async whether the copy should be performed asynchronously wrt. the host
Jens Petit's avatar
Jens Petit committed
99
         *
100
101
102
         * \param hostData pointer to host data
         * \param gpuData pointer to gpu data
         * \param[in] extent specifies the amount of data to be copied
Jens Petit's avatar
Jens Petit committed
103
104
105
106
107
         *
         * Note that hostData is expected to be a pointer to a linear memory region with no padding
         * between dimensions - e.g. the data in DataContainer is stored as a vector with no extra
         * padding, and the pointer to the start of the memory region can be retrieved as follows:
         *
108
109
110
         * DataContainer x;
         * void* hostData = (void*)&x[0];
         */
Jens Petit's avatar
Jens Petit committed
111
112
113
        template <cudaMemcpyKind direction, bool async = true>
        void copy3DDataContainerGPU(void* hostData, const cudaPitchedPtr& gpuData,
                                    const cudaExtent& extent) const;
Nikola Dinev's avatar
Nikola Dinev committed
114
115
116
117
118

        /// lift from base class
        using LinearOperator<data_t>::_domainDescriptor;
        using LinearOperator<data_t>::_rangeDescriptor;
    };
Jens Petit's avatar
Jens Petit committed
119
} // namespace elsa