Line data    Source code

       1             : #include "VoxelProjector.h"
       2             : #include "Timer.h"
       3             : #include "Assertions.h"
       4             : 
       5             : namespace elsa
       6             : {
       7             :     template <typename data_t>
       8             :     BlobVoxelProjector<data_t>::BlobVoxelProjector(const VolumeDescriptor& domainDescriptor,
       9             :                                                    const DetectorDescriptor& rangeDescriptor,
      10             :                                                    data_t radius, data_t alpha, index_t order)
      11             :         : LinearOperator<data_t>(domainDescriptor, rangeDescriptor),
      12             :           blob(radius, alpha, order),
      13             :           _dim(domainDescriptor.getNumberOfDimensions())
      14        2319 :     {
      15             :         // sanity checks
      16        2319 :         if (_dim < 2 || _dim > 3) {
      17           0 :             throw InvalidArgumentError("BlobVoxelProjector: only supporting 2d/3d operations");
      18           0 :         }
      19             : 
      20        2319 :         if (_dim != rangeDescriptor.getNumberOfDimensions()) {
      21           0 :             throw InvalidArgumentError(
      22           0 :                 "BlobVoxelProjector: domain and range dimension need to match");
      23           0 :         }
      24             : 
      25        2319 :         if (rangeDescriptor.getNumberOfGeometryPoses() == 0) {
      26           0 :             throw InvalidArgumentError("BlobVoxelProjector: rangeDescriptor without any geometry");
      27           0 :         }
      28        2319 :     }
      29             : 
      30             :     template <typename data_t>
      31             :     void BlobVoxelProjector<data_t>::applyImpl(const elsa::DataContainer<data_t>& x,
      32             :                                                elsa::DataContainer<data_t>& Ax) const
      33        2012 :     {
      34        2012 :         if (_dim == 2)
      35        1472 :             voxel::forwardVoxel<2>(x, Ax, blob.get_lut(),
      36        1472 :                                    voxel::classic_weight_function<1, data_t>);
      37         540 :         else
      38         540 :             voxel::forwardVoxel<3>(x, Ax, blob.get_lut(),
      39         540 :                                    voxel::classic_weight_function<2, data_t>);
      40        2012 :     }
      41             : 
      42             :     template <typename data_t>
      43             :     void BlobVoxelProjector<data_t>::applyAdjointImpl(const elsa::DataContainer<data_t>& y,
      44             :                                                       elsa::DataContainer<data_t>& Aty) const
      45          17 :     {
      46          17 :         if (_dim == 2)
      47          17 :             voxel::backwardVoxel<2>(y, Aty, blob.get_lut(),
      48          17 :                                     voxel::classic_weight_function<1, data_t>);
      49           0 :         else
      50           0 :             voxel::backwardVoxel<3>(y, Aty, blob.get_lut(),
      51           0 :                                     voxel::classic_weight_function<2, data_t>);
      52          17 :     }
      53             : 
      54             :     template <typename data_t>
      55             :     bool BlobVoxelProjector<data_t>::isEqual(const elsa::LinearOperator<data_t>& other) const
      56           0 :     {
      57           0 :         return LinearOperator<data_t>::isEqual(other);
      58           0 :     }
      59             : 
      60             :     template <typename data_t>
      61             :     BlobVoxelProjector<data_t>* BlobVoxelProjector<data_t>::cloneImpl() const
      62           0 :     {
      63           0 :         return new BlobVoxelProjector(downcast<VolumeDescriptor>(*this->_domainDescriptor),
      64           0 :                                       downcast<DetectorDescriptor>(*this->_rangeDescriptor),
      65           0 :                                       blob.radius(), blob.alpha(), blob.order());
      66           0 :     }
      67             : 
      68             :     template <typename data_t>
      69             :     BSplineVoxelProjector<data_t>::BSplineVoxelProjector(const VolumeDescriptor& domainDescriptor,
      70             :                                                          const DetectorDescriptor& rangeDescriptor,
      71             :                                                          const index_t order)
      72             :         : LinearOperator<data_t>(domainDescriptor, rangeDescriptor),
      73             :           _dim(domainDescriptor.getNumberOfDimensions()),
      74             :           bspline(_dim, order)
      75        1612 :     {
      76             :         // sanity checks
      77        1612 :         if (_dim < 2 || _dim > 3) {
      78           0 :             throw InvalidArgumentError("BSplineVoxelProjector: only supporting 2d/3d operations");
      79           0 :         }
      80             : 
      81        1612 :         if (_dim != rangeDescriptor.getNumberOfDimensions()) {
      82           0 :             throw InvalidArgumentError(
      83           0 :                 "BSplineVoxelProjector: domain and range dimension need to match");
      84           0 :         }
      85             : 
      86        1612 :         if (rangeDescriptor.getNumberOfGeometryPoses() == 0) {
      87           0 :             throw InvalidArgumentError(
      88           0 :                 "BSplineVoxelProjector: rangeDescriptor without any geometry");
      89           0 :         }
      90        1612 :     }
      91             : 
      92             :     template <typename data_t>
      93             :     void BSplineVoxelProjector<data_t>::applyImpl(const elsa::DataContainer<data_t>& x,
      94             :                                                   elsa::DataContainer<data_t>& Ax) const
      95        1595 :     {
      96        1595 :         if (_dim == 2)
      97        1055 :             voxel::forwardVoxel<2>(x, Ax, bspline.get_lut(),
      98        1055 :                                    voxel::classic_weight_function<1, data_t>);
      99         540 :         else
     100         540 :             voxel::forwardVoxel<3>(x, Ax, bspline.get_lut(),
     101         540 :                                    voxel::classic_weight_function<2, data_t>);
     102        1595 :     }
     103             : 
     104             :     template <typename data_t>
     105             :     void BSplineVoxelProjector<data_t>::applyAdjointImpl(const elsa::DataContainer<data_t>& y,
     106             :                                                          elsa::DataContainer<data_t>& Aty) const
     107          17 :     {
     108          17 :         if (_dim == 2)
     109          17 :             voxel::backwardVoxel<2>(y, Aty, bspline.get_lut(),
     110          17 :                                     voxel::classic_weight_function<1, data_t>);
     111           0 :         else
     112           0 :             voxel::backwardVoxel<3>(y, Aty, bspline.get_lut(),
     113           0 :                                     voxel::classic_weight_function<2, data_t>);
     114          17 :     }
     115             : 
     116             :     template <typename data_t>
     117             :     bool BSplineVoxelProjector<data_t>::isEqual(const elsa::LinearOperator<data_t>& other) const
     118           0 :     {
     119           0 :         return LinearOperator<data_t>::isEqual(other);
     120           0 :     }
     121             : 
     122             :     template <typename data_t>
     123             :     BSplineVoxelProjector<data_t>* BSplineVoxelProjector<data_t>::cloneImpl() const
     124           0 :     {
     125           0 :         return new BSplineVoxelProjector(downcast<VolumeDescriptor>(*this->_domainDescriptor),
     126           0 :                                          downcast<DetectorDescriptor>(*this->_rangeDescriptor),
     127           0 :                                          bspline.order());
     128           0 :     }
     129             : 
     130             :     template class BlobVoxelProjector<float>;
     131             :     template class BlobVoxelProjector<double>;
     132             :     template class BSplineVoxelProjector<float>;
     133             :     template class BSplineVoxelProjector<double>;
     134             : }; // namespace elsa