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      652107 :         {
      23      652107 :             return lut(distance.norm());
      24      652107 :         }
      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    14758813 :             for (index_t i = min_corner[0]; i <= max_corner[0]; i++) {
      91     1291743 :                 if constexpr (dim == 2) {
      92             :                     // traverse detector pixel in voxel footprint
      93     1190223 :                     const Eigen::Matrix<real_t, 1, 1> distance((center[0] - i) / scaling);
      94     1190223 :                     apply(detectorZeroIndex + i, distance);
      95     1190223 :                 } 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     1291743 :             }
     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         448 :                 for (index_t geomIndex = 0; geomIndex < detectorDesc.getNumberOfGeometryPoses();
     172         448 :                      geomIndex++) {
     173             : 
     174         448 :                     auto& geometry = detectorDesc.getGeometry()[asUnsigned(geomIndex)];
     175             : 
     176         448 :                     visitDetector<dim>(
     177         448 :                         domainIndex, geomIndex, geometry, lut.support(), detectorDims,
     178         448 :                         volumeStrides,
     179         528 :                         [&](const index_t index, Eigen::Matrix<real_t, dim - 1, 1> distance) {
     180         528 :                             auto wght = basis_function(lut, distance);
     181         528 :                             Aty[domainIndex] += wght * y[index];
     182         528 :                         });
     183         448 :                 }
     184           0 :             }
     185          34 :         }
     186             :     }; // namespace voxel
     187             : };     // namespace elsa