Line data    Source code

       1             : #include "DetectorDescriptor.h"
       2             : 
       3             : namespace elsa
       4             : {
       5             :     DetectorDescriptor::DetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
       6             :                                            const std::vector<Geometry>& geometryList)
       7             :         : DataDescriptor(numOfCoeffsPerDim), _geometry(geometryList)
       8        5779 :     {
       9             :         // TODO Clarify: What about empty geometryList? Do we want to support it, or throw an
      10             :         // exception?
      11        5779 :     }
      12             : 
      13             :     DetectorDescriptor::DetectorDescriptor(const IndexVector_t& numOfCoeffsPerDim,
      14             :                                            const RealVector_t& spacingPerDim,
      15             :                                            const std::vector<Geometry>& geometryList)
      16             :         : DataDescriptor(numOfCoeffsPerDim, spacingPerDim), _geometry(geometryList)
      17       50573 :     {
      18       50573 :     }
      19             : 
      20             :     RealRay_t DetectorDescriptor::computeRayFromDetectorCoord(const index_t detectorIndex) const
      21      226010 :     {
      22             : 
      23             :         // Return empty, if access out of bounds
      24      226010 :         assert(detectorIndex < getNumberOfCoefficients()
      25      226010 :                && "PlanarDetectorDescriptor::computeRayToDetector(index_t): Assumption "
      26      226010 :                   "detectorIndex smaller than number of coeffs, broken");
      27             : 
      28      226010 :         auto coord = getCoordinateFromIndex(detectorIndex);
      29      226010 :         return computeRayFromDetectorCoord(coord);
      30      226010 :     }
      31             : 
      32             :     RealRay_t DetectorDescriptor::computeRayFromDetectorCoord(const IndexVector_t coord) const
      33      229933 :     {
      34             :         // Assume all of the coordinates are inside of the volume
      35             :         // auto tmp = (coord.array() < getNumberOfCoefficientsPerDimension().array());
      36             :         // assert(tmp.all()
      37             :         // && "DetectorDescriptor::computeRayToDetector(IndexVector_t): Assumption coord "
      38             :         // "in bound wrong");
      39             : 
      40      229933 :         auto dim = getNumberOfDimensions();
      41             : 
      42             :         // Assume dimension of coord is equal to dimension of descriptor
      43      229933 :         assert(dim == coord.size());
      44             : 
      45             :         // Cast to real_t and shift to center of pixel
      46      229933 :         auto detectorCoord = coord.head(dim - 1).template cast<real_t>().array() + 0.5;
      47             : 
      48             :         // Last dimension is always the pose index
      49      229933 :         auto poseIndex = coord[dim - 1];
      50             : 
      51      229933 :         return computeRayFromDetectorCoord(detectorCoord, poseIndex);
      52      229933 :     }
      53             : 
      54             :     const std::vector<Geometry>& DetectorDescriptor::getGeometry() const
      55       40695 :     {
      56       40695 :         return _geometry;
      57       40695 :     }
      58             : 
      59             :     index_t DetectorDescriptor::getNumberOfGeometryPoses() const
      60       61946 :     {
      61       61946 :         return static_cast<index_t>(_geometry.size());
      62       61946 :     }
      63             : 
      64             :     Geometry DetectorDescriptor::getGeometryAt(const index_t index) const
      65       10857 :     {
      66       10857 :         return _geometry.at(asUnsigned(index));
      67       10857 :     }
      68             : 
      69             :     bool DetectorDescriptor::isEqual(const DataDescriptor& other) const
      70          70 :     {
      71          70 :         if (!DataDescriptor::isEqual(other))
      72           0 :             return false;
      73             : 
      74             :         // static cast as type checked in base comparison
      75          70 :         const auto& otherBlock = downcast<const DetectorDescriptor>(other);
      76             : 
      77          70 :         return _geometry == otherBlock._geometry;
      78          70 :     }
      79             : 
      80             :     RealRay_t DetectorDescriptor::computeRayFromDetectorCoord(const RealVector_t& detectorCoord,
      81             :                                                               const index_t poseIndex) const
      82      234405 :     {
      83             :         // Assert that for all dimension of detectorCoord is in bounds and poseIndex can
      84             :         // be index in the _geometry. If not the calculation will not be correct, but
      85             :         // as this is the hot path, I don't want exceptions and unpacking everything
      86             :         // We'll just have to ensure, that we don't mess up in our hot path! :-)
      87      234405 :         assert((detectorCoord.block(0, 0, getNumberOfDimensions() - 1, 0).array()
      88      234405 :                 < getNumberOfCoefficientsPerDimension()
      89      234405 :                       .block(0, 0, getNumberOfDimensions() - 1, 0)
      90      234405 :                       .template cast<real_t>()
      91      234405 :                       .array())
      92      234405 :                    .all()
      93      234405 :                && "PlanarDetectorDescriptor::computeRayToDetector: Assumption detectorCoord in "
      94      234405 :                   "bounds, wrong");
      95      234405 :         assert(asUnsigned(poseIndex) < _geometry.size()
      96      234405 :                && "PlanarDetectorDescriptor::computeRayToDetector: Assumption poseIndex smaller "
      97      234405 :                   "than number of poses, wrong");
      98             : 
      99      234405 :         auto dim = getNumberOfDimensions();
     100             : 
     101             :         // get the pose of trajectory
     102      234405 :         auto geometry = _geometry[asUnsigned(poseIndex)];
     103             : 
     104      234405 :         auto projInvMatrix = geometry.getInverseProjectionMatrix();
     105             : 
     106             :         // homogeneous coordinates [p;1], with p in detector space
     107      234405 :         RealVector_t homogeneousPixelCoord(dim);
     108      234405 :         homogeneousPixelCoord << detectorCoord, 1;
     109             : 
     110             :         // Camera center is always the ray origin
     111      234405 :         auto ro = geometry.getCameraCenter();
     112             : 
     113      234405 :         auto rd = (projInvMatrix * homogeneousPixelCoord) // Matrix-Vector multiplication
     114      234405 :                       .head(dim)                          // Transform to non-homogeneous
     115      234405 :                       .normalized();                      // normalize vector
     116             : 
     117      234405 :         return RealRay_t(ro, rd);
     118      234405 :     }
     119             : 
     120             :     std::pair<RealVector_t, real_t>
     121             :         DetectorDescriptor::projectAndScaleVoxelOnDetector(const RealVector_t& voxelCoord,
     122             :                                                            const index_t poseIndex) const
     123         504 :     {
     124         504 :         assert(asUnsigned(poseIndex) < _geometry.size()
     125         504 :                && "PlanarDetectorDescriptor::computeRayToDetector: Assumption poseIndex smaller "
     126         504 :                   "than number of poses, wrong");
     127             : 
     128             :         // get the pose of trajectory
     129         504 :         const auto& geometry = _geometry[asUnsigned(poseIndex)];
     130             : 
     131         504 :         auto homogeneousVoxelCoord = voxelCoord.homogeneous();
     132             : 
     133             :         // Project voxel center onto detector
     134         504 :         RealVector_t center =
     135         504 :             (geometry.getProjectionMatrix() * homogeneousVoxelCoord).hnormalized();
     136             : 
     137         504 :         auto scaling = geometry.getSourceDetectorDistance()
     138         504 :                        / (geometry.getExtrinsicMatrix() * homogeneousVoxelCoord).norm();
     139             : 
     140         504 :         return {std::move(center), scaling};
     141         504 :     }
     142             : } // namespace elsa