Line data    Source code

       1             : #include "CurvedDetectorDescriptor.h"
       2             : #include "TypeCasts.hpp"
       3             : 
       4             : #include <vector>
       5             : 
       6             : namespace elsa
       7             : {
       8             :     CurvedDetectorDescriptor::CurvedDetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
       9             :                                                        const RealVector_t& spacingPerDim,
      10             :                                                        const std::vector<Geometry>& geometryList,
      11             :                                                        const geometry::Radian angle,
      12             :                                                        const real_t s2d)
      13             :         : DetectorDescriptor(numOfCoeffsPerDim, spacingPerDim, geometryList),
      14             :           _angle{angle},
      15             :           _radius{static_cast<real_t>(numOfCoeffsPerDim[0]) / static_cast<real_t>(angle)},
      16             :           _angleDelta{_angle / static_cast<real_t>(numOfCoeffsPerDim[0])},
      17             :           _s2d{s2d}
      18          11 :     {
      19          11 :         setup();
      20          11 :     }
      21             :     CurvedDetectorDescriptor::CurvedDetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
      22             :                                                        const std::vector<Geometry>& geometryList,
      23             :                                                        const geometry::Radian angle,
      24             :                                                        const real_t s2d)
      25             :         : DetectorDescriptor(numOfCoeffsPerDim, geometryList),
      26             :           _angle{angle},
      27             :           _radius{static_cast<real_t>(numOfCoeffsPerDim[0]) / static_cast<real_t>(angle)},
      28             :           _angleDelta{_angle / static_cast<real_t>(numOfCoeffsPerDim[0])},
      29             :           _s2d{s2d}
      30          20 :     {
      31          20 :         setup();
      32          20 :     }
      33             : 
      34             :     void CurvedDetectorDescriptor::setup()
      35          31 :     {
      36          31 :         index_t detectorLen{_numberOfCoefficientsPerDimension[0]};
      37          31 :         index_t detectorHeight{_numberOfDimensions == 3 ? _numberOfCoefficientsPerDimension[1] : 1};
      38             : 
      39          31 :         _detectorCenter = RealVector_t(2);
      40          31 :         _detectorCenter << static_cast<real_t>(detectorLen) / static_cast<real_t>(2.0),
      41          31 :             static_cast<real_t>(detectorHeight) / static_cast<real_t>(2.0);
      42             : 
      43          31 :         _centerOfCircle = RealVector_t(2);
      44          31 :         _centerOfCircle << _detectorCenter[0], _s2d - _radius;
      45             : 
      46          31 :         _planarCoords.reserve(static_cast<size_t>(detectorLen * detectorHeight));
      47             : 
      48          81 :         for (index_t coordH{0}; coordH < detectorHeight; coordH++) {
      49         309 :             for (index_t coordL{0}; coordL < detectorLen; coordL++) {
      50             :                 // shift by 0.5 to get the midpoint of the detector pixel
      51         259 :                 RealVector_t pixelCenter(_numberOfDimensions - 1);
      52         259 :                 if (_numberOfDimensions == 2) {
      53         136 :                     pixelCenter << static_cast<real_t>(coordL) + static_cast<real_t>(0.5);
      54         136 :                 } else {
      55         123 :                     pixelCenter << static_cast<real_t>(coordL) + static_cast<real_t>(0.5),
      56         123 :                         static_cast<real_t>(coordH) + static_cast<real_t>(0.5);
      57         123 :                 }
      58         259 :                 _planarCoords.push_back(mapCurvedCoordToPlanarCoord(pixelCenter));
      59         259 :             }
      60          50 :         }
      61          31 :     }
      62             : 
      63             :     inline RealVector_t
      64             :         CurvedDetectorDescriptor::mapCurvedCoordToPlanarCoord(const RealVector_t& coord) const
      65         259 :     {
      66             :         // angle between the principal ray and the current ray we want to compute
      67         259 :         real_t rayAngle{(coord[0] - _detectorCenter[0]) * _angleDelta};
      68             : 
      69             :         // deltas between the center and the point on the curved detector
      70         259 :         real_t delta_x{static_cast<real_t>(sin(rayAngle)) * _radius};
      71         259 :         real_t delta_depth{static_cast<real_t>(cos(rayAngle)) * _radius};
      72             : 
      73             :         // the coordinate of the point on the curved detector in detector space
      74         259 :         RealVector_t pointOnCurved(2);
      75         259 :         pointOnCurved << _centerOfCircle[0] + delta_x, _centerOfCircle[1] + delta_depth;
      76             :         // if the curvature is very small and the points are on top of each other, the factor is
      77             :         // zero
      78         259 :         real_t factor{pointOnCurved[1] == 0 ? 1 : _s2d / pointOnCurved[1]};
      79             : 
      80             :         // return new coordinate by scaling the vector onto the flat detector
      81         259 :         RealVector_t pointOnFlat(_numberOfDimensions - 1);
      82         259 :         if (_numberOfDimensions == 3) {
      83         123 :             pointOnFlat << (pointOnCurved[0] - _detectorCenter[0]) * factor + _detectorCenter[0],
      84         123 :                 (coord[1] - _detectorCenter[1]) * factor + _detectorCenter[1];
      85         136 :         } else {
      86         136 :             pointOnFlat << (pointOnCurved[0] - _detectorCenter[0]) * factor + _detectorCenter[0];
      87         136 :         }
      88             : 
      89         259 :         return pointOnFlat;
      90         259 :     }
      91             : 
      92             :     RealRay_t
      93             :         CurvedDetectorDescriptor::computeRayFromDetectorCoord(const RealVector_t& detectorCoord,
      94             :                                                               const index_t poseIndex) const
      95          48 :     {
      96          48 :         size_t pixelIndex{static_cast<size_t>(detectorCoord[0])};
      97             : 
      98             :         // instead of multiplying the detectorCoord with the projection matrix we retrieve the
      99             :         // actual coordinate in detector space by treating detectorCoord as an Index into
     100             :         // _planarCoords
     101          48 :         bool isPixelCenter{
     102          48 :             std::abs(detectorCoord[0] - static_cast<real_t>(pixelIndex) - static_cast<real_t>(0.5))
     103          48 :             <= static_cast<real_t>(0.0001)};
     104             : 
     105          48 :         if (_numberOfDimensions == 3) {
     106           9 :             size_t pixelIndexWidth{static_cast<size_t>(detectorCoord[1])};
     107           9 :             isPixelCenter = isPixelCenter
     108           9 :                             && std::abs(detectorCoord[1] - static_cast<real_t>(pixelIndexWidth)
     109           9 :                                         - static_cast<real_t>(0.5))
     110           9 :                                    <= static_cast<real_t>(0.0001);
     111           9 :             pixelIndex +=
     112           9 :                 pixelIndexWidth * static_cast<size_t>(_numberOfCoefficientsPerDimension[0]);
     113           9 :         }
     114             : 
     115          48 :         RealVector_t curvedDetectorCoord{
     116          48 :             isPixelCenter ? _planarCoords[pixelIndex] : mapCurvedCoordToPlanarCoord(detectorCoord)};
     117             : 
     118          48 :         return DetectorDescriptor::computeRayFromDetectorCoord(curvedDetectorCoord, poseIndex);
     119          48 :     }
     120             : 
     121             :     bool CurvedDetectorDescriptor::isEqual(const DataDescriptor& other) const
     122           1 :     {
     123           1 :         const CurvedDetectorDescriptor* otherPtr{
     124           1 :             downcast_safe<const CurvedDetectorDescriptor, const DataDescriptor>(&other)};
     125             : 
     126           1 :         return (otherPtr != nullptr && _angle == otherPtr->_angle && _radius == otherPtr->_radius
     127           1 :                 && _s2d == otherPtr->_s2d && DetectorDescriptor::isEqual(other));
     128           1 :     }
     129             : 
     130             :     CurvedDetectorDescriptor* CurvedDetectorDescriptor::cloneImpl() const
     131           9 :     {
     132           9 :         return new CurvedDetectorDescriptor(getNumberOfCoefficientsPerDimension(),
     133           9 :                                             getSpacingPerDimension(), _geometry, _angle, _s2d);
     134           9 :     }
     135             : 
     136             :     const std::vector<RealVector_t>& CurvedDetectorDescriptor::getPlanarCoords() const
     137           4 :     {
     138           4 :         return _planarCoords;
     139           4 :     }
     140             : 
     141             :     real_t CurvedDetectorDescriptor::getRadius() const
     142           4 :     {
     143           4 :         return _radius;
     144           4 :     }
     145             :     std::unique_ptr<DetectorDescriptor>
     146             :         CurvedDetectorDescriptor::cloneWithGeometry(std::vector<Geometry> geometries) const
     147           2 :     {
     148           2 :         auto coeff = getNumberOfCoefficientsPerDimension();
     149           2 :         auto dim = coeff.size();
     150           2 :         coeff[dim - 1] = geometries.size();
     151           2 :         return std::make_unique<CurvedDetectorDescriptor>(coeff, getSpacingPerDimension(),
     152           2 :                                                           geometries, _angle, _s2d);
     153           2 :     }
     154             : 
     155             : } // namespace elsa