tuned ray-triangle intersection, now there are three algorithms to choose from:
Plucker, Barycentric and Barycentric with preprocessing (Wald)
methods in Vector and Shape (and derivates) made const
/* * C++ RayTracer * file: scene.cc * * Radek Brich, 2006 */#include <math.h>#include "common.h"#include "scene.h"void Camera::rotate(const Quaternion &q){ /* //non-optimized Quaternion res; res = q * Quaternion(u) * conjugate(q); u = res.toVector(); res = q * Quaternion(v) * conjugate(q); v = res.toVector(); res = q * Quaternion(p) * conjugate(q); p = res.toVector(); */ // optimized Float t2 = q.a*q.b; Float t3 = q.a*q.c; Float t4 = q.a*q.d; Float t5 = -q.b*q.b; Float t6 = q.b*q.c; Float t7 = q.b*q.d; Float t8 = -q.c*q.c; Float t9 = q.c*q.d; Float t10 = -q.d*q.d; Float x,y,z; x = 2*( (t8 + t10)*p.x + (t6 - t4)*p.y + (t3 + t7)*p.z ) + p.x; y = 2*( (t4 + t6)*p.x + (t5 + t10)*p.y + (t9 - t2)*p.z ) + p.y; z = 2*( (t7 - t3)*p.x + (t2 + t9)*p.y + (t5 + t8)*p.z ) + p.z; p = Vector3(x,y,z); x = 2*( (t8 + t10)*u.x + (t6 - t4)*u.y + (t3 + t7)*u.z ) + u.x; y = 2*( (t4 + t6)*u.x + (t5 + t10)*u.y + (t9 - t2)*u.z ) + u.y; z = 2*( (t7 - t3)*u.x + (t2 + t9)*u.y + (t5 + t8)*u.z ) + u.z; u = Vector3(x,y,z); x = 2*( (t8 + t10)*v.x + (t6 - t4)*v.y + (t3 + t7)*v.z ) + v.x; y = 2*( (t4 + t6)*v.x + (t5 + t10)*v.y + (t9 - t2)*v.z ) + v.y; z = 2*( (t7 - t3)*v.x + (t2 + t9)*v.y + (t5 + t8)*v.z ) + v.z; v = Vector3(x,y,z); p.normalize(); u.normalize(); v.normalize();}void Camera::move(const Float fw, const Float left, const Float up){ eye = eye + fw*p + left*u + up*v;}/* http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter3.htm */bool BBox::intersect(const Ray &ray, Float &a, Float &b){ Float tnear = -FLT_MAX; Float tfar = FLT_MAX; Float t1, t2; for (int i = 0; i < 3; i++) { if (ray.dir.cell[i] == 0) { /* ray is parallel to these planes */ if (ray.o.cell[i] < L.cell[i] || ray.o.cell[i] > H.cell[i]) return false; } else { /* compute the intersection distance of the planes */ t1 = (L.cell[i] - ray.o.cell[i]) / ray.dir.cell[i]; t2 = (H.cell[i] - ray.o.cell[i]) / ray.dir.cell[i]; if (t1 > t2) swap(t1, t2); if (t1 > tnear) tnear = t1; /* want largest Tnear */ if (t2 < tfar) tfar = t2; /* want smallest Tfar */ if (tnear > tfar) return false; /* box missed */ if (tfar < 0) return false; /* box is behind ray */ } } a = tnear; b = tfar; return true;}bool Sphere::intersect(const Ray &ray, Float &dist) const{ Vector3 V = ray.o - center; register Float d = -dot(V, ray.dir); register Float Det = d * d - (dot(V,V) - sqr_radius); if (Det > 0) { d -= sqrtf(Det); if (d > 0 && d < dist) { dist = d; return true; } } return false;}bool Sphere::intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const{ //allts = new vector<Float>(); Vector3 V = ((Ray)ray).o - center; Float Vd = - dot(V, ray.dir); Float Det = Vd * Vd - (dot(V,V) - sqr_radius); if (Det > 0) { Det = sqrtf(Det); Float t1 = Vd - Det; Float t2 = Vd + Det; if (t1 < 0) { if (t2 > 0) { // ray from inside of the sphere allts.push_back(0.0); allts.push_back(t2); return true; } else return false; } else { allts.push_back(t1); allts.push_back(t2); return true; } } return false;}BBox Sphere::get_bbox() const{ BBox bbox = BBox(); bbox.L = center - radius; bbox.H = center + radius; return bbox;}bool Box::intersect(const Ray &ray, Float &dist) const{ Float a,b; bool res = get_bbox().intersect(ray, a, b); if (res && a < dist) { dist = a; return true; } else return false;}Vector3 Box::normal(Vector3 &P) const{ Vector3 N; for (int i = 0; i < 3; i++) { if (P.cell[i] >= L.cell[i]-Eps && P.cell[i] <= L.cell[i]+Eps) //if (P.cell[i] == L.cell[i]) { N.cell[i] = -1.0; break; } if (P.cell[i] >= H.cell[i]-Eps && P.cell[i] <= H.cell[i]+Eps) //if (P.cell[i] == H.cell[i]) { N.cell[i] = +1.0; break; } } return N;}#ifdef TRI_PLUCKERinline void Plucker(const Vector3 &p, const Vector3 &q, Float* pl){ pl[0] = p.x*q.y - q.x*p.y; pl[1] = p.x*q.z - q.x*p.z; pl[2] = p.x - q.x; pl[3] = p.y*q.z - q.y*p.z; pl[4] = p.z - q.z; pl[5] = q.y - p.y;}inline Float Side(const Float* pla, const Float* plb){ return pla[0]*plb[4] + pla[1]*plb[5] + pla[2]*plb[3] + pla[4]*plb[0] + pla[5]*plb[1] + pla[3]*plb[2];}#endifTriangle::Triangle(const Vector3 &aA, const Vector3 &aB, const Vector3 &aC, Material *amaterial) : A(aA), B(aB), C(aC){ material = amaterial; material->reflection = 0; const Vector3 c = B - A; const Vector3 b = C - A; N = -cross(c, b); N.normalize();#ifdef TRI_PLUCKER Plucker(B,C,pla); Plucker(C,A,plb); Plucker(A,B,plc);#endif#if defined(TRI_BARI) || defined(TRI_BARI_PRE) if (fabsf(N.x) > fabsf(N.y)) { if (fabsf(N.x) > fabsf(N.z)) k = 0; else k = 2; } else { if (fabsf(N.y) > fabsf(N.z)) k = 1; else k = 2; }#endif#ifdef TRI_BARI_PRE int u = (k + 1) % 3; int v = (k + 2) % 3; Float krec = 1.0 / N[k]; nu = N[u] * krec; nv = N[v] * krec; nd = dot(N, A) * krec; // first line equation Float reci = 1.0f / (b[u] * c[v] - b[v] * c[u]); bnu = b[u] * reci; bnv = -b[v] * reci; // second line equation cnu = -c[u] * reci; cnv = c[v] * reci;#endif}bool Triangle::intersect(const Ray &ray, Float &dist) const{#ifdef TRI_PLUCKER Float plr[6]; Plucker(ray.o, ray.o+ray.dir, plr); const bool side0 = Side(plr, pla) >= 0.0; const bool side1 = Side(plr, plb) >= 0.0; if (side0 != side1) return false; const bool side2 = Side(plr, plc) >= 0.0; if (side0 != side2) return false; const Float t = - dot( (ray.o-A), N) / dot(ray.dir,N); if(t <= Eps || t >= dist) return false; dist = t; return true;#endif#if defined(TRI_BARI) || defined(TRI_BARI_PRE) static const int modulo3[5] = {0,1,2,0,1}; const Vector3 &O = ray.o; const Vector3 &D = ray.dir; register const int u = modulo3[k+1]; register const int v = modulo3[k+2];#endif#ifdef TRI_BARI_PRE const Float t = (nd - O[k] - nu * O[u] - nv * O[v]) / (D[k] + nu * D[u] + nv * D[v]); if (t >= dist || t < Eps) return false; const Float hu = O[u] + t * D[u] - A[u]; const Float hv = O[v] + t * D[v] - A[v]; const Float beta = hv * bnu + hu * bnv; if (beta < 0.) return false; const Float gamma = hu * cnv + hv * cnu; if (gamma < 0. || beta + gamma > 1) return false; dist = t; return true;#endif#ifdef TRI_BARI // original barycentric coordinates based intesection // not optimized, just for reference const Vector3 c = B - A; const Vector3 b = C - A; // distance test const Float t = - dot( (O-A), N) / dot(D,N); if (t < Eps || t > dist) return false; // calc hitpoint const Float Hu = O[u] + t * D[u] - A[u]; const Float Hv = O[v] + t * D[v] - A[v]; const Float beta = (b[u] * Hv - b[v] * Hu) / (b[u] * c[v] - b[v] * c[u]); if (beta < 0) return false; const Float gamma = (c[v] * Hu - c[u] * Hv) / (b[u] * c[v] - b[v] * c[u]); if (gamma < 0) return false; if (beta+gamma > 1) return false; dist = t; return true;#endif}BBox Triangle::get_bbox() const{ BBox bbox = BBox(); bbox.L = A; if (B.x < bbox.L.x) bbox.L.x = B.x; if (C.x < bbox.L.x) bbox.L.x = C.x; if (B.y < bbox.L.y) bbox.L.y = B.y; if (C.y < bbox.L.y) bbox.L.y = C.y; if (B.z < bbox.L.z) bbox.L.z = B.z; if (C.z < bbox.L.z) bbox.L.z = C.z; bbox.H = A; if (B.x > bbox.H.x) bbox.H.x = B.x; if (C.x > bbox.H.x) bbox.H.x = C.x; if (B.y > bbox.H.y) bbox.H.y = B.y; if (C.y > bbox.H.y) bbox.H.y = C.y; if (B.z > bbox.H.z) bbox.H.z = B.z; if (C.z > bbox.H.z) bbox.H.z = C.z; return bbox;};