C++ demos: prepare infrastructure, add spheres_shadow.cc
rename Ray::a to Ray::o
KdNode::shapes changed to pointer and added to union together with *children (to save memory)
/*
* C++ RayTracer
* file: scene.cc
*
* Radek Brich, 2006
*/
#include <math.h>
#include "common.h"
#include "scene.h"
/* 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)
{
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;
if (t1 > 0)
{
if (t1 < dist) {
dist = t1;
return true;
}
} else {
float t2 = Vd + Det;
if (t2 > 0)
{
// ray from inside of the sphere
dist = t2;
return true;
}
}
}
return false;
}
bool Sphere::intersect_all(const Ray &ray, float dist, vector<float> &allts)
{
//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()
{
BBox bbox = BBox();
bbox.L = center - radius;
bbox.H = center + radius;
return bbox;
}
bool Box::intersect(const Ray &ray, float &dist)
{
float b;
return get_bbox().intersect(ray, dist, b);
}
Vector3 Box::normal(Vector3 &P)
{
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;
}
// this initialization and following intersection methods implements
// Fast Triangle Intersection algorithm from
// http://www.mpi-inf.mpg.de/~wald/PhD/
Triangle::Triangle(const Vector3 &aA, const Vector3 &aB, const Vector3 &aC, Material *amaterial)
: A(aA), B(aB), C(aC)
{
material = amaterial;
Vector3 c = B - A;
Vector3 b = C - A;
N = cross(c, b);
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;
}
int u = (k + 1) % 3;
int v = (k + 2) % 3;
float krec = 1.0f / 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[v] * reci;
cnv = -c[u] * reci;
// finalize normal
N.normalize();
}
// see comment for previous method
bool Triangle::intersect(const Ray &ray, float &dist)
{
Vector3 O = ray.o;
Vector3 D = ray.dir;
const int modulo3[5] = {0,1,2,0,1};
const int ku = modulo3[k+1];
const int kv = modulo3[k+2];
const float lnd = 1.0f / (D[k] + nu * D[ku] + nv * D[kv]);
const float t = (nd - O[k] - nu * O[ku] - nv * O[kv]) * lnd;
if (!(t < dist && t > 0))
return false;
float hu = O[ku] + t * D[ku] - A[ku];
float hv = O[kv] + t * D[kv] - A[kv];
float beta = hv * bnu + hu * bnv;
if (beta < 0)
return false;
float gamma = hu * cnu + hv * cnv;
if (gamma < 0)
return false;
if ((beta + gamma) > 1)
return false;
dist = t;
return true;
}
BBox Triangle::get_bbox()
{
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;
};