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             :     std::tuple<real_t, real_t, real_t, real_t> intersect(const BoundingBox& aabb,
       8             :                                                          const RealRay_t& r)
       9      312880 :     {
      10      312880 :         real_t invDir = 1 / r.direction()(0);
      11             : 
      12      312880 :         real_t t1 = (aabb._min(0) - r.origin()(0)) * invDir;
      13      312880 :         real_t t2 = (aabb._max(0) - r.origin()(0)) * invDir;
      14             : 
      15      312880 :         real_t tmin = invDir >= 0 ? t1 : t2;
      16      312880 :         real_t tmax = invDir >= 0 ? t2 : t1;
      17      312880 :         const auto txmin = tmin;
      18      312880 :         const auto txmax = tmax;
      19             : 
      20      637236 :         for (int i = 1; i < aabb._min.rows(); ++i) {
      21      324356 :             invDir = 1 / r.direction()(i);
      22             : 
      23      324356 :             t1 = (aabb._min(i) - r.origin()(i)) * invDir;
      24      324356 :             t2 = (aabb._max(i) - r.origin()(i)) * invDir;
      25             : 
      26      324356 :             tmin = maxNum(tmin, invDir >= 0 ? t1 : t2);
      27      324356 :             tmax = minNum(tmax, invDir >= 0 ? t2 : t1);
      28      324356 :         }
      29             : 
      30      312880 :         return {tmin, tmax, txmin, txmax};
      31      312880 :     }
      32             :     std::optional<IntersectionResult> Intersection::withRay(const BoundingBox& aabb,
      33             :                                                             const RealRay_t& r)
      34      305746 :     {
      35      305746 :         const auto [tmin, tmax, txmin, txmax] = intersect(aabb, r);
      36             : 
      37      305746 :         if (tmax == 0.0 && tmin == 0.0)
      38           0 :             return std::nullopt;
      39      305746 :         if (tmax >= maxNum(tmin, 0.0f)) // hit
      40      283865 :             return std::make_optional<IntersectionResult>(tmin, tmax);
      41       21881 :         return std::nullopt;
      42       21881 :     }
      43             : 
      44             :     std::optional<IntersectionResult> Intersection::xPlanesWithRay(BoundingBox aabb,
      45             :                                                                    const RealRay_t& r)
      46        6852 :     {
      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        6852 :         aabb._min[0] += 0.5f;
      51        6852 :         aabb._max[0] -= 0.5f;
      52             : 
      53        6852 :         auto [tmin, tmax, txmin, txmax] = intersect(aabb, r);
      54             : 
      55        6852 :         auto dist_to_integer = [](auto real) {
      56         384 :             const auto aabbmin = static_cast<real_t>(static_cast<int>(std::round(real)));
      57         384 :             return std::abs(real - aabbmin);
      58         384 :         };
      59             : 
      60        6852 :         auto fractional = [](auto real) {
      61         383 :             real_t trash = 0;
      62         383 :             return std::modf(real, &trash);
      63         383 :         };
      64             : 
      65        6852 :         auto advance_to_next_voxel = [&](auto aabb, auto& tmin) -> real_t {
      66             :             // the x-axis coord for voxel centers
      67         304 :             const auto xvoxelcenter = dist_to_integer(aabb._min[0]);
      68             : 
      69         304 :             RealVector_t entry = r.pointAt(tmin);
      70             : 
      71             :             // Calculate the distance from entry.x to the x coord of the next voxel
      72         304 :             const auto tmp = std::abs(fractional(entry[0] - xvoxelcenter));
      73         304 :             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         304 :             return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
      77         304 :         };
      78             : 
      79        6852 :         auto retreat_to_last_voxel = [&](auto aabb, auto& tmax) -> real_t {
      80          80 :             const auto xvoxelcenter = dist_to_integer(aabb._min[0]);
      81             : 
      82          80 :             RealVector_t exit = r.pointAt(tmax);
      83             : 
      84             :             // calculate the distance from entry.x to the x coord of the previous voxel
      85          80 :             const auto tmp = std::abs(fractional(exit[0] - xvoxelcenter));
      86          80 :             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          80 :             return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
      90          80 :         };
      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        6852 :         if (txmin < tmin) {
      95         304 :             tmin += advance_to_next_voxel(aabb, tmin);
      96         304 :         }
      97             : 
      98             :         // This can also happen if we leave at the top of bottom
      99        6852 :         if (txmax > tmax) {
     100          80 :             tmax -= retreat_to_last_voxel(aabb, tmax);
     101          80 :         }
     102             : 
     103        6852 :         if (tmax == 0.0 && tmin == 0.0)
     104           0 :             return std::nullopt;
     105             : 
     106        6852 :         if (tmax - maxNum(tmin, 0.0f) > -0.000001) // hit
     107        6856 :             return std::make_optional<IntersectionResult>(tmin, tmax);
     108             : 
     109 >1844*10^16 :         return std::nullopt;
     110 >1844*10^16 :     }
     111             : } // namespace elsa