JosephsMethodCUDA.h 5.95 KB
Newer Older
1
2
3
4
5
6
#pragma once

#include <cstring>

#include <cuda_runtime.h>

7
#include "elsaDefines.h"
8
9
10
11
#include "LinearOperator.h"
#include "Geometry.h"
#include "BoundingBox.h"

12
#include "TraverseJosephsCUDA.cuh"
13
14

namespace elsa
Jens Petit's avatar
Jens Petit committed
15
{
16
    /**
Jens Petit's avatar
Jens Petit committed
17
18
     * \brief GPU-operator representing the discretized X-ray transform in 2d/3d using Joseph's
     * method.
19
20
     *
     * \author Nikola Dinev
Jens Petit's avatar
Jens Petit committed
21
     *
22
     * \tparam data_t data type for the domain and range of the operator, defaulting to real_t
Jens Petit's avatar
Jens Petit committed
23
     *
24
25
26
27
     * 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
28
29
30
31
     *
     * The geometry is represented as a list of projection matrices (see class Geometry), one for
     * each acquisition pose.
     *
32
     * Forward projection is accomplished using apply(), backward projection using applyAdjoint().
Jens Petit's avatar
Jens Petit committed
33
34
35
36
37
38
     * The projector provides two implementations for the backward projection. The slow version is
     * matched, while the fast one is not.
     *
     * Currently only utilizes a single GPU. Volume and images should both fit in device memory at
     * the same time.
     *
39
     * \warning Hardware interpolation is only supported for JosephsMethodCUDA<float>
Jens Petit's avatar
Jens Petit committed
40
41
     * \warning Hardware interpolation is significantly less accurate than the software
     * interpolation
42
     */
Jens Petit's avatar
Jens Petit committed
43
    template <typename data_t = real_t>
44
45
46
47
    class JosephsMethodCUDA : public LinearOperator<data_t>
    {
    public:
        /**
48
         * \brief Constructor for Joseph's traversal method.
Jens Petit's avatar
Jens Petit committed
49
         *
50
51
52
         * \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
Jens Petit's avatar
Jens Petit committed
53
54
55
         * \param[in] fast performs fast backward projection if set, otherwise matched; forward
         * projection is unaffected
         *
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).
58
         */
Jens Petit's avatar
Jens Petit committed
59
60
61
        JosephsMethodCUDA(const DataDescriptor& domainDescriptor,
                          const DataDescriptor& rangeDescriptor,
                          const std::vector<Geometry>& geometryList, bool fast = true);
62

63
        /// destructor
64
65
66
        virtual ~JosephsMethodCUDA();

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

70
        /// apply Joseph's method (i.e. forward projection)
71
        void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
72

73
        /// apply the adjoint of Joseph's method (i.e. backward projection)
74
75
        void applyAdjointImpl(const DataContainer<data_t>& y,
                              DataContainer<data_t>& Aty) const override;
76
77
78
79
80
81

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

        /// implement the polymorphic comparison operation
        bool isEqual(const LinearOperator<data_t>& other) const override;
Jens Petit's avatar
Jens Petit committed
82

83
    private:
84
85
86
87
88
89
90
        /// the bounding box of the volume
        BoundingBox _boundingBox;

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

        /// threads per block used in the kernel execution configuration
91
        const int _threadsPerBlock = TraverseJosephsCUDA<data_t>::MAX_THREADS_PER_BLOCK;
92
93

        /// flag specifying which version of the backward projection should be used
94
95
        const bool _fast;

96
        /// inverse of of projection matrices; stored column-wise on GPU
97
98
        cudaPitchedPtr _projInvMatrices;

99
        /// projection matrices; stored column-wise on GPU
100
101
        cudaPitchedPtr _projMatrices;

102
        /// ray origins for each acquisition angle
103
104
        cudaPitchedPtr _rayOrigins;

105
        /// convenience typedef for cuda array flags
106
107
108
        using cudaArrayFlags = unsigned int;

        /**
109
         * \brief Copies contents of a 3D data container between GPU and host memory
Jens Petit's avatar
Jens Petit committed
110
         *
111
112
         * \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
113
         *
114
115
116
         * \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
117
118
119
120
121
         *
         * 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:
         *
122
123
         * DataContainer x;
         * void* hostData = (void*)&x[0];
124
         */
Jens Petit's avatar
Jens Petit committed
125
126
127
        template <cudaMemcpyKind direction, bool async = true>
        void copy3DDataContainer(void* hostData, const cudaPitchedPtr& gpuData,
                                 const cudaExtent& extent) const;
128
129

        /**
130
         * \brief Copies the entire contents of DataContainer to the GPU texture memory
Jens Petit's avatar
Jens Petit committed
131
132
133
134
         *
         * \tparam cudaArrayFlags flags used for the creation of the cudaArray which will contain
         * the data
         *
135
         * \param[in] hostData the host data container
Jens Petit's avatar
Jens Petit committed
136
         *
137
         * \returns a pair of the created texture object and its associated cudaArray
138
139
         */
        template <cudaArrayFlags flags = 0U>
Jens Petit's avatar
Jens Petit committed
140
141
        std::pair<cudaTextureObject_t, cudaArray*>
            copyTextureToGPU(const DataContainer<data_t>& hostData) const;
142
143
144
145
146

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