LCOV - code coverage report
Current view: top level - projectors - SliceTraversal.cpp (source / functions) Hit Total Coverage
Test: test_coverage.info.cleaned Lines: 0 124 0.0 %
Date: 2022-08-04 03:43:28 Functions: 0 32 0.0 %

          Line data    Source code
       1             : #include "SliceTraversal.h"
       2             : #include "Error.h"
       3             : #include "elsaDefines.h"
       4             : #include "spdlog/fmt/fmt.h"
       5             : #include "spdlog/fmt/ostr.h"
       6             : 
       7             : namespace elsa
       8             : {
       9           0 :     RealMatrix_t TransformToTraversal::create2DRotation(const real_t leadingCoeff,
      10             :                                                         const index_t leadingAxisIndex)
      11             :     {
      12           0 :         auto createRotationMatrix = [](real_t radian) {
      13           0 :             real_t c = std::cos(radian);
      14           0 :             real_t s = std::sin(radian);
      15           0 :             return RealMatrix_t({{c, -s}, {s, c}});
      16             :         };
      17             : 
      18             :         // TODO: Does a closed form solution exist?
      19             :         // Use conversion to polar coordinates
      20           0 :         if (leadingAxisIndex == 0 && leadingCoeff >= 0) {
      21             :             // Already identity, so do nothing
      22           0 :             return createRotationMatrix(0);
      23           0 :         } else if (leadingAxisIndex == 0 && leadingCoeff < 0) {
      24             :             // create a 2D 180° rotation matrix
      25           0 :             return createRotationMatrix(pi_t);
      26           0 :         } else if (leadingAxisIndex == 1 && leadingCoeff >= 0) {
      27             :             // create a 2D 270° rotation matrix counter clockwise
      28           0 :             return createRotationMatrix(3 * pi_t / 2);
      29           0 :         } else if (leadingAxisIndex == 1 && leadingCoeff <= 0) {
      30             :             // create a 2D 90° rotation matrix counter clockwise
      31           0 :             return createRotationMatrix(pi_t / 2);
      32             :         }
      33           0 :         return createRotationMatrix(0);
      34             :     }
      35             : 
      36           0 :     RealMatrix_t TransformToTraversal::create3DRotation(const real_t leadingCoeff,
      37             :                                                         const index_t leadingAxisIndex)
      38             :     {
      39           0 :         auto identity = []() { return RealMatrix_t({{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}); };
      40             : 
      41             :         auto rotationX = [](real_t radian) {
      42             :             real_t c = std::cos(radian);
      43             :             real_t s = std::sin(radian);
      44             :             return RealMatrix_t({{1, 0, 0}, {0, c, -s}, {0, s, c}});
      45             :         };
      46             : 
      47           0 :         auto rotationY = [](real_t radian) {
      48           0 :             real_t c = std::cos(radian);
      49           0 :             real_t s = std::sin(radian);
      50           0 :             return RealMatrix_t({{c, 0, s}, {0, 1, 0}, {-s, 0, c}});
      51             :         };
      52             : 
      53           0 :         auto rotationZ = [](real_t radian) {
      54           0 :             real_t c = std::cos(radian);
      55           0 :             real_t s = std::sin(radian);
      56           0 :             return RealMatrix_t({{c, -s, 0}, {s, c, 0}, {0, 0, 1}});
      57             :         };
      58             : 
      59             :         // TODO: Does a closed form solution exist?
      60             :         // Use conversion to polar coordinates
      61           0 :         if (leadingAxisIndex == 0 && leadingCoeff >= 0) {
      62             :             // Already identity, so do nothing
      63           0 :             return identity();
      64           0 :         } else if (leadingAxisIndex == 0 && leadingCoeff < 0) {
      65           0 :             return rotationY(pi_t);
      66           0 :         } else if (leadingAxisIndex == 1 && leadingCoeff >= 0) {
      67           0 :             return rotationZ(-0.5 * pi_t);
      68           0 :         } else if (leadingAxisIndex == 1 && leadingCoeff <= 0) {
      69           0 :             return rotationZ(0.5 * pi_t);
      70           0 :         } else if (leadingAxisIndex == 2 && leadingCoeff >= 0) {
      71           0 :             return rotationY(0.5 * pi_t);
      72           0 :         } else if (leadingAxisIndex == 2 && leadingCoeff <= 0) {
      73           0 :             return rotationY(-0.5 * pi_t);
      74             :         }
      75             : 
      76           0 :         return identity();
      77             :     }
      78             : 
      79           0 :     RealMatrix_t TransformToTraversal::createRotation(RealRay_t ray)
      80             :     {
      81             :         // Get the leading axis, by absolute value
      82           0 :         index_t leadingAxisIndex = 0;
      83           0 :         ray.direction().array().abs().maxCoeff(&leadingAxisIndex);
      84             : 
      85             :         // Get the signed leading coefficient
      86           0 :         const auto leadingCoeff = ray.direction()(leadingAxisIndex);
      87             : 
      88           0 :         if (ray.dim() == 2) {
      89           0 :             return create2DRotation(leadingCoeff, leadingAxisIndex);
      90           0 :         } else if (ray.dim() == 3) {
      91           0 :             return create3DRotation(leadingCoeff, leadingAxisIndex);
      92             :         }
      93           0 :         throw Error("Can not create a {}-dimensional transformation for the slice traversal",
      94           0 :                     ray.dim());
      95             :     }
      96             : 
      97           0 :     RealMatrix_t TransformToTraversal::createTransformation(const RealRay_t& ray,
      98             :                                                             const RealVector_t& centerOfRotation)
      99             :     {
     100           0 :         const auto dim = ray.dim();
     101             : 
     102             :         // Setup translation
     103           0 :         RealMatrix_t translation(dim + 1, dim + 1);
     104           0 :         translation.setIdentity();
     105           0 :         translation.block(0, dim, dim, 1) = -centerOfRotation;
     106             : 
     107             :         // Setup rotation
     108           0 :         RealMatrix_t rotation(dim + 1, dim + 1);
     109           0 :         rotation.setIdentity();
     110           0 :         rotation.block(0, 0, dim, dim) = createRotation(ray);
     111             : 
     112           0 :         return rotation * translation;
     113           0 :     }
     114             : 
     115           0 :     TransformToTraversal::TransformToTraversal(const RealRay_t& ray,
     116           0 :                                                const RealVector_t& centerOfRotation)
     117           0 :         : transformation_(createTransformation(ray, centerOfRotation)),
     118           0 :           translation_(-centerOfRotation)
     119             :     {
     120           0 :     }
     121             : 
     122           0 :     RealRay_t TransformToTraversal::toTraversalCoordinates(RealRay_t ray) const
     123             :     {
     124           0 :         ray.origin() = *this * Point(ray.origin());
     125           0 :         ray.direction() = *this * Vec(ray.direction());
     126           0 :         return ray;
     127             :     }
     128             : 
     129           0 :     BoundingBox TransformToTraversal::toTraversalCoordinates(BoundingBox aabb) const
     130             :     {
     131             :         // Only translate, as the rotations are always 90, 180, 270 degrees,
     132             :         // TODO: this might be wrong, idk, what about non square things
     133           0 :         aabb._min = *this * Point(aabb._min);
     134           0 :         aabb._max = *this * Point(aabb._max);
     135           0 :         aabb.recomputeBounds();
     136           0 :         return aabb;
     137             :     }
     138             : 
     139           0 :     RealMatrix_t TransformToTraversal::transformation() const { return transformation_; }
     140             : 
     141           0 :     RealMatrix_t TransformToTraversal::invTransformation() const
     142             :     {
     143           0 :         return transformation().inverse();
     144             :     }
     145             : 
     146           0 :     RealMatrix_t TransformToTraversal::rotation() const
     147             :     {
     148           0 :         const auto dim = transformation_.rows();
     149           0 :         return transformation_.block(0, 0, dim - 1, dim - 1);
     150             :     }
     151             : 
     152           0 :     RealMatrix_t TransformToTraversal::invRotation() const { return rotation().transpose(); }
     153             : 
     154           0 :     RealVector_t TransformToTraversal::translation() const { return translation_; }
     155             : 
     156           0 :     RealMatrix_t TransformToTraversal::linear() const { return rotation(); }
     157             : 
     158           0 :     RealVector_t TransformToTraversal::operator*(const Point<real_t>& point) const
     159             :     {
     160           0 :         return (transformation() * point.point_.homogeneous()).hnormalized();
     161             :     }
     162             : 
     163           0 :     RealVector_t TransformToTraversal::operator*(const Vec<real_t>& vec) const
     164             :     {
     165           0 :         return linear() * vec.vec_;
     166             :     }
     167             : 
     168           0 :     SliceTraversal::SliceTraversal(BoundingBox aabb, RealRay_t ray)
     169           0 :         : transformation_(ray, aabb._max.array() / 2), ray_(ray)
     170             :     {
     171             :         // Transform ray and aabb to traversal coordinate space
     172           0 :         ray = transformation_.toTraversalCoordinates(ray);
     173           0 :         aabb = transformation_.toTraversalCoordinates(aabb);
     174             : 
     175             :         // TODO: We only want to shift in x direction by 0.5, not in y
     176           0 :         auto hit = Intersection::xPlanesWithRay(aabb, ray);
     177             : 
     178             :         // Default init is sufficient if we don't hit,  TODO: test that!
     179           0 :         if (hit) {
     180           0 :             const auto firstVoxel = ray.pointAt(hit->_tmin);
     181           0 :             const auto lastVoxel = ray.pointAt(hit->_tmax);
     182             : 
     183           0 :             endIndex_ = std::max<index_t>(
     184           0 :                 1, static_cast<index_t>(std::ceil(lastVoxel[0] - firstVoxel[0] + 0.5f)));
     185             : 
     186           0 :             tDelta_ = 1 / ray.direction()[0];
     187             : 
     188           0 :             t_ = hit->_tmin;
     189           0 :         }
     190           0 :     }
     191             : 
     192           0 :     index_t SliceTraversal::leadingDirection() const
     193             :     {
     194           0 :         index_t idx = 0;
     195           0 :         ray_.direction().array().cwiseAbs().maxCoeff(&idx);
     196           0 :         return idx;
     197             :     }
     198             : 
     199           0 :     real_t SliceTraversal::t() const { return t_; }
     200             : 
     201           0 :     real_t SliceTraversal::tDelta() const { return tDelta_; }
     202             : 
     203           0 :     SliceTraversal::Iter SliceTraversal::begin() const
     204             :     {
     205           0 :         return {startIndex_, ray_.pointAt(t_), ray_.direction() * tDelta_};
     206             :     }
     207             : 
     208           0 :     SliceTraversal::Iter SliceTraversal::end() const { return {endIndex_}; }
     209             : 
     210           0 :     index_t SliceTraversal::startIndex() const { return startIndex_; }
     211             : 
     212           0 :     index_t SliceTraversal::endIndex() const { return endIndex_; }
     213             : 
     214           0 :     SliceTraversal::Iter::value_type SliceTraversal::Iter::operator*() const
     215             :     {
     216           0 :         return cur_.template cast<index_t>();
     217             :     }
     218             : 
     219           0 :     SliceTraversal::Iter& SliceTraversal::Iter::operator++()
     220             :     {
     221           0 :         ++pos_;
     222           0 :         cur_ += dir_;
     223           0 :         return *this;
     224             :     }
     225             : 
     226           0 :     SliceTraversal::Iter SliceTraversal::Iter::operator++(int)
     227             :     {
     228           0 :         auto copy = *this;
     229           0 :         ++(*this);
     230           0 :         return copy;
     231           0 :     }
     232             : 
     233           0 :     bool operator==(const SliceTraversal::Iter& lhs, const SliceTraversal::Iter& rhs)
     234             :     {
     235           0 :         return lhs.pos_ == rhs.pos_;
     236             :     }
     237             : 
     238           0 :     bool operator!=(const SliceTraversal::Iter& lhs, const SliceTraversal::Iter& rhs)
     239             :     {
     240           0 :         return !(lhs == rhs);
     241             :     }
     242             : 
     243             : } // namespace elsa

Generated by: LCOV version 1.14