Skip to content
Snippets Groups Projects
Commit 9ac65116 authored by Lukas Werner's avatar Lukas Werner
Browse files

Added ray - capsule intersection test

parent 08b57b88
1 merge request!89Add image generation using a raytracer
Hide whitespace changes
Inline Side-by-side
src/pe/raytracing/Intersects.h
+ 156
− 20
View file @ 9ac65116
......@@ -40,13 +40,14 @@ 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, 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);
inline bool intersects(const CapsuleID capsule, const Ray& ray, real_t& t, Vec3& n);
inline bool intersects(const AABB& aabb, const Ray& ray, real_t& t, real_t padding = real_t(0.0));
inline bool intersectsSphere(const Vec3& gpos, real_t radius, const Ray& ray, real_t& t0, real_t& t1);
struct IntersectsFunctor
{
const Ray& ray_;
......@@ -63,26 +64,13 @@ struct IntersectsFunctor
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();
real_t a = direction * direction;
real_t b = real_t(2.) * (displacement * direction);
real_t c = (displacement * displacement) - (sphere->getRadius() * sphere->getRadius());
real_t discriminant = b*b - real_t(4.)*a*c;
if (discriminant < EPSILON) {
// with discriminant smaller than 0, sphere is not hit by ray
// (no solution for quadratic equation)
// with discriminant being 0, sphere only tangentially hits the ray (not enough)
t = inf;
return false;
}
real_t root = real_t(std::sqrt(discriminant));
real_t t0 = (-b - root) / (real_t(2.) * a); // point where the ray enters the sphere
real_t t1 = (-b + root) / (real_t(2.) * a); // point where the ray leaves the sphere
if (t0 < 0 && t1 < 0) {
real_t t0, t1;
if (!intersectsSphere(sphere->getPosition(), sphere->getRadius(), ray, t0, t1)) {
t = inf;
return false;
}
t = (t0 < t1) ? t0 : t1; // assign the closest hit point to t
if (t < 0) {
// at least one of the calculated hit points is behind the rays origin
......@@ -95,7 +83,8 @@ inline bool intersects(const SphereID sphere, const Ray& ray, real_t& t, Vec3& n
t = t1;
}
}
Vec3 intersectionPoint = ray.getOrigin() + direction*t;
Vec3 intersectionPoint = ray.getOrigin() + ray.getDirection()*t;
n = (intersectionPoint - sphere->getPosition()).getNormalized();
return true;
......@@ -196,6 +185,153 @@ inline bool intersects(const BoxID box, const Ray& ray, real_t& t, Vec3& n) {
t = t_;
return true;
}
inline bool intersects(const CapsuleID capsule, const Ray& ray, real_t& t, Vec3& n) {
Ray transformedRay = ray.transformedToBF(capsule);
Vec3 direction = transformedRay.getDirection();
Vec3 origin = transformedRay.getOrigin();
real_t halfLength = capsule->getLength()/real_t(2);
real_t inf = std::numeric_limits<real_t>::max();
t = inf;
bool t0hit = false, t1hit = false;
size_t intersectedPrimitive = 0; // 1 for capsule, 2 for left half sphere, 3 for right half sphere
real_t a = direction[2]*direction[2] + direction[1]*direction[1];
real_t b = real_t(2)*origin[2]*direction[2] + real_t(2)*origin[1]*direction[1];
real_t c = origin[2]*origin[2] + origin[1]*origin[1] - capsule->getRadius();
real_t discriminant = b*b - real_t(4.)*a*c;
if (discriminant >= EPSILON) {
// With discriminant smaller than 0, cylinder is not hit by ray (no solution for quadratic equation).
// Thus only enter this section if the equation is actually solvable.
real_t root = real_t(std::sqrt(discriminant));
real_t t0 = (-b - root) / (real_t(2) * a); // point where the ray enters the cylinder
real_t t1 = (-b + root) / (real_t(2) * a); // point where the ray leaves the cylinder
real_t tx0 = origin[0] + direction[0]*t0;
real_t tx1 = origin[0] + direction[0]*t1;
if (t0 > 0 && tx0 >= -halfLength && tx0 <= halfLength) {
t0hit = true;
intersectedPrimitive = 1;
t = t0;
}
if (t1 > 0 && tx1 >= -halfLength && tx1 <= halfLength && t1 < t) {
t1hit = true;
if (t1 < t) {
intersectedPrimitive = 1;
t = t1;
}
}
}
// check now for end capping half spheres.
if (!t0hit || !t1hit) {
// only check capping half spheres if the ray didnt both enter and leave the cylinder part of the capsule already.
real_t t0_left, t1_left;
Vec3 leftSpherePos(-halfLength, 0, 0);
if (intersectsSphere(leftSpherePos, capsule->getRadius(), transformedRay, t0_left, t1_left)) {
// at least one of t0_left and t1_left are not behind the rays origin
real_t t0x_left = origin[0] + direction[0]*t0_left;
real_t t1x_left = origin[0] + direction[0]*t1_left;
real_t t_left;
if (t0_left > 0 && t0x_left < -halfLength) {
t_left = t0_left;
}
if (t1_left > 0 && t1x_left < -halfLength && t1_left < t_left) {
t_left = t1_left;
}
if (t_left < t) {
intersectedPrimitive = 2;
t = t_left;
}
}
real_t t0_right, t1_right;
Vec3 rightSpherePos(halfLength, 0, 0);
if (intersectsSphere(rightSpherePos, capsule->getRadius(), transformedRay, t0_right, t1_right)) {
// at least one of t0_right and t1_right are not behind the rays origin
real_t t0x_right = origin[0] + direction[0]*t0_right;
real_t t1x_right = origin[0] + direction[0]*t1_right;
real_t t_right;
if (t0_right > 0 && t0x_right > halfLength) {
t_right = t0_right;
}
if (t1_right > 0 && t1x_right > halfLength && t1_right < t_right) {
t_right = t1_right;
}
if (t_right < t) {
intersectedPrimitive = 3;
t = t_right;
}
}
if (realIsIdentical(t, inf)) {
return false;
}
Vec3 intersectionPoint = origin + direction*t;
if (intersectedPrimitive == 2) {
n = (intersectionPoint - leftSpherePos).getNormalized();
} else if (intersectedPrimitive == 3) {
n = (intersectionPoint - rightSpherePos).getNormalized();
}
}
WALBERLA_ASSERT(intersectedPrimitive != 0);
if (intersectedPrimitive == 1) {
Vec3 intersectionPoint = origin + direction*t;
Vec3 intersectionPointOnXAxis(intersectionPoint[0], real_t(0), real_t(0));
n = (intersectionPoint - intersectionPointOnXAxis).getNormalized();
}
n = capsule->vectorFromBFtoWF(n);
return true;
}
inline bool intersectsSphere(const Vec3& gpos, real_t radius, const Ray& ray, real_t& t0, real_t& t1) {
real_t inf = std::numeric_limits<real_t>::max();
const Vec3& direction = ray.getDirection();
Vec3 displacement = ray.getOrigin() - gpos;
real_t a = direction * direction;
real_t b = real_t(2.) * (displacement * direction);
real_t c = (displacement * displacement) - (radius * radius);
real_t discriminant = b*b - real_t(4.)*a*c;
if (discriminant < 0) {
// with discriminant smaller than 0, sphere is not hit by ray
// (no solution for quadratic equation)
// with discriminant being 0, sphere only tangentially hits the ray (not enough)
t0 = inf;
t1 = inf;
return false;
}
real_t root = real_t(std::sqrt(discriminant));
t0 = (-b - root) / (real_t(2.) * a); // distance to point where the ray enters the sphere
t1 = (-b + root) / (real_t(2.) * a); // distance to point where the ray leaves the sphere
if (t0 < 0 && t1 < 0) {
return false;
}
real_t t = (t0 < t1) ? t0 : t1; // assign the closest distance to t
if (t < 0) {
// at least one of the calculated distances is behind the rays origin
if (t1 < 0) {
// both of the points are behind the origin (ray does not hit sphere)
return false;
}
}
return true;
}
inline bool intersects(const AABB& aabb, const Ray& ray, real_t& t, real_t padding) {
// An Efficient and Robust Ray–Box Intersection Algorithm: http://people.csail.mit.edu/amy/papers/box-jgt.pdf
......
tests/pe/Raytracing.cpp
+ 53
− 20
View file @ 9ac65116
......@@ -24,7 +24,7 @@ using namespace walberla;
using namespace walberla::pe;
using namespace walberla::pe::raytracing;
typedef boost::tuple<Box, Plane, Sphere> BodyTuple ;
typedef boost::tuple<Box, Plane, Sphere, Capsule> BodyTuple ;
void SphereIntersectsTest()
{
......@@ -183,19 +183,44 @@ void AABBIntersectsTest() {
WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(4));
}
void CapsuleIntersectsTest() {
MaterialID iron = Material::find("iron");
real_t t;
Vec3 n;
Capsule cp1(0, 0, Vec3(2,3,3), Vec3(0,0,0), Quat(), real_t(2), real_t(2), iron, false, true, false);
// ray through the center
Ray ray1(Vec3(3,-5,3), Vec3(0,1,0));
WALBERLA_LOG_INFO("RAY -> CAPSULE");
WALBERLA_CHECK(intersects(&cp1, 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 ray2(Vec3(-5,3,3), Vec3(1,0,0));
WALBERLA_CHECK(intersects(&cp1, ray2, t, n));
WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(4));
WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1));
WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0));
WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0));
}
void RaytracerTest() {
WALBERLA_LOG_INFO("Raytracer");
shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
shared_ptr<BlockForest> forest = createBlockForest(AABB(0,-5,-5,10,5,5), Vec3(1,1,1), Vec3(false, false, false));
shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,10,10,10), Vec3(1,1,1), Vec3(false, false, false));
auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
Lighting lighting(Vec3(0, 3, 3),
Lighting lighting(Vec3(0, 5, 8), // 8, 5, 9.5 gut für ebenen, 0,5,8
Vec3(1, 1, 1), //diffuse
Vec3(1, 1, 1), //specular
Vec3(0.4, 0.4, 0.4)); //ambient
Raytracer raytracer(forest, storageID, globalBodyStorage,
size_t(640), size_t(480),
49.13,
Vec3(-5,0,0), Vec3(-1,0,0), Vec3(0,0,1),
Vec3(-5,5,5), Vec3(-1,5,5), Vec3(0,0,1),
lighting);
MaterialID iron = Material::find("iron");
......@@ -206,26 +231,33 @@ void RaytracerTest() {
//PlaneID xNegPlaneClose = createPlane(*globalBodyStorage, 0, Vec3(-1,0,0), Vec3(1,0,0), iron);
// Test Scene v1 - Spheres, (rotated) boxes, confining walls, tilted plane in right bottom back corner
createPlane(*globalBodyStorage, 0, Vec3(0,-1,0), Vec3(0,5,0), iron); // left wall
createPlane(*globalBodyStorage, 0, Vec3(0,1,0), Vec3(0,-5,0), iron); // right wall
createPlane(*globalBodyStorage, 0, Vec3(0,0,1), Vec3(0,0,-5), iron); // floor
createPlane(*globalBodyStorage, 0, Vec3(0,0,-1), Vec3(0,0,5), iron); // ceiling
createPlane(*globalBodyStorage, 0, Vec3(0,-1,0), Vec3(0,10,0), iron); // left wall
createPlane(*globalBodyStorage, 0, Vec3(0,1,0), Vec3(0,0,0), iron); // right wall
createPlane(*globalBodyStorage, 0, Vec3(0,0,1), Vec3(0,0,0), iron); // floor
createPlane(*globalBodyStorage, 0, Vec3(0,0,-1), Vec3(0,0,10), iron); // ceiling
createPlane(*globalBodyStorage, 0, Vec3(-1,0,0), Vec3(10,0,0), iron); // back wall
createPlane(*globalBodyStorage, 0, Vec3(1,0,0), Vec3(0,0,0), iron); // front wall, should not get rendered
createPlane(*globalBodyStorage, 0, Vec3(-1,1,1), Vec3(8,-3,-3), iron); // tilted plane in right bottom back corner
createPlane(*globalBodyStorage, 0, Vec3(-1,1,1), Vec3(8,2,2), iron); // tilted plane in right bottom back corner
createSphere(*globalBodyStorage, *forest, storageID, 2, Vec3(6,4.5,4.5), real_t(0.5));
createSphere(*globalBodyStorage, *forest, storageID, 3, Vec3(4,0.5,0), real_t(1));
createSphere(*globalBodyStorage, *forest, storageID, 6, Vec3(3,3.5,0), real_t(1));
BoxID box = createBox(*globalBodyStorage, *forest, storageID, 7, Vec3(5,1.5,0), Vec3(2,4,3));
createSphere(*globalBodyStorage, *forest, storageID, 2, Vec3(6,9.5,9.5), real_t(0.5));
createSphere(*globalBodyStorage, *forest, storageID, 3, Vec3(4,5.5,5), real_t(1));
createSphere(*globalBodyStorage, *forest, storageID, 6, Vec3(3,8.5,5), real_t(1));
BoxID box = createBox(*globalBodyStorage, *forest, storageID, 7, Vec3(5,6.5,5), Vec3(2,4,3));
box->rotate(0,math::M_PI/4,math::M_PI/4);
createBox(*globalBodyStorage, *forest, storageID, 8, Vec3(5,-4,3), Vec3(2,2,2));
createBox(*globalBodyStorage, *forest, storageID, 8, Vec3(5,1,8), Vec3(2,2,2));
// Test scene v1 end
raytracer.setTBufferOutputDirectory("/Users/ng/Desktop/walberla");
// Test scene v2 additions start
createBox(*globalBodyStorage, *forest, storageID, 7, Vec3(9,9,5), Vec3(1,1,10));
createCapsule(*globalBodyStorage, *forest, storageID, 9, Vec3(3, 9, 1), real_t(0.5), real_t(7), iron);
CapsuleID capsule2 = createCapsule(*globalBodyStorage, *forest, storageID, 9, Vec3(7, 3.5, 7.5), real_t(1), real_t(2), iron);
capsule2->rotate(0,math::M_PI/4,math::M_PI/4-math::M_PI/8);
// Test scene v2 end
raytracer.setTBufferOutputDirectory("tbuffer");
raytracer.setTBufferOutputEnabled(true);
raytracer.setImageOutputDirectory("/Users/ng/Desktop/walberla");
raytracer.setImageOutputDirectory("image");
raytracer.setImageOutputEnabled(true);
raytracer.rayTrace<BodyTuple>(0);
......@@ -238,10 +270,11 @@ int main( int argc, char** argv )
SetBodyTypeIDs<BodyTuple>::execute();
//SphereIntersectsTest();
//PlaneIntersectsTest();
//BoxIntersectsTest();
//AABBIntersectsTest();
SphereIntersectsTest();
PlaneIntersectsTest();
BoxIntersectsTest();
AABBIntersectsTest();
CapsuleIntersectsTest();
RaytracerTest();
return EXIT_SUCCESS;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment