diff --git a/src/pe/raytracing/Intersects.h b/src/pe/raytracing/Intersects.h
index ae045d6d3fc20e91ab4e0f1a91d8137f0aba38a0..d8ffff061bc7a9d02ea9a506c01a799c96e65685 100644
--- a/src/pe/raytracing/Intersects.h
+++ b/src/pe/raytracing/Intersects.h
@@ -28,34 +28,39 @@
#include "pe/rigidbody/Sphere.h"
#include "pe/rigidbody/Union.h"
#include "pe/utility/BodyCast.h"
+#include <boost/math/special_functions/sign.hpp>
+#include <boost/tuple/tuple.hpp>
#include <pe/raytracing/Ray.h>
-#include <boost/tuple/tuple.hpp>
-
#define EPSILON real_t(1e-4)
+using namespace boost::math;
+
namespace walberla {
namespace pe {
namespace raytracing {
inline bool intersects(const AABB& aabb, const Ray& ray, real_t& t, real_t padding = real_t(0.0));
-inline bool intersects(const SphereID sphere, const Ray& ray, real_t& t);
-inline bool intersects(const PlaneID plane, const Ray& ray, real_t& t);
-inline bool intersects(const BoxID box, const Ray& ray, real_t& t);
+inline bool intersects(const SphereID sphere, const Ray& ray, real_t& t, Vec3& n);
+inline bool intersects(const PlaneID plane, const Ray& ray, real_t& t, Vec3& n);
+inline bool intersects(const BoxID box, const Ray& ray, real_t& t, Vec3& n);
struct IntersectsFunctor
{
const Ray& ray_;
real_t& t_;
+ Vec3& n_;
- IntersectsFunctor(const Ray& ray, real_t& t) : ray_(ray), t_(t) {}
+ IntersectsFunctor(const Ray& ray, real_t& t, Vec3& n) : ray_(ray), t_(t), n_(n) {}
template< typename BodyType >
- bool operator()( BodyType* bd1 ) { return intersects( bd1, ray_, t_); }
+ bool operator()( BodyType* bd1 ) {
+ return intersects( bd1, ray_, t_, n_);
+ }
};
-inline bool intersects(const SphereID sphere, const Ray& ray, real_t& t) {
+inline bool intersects(const SphereID sphere, const Ray& ray, real_t& t, Vec3& n) {
real_t inf = std::numeric_limits<real_t>::max();
const Vec3& direction = ray.getDirection();
Vec3 displacement = ray.getOrigin() - sphere->getPosition();
@@ -89,22 +94,24 @@ inline bool intersects(const SphereID sphere, const Ray& ray, real_t& t) {
t = t1;
}
}
+ Vec3 intersectionPoint = ray.getOrigin() + direction*t;
+ n = (intersectionPoint - sphere->getPosition()).getNormalized();
+
return true;
}
-inline bool intersects(const PlaneID plane, const Ray& ray, real_t& t) {
+inline bool intersects(const PlaneID plane, const Ray& ray, real_t& t, Vec3& n) {
real_t inf = std::numeric_limits<real_t>::max();
const Vec3& direction = ray.getDirection();
const Vec3& origin = ray.getOrigin();
- real_t denominator = plane->getNormal() * direction;
+ const Vec3& planeNormal = plane->getNormal();
+ real_t denominator = planeNormal * direction;
if (std::abs(denominator) > EPSILON) {
real_t t_;
- t_ = ((plane->getPosition() - origin) * plane->getNormal()) / denominator;
-
+ t_ = ((plane->getPosition() - origin) * planeNormal) / denominator;
if (t_ > EPSILON) {
- Vec3 intersectionPoint = (t_*direction + origin);
- Vec3 originToIntersection = origin - intersectionPoint;
- t = originToIntersection.length();
+ t = t_;
+ n = planeNormal * sign(-denominator);
return true;
} else {
t = inf;
@@ -114,7 +121,7 @@ inline bool intersects(const PlaneID plane, const Ray& ray, real_t& t) {
return false;
}
-inline bool intersects(const BoxID box, const Ray& ray, real_t& t) {
+inline bool intersects(const BoxID box, const Ray& ray, real_t& t, Vec3& n) {
Ray transformedRay(box->pointFromWFtoBF(ray.getOrigin()), box->vectorFromWFtoBF(ray.getDirection()));
const Vec3& lengths = box->getLengths();
@@ -131,6 +138,7 @@ inline bool intersects(const BoxID box, const Ray& ray, real_t& t) {
real_t inf = std::numeric_limits<real_t>::max();
+ size_t tminAxis = 0, tmaxAxis = 0;
real_t txmin, txmax;
real_t tmin = txmin = (bounds[sign[0]][0] - origin[0]) * invDirection[0];
real_t tmax = txmax = (bounds[1-sign[0]][0] - origin[0]) * invDirection[0];
@@ -141,9 +149,11 @@ inline bool intersects(const BoxID box, const Ray& ray, real_t& t) {
return false;
}
if (tymin > tmin) {
+ tminAxis = 1;
tmin = tymin;
}
if (tymax < tmax) {
+ tmaxAxis = 1;
tmax = tymax;
}
real_t tzmin = (bounds[sign[2]][2] - origin[2]) * invDirection[2];
@@ -153,16 +163,20 @@ inline bool intersects(const BoxID box, const Ray& ray, real_t& t) {
return false;
}
if (tzmin > tmin) {
+ tminAxis = 2;
tmin = tzmin;
}
if (tzmax < tmax) {
+ tmaxAxis = 2;
tmax = tzmax;
}
+ n[0] = n[1] = n[2] = real_t(0);
real_t t_;
if (tmin > 0) {
// ray hit box from outside
t_ = tmin;
+ n[tminAxis] = real_t(1);
} else if (tmax < 0) {
// tmin and tmax are smaller than 0 -> box is in rays negative direction
t = inf;
@@ -170,8 +184,16 @@ inline bool intersects(const BoxID box, const Ray& ray, real_t& t) {
} else {
// ray origin within box
t_ = tmax;
+ n[tmaxAxis] = real_t(1);
+ }
+
+ if (transformedRay.getDirection() * n > 0) {
+ n = -n;
}
+ n = box->vectorFromBFtoWF(n);
+ WALBERLA_LOG_INFO("t_: " << t_ << ", n: " << n);
+
t = t_;
return true;
}
diff --git a/src/pe/raytracing/Raytracer.h b/src/pe/raytracing/Raytracer.h
index 3c2cc9959b336ef34e89877060bcfb69eaa814d4..1a1c239819ea214f22d091dae5fef74deb91a4f4 100644
--- a/src/pe/raytracing/Raytracer.h
+++ b/src/pe/raytracing/Raytracer.h
@@ -254,10 +254,10 @@ void Raytracer::rayTrace(const size_t timestep) {
std::vector<BodyIntersectionInfo> intersections; // contains for each pixel information about an intersection, if existent
real_t t, t_closest;
- walberla::id_t id_closest;
+ Vec3 n;
RigidBody* body_closest = NULL;
Ray ray(cameraPosition_, Vec3(1,0,0));
- IntersectsFunctor func(ray, t);
+ IntersectsFunctor func(ray, t, n);
tp_["Raytracing"].start();
for (size_t x = 0; x < pixelsHorizontal_; x++) {
for (size_t y = 0; y < pixelsVertical_; y++) {
@@ -265,8 +265,8 @@ void Raytracer::rayTrace(const size_t timestep) {
Vec3 direction = (pixelLocation - cameraPosition_).getNormalized();
ray.setDirection(direction);
+ n[0] = n[1] = n[2] = real_t(0);
t_closest = inf;
- id_closest = 0;
body_closest = NULL;
for (auto blockIt = forest_->begin(); blockIt != forest_->end(); ++blockIt) {
#ifndef DISABLE_BLOCK_AABB_INTERSECTION_PRECHECK
@@ -288,7 +288,6 @@ void Raytracer::rayTrace(const size_t timestep) {
if (intersects && t < t_closest) {
// body was shot by ray and currently closest to camera
t_closest = t;
- id_closest = bodyIt->getID();
body_closest = *bodyIt;
}
}
@@ -310,7 +309,6 @@ void Raytracer::rayTrace(const size_t timestep) {
if (intersects && t < t_closest) {
// body was shot by ray and currently closest to camera
t_closest = t;
- id_closest = bodyIt->getID();
body_closest = *bodyIt;
}
}
diff --git a/tests/pe/Raytracing.cpp b/tests/pe/Raytracing.cpp
index 7e1e5f22574774d9f281dc11a7b057ec337aceb1..3722f69f8f9ef131cc46e7cf0fb91ba19bb45275 100644
--- a/tests/pe/Raytracing.cpp
+++ b/tests/pe/Raytracing.cpp
@@ -16,6 +16,8 @@
#include "core/DataTypes.h"
#include "core/math/Vector3.h"
+#include <pe/raytracing/Ray.h>
+#include <pe/raytracing/Intersects.h>
#include <pe/raytracing/Raytracer.h>
using namespace walberla;
@@ -29,25 +31,29 @@ void SphereIntersectsTest()
MaterialID iron = Material::find("iron");
Sphere sp1(123, 1, Vec3(3,3,3), Vec3(0,0,0), Quat(), 2, iron, false, true, false);
real_t t;
+ Vec3 n;
// ray through the center
Ray ray1(Vec3(3,-5,3), Vec3(0,1,0));
WALBERLA_LOG_INFO("RAY -> SPHERE");
- WALBERLA_CHECK(intersects(&sp1, ray1, t));
+ WALBERLA_CHECK(intersects(&sp1, ray1, t, n));
WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(6));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(-1));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0));
// ray tangential
Ray ray2(Vec3(3,-5,3), Vec3(0,7.5,real_t(std::sqrt(real_t(15))/real_t(2))).getNormalized());
- WALBERLA_CHECK(!intersects(&sp1, ray2, t));
+ WALBERLA_CHECK(!intersects(&sp1, ray2, t, n));
// sphere behind ray origin
Sphere sp2(123, 1, Vec3(3,-8,3), Vec3(0,0,0), Quat(), 2, iron, false, true, false);
- WALBERLA_CHECK(!intersects(&sp2, ray1, t));
+ WALBERLA_CHECK(!intersects(&sp2, ray1, t, n));
// sphere around ray origin
Sphere sp3(123, 1, Vec3(3,-5,3), Vec3(0,0,0), Quat(), 2, iron, false, true, false);
- WALBERLA_CHECK(intersects(&sp3, ray1, t));
+ WALBERLA_CHECK(intersects(&sp3, ray1, t, n));
WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(2));
}
@@ -58,22 +64,40 @@ void PlaneIntersectsTest() {
Ray ray1(Vec3(-5,3,3), Vec3(1,0,0));
real_t t;
+ Vec3 n;
WALBERLA_LOG_INFO("RAY -> PLANE");
- WALBERLA_CHECK(intersects(&pl1, ray1, t), "ray through center did not hit");
+ WALBERLA_CHECK(intersects(&pl1, ray1, t, n), "ray through center did not hit");
WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(8), "distance between ray and plane is incorrect");
Ray ray2(Vec3(-5,3,3), Vec3(1,0,-1).getNormalized());
- WALBERLA_CHECK(intersects(&pl1, ray2, t), "ray towards random point on plane didn't hit");
+ WALBERLA_CHECK(intersects(&pl1, ray2, t, n), "ray towards random point on plane didn't hit");
WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(sqrt(real_t(128))), "distance between ray and plane is incorrect");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+
+ Plane pl1neg(1, 1, Vec3(3, 3, 3), Vec3(-1, 0, 0), real_t(1.0), iron);
+ WALBERLA_CHECK(intersects(&pl1neg, ray2, t, n), "ray towards random point on plane didn't hit");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+
+ Ray ray3(Vec3(-5,3,3), Vec3(-1,0,0).getNormalized());
+ Plane pl5(1, 1, Vec3(-7, 3, 3), Vec3(1, 0, 0), real_t(1.0), iron);
+ WALBERLA_CHECK(intersects(&pl5, ray3, t, n), "ray towards random point on plane didn't hit");
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(2), "distance between ray and plane is incorrect");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
// plane with center 3,3,3 and parallel to x-z plane
Plane pl2(1, 1, Vec3(3, 3, 3), Vec3(0, 1, 0), real_t(1.0), iron);
- WALBERLA_CHECK(!intersects(&pl2, ray1, t), "ray parallel to plane shouldnt hit");
+ WALBERLA_CHECK(!intersects(&pl2, ray1, t, n), "ray parallel to plane shouldnt hit");
// plane with center -10,3,3 and parallel to y-z plane
Plane pl4(1, 1, Vec3(-10, 3, 3), Vec3(1, 0, 0), real_t(1.0), iron);
- WALBERLA_CHECK(!intersects(&pl4, ray1, t), "ray hit plane behind origin");
+ WALBERLA_CHECK(!intersects(&pl4, ray1, t, n), "ray hit plane behind origin");
}
void BoxIntersectsTest() {
@@ -81,37 +105,55 @@ void BoxIntersectsTest() {
MaterialID iron = Material::find("iron");
real_t t;
+ Vec3 n;
Box box1(127, 5, Vec3(0, -15, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
Ray ray1(Vec3(3,-5,3), Vec3(0,1,0));
- WALBERLA_CHECK(!intersects(&box1, ray1, t));
+ WALBERLA_CHECK(!intersects(&box1, ray1, t, n));
Box box2(128, 5, Vec3(0, -2, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
- WALBERLA_CHECK(intersects(&box2, ray1, t));
+ WALBERLA_CHECK(intersects(&box2, ray1, t, n));
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(8), real_t(1e-7));
Box box3(128, 5, Vec3(0, 5, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
- WALBERLA_CHECK(intersects(&box3, ray1, t));
+ WALBERLA_CHECK(intersects(&box3, ray1, t, n));
WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(5));
+ Ray ray6(Vec3(-8,5,0), Vec3(1,0,0));
+ WALBERLA_CHECK(intersects(&box3, ray6, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(3));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+
+ Ray ray7(Vec3(8,5,0), Vec3(-1,0,0));
+ WALBERLA_CHECK(intersects(&box3, ray7, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(3));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+
// ray origin within box
Ray ray2(Vec3(-2,0,0), Vec3(1,0,1).getNormalized());
- WALBERLA_CHECK(intersects(&box3, ray2, t));
+ WALBERLA_CHECK(intersects(&box3, ray2, t, n));
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(7.0710), real_t(1e-4));
Ray ray3(Vec3(3,-5,3), Vec3(2, -1.5, 0.5).getNormalized());
Box box4(128, 5, Vec3(0, 8, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
- WALBERLA_CHECK(!intersects(&box4, ray3, t));
+ WALBERLA_CHECK(!intersects(&box4, ray3, t, n));
Ray ray4(Vec3(3,-5,3), Vec3(-2, 3, 0.5).getNormalized());
- WALBERLA_CHECK(intersects(&box4, ray4, t));
+ WALBERLA_CHECK(intersects(&box4, ray4, t, n));
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(9.7068), real_t(1e-4));
Box box5(128, 5, Vec3(4, 0, 0), Vec3(0, 0, 0), Quat(), Vec3(4, 4, 4), iron, false, true, false);
box5.rotate(0,0,math::M_PI/4);
Ray ray5(Vec3(0,1.5,0), Vec3(1,0,0));
- WALBERLA_CHECK(intersects(&box5, ray5, t));
+ WALBERLA_CHECK(intersects(&box5, ray5, t, n));
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(2.67157), real_t(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(n[0], real_t(-0.707107), real_t(1e-5), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(n[1], real_t(0.707107), real_t(1e-5), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
}
void AABBIntersectsTest() {