LCOV - code coverage report
Current view: top level - projectors - Intersection.cpp (source / functions) Hit Total Coverage
Test: test_coverage.info.cleaned Lines: 52 55 94.5 %
Date: 2022-05-27 02:48:28 Functions: 7 7 100.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      329710 :     std::tuple<real_t, real_t, real_t, real_t> intersect(const BoundingBox& aabb,
       8             :                                                          const RealRay_t& r)
       9             :     {
      10      329710 :         real_t invDir = 1 / r.direction()(0);
      11             : 
      12      329710 :         real_t t1 = (aabb._min(0) - r.origin()(0)) * invDir;
      13      329710 :         real_t t2 = (aabb._max(0) - r.origin()(0)) * invDir;
      14             : 
      15      329710 :         real_t tmin = invDir >= 0 ? t1 : t2;
      16      329710 :         real_t tmax = invDir >= 0 ? t2 : t1;
      17      329710 :         const auto txmin = tmin;
      18      329710 :         const auto txmax = tmax;
      19             : 
      20      683321 :         for (int i = 1; i < aabb._min.rows(); ++i) {
      21      353611 :             invDir = 1 / r.direction()(i);
      22             : 
      23      353611 :             t1 = (aabb._min(i) - r.origin()(i)) * invDir;
      24      353611 :             t2 = (aabb._max(i) - r.origin()(i)) * invDir;
      25             : 
      26      353611 :             tmin = maxNum(tmin, invDir >= 0 ? t1 : t2);
      27      353611 :             tmax = minNum(tmax, invDir >= 0 ? t2 : t1);
      28             :         }
      29             : 
      30      659420 :         return {tmin, tmax, txmin, txmax};
      31             :     }
      32      309350 :     std::optional<IntersectionResult> Intersection::withRay(const BoundingBox& aabb,
      33             :                                                             const RealRay_t& r)
      34             :     {
      35      309350 :         const auto [tmin, tmax, txmin, txmax] = intersect(aabb, r);
      36             : 
      37      309350 :         if (tmax == 0.0 && tmin == 0.0)
      38           0 :             return std::nullopt;
      39      309350 :         if (tmax >= maxNum(tmin, 0.0f)) // hit
      40      285074 :             return std::make_optional<IntersectionResult>(tmin, tmax);
      41       24276 :         return std::nullopt;
      42             :     }
      43             : 
      44       20360 :     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       20360 :         aabb._min[0] += 0.5f;
      51       20360 :         aabb._max[0] -= 0.5f;
      52             : 
      53       20360 :         auto [tmin, tmax, txmin, txmax] = intersect(aabb, r);
      54             : 
      55         212 :         auto dist_to_integer = [](auto real) {
      56         212 :             const auto aabbmin = static_cast<real_t>(static_cast<int>(std::round(real)));
      57         212 :             return std::abs(real - aabbmin);
      58             :         };
      59             : 
      60         212 :         auto fractional = [](auto real) {
      61         212 :             real_t trash = 0;
      62         212 :             return std::modf(real, &trash);
      63             :         };
      64             : 
      65         204 :         auto advance_to_next_voxel = [&](auto aabb, auto& tmin) -> real_t {
      66             :             // the x-axis coord for voxel centers
      67         204 :             const auto xvoxelcenter = dist_to_integer(aabb._min[0]);
      68             : 
      69         204 :             RealVector_t entry = r.pointAt(tmin);
      70             : 
      71             :             // Calculate the distance from entry.x to the x coord of the next voxel
      72         204 :             const auto tmp = std::abs(fractional(entry[0] - xvoxelcenter));
      73         204 :             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         408 :             return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
      77       20360 :         };
      78             : 
      79           8 :         auto retreat_to_last_voxel = [&](auto aabb, auto& tmax) -> real_t {
      80           8 :             const auto xvoxelcenter = dist_to_integer(aabb._min[0]);
      81             : 
      82           8 :             RealVector_t exit = r.pointAt(tmax);
      83             : 
      84             :             // calculate the distance from entry.x to the x coord of the previous voxel
      85           8 :             const auto tmp = std::abs(fractional(exit[0] - xvoxelcenter));
      86           8 :             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          16 :             return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
      90       20360 :         };
      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       20360 :         if (txmin < tmin) {
      95         204 :             tmin += advance_to_next_voxel(aabb, tmin);
      96             :         }
      97             : 
      98             :         // This can also happen if we leave at the top of bottom
      99       20360 :         if (txmax > tmax) {
     100           8 :             tmax -= retreat_to_last_voxel(aabb, tmax);
     101             :         }
     102             : 
     103       20360 :         if (tmax == 0.0 && tmin == 0.0)
     104           0 :             return std::nullopt;
     105             : 
     106       20360 :         if (tmax - maxNum(tmin, 0.0f) > -0.000001) // hit
     107       20360 :             return std::make_optional<IntersectionResult>(tmin, tmax);
     108             : 
     109           0 :         return std::nullopt;
     110             :     }
     111             : } // namespace elsa

Generated by: LCOV version 1.15