Adding point in volume implementation and tests

include/xdg/mesh_manager_interface.h

@@ -58,6 +58,8 @@ class MeshManager {
 
   virtual void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) = 0;
 
+  MeshID next_volume(MeshID surface, MeshID current_volume) const;
+
   // Methods
   MeshID next_volume_id() const;
 

include/xdg/ray.h

@@ -80,6 +80,11 @@ struct RTCDRay: RTCRay {
 
 /*! Structure extending Embree's RayHit to include a double precision version of the primitive normal */
 struct RTCDHit : RTCHit {
+  RTCDHit() {
+    this->geomID = RTC_INVALID_GEOMETRY_ID;
+    this->primID = RTC_INVALID_GEOMETRY_ID;
+  }
+
   // data members
   const TriangleRef* tri_ref {nullptr}; //!< Pointer to the triangle reference for this hit
   Vec3da dNg; //!< Double precision version of the primitive normal

include/xdg/ray_tracing_interface.h

@@ -31,6 +31,14 @@ class RayTracer {
     }
   }
 
+
+  // Query Methods
+
+  bool point_in_volume(MeshID volume,
+                       const Position& point,
+                       const Direction* direction = nullptr,
+                       const std::vector<MeshID>* exclude_primitives = nullptr);
+
   void ray_fire(MeshID volume,
                 const Position& origin,
                 const Direction& direction,
@@ -49,15 +57,13 @@ class RayTracer {
 // Accessors
   int num_registered_volumes() const { return volume_map_.size(); }
 
-
 // Data members
 private:
   // Embree members
   RTCDevice device_;
-
   // Mesh-to-Scene maps
-  std::map<MeshID, RTCScene> volume_map_;
-  std::map<MeshID, RTCScene> surface_map_;
+  std::map<MeshID, RTCScene> volume_map_;  //<! Map from mesh volume to embree scene
+  std::map<MeshID, RTCScene> surface_map_; //<! Map from mesh surface to embree scnee
   RTCScene gloabal_scene_;
 
   // Internal Embree Mappings

src/mesh_manager_interface.cpp

@@ -2,6 +2,8 @@
 
 #include <algorithm>
 
+#include "xdg/error.h"
+
 namespace xdg {
 
 void
@@ -38,4 +40,16 @@ MeshID MeshManager::next_surface_id() const
     return *std::max_element(surfaces().begin(), surfaces().end()) + 1;
 }
 
-}
+MeshID MeshManager::next_volume(MeshID surface, MeshID current_volume) const
+{
+  auto parent_vols = this->get_parent_volumes(surface);
+
+  if (parent_vols.first == current_volume)
+    return parent_vols.second;
+  else if (parent_vols.second == current_volume)
+    return parent_vols.first;
+  else
+    fatal_error("Volume {} is not a parent of surface {}", current_volume, surface);
+}
+
+} // namespace xdg

src/ray_tracing_interface.cpp

@@ -107,26 +107,54 @@ RayTracer::register_volume(const std::shared_ptr<MeshManager> mesh_manager,
   rtcCommitScene(volume_scene);
 }
 
+bool RayTracer::point_in_volume(MeshID volume,
+                                const Position& point,
+                                const Direction* direction,
+                                const std::vector<MeshID>* exclude_primitives)
+{
+  RTCScene scene = volume_map_[volume];
+
+  RTCDRayHit rayhit;
+  rayhit.ray.set_org(point);
+  if (direction != nullptr) rayhit.ray.set_dir(*direction);
+  else rayhit.ray.set_dir({1. / std::sqrt(2.0), 1 / std::sqrt(2.0), 0.0});
+  rayhit.ray.rf_type = RayFireType::VOLUME;
+  rayhit.ray.orientation = HitOrientation::ANY;
+  rayhit.ray.set_tfar(INFTY);
+  rayhit.ray.set_tnear(0.0);
+  if (exclude_primitives != nullptr) rayhit.ray.exclude_primitives = exclude_primitives;
+
+  {
+    rtcIntersect1(scene, (RTCRayHit*)&rayhit);
+  }
+
+  // if the ray hit nothing, the point is outside of the volume
+  if (rayhit.hit.geomID == RTC_INVALID_GEOMETRY_ID) return false;
+
+  // use the hit triangle normal to determine if the intersection is
+  // exiting or entering
+  return rayhit.ray.ddir.dot(rayhit.hit.dNg) > 0.0;
+}
+
 void
 RayTracer::ray_fire(MeshID volume,
-                              const Position& origin,
-                              const Direction& direction,
-                              double& distance,
-                              const std::vector<MeshID>* exclude_primitves)
+                    const Position& origin,
+                    const Direction& direction,
+                    double& distance,
+                    const std::vector<MeshID>* exclude_primitves)
 {
   RTCScene scene = volume_map_[volume];
 
   RTCDRayHit rayhit;
+  // set ray data
   rayhit.ray.set_org(origin);
   rayhit.ray.set_dir(direction);
   rayhit.ray.set_tfar(INFTY);
   rayhit.ray.set_tnear(0.0);
   rayhit.ray.rf_type = RayFireType::VOLUME;
   rayhit.ray.orientation = HitOrientation::EXITING;
   rayhit.ray.mask = -1; // no mask
-
-  rayhit.hit.geomID = RTC_INVALID_GEOMETRY_ID;
-  rayhit.hit.primID = RTC_INVALID_GEOMETRY_ID;
+  if (exclude_primitves != nullptr) rayhit.ray.exclude_primitives = exclude_primitves;
 
   // fire the ray
   {
@@ -137,11 +165,10 @@ RayTracer::ray_fire(MeshID volume,
     rayhit.hit.Ng_z *= -1.0;
   }
 
-  if (rayhit.hit.geomID == RTC_INVALID_GEOMETRY_ID) {
+  if (rayhit.hit.geomID == RTC_INVALID_GEOMETRY_ID)
     distance = INFTY;
-  } else {
+  else
     distance = rayhit.ray.dtfar;
-  }
 }
 
 void RayTracer::closest(MeshID volume,

src/triangle_ref.cpp

@@ -11,12 +11,14 @@ namespace xdg
 
 bool orientation_cull(const Direction& ray_dir, const Direction& normal, HitOrientation orientation)
 {
+  if (orientation == HitOrientation::ANY) return false;
+
   double dot_prod = ray_dir.dot(normal);
-  if (orientation == HitOrientation::EXITING) {
-    if (dot_prod < 0.0) return true;
+  if (orientation == HitOrientation::EXITING && dot_prod < 0.0) {
+    return true;
   }
-  else if (orientation == HitOrientation::ENTERING) {
-    if (dot_prod >= 0.0) return true;
+  else if (orientation == HitOrientation::ENTERING && dot_prod >= 0.0) {
+    return true;
   }
   return false;
 }
@@ -47,7 +49,6 @@ void TriangleBoundsFunc(RTCBoundsFunctionArguments* args)
   args->bounds_o->upper_z = bounds.max_z + user_data->box_bump;
 }
 
-
 void TriangleIntersectionFunc(RTCIntersectFunctionNArguments* args) {
   const GeometryUserData* user_data = (const GeometryUserData*)args->geometryUserPtr;
   const MeshManager* mesh_manager = user_data->mesh_manager;
@@ -91,9 +92,7 @@ void TriangleIntersectionFunc(RTCIntersectFunctionNArguments* args) {
   rayhit->hit.primID = args->primID;
   rayhit->hit.tri_ref = &tri_ref;
 
-  rayhit->hit.dNg[0] = normal[0];
-  rayhit->hit.dNg[1] = normal[1];
-  rayhit->hit.dNg[2] = normal[2];
+  rayhit->hit.dNg = normal;
 }
 
 bool TriangleClosestFunc(RTCPointQueryFunctionArguments* args) {

tests/CMakeLists.txt

@@ -5,9 +5,10 @@ TEST_NAMES
 test_bbox
 test_bvh
 test_closest
+test_mesh_internal
 test_occluded
 test_ray_fire
-test_mesh_internal
+test_point_in_volume
 test_xdg_interface
 )
 

tests/test_point_in_volume.cpp

@@ -0,0 +1,61 @@
+// for testing
+#include <catch2/catch_test_macros.hpp>
+
+// xdg includes
+#include "xdg/mesh_manager_interface.h"
+#include "xdg/ray_tracing_interface.h"
+
+#include "mesh_mock.h"
+
+using namespace xdg;
+
+TEST_CASE("Test Point in Volume")
+{
+  std::shared_ptr<MeshManager> mm = std::make_shared<MeshMock>();
+  mm->init(); // this should do nothing, just good practice to call it
+  REQUIRE(mm->mesh_library() == MeshLibrary::INTERANAL);
+  std::shared_ptr<RayTracer> rti = std::make_shared<RayTracer>();
+
+  rti->register_all_volumes(mm);
+
+  MeshID volume = mm->volumes()[0];
+
+  Position point {0.0, 0.0, 0.0};
+  bool result = rti->point_in_volume(volume, point);
+  REQUIRE(result == true);
+
+  point = {0.0, 0.0, 1000.0};
+  result = rti->point_in_volume(volume, point);
+  REQUIRE(result == false);
+
+  // test a point just inside the positive x boundary
+  point = {5.0 - 1e-06, 0.0, 0.0};
+  result = rti->point_in_volume(volume, point);
+  REQUIRE(result == true);
+
+  // test a point on the positive x boundary
+  point = {5.0, 0.0, 0,0};
+  result = rti->point_in_volume(volume, point);
+  REQUIRE(result == false);
+
+  // test a point just outside the positive x boundary
+  // and provide a direction
+  point = {5.0, 0.0, 0.0};
+  Direction dir = {1.0, 0.0, 0.0};
+  result = rti->point_in_volume(volume, point, &dir);
+  REQUIRE(result == true);
+
+  // test a point just outside the positive x boundary
+  // and provide a direction
+  point = {5.1, 0.0, 0.0};
+  dir = {1.0, 0.0, 0.0};
+  result = rti->point_in_volume(volume, point, &dir);
+  REQUIRE(result == false);
+
+  // test a point just outside the positive x boundary,
+  // flip the direction
+  point = {5.1, 0.0, 0.0};
+  dir = {-1.0, 0.0, 0.0};
+  result = rti->point_in_volume(volume, point, &dir);
+  REQUIRE(result == false);
+}