Line data    Source code

       1             : #include "EllipCylinderFree.h"
       2             : 
       3             : namespace elsa::phantoms
       4             : {
       5             : 
       6             :     template <typename data_t>
       7             :     EllipCylinderFree<data_t>::EllipCylinderFree(data_t amplit, elsa::phantoms::Vec3i center,
       8             :                                                  Vec2X<data_t> halfAxis, data_t length,
       9             :                                                  Vec3X<data_t> eulers)
      10             :         : _amplit{amplit}, _center{center}, _halfAxis{halfAxis}, _length{length}, _eulers{eulers}
      11         194 :     {
      12             : 
      13         194 :         aSqr = _halfAxis[INDEX_A] * _halfAxis[INDEX_A];
      14         194 :         bSqr = _halfAxis[INDEX_B] * _halfAxis[INDEX_B];
      15         194 :         aSqrbSqr = aSqr * bSqr;
      16             : 
      17         194 :         _centerX = _center.cast<data_t>();
      18             : 
      19         194 :         fillRotationMatrix(_eulers, rot);
      20         194 :         rot.transposeInPlace();
      21         194 :     };
      22             : 
      23             :     template <typename data_t>
      24             :     bool EllipCylinderFree<data_t>::isInEllipCylinderFree(const Vec3X<data_t>& idx,
      25             :                                                           index_t halfLength) const
      26      710853 :     {
      27             :         // check length on z axis
      28      710853 :         Vec3X<data_t> shifted{idx - _centerX};
      29      710853 :         Vec3X<data_t> rotated{getInvRotationMatrix() * shifted};
      30      710853 :         if (std::abs(rotated[INDEX_Z]) > data_t(halfLength)) {
      31      125170 :             return false;
      32      585683 :         } else {
      33      585683 :             return (double(rotated[INDEX_X] * rotated[INDEX_X]) * bSqr
      34      585683 :                     + double(rotated[INDEX_Y] * rotated[INDEX_Y]) * aSqr)
      35      585683 :                    <= aSqrbSqr;
      36      585683 :         }
      37      710853 :     };
      38             : 
      39             :     template <typename data_t>
      40             :     index_t getHalfDiagonal(EllipCylinderFree<data_t>& el, const index_t centerHalfLength)
      41         194 :     {
      42         194 :         auto maxAxis = std::max(el.getHalfAxis()[INDEX_A], el.getHalfAxis()[INDEX_A]);
      43             :         // Diagonal from center to edge of ellipse on the end of the cylinder with theorem of
      44             :         // pythagoras
      45         194 :         return std::lround(std::ceil(
      46         194 :             std::sqrt(double(centerHalfLength) * double(centerHalfLength) + maxAxis * maxAxis)));
      47         194 :     }
      48             : 
      49             :     template <typename data_t>
      50             :     std::array<data_t, 6>
      51             :         boundingBox(const index_t halfLength, const Vec3i& _center,
      52             :                     const IndexVector_t& ncpd /*numberOfCoefficientsPerDimension zero based*/)
      53         194 :     {
      54             : 
      55         194 :         index_t minX, minY, minZ, maxX, maxY, maxZ;
      56             : 
      57         194 :         minX = _center[INDEX_X] - halfLength > 0 ? _center[INDEX_X] - halfLength : 0;
      58         194 :         minY = _center[INDEX_Y] - halfLength > 0 ? _center[INDEX_Y] - halfLength : 0;
      59         194 :         minZ = _center[INDEX_Z] - halfLength > 0 ? _center[INDEX_Z] - halfLength : 0;
      60             : 
      61         194 :         maxX = _center[INDEX_X] + halfLength > ncpd[INDEX_X] ? ncpd[INDEX_X]
      62         194 :                                                              : _center[INDEX_X] + halfLength;
      63         194 :         maxY = _center[INDEX_Y] + halfLength > ncpd[INDEX_Y] ? ncpd[INDEX_Y]
      64         194 :                                                              : _center[INDEX_Y] + halfLength;
      65         194 :         maxZ = _center[INDEX_Z] + halfLength > ncpd[INDEX_Z] ? ncpd[INDEX_Z]
      66         194 :                                                              : _center[INDEX_Z] + halfLength;
      67             : 
      68         194 :         return {data_t(minX), data_t(minY), data_t(minZ), data_t(maxX), data_t(maxY), data_t(maxZ)};
      69         194 :     };
      70             : 
      71             :     template <Blending b, typename data_t>
      72             :     void rasterize(EllipCylinderFree<data_t>& el, VolumeDescriptor& dd, DataContainer<data_t>& dc)
      73         194 :     {
      74         194 :         auto strides = dd.getProductOfCoefficientsPerDimension();
      75             :         // global offests for a fast memory reuse in the for loops
      76         194 :         index_t xOffset = 0;
      77         194 :         index_t yOffset = 0;
      78         194 :         index_t zOffset = 0;
      79             : 
      80         194 :         Vec3i _center = el.getCenter();
      81         194 :         data_t _amplit = el.getAmplitude();
      82             : 
      83         194 :         index_t halfLength = std::lround(el.getLength() / 2);
      84             : 
      85         194 :         auto [minX, minY, minZ, maxX, maxY, maxZ] = boundingBox<data_t>(
      86         194 :             getHalfDiagonal(el, halfLength), _center, dd.getNumberOfCoefficientsPerDimension());
      87             : 
      88         194 :         Vec3X<data_t> idx(3);
      89         194 :         idx << minX, minY, minZ;
      90             : 
      91        3134 :         for (index_t z = index_t(minZ); z <= index_t(maxZ); z++, idx[INDEX_Z]++) {
      92        2940 :             zOffset = z * strides[INDEX_Z];
      93             : 
      94       48174 :             for (index_t y = index_t(minY); y <= index_t(maxY); y++, idx[INDEX_Y]++) {
      95       45234 :                 yOffset = y * strides[INDEX_Y];
      96             : 
      97      756087 :                 for (index_t x = index_t(minX); x <= index_t(maxX); x++, idx[INDEX_X]++) {
      98      710853 :                     xOffset = x * strides[INDEX_X];
      99      710853 :                     if (el.isInEllipCylinderFree(idx, halfLength)) {
     100        2682 :                         blend<b>(dc, zOffset + yOffset + xOffset, _amplit);
     101        2682 :                     }
     102      710853 :                 }
     103             : 
     104       45234 :                 idx[INDEX_X] = minX;
     105       45234 :             }
     106        2940 :             idx[INDEX_Y] = minY;
     107        2940 :         }
     108         194 :     };
     109             : 
     110             :     // ------------------------------------------
     111             :     // explicit template instantiation
     112             :     template class EllipCylinderFree<float>;
     113             :     template class EllipCylinderFree<double>;
     114             : 
     115             :     template void rasterize<Blending::ADDITION, float>(EllipCylinderFree<float>& el,
     116             :                                                        VolumeDescriptor& dd,
     117             :                                                        DataContainer<float>& dc);
     118             :     template void rasterize<Blending::ADDITION, double>(EllipCylinderFree<double>& el,
     119             :                                                         VolumeDescriptor& dd,
     120             :                                                         DataContainer<double>& dc);
     121             : 
     122             :     template void rasterize<Blending::OVERWRITE, float>(EllipCylinderFree<float>& el,
     123             :                                                         VolumeDescriptor& dd,
     124             :                                                         DataContainer<float>& dc);
     125             :     template void rasterize<Blending::OVERWRITE, double>(EllipCylinderFree<double>& el,
     126             :                                                          VolumeDescriptor& dd,
     127             :                                                          DataContainer<double>& dc);
     128             : 
     129             : } // namespace elsa::phantoms