Line data Source code
1 : #pragma once
2 :
3 : #include "elsaDefines.h"
4 : #include "Math.hpp"
5 : #include "Geometry.h"
6 : #include "Blobs.h"
7 :
8 : namespace elsa
9 : {
10 :
11 : namespace voxel
12 : {
13 :
14 : using RealVector2D_t = Eigen::Matrix<real_t, 2, 1>;
15 : using RealVector3D_t = Eigen::Matrix<real_t, 3, 1>;
16 : using RealVector4D_t = Eigen::Matrix<real_t, 4, 1>;
17 : using IndexVector2D_t = Eigen::Matrix<index_t, 2, 1>;
18 : using IndexVector3D_t = Eigen::Matrix<index_t, 3, 1>;
19 :
20 : template <index_t dim, typename data_t, size_t N = DEFAULT_PROJECTOR_LUT_SIZE>
21 : data_t classic_weight_function(Lut<data_t, N> lut, Eigen::Matrix<real_t, dim, 1> distance)
22 652111 : {
23 652111 : return lut(distance.norm());
24 652111 : }
25 :
26 : template <typename data_t, size_t N = DEFAULT_PROJECTOR_LUT_SIZE>
27 : data_t differential_weight_function_2D(Lut<data_t, N> lut,
28 : Eigen::Matrix<real_t, 1, 1> distance)
29 612928 : {
30 612928 : return lut(distance.norm()) * math::sgn(distance[0]);
31 612928 : }
32 :
33 : template <typename data_t, size_t N = DEFAULT_PROJECTOR_LUT_SIZE>
34 : data_t differential_weight_function_3D(Lut<data_t, N> lut,
35 : Eigen::Matrix<real_t, 2, 1> distance)
36 0 : {
37 0 : return lut(distance.norm()) * distance[0];
38 0 : }
39 :
40 : template <index_t dim, class Fn>
41 : void visitDetector(index_t domainIndex, index_t geomIndex, Geometry geometry, real_t radius,
42 : Eigen::Matrix<index_t, dim - 1, 1> detectorDims,
43 : Eigen::Matrix<index_t, dim, 1> volumeStrides, Fn apply)
44 13467070 : {
45 13467070 : using IndexVectorVolume_t = Eigen::Matrix<index_t, dim, 1>;
46 13467070 : using IndexVectorDetector_t = Eigen::Matrix<index_t, dim - 1, 1>;
47 13467070 : using RealVectorVolume_t = Eigen::Matrix<real_t, dim, 1>;
48 13467070 : using RealVectorDetector_t = Eigen::Matrix<real_t, dim - 1, 1>;
49 13467070 : using RealVectorHom_t = Eigen::Matrix<real_t, dim + 1, 1>;
50 :
51 : // compute coordinate from index
52 13467070 : IndexVectorVolume_t coordinate = detail::idx2Coord(domainIndex, volumeStrides);
53 :
54 : // Cast to real_t and shift to center of voxel according to origin
55 13467070 : RealVectorVolume_t coordinateShifted = coordinate.template cast<real_t>().array() + 0.5;
56 13467070 : RealVectorHom_t homogeneousVoxelCoord = coordinateShifted.homogeneous();
57 :
58 : // Project voxel center onto detector
59 13467070 : RealVectorDetector_t center =
60 13467070 : (geometry.getProjectionMatrix() * homogeneousVoxelCoord).hnormalized();
61 13467070 : center = center.array() - 0.5;
62 :
63 : // transform voxel into camera space (camera at 0,0) and use the norm of the vector to
64 : // compute the source voxel distance
65 13467070 : auto scaling = geometry.getSourceDetectorDistance()
66 13467070 : / (geometry.getExtrinsicMatrix() * homogeneousVoxelCoord).norm();
67 :
68 13467070 : auto radiusOnDetector = scaling * radius;
69 13467070 : IndexVectorDetector_t detector_max = detectorDims.array() - 1;
70 13467070 : index_t detectorZeroIndex = (detector_max[0] + 1) * geomIndex;
71 :
72 13467070 : IndexVectorDetector_t min_corner =
73 13467070 : (center.array() - radiusOnDetector).ceil().template cast<index_t>();
74 13467070 : min_corner = min_corner.cwiseMax(IndexVectorDetector_t::Zero());
75 13467070 : IndexVectorDetector_t max_corner =
76 13467070 : (center.array() + radiusOnDetector).floor().template cast<index_t>();
77 13467070 : max_corner = max_corner.cwiseMin(detector_max);
78 :
79 13467070 : RealVectorDetector_t current = min_corner.template cast<real_t>();
80 13467070 : index_t currentIndex{0}, iStride{0}, jStride{0};
81 13467070 : if constexpr (dim == 2) {
82 135000 : currentIndex = detectorZeroIndex + min_corner[0];
83 135000 : } else {
84 135000 : currentIndex = detectorZeroIndex * (detector_max[1] + 1)
85 135000 : + min_corner[1] * (detector_max[0] + 1) + min_corner[0];
86 :
87 135000 : iStride = max_corner[0] - min_corner[0] + 1;
88 135000 : jStride = (detector_max[0] + 1) - iStride;
89 135000 : }
90 14758809 : for (index_t i = min_corner[0]; i <= max_corner[0]; i++) {
91 1291739 : if constexpr (dim == 2) {
92 : // traverse detector pixel in voxel footprint
93 1190219 : const Eigen::Matrix<real_t, 1, 1> distance((center[0] - i) / scaling);
94 1190219 : apply(detectorZeroIndex + i, distance);
95 1190219 : } else {
96 176312 : for (index_t j = min_corner[1]; j <= max_corner[1]; j++) {
97 74792 : const RealVector2D_t distanceVec = (center - current) / scaling;
98 74792 : apply(currentIndex, distanceVec);
99 74792 : currentIndex += 1;
100 74792 : current[0] += 1;
101 74792 : }
102 101520 : currentIndex += jStride;
103 101520 : current[0] -= iStride;
104 101520 : current[1] += 1;
105 101520 : }
106 1291739 : }
107 13467070 : }
108 :
109 : template <int dim, typename data_t, typename basis_function_t,
110 : size_t N = DEFAULT_PROJECTOR_LUT_SIZE>
111 : void forwardVoxel(const DataContainer<data_t>& x, DataContainer<data_t>& Ax,
112 : const Lut<data_t, N>& lut, basis_function_t& basis_function)
113 4099 : {
114 :
115 : // Just to be sure, zero out the result
116 4099 : Ax = 0;
117 :
118 4099 : const DetectorDescriptor& detectorDesc =
119 4099 : downcast<DetectorDescriptor>(Ax.getDataDescriptor());
120 :
121 4099 : auto& volume = x.getDataDescriptor();
122 4099 : const Eigen::Matrix<index_t, dim, 1>& volumeStrides =
123 4099 : volume.getProductOfCoefficientsPerDimension();
124 4099 : const Eigen::Matrix<index_t, dim - 1, 1>& detectorDims =
125 4099 : detectorDesc.getNumberOfCoefficientsPerDimension().head(dim - 1);
126 :
127 : // loop over geometries/poses in parallel
128 4099 : #pragma omp parallel for
129 4099 : for (index_t geomIndex = 0; geomIndex < detectorDesc.getNumberOfGeometryPoses();
130 0 : geomIndex++) {
131 0 : auto& geometry = detectorDesc.getGeometry()[asUnsigned(geomIndex)];
132 : // loop over voxels
133 13466395 : for (index_t domainIndex = 0; domainIndex < volume.getNumberOfCoefficients();
134 13466395 : ++domainIndex) {
135 13466395 : auto voxelWeight = x[domainIndex];
136 :
137 13466395 : visitDetector<dim>(
138 13466395 : domainIndex, geomIndex, geometry, lut.support(), detectorDims,
139 13466395 : volumeStrides,
140 13466395 : [&](const index_t index, Eigen::Matrix<real_t, dim - 1, 1> distance) {
141 1264489 : auto wght = basis_function(lut, distance);
142 1264489 : Ax[index] += voxelWeight * wght;
143 1264489 : });
144 13466395 : }
145 0 : }
146 4099 : }
147 :
148 : template <int dim, typename data_t, typename basis_function_t,
149 : size_t N = DEFAULT_PROJECTOR_LUT_SIZE>
150 : void backwardVoxel(const DataContainer<data_t>& y, DataContainer<data_t>& Aty,
151 : const Lut<data_t, N>& lut, basis_function_t basis_function)
152 34 : {
153 :
154 : // Just to be sure, zero out the result
155 34 : Aty = 0;
156 :
157 34 : const DetectorDescriptor& detectorDesc =
158 34 : downcast<DetectorDescriptor>(y.getDataDescriptor());
159 :
160 34 : auto& volume = Aty.getDataDescriptor();
161 34 : const Eigen::Matrix<index_t, dim, 1>& volumeStrides =
162 34 : volume.getProductOfCoefficientsPerDimension();
163 34 : const Eigen::Matrix<index_t, dim - 1, 1>& detectorDims =
164 34 : detectorDesc.getNumberOfCoefficientsPerDimension().head(dim - 1);
165 :
166 34 : #pragma omp parallel for
167 : // loop over voxels in parallel
168 34 : for (index_t domainIndex = 0; domainIndex < volume.getNumberOfCoefficients();
169 0 : ++domainIndex) {
170 : // loop over geometries/poses in parallel
171 411 : for (index_t geomIndex = 0; geomIndex < detectorDesc.getNumberOfGeometryPoses();
172 411 : geomIndex++) {
173 :
174 411 : auto& geometry = detectorDesc.getGeometry()[asUnsigned(geomIndex)];
175 :
176 411 : visitDetector<dim>(
177 411 : domainIndex, geomIndex, geometry, lut.support(), detectorDims,
178 411 : volumeStrides,
179 516 : [&](const index_t index, Eigen::Matrix<real_t, dim - 1, 1> distance) {
180 516 : auto wght = basis_function(lut, distance);
181 516 : Aty[domainIndex] += wght * y[index];
182 516 : });
183 411 : }
184 0 : }
185 34 : }
186 : }; // namespace voxel
187 : }; // namespace elsa
|