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

          Line data    Source code
       1             : #include "Intersection.h"
       2             : #include "spdlog/fmt/fmt.h"
       3             : #include "spdlog/fmt/ostr.h"
       4             : 
       5             : namespace elsa
       6             : {
       7           0 :     std::tuple<real_t, real_t, real_t, real_t> intersect(const BoundingBox& aabb,
       8             :                                                          const RealRay_t& r)
       9             :     {
      10           0 :         real_t invDir = 1 / r.direction()(0);
      11             : 
      12           0 :         real_t t1 = (aabb._min(0) - r.origin()(0)) * invDir;
      13           0 :         real_t t2 = (aabb._max(0) - r.origin()(0)) * invDir;
      14             : 
      15           0 :         real_t tmin = invDir >= 0 ? t1 : t2;
      16           0 :         real_t tmax = invDir >= 0 ? t2 : t1;
      17           0 :         const auto txmin = tmin;
      18           0 :         const auto txmax = tmax;
      19             : 
      20           0 :         for (int i = 1; i < aabb._min.rows(); ++i) {
      21           0 :             invDir = 1 / r.direction()(i);
      22             : 
      23           0 :             t1 = (aabb._min(i) - r.origin()(i)) * invDir;
      24           0 :             t2 = (aabb._max(i) - r.origin()(i)) * invDir;
      25             : 
      26           0 :             tmin = maxNum(tmin, invDir >= 0 ? t1 : t2);
      27           0 :             tmax = minNum(tmax, invDir >= 0 ? t2 : t1);
      28             :         }
      29             : 
      30           0 :         return {tmin, tmax, txmin, txmax};
      31             :     }
      32           0 :     std::optional<IntersectionResult> Intersection::withRay(const BoundingBox& aabb,
      33             :                                                             const RealRay_t& r)
      34             :     {
      35           0 :         const auto [tmin, tmax, txmin, txmax] = intersect(aabb, r);
      36             : 
      37           0 :         if (tmax == 0.0 && tmin == 0.0)
      38           0 :             return std::nullopt;
      39           0 :         if (tmax >= maxNum(tmin, 0.0f)) // hit
      40           0 :             return std::make_optional<IntersectionResult>(tmin, tmax);
      41           0 :         return std::nullopt;
      42             :     }
      43             : 
      44           0 :     std::optional<IntersectionResult> Intersection::xPlanesWithRay(BoundingBox aabb,
      45             :                                                                    const RealRay_t& r)
      46             :     {
      47             :         // Slightly increase the size of the bounding box, such that center of voxels are actually
      48             :         // inside the bounding box, parallel rays which directly go through the center of the
      49             :         // voxel, will be recognized with this cheapo hack
      50           0 :         aabb._min[0] += 0.5f;
      51           0 :         aabb._max[0] -= 0.5f;
      52             : 
      53           0 :         auto [tmin, tmax, txmin, txmax] = intersect(aabb, r);
      54             : 
      55           0 :         auto dist_to_integer = [](auto real) {
      56           0 :             const auto aabbmin = static_cast<real_t>(static_cast<int>(std::round(real)));
      57           0 :             return std::abs(real - aabbmin);
      58             :         };
      59             : 
      60           0 :         auto fractional = [](auto real) {
      61           0 :             real_t trash = 0;
      62           0 :             return std::modf(real, &trash);
      63             :         };
      64             : 
      65           0 :         auto advance_to_next_voxel = [&](auto aabb, auto& tmin) -> real_t {
      66             :             // the x-axis coord for voxel centers
      67           0 :             const auto xvoxelcenter = dist_to_integer(aabb._min[0]);
      68             : 
      69           0 :             RealVector_t entry = r.pointAt(tmin);
      70             : 
      71             :             // Calculate the distance from entry.x to the x coord of the next voxel
      72           0 :             const auto tmp = std::abs(fractional(entry[0] - xvoxelcenter));
      73           0 :             const auto frac = entry[0] < xvoxelcenter ? tmp : 1 - tmp;
      74             : 
      75             :             // Calculate distance in t, but don't advance if it's too small
      76           0 :             return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
      77           0 :         };
      78             : 
      79           0 :         auto retreat_to_last_voxel = [&](auto aabb, auto& tmax) -> real_t {
      80           0 :             const auto xvoxelcenter = dist_to_integer(aabb._min[0]);
      81             : 
      82           0 :             RealVector_t exit = r.pointAt(tmax);
      83             : 
      84             :             // calculate the distance from entry.x to the x coord of the previous voxel
      85           0 :             const auto tmp = std::abs(fractional(exit[0] - xvoxelcenter));
      86           0 :             const auto frac = exit[0] < xvoxelcenter ? 1 - tmp : tmp;
      87             : 
      88             :             // Calculate distance in t, but don't advance if it's too small
      89           0 :             return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
      90           0 :         };
      91             : 
      92             :         // Sometimes we hit the y axis first, not the x axis, but we don't want
      93             :         // that, we always want the intersection with the x-plane
      94           0 :         if (txmin < tmin) {
      95           0 :             tmin += advance_to_next_voxel(aabb, tmin);
      96             :         }
      97             : 
      98             :         // This can also happen if we leave at the top of bottom
      99           0 :         if (txmax > tmax) {
     100           0 :             tmax -= retreat_to_last_voxel(aabb, tmax);
     101             :         }
     102             : 
     103           0 :         if (tmax == 0.0 && tmin == 0.0)
     104           0 :             return std::nullopt;
     105             : 
     106           0 :         if (tmax - maxNum(tmin, 0.0f) > -0.000001) // hit
     107           0 :             return std::make_optional<IntersectionResult>(tmin, tmax);
     108             : 
     109           0 :         return std::nullopt;
     110             :     }
     111             : } // namespace elsa

Generated by: LCOV version 1.14