intersect function additions to scene

src/core/scene.cpp

@@ -39,6 +39,7 @@
 #include "parallel.h"
 #include "progressreporter.h"
 #include "stats.h"
+#include <boost/numeric/odeint.hpp>
 
 namespace pbrt {
 
@@ -74,4 +75,71 @@ bool Scene::IntersectTr(Ray ray, Sampler &sampler, SurfaceInteraction *isect,
     }
 }
 
+int Scene::forward(state_type &prev, state_type &state, SurfaceInteraction *isect) const{
+   Point3f vc1 = Point3f(prev[0], prev[1], prev[2]);
+   Point3f vc2 = Point3f(state[0],state[1],state[2]);
+
+   return aggregate->IntersectCu(vc1, vc2, isect);
+   return -1;
+}
+
+bool Scene::IntersectCu(const Ray &ray, SurfaceInteraction *isect) const {
+    ++nIntersectionTests;
+    DCHECK_NE(ray.d, Vector3f(0,0,0));
+
+    state_type x(6);
+    x[0] = ray.o.x;
+    x[1] = ray.o.y;
+    x[2] = ray.o.z;
+
+    x[3] = ray.d.x;
+    x[4] = ray.d.y;
+    x[5] = ray.d.z;
+
+    boost::numeric::odeint::runge_kutta4< state_type > rk;
+
+    double dt = 0.5;
+    double t = 0.0;
+
+    state_type prev = x;
+    for(int i = 0; i<ITERATIONS; ++i, t+=dt ){
+        prev[0] = x[0];
+        prev[1] = x[1];
+        prev[2] = x[2];
+        prev[3] = x[3];
+        prev[4] = x[4];
+        prev[5] = x[5];
+        rk.do_step(ray, x, t, dt);
+        int stage = forward(prev, x, isect);
+        /*if(stage==1)break;
+        else if(stage == 0)
+        {
+            x[0] = prev[0];
+            x[1] = prev[1];
+            x[2] = prev[2];
+            x[3] = prev[3];
+            x[4] = prev[4];
+            x[5] = prev[5];
+            dt = dt/2.0;
+            t -= dt;
+        }
+    }
+    if(id==-1)
+      return false;
+    else{
+      o.x = x[0];
+      o.y = x[1];
+      o.z = x[2];
+
+      d.x = x[3];
+      d.y = x[4];
+      d.z = x[5];
+      d = d.norm();
+      return true;
+    }*/
+    }
+
+    return aggregate->Intersect(ray, isect);
+}
+
 }  // namespace pbrt

src/core/scene.h

@@ -65,6 +65,10 @@ class Scene {
     bool IntersectP(const Ray &ray) const;
     bool IntersectTr(Ray ray, Sampler &sampler, SurfaceInteraction *isect,
                      Spectrum *transmittance) const;
+    // curved spacetime intersection
+    bool IntersectCu(const Ray &ray, SurfaceInteraction *isect) const;
+
+    int forward(state_type &prev, state_type &state, SurfaceInteraction *isect) const;
 
     // Scene Public Data
     std::vector<std::shared_ptr<Light>> lights;