https://github.com/ialhashim/topo-blend
Tip revision: 39b13612ebd645a65eda854771b517371f2f858a authored by ennetws on 13 March 2015, 18:17:18 UTC
Create README.md
Create README.md
Tip revision: 39b1361
BoundingBox.h
#pragma once
#include <float.h>
/* AABB-triangle overlap test code */
/* by Tomas Akenine-Möller */
#define X 0
#define Y 1
#define Z 2
#define FINDMINMAX(x0,x1,x2,min,max) \
min = max = x0; \
if(x1<min) min=x1;\
if(x1>max) max=x1;\
if(x2<min) min=x2;\
if(x2>max) max=x2;
static inline int planeBoxOverlap(const Vec3d& normal, const Vec3d& vert, const Vec3d& maxbox)
{
Vec3d vmin,vmax;
for(int q=X; q<=Z; q++){
double v = vert[q];
if(normal[q] > 0.0){
vmin[q]=-maxbox[q] - v;
vmax[q]= maxbox[q] - v;
}
else{
vmin[q]= maxbox[q] - v;
vmax[q]=-maxbox[q] - v;
}
}
if(dot(normal , vmin) > 0.0) return 0;
if(dot(normal , vmax) >= 0.0) return 1;
return 0;
}
/*======================== X-tests ========================*/
#define AXISTEST_X01(a, b, fa, fb) \
p0 = a*v0[Y] - b*v0[Z]; \
p2 = a*v2[Y] - b*v2[Z]; \
if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z]; \
if(min>rad || max<-rad) return 0;
#define AXISTEST_X2(a, b, fa, fb) \
p0 = a*v0[Y] - b*v0[Z]; \
p1 = a*v1[Y] - b*v1[Z]; \
if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z]; \
if(min>rad || max<-rad) return 0;
/*======================== Y-tests ========================*/
#define AXISTEST_Y02(a, b, fa, fb) \
p0 = -a*v0[X] + b*v0[Z]; \
p2 = -a*v2[X] + b*v2[Z]; \
if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z]; \
if(min>rad || max<-rad) return 0;
#define AXISTEST_Y1(a, b, fa, fb) \
p0 = -a*v0[X] + b*v0[Z]; \
p1 = -a*v1[X] + b*v1[Z]; \
if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z]; \
if(min>rad || max<-rad) return 0;
/*======================== Z-tests ========================*/
#define AXISTEST_Z12(a, b, fa, fb) \
p1 = a*v1[X] - b*v1[Y]; \
p2 = a*v2[X] - b*v2[Y]; \
if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;} \
rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y]; \
if(min>rad || max<-rad) return 0;
#define AXISTEST_Z0(a, b, fa, fb) \
p0 = a*v0[X] - b*v0[Y]; \
p1 = a*v1[X] - b*v1[Y]; \
if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y]; \
if(min>rad || max<-rad) return 0;
#undef X
#undef Y
#undef Z
struct Ray
{
Vec3d origin;
Vec3d direction;
int index;
double thickness;
Ray(const Vec3d & Origin = Vec3d(), const Vec3d & Direction = Vec3d(), double Thickness = 0.0, int Index = -1) : origin(Origin), index(Index){
direction = Direction.normalized();
thickness = Thickness;
}
inline Ray inverse() const { return Ray(origin, -direction); }
Ray& operator= (const Ray& other){
this->origin = other.origin;
this->direction = other.direction;
this->index = other.index;
this->thickness = other.thickness;
return *this;
}
};
class BoundingBox
{
public:
Vec3d center;
Vec3d vmax, vmin;
double xExtent, yExtent, zExtent;
BoundingBox()
{
this->center = Vec3d(0,0,0);
this->xExtent = 0;
this->yExtent = 0;
this->zExtent = 0;
}
BoundingBox( const Vec3d& c, double x, double y, double z )
{
this->center = c;
this->xExtent = x;
this->yExtent = y;
this->zExtent = z;
Vec3d corner(x, y, z);
vmin = center - corner;
vmax = center + corner;
}
BoundingBox& operator=( const BoundingBox& other )
{
this->center = other.center;
this->xExtent = other.xExtent;
this->yExtent = other.yExtent;
this->zExtent = other.zExtent;
this->vmax = other.vmax;
this->vmin = other.vmin;
return *this;
}
void computeFromTris( const std::vector< std::vector<Vec3d> > & tris )
{
vmin = Vec3d(DBL_MAX, DBL_MAX, DBL_MAX);
vmax = -vmin;
double minx = 0, miny = 0, minz = 0;
double maxx = 0, maxy = 0, maxz = 0;
minx = maxx = tris[0][0].x();
miny = maxy = tris[0][1].y();
minz = maxz = tris[0][2].z();
for (int i = 0; i < (int)tris.size(); i++)
{
for(int v = 0; v < 3; v++)
{
Vec3d vec = tris[i][v];
if (vec.x() < minx) minx = vec.x();
if (vec.x() > maxx) maxx = vec.x();
if (vec.y() < miny) miny = vec.y();
if (vec.y() > maxy) maxy = vec.y();
if (vec.z() < minz) minz = vec.z();
if (vec.z() > maxz) maxz = vec.z();
}
}
vmax = Vec3d(maxx, maxy, maxz);
vmin = Vec3d(minx, miny, minz);
this->center = (vmin + vmax) / 2.0;
this->xExtent = abs(vmax.x() - center.x());
this->yExtent = abs(vmax.y() - center.y());
this->zExtent = abs(vmax.z() - center.z());
}
std::vector<Vec3d> getCorners()
{
std::vector<Vec3d> corners;
Vec3d x = (Vec3d(1,0,0) * xExtent);
Vec3d y = (Vec3d(0,1,0) * yExtent);
Vec3d z = (Vec3d(0,0,1) * zExtent);
Vec3d c = center + x + y + z;
corners.push_back(c);
corners.push_back(c - (x*2));
corners.push_back(c - (y*2));
corners.push_back(c - (x*2) - (y*2));
corners.push_back(corners[0] - (z*2));
corners.push_back(corners[1] - (z*2));
corners.push_back(corners[2] - (z*2));
corners.push_back(corners[3] - (z*2));
return corners;
}
bool intersects( const Ray& ray ) const
{
// r.dir is unit direction vector of ray
Vec3d dirfrac;
dirfrac.x() = 1.0f / ray.direction.x();
dirfrac.y() = 1.0f / ray.direction.y();
dirfrac.z() = 1.0f / ray.direction.z();
double overlap = ray.thickness;
Vec3d minv = vmin + ((vmin - center).normalized() * overlap);
Vec3d maxv = vmax + ((vmax - center).normalized() * overlap);
// lb is the corner of AABB with minimal coordinates - left bottom, rt is maximal corner
// ray.origin is origin of ray
float t1 = (minv.x() - ray.origin.x())*dirfrac.x();
float t2 = (maxv.x() - ray.origin.x())*dirfrac.x();
float t3 = (minv.y() - ray.origin.y())*dirfrac.y();
float t4 = (maxv.y() - ray.origin.y())*dirfrac.y();
float t5 = (minv.z() - ray.origin.z())*dirfrac.z();
float t6 = (maxv.z() - ray.origin.z())*dirfrac.z();
float tmin = qMax(qMax(qMin(t1, t2), qMin(t3, t4)), qMin(t5, t6));
float tmax = qMin(qMin(qMax(t1, t2), qMax(t3, t4)), qMax(t5, t6));
double t = 0;
// if tmax < 0, ray (line) is intersecting AABB, but whole AABB is behing us
if (tmax < 0)
{
t = tmax;
return false;
}
// if tmin > tmax, ray doesn't intersect AABB
if (tmin > tmax)
{
t = tmax;
return false;
}
t = tmin;
return true;
}
bool intersectsWorking( const Ray& ray )
{
Vec3d T_1, T_2; // vectors to hold the T-values for every direction
double t_near = -DBL_MAX; // maximums defined in float.h
double t_far = DBL_MAX;
//vmin = center - Vec3d(xExtent);
//vmax = center + Vec3d(xExtent);
for (int i = 0; i < 3; i++){ //we test slabs in every direction
if (ray.direction[i] == 0){ // ray parallel to planes in this direction
if ((ray.origin[i] < vmin[i]) || (ray.origin[i] > vmax[i])) {
return false; // parallel AND outside box : no intersection possible
}
} else { // ray not parallel to planes in this direction
T_1[i] = (vmin[i] - ray.origin[i]) / ray.direction[i];
T_2[i] = (vmax[i] - ray.origin[i]) / ray.direction[i];
if(T_1[i] > T_2[i]){ // we want T_1 to hold values for intersection with near plane
std::swap(T_1,T_2);
}
if (T_1[i] > t_near){
t_near = T_1[i];
}
if (T_2[i] < t_far){
t_far = T_2[i];
}
if( (t_near > t_far) || (t_far < 0) ){
return false;
}
}
}
double tnear = t_near, tfar = t_far; // put return values in place
Q_UNUSED(tnear);
Q_UNUSED(tfar);
return true; // if we made it here, there was an intersection - YAY
}
bool intersectsOld( const Ray& ray ) const
{
double rhs;
double fWdU[3];
double fAWdU[3];
double fDdU[3];
double fADdU[3];
double fAWxDdU[3];
Vec3d UNIT_X(1.0, 0.0, 0.0);
Vec3d UNIT_Y(0.0, 1.0, 0.0);
Vec3d UNIT_Z(0.0, 0.0, 1.0);
Vec3d diff = ray.origin - center;
Vec3d wCrossD = cross(ray.direction , diff);
fWdU[0] = dot(ray.direction , UNIT_X);
fAWdU[0] = abs(fWdU[0]);
fDdU[0] = dot(diff , UNIT_X);
fADdU[0] = abs(fDdU[0]);
if (fADdU[0] > xExtent && fDdU[0] * fWdU[0] >= 0.0) return false;
fWdU[1] = dot(ray.direction , UNIT_Y);
fAWdU[1] = abs(fWdU[1]);
fDdU[1] = dot(diff , UNIT_Y);
fADdU[1] = abs(fDdU[1]);
if (fADdU[1] > yExtent && fDdU[1] * fWdU[1] >= 0.0) return false;
fWdU[2] = dot(ray.direction , UNIT_Z);
fAWdU[2] = abs(fWdU[2]);
fDdU[2] = dot(diff , UNIT_Z);
fADdU[2] = abs(fDdU[2]);
if (fADdU[2] > zExtent && fDdU[2] * fWdU[2] >= 0.0) return false;
fAWxDdU[0] = abs(dot(wCrossD , UNIT_X));
rhs = yExtent , fAWdU[2] + zExtent * fAWdU[1];
if (fAWxDdU[0] > rhs) return false;
fAWxDdU[1] = abs(dot(wCrossD , UNIT_Y));
rhs = xExtent * fAWdU[2] + zExtent * fAWdU[0];
if (fAWxDdU[1] > rhs) return false;
fAWxDdU[2] = abs(dot(wCrossD , UNIT_Z));
rhs = xExtent * fAWdU[1] + yExtent * fAWdU[0];
if (fAWxDdU[2] > rhs) return false;
return true;
}
/* AABB-triangle overlap test code */
/* by Tomas Akenine-Möller */
bool containsTriangle( const Vec3d& tv0, const Vec3d& tv1, const Vec3d& tv2 ) const
{
Vec3d boxcenter(center);
Vec3d boxhalfsize(xExtent, yExtent, zExtent);
int X = 0, Y = 1, Z = 2;
/* use separating axis theorem to test overlap between triangle and box */
/* need to test for overlap in these directions: */
/* 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle */
/* we do not even need to test these) */
/* 2) normal of the triangle */
/* 3) crossproduct(edge from tri, {x,y,z}-directin) */
/* this gives 3x3=9 more tests */
Vec3d v0,v1,v2;
double min,max,p0,p1,p2,rad,fex,fey,fez;
Vec3d normal,e0,e1,e2;
/* This is the fastest branch on Sun */
/* move everything so that the box center is in (0,0,0) */
v0=tv0-boxcenter;
v1=tv1-boxcenter;
v2=tv2-boxcenter;
/* compute triangle edges */
e0=v1-v0; /* tri edge 0 */
e1=v2-v1; /* tri edge 1 */
e2=v0-v2; /* tri edge 2 */
/* Bullet 3: */
/* test the 9 tests first (this was faster) */
fex = fabsf(e0[X]);
fey = fabsf(e0[Y]);
fez = fabsf(e0[Z]);
AXISTEST_X01(e0[Z], e0[Y], fez, fey);
AXISTEST_Y02(e0[Z], e0[X], fez, fex);
AXISTEST_Z12(e0[Y], e0[X], fey, fex);
fex = fabsf(e1[X]);
fey = fabsf(e1[Y]);
fez = fabsf(e1[Z]);
AXISTEST_X01(e1[Z], e1[Y], fez, fey);
AXISTEST_Y02(e1[Z], e1[X], fez, fex);
AXISTEST_Z0(e1[Y], e1[X], fey, fex);
fex = fabsf(e2[X]);
fey = fabsf(e2[Y]);
fez = fabsf(e2[Z]);
AXISTEST_X2(e2[Z], e2[Y], fez, fey);
AXISTEST_Y1(e2[Z], e2[X], fez, fex);
AXISTEST_Z12(e2[Y], e2[X], fey, fex);
/* Bullet 1: */
/* first test overlap in the {x,y,z}-directions */
/* find min, max of the triangle each direction, and test for overlap in */
/* that direction -- this is equivalent to testing a minimal AABB around */
/* the triangle against the AABB */
/* test in X-direction */
FINDMINMAX(v0[X],v1[X],v2[X],min,max);
if(min>boxhalfsize[X] || max<-boxhalfsize[X]) return 0;
/* test in Y-direction */
FINDMINMAX(v0[Y],v1[Y],v2[Y],min,max);
if(min>boxhalfsize[Y] || max<-boxhalfsize[Y]) return 0;
/* test in Z-direction */
FINDMINMAX(v0[Z],v1[Z],v2[Z],min,max);
if(min>boxhalfsize[Z] || max<-boxhalfsize[Z]) return 0;
/* Bullet 2: */
/* test if the box intersects the plane of the triangle */
/* compute plane equation of triangle: normal*x+d=0 */
normal = cross(e0, e1);
if(!planeBoxOverlap(normal,v0,boxhalfsize)) return 0;
return 1; /* box and triangle overlaps */
}
bool intersectsBoundingBox( const BoundingBox& bb ) const
{
if (center.x() + xExtent < bb.center.x() - bb.xExtent || center.x() - xExtent > bb.center.x() + bb.xExtent)
return false;
else if (center.y() + yExtent < bb.center.y() - bb.yExtent || center.y() - yExtent > bb.center.y() + bb.yExtent)
return false;
else if (center.z() + zExtent < bb.center.z() - bb.zExtent || center.z() - zExtent > bb.center.z() + bb.zExtent)
return false;
else
return true;
}
bool intersectsSphere( const Vec3d& sphere_center, double radius )
{
if (abs(center.x() - sphere_center.x()) < radius + xExtent
&& abs(center.y() - sphere_center.y()) < radius + yExtent
&& abs(center.z() - sphere_center.z()) < radius + zExtent)
return true;
return false;
}
bool contains( const Vec3d& point ) const
{
return abs(center.x() - point.x()) < xExtent
&& abs(center.y() - point.y()) < yExtent
&& abs(center.z() - point.z()) < zExtent;
}
Vec3d Center()
{
return center;
}
double Diag()
{
return (vmin - vmax).norm();
}
};
struct HitResult{
bool hit;
double distance;
double u;
double v;
int index;
HitResult(bool isHit = false, double hitDistance = DBL_MAX) : hit(isHit), distance(hitDistance)
{
u = -1.0;
v = -1.0;
index = -1;
}
HitResult(const HitResult& other)
{
this->hit = other.hit;
this->distance = other.distance;
this->u = other.u;
this->v = other.v;
this->index = other.index;
}
HitResult& operator= (const HitResult& other)
{
this->hit = other.hit;
this->distance = other.distance;
this->u = other.u;
this->v = other.v;
this->index = other.index;
return *this;
}
};