489 lines
12 KiB
C
489 lines
12 KiB
C
/*
|
|
Copyright (c) 2009 David Bucciarelli (davibu@interfree.it)
|
|
|
|
Permission is hereby granted, free of charge, to any person obtaining
|
|
a copy of this software and associated documentation files (the
|
|
"Software"), to deal in the Software without restriction, including
|
|
without limitation the rights to use, copy, modify, merge, publish,
|
|
distribute, sublicense, and/or sell copies of the Software, and to
|
|
permit persons to whom the Software is furnished to do so, subject to
|
|
the following conditions:
|
|
|
|
The above copyright notice and this permission notice shall be included
|
|
in all copies or substantial portions of the Software.
|
|
|
|
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
|
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
|
|
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
|
|
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
|
|
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
|
|
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
|
|
*/
|
|
|
|
#ifndef _GEOMFUNC_H
|
|
#define _GEOMFUNC_H
|
|
|
|
#include "geom.h"
|
|
#include "simplernd.h"
|
|
|
|
#ifndef SMALLPT_GPU
|
|
|
|
static float SphereIntersect(
|
|
#ifdef GPU_KERNEL
|
|
OCL_CONSTANT_BUFFER
|
|
#endif
|
|
const Sphere *s,
|
|
const Ray *r) { /* returns distance, 0 if nohit */
|
|
Vec op; /* Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0 */
|
|
vsub(op, s->p, r->o);
|
|
|
|
float b = vdot(op, r->d);
|
|
float det = b * b - vdot(op, op) + s->rad * s->rad;
|
|
if (det < 0.f)
|
|
return 0.f;
|
|
else
|
|
det = sqrt(det);
|
|
|
|
float t = b - det;
|
|
if (t > EPSILON)
|
|
return t;
|
|
else {
|
|
t = b + det;
|
|
|
|
if (t > EPSILON)
|
|
return t;
|
|
else
|
|
return 0.f;
|
|
}
|
|
}
|
|
|
|
static void UniformSampleSphere(const float u1, const float u2, Vec *v) {
|
|
const float zz = 1.f - 2.f * u1;
|
|
const float r = sqrt(max(0.f, 1.f - zz * zz));
|
|
const float phi = 2.f * FLOAT_PI * u2;
|
|
const float xx = r * cos(phi);
|
|
const float yy = r * sin(phi);
|
|
|
|
vinit(*v, xx, yy, zz);
|
|
}
|
|
|
|
static int Intersect(
|
|
#ifdef GPU_KERNEL
|
|
OCL_CONSTANT_BUFFER
|
|
#endif
|
|
const Sphere *spheres,
|
|
const unsigned int sphereCount,
|
|
const Ray *r,
|
|
float *t,
|
|
unsigned int *id) {
|
|
float inf = (*t) = 1e20f;
|
|
|
|
unsigned int i = sphereCount;
|
|
for (; i--;) {
|
|
const float d = SphereIntersect(&spheres[i], r);
|
|
if ((d != 0.f) && (d < *t)) {
|
|
*t = d;
|
|
*id = i;
|
|
}
|
|
}
|
|
|
|
return (*t < inf);
|
|
}
|
|
|
|
static int IntersectP(
|
|
#ifdef GPU_KERNEL
|
|
OCL_CONSTANT_BUFFER
|
|
#endif
|
|
const Sphere *spheres,
|
|
const unsigned int sphereCount,
|
|
const Ray *r,
|
|
const float maxt) {
|
|
unsigned int i = sphereCount;
|
|
for (; i--;) {
|
|
const float d = SphereIntersect(&spheres[i], r);
|
|
if ((d != 0.f) && (d < maxt))
|
|
return 1;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
static void SampleLights(
|
|
#ifdef GPU_KERNEL
|
|
OCL_CONSTANT_BUFFER
|
|
#endif
|
|
const Sphere *spheres,
|
|
const unsigned int sphereCount,
|
|
unsigned int *seed0, unsigned int *seed1,
|
|
const Vec *hitPoint,
|
|
const Vec *normal,
|
|
Vec *result) {
|
|
vclr(*result);
|
|
|
|
/* For each light */
|
|
unsigned int i;
|
|
for (i = 0; i < sphereCount; i++) {
|
|
#ifdef GPU_KERNEL
|
|
OCL_CONSTANT_BUFFER
|
|
#endif
|
|
const Sphere *light = &spheres[i];
|
|
if (!viszero(light->e)) {
|
|
/* It is a light source */
|
|
Ray shadowRay;
|
|
shadowRay.o = *hitPoint;
|
|
|
|
/* Choose a point over the light source */
|
|
Vec unitSpherePoint;
|
|
UniformSampleSphere(GetRandom(seed0, seed1), GetRandom(seed0, seed1), &unitSpherePoint);
|
|
Vec spherePoint;
|
|
vsmul(spherePoint, light->rad, unitSpherePoint);
|
|
vadd(spherePoint, spherePoint, light->p);
|
|
|
|
/* Build the shadow ray direction */
|
|
vsub(shadowRay.d, spherePoint, *hitPoint);
|
|
const float len = sqrt(vdot(shadowRay.d, shadowRay.d));
|
|
vsmul(shadowRay.d, 1.f / len, shadowRay.d);
|
|
|
|
float wo = vdot(shadowRay.d, unitSpherePoint);
|
|
if (wo > 0.f) {
|
|
/* It is on the other half of the sphere */
|
|
continue;
|
|
} else
|
|
wo = -wo;
|
|
|
|
/* Check if the light is visible */
|
|
const float wi = vdot(shadowRay.d, *normal);
|
|
if ((wi > 0.f) && (!IntersectP(spheres, sphereCount, &shadowRay, len - EPSILON))) {
|
|
Vec c; vassign(c, light->e);
|
|
const float s = (4.f * FLOAT_PI * light->rad * light->rad) * wi * wo / (len *len);
|
|
vsmul(c, s, c);
|
|
vadd(*result, *result, c);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
static void RadiancePathTracing(
|
|
#ifdef GPU_KERNEL
|
|
OCL_CONSTANT_BUFFER
|
|
#endif
|
|
const Sphere *spheres,
|
|
const unsigned int sphereCount,
|
|
const Ray *startRay,
|
|
unsigned int *seed0, unsigned int *seed1,
|
|
Vec *result) {
|
|
Ray currentRay; rassign(currentRay, *startRay);
|
|
Vec rad; vinit(rad, 0.f, 0.f, 0.f);
|
|
Vec throughput; vinit(throughput, 1.f, 1.f, 1.f);
|
|
|
|
unsigned int depth = 0;
|
|
int specularBounce = 1;
|
|
for (;; ++depth) {
|
|
// Removed Russian Roulette in order to improve execution on SIMT
|
|
if (depth > 6) {
|
|
*result = rad;
|
|
return;
|
|
}
|
|
|
|
float t; /* distance to intersection */
|
|
unsigned int id = 0; /* id of intersected object */
|
|
if (!Intersect(spheres, sphereCount, ¤tRay, &t, &id)) {
|
|
*result = rad; /* if miss, return */
|
|
return;
|
|
}
|
|
|
|
#ifdef GPU_KERNEL
|
|
OCL_CONSTANT_BUFFER
|
|
#endif
|
|
const Sphere *obj = &spheres[id]; /* the hit object */
|
|
|
|
Vec hitPoint;
|
|
vsmul(hitPoint, t, currentRay.d);
|
|
vadd(hitPoint, currentRay.o, hitPoint);
|
|
|
|
Vec normal;
|
|
vsub(normal, hitPoint, obj->p);
|
|
vnorm(normal);
|
|
|
|
const float dp = vdot(normal, currentRay.d);
|
|
|
|
Vec nl;
|
|
// SIMT optimization
|
|
const float invSignDP = -1.f * sign(dp);
|
|
vsmul(nl, invSignDP, normal);
|
|
|
|
/* Add emitted light */
|
|
Vec eCol; vassign(eCol, obj->e);
|
|
if (!viszero(eCol)) {
|
|
if (specularBounce) {
|
|
vsmul(eCol, fabs(dp), eCol);
|
|
vmul(eCol, throughput, eCol);
|
|
vadd(rad, rad, eCol);
|
|
}
|
|
|
|
*result = rad;
|
|
return;
|
|
}
|
|
|
|
if (obj->refl == DIFF) { /* Ideal DIFFUSE reflection */
|
|
specularBounce = 0;
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
/* Direct lighting component */
|
|
|
|
Vec Ld;
|
|
SampleLights(spheres, sphereCount, seed0, seed1, &hitPoint, &nl, &Ld);
|
|
vmul(Ld, throughput, Ld);
|
|
vadd(rad, rad, Ld);
|
|
|
|
/* Diffuse component */
|
|
|
|
float r1 = 2.f * FLOAT_PI * GetRandom(seed0, seed1);
|
|
float r2 = GetRandom(seed0, seed1);
|
|
float r2s = sqrt(r2);
|
|
|
|
Vec w; vassign(w, nl);
|
|
|
|
Vec u, a;
|
|
if (fabs(w.x) > .1f) {
|
|
vinit(a, 0.f, 1.f, 0.f);
|
|
} else {
|
|
vinit(a, 1.f, 0.f, 0.f);
|
|
}
|
|
vxcross(u, a, w);
|
|
vnorm(u);
|
|
|
|
Vec v;
|
|
vxcross(v, w, u);
|
|
|
|
Vec newDir;
|
|
vsmul(u, cos(r1) * r2s, u);
|
|
vsmul(v, sin(r1) * r2s, v);
|
|
vadd(newDir, u, v);
|
|
vsmul(w, sqrt(1 - r2), w);
|
|
vadd(newDir, newDir, w);
|
|
|
|
currentRay.o = hitPoint;
|
|
currentRay.d = newDir;
|
|
continue;
|
|
} else if (obj->refl == SPEC) { /* Ideal SPECULAR reflection */
|
|
specularBounce = 1;
|
|
|
|
Vec newDir;
|
|
vsmul(newDir, 2.f * vdot(normal, currentRay.d), normal);
|
|
vsub(newDir, currentRay.d, newDir);
|
|
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
rinit(currentRay, hitPoint, newDir);
|
|
continue;
|
|
} else {
|
|
specularBounce = 1;
|
|
|
|
Vec newDir;
|
|
vsmul(newDir, 2.f * vdot(normal, currentRay.d), normal);
|
|
vsub(newDir, currentRay.d, newDir);
|
|
|
|
Ray reflRay; rinit(reflRay, hitPoint, newDir); /* Ideal dielectric REFRACTION */
|
|
int into = (vdot(normal, nl) > 0); /* Ray from outside going in? */
|
|
|
|
float nc = 1.f;
|
|
float nt = 1.5f;
|
|
float nnt = into ? nc / nt : nt / nc;
|
|
float ddn = vdot(currentRay.d, nl);
|
|
float cos2t = 1.f - nnt * nnt * (1.f - ddn * ddn);
|
|
|
|
if (cos2t < 0.f) { /* Total internal reflection */
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
rassign(currentRay, reflRay);
|
|
continue;
|
|
}
|
|
|
|
float kk = (into ? 1 : -1) * (ddn * nnt + sqrt(cos2t));
|
|
Vec nkk;
|
|
vsmul(nkk, kk, normal);
|
|
Vec transDir;
|
|
vsmul(transDir, nnt, currentRay.d);
|
|
vsub(transDir, transDir, nkk);
|
|
vnorm(transDir);
|
|
|
|
float a = nt - nc;
|
|
float b = nt + nc;
|
|
float R0 = a * a / (b * b);
|
|
float c = 1 - (into ? -ddn : vdot(transDir, normal));
|
|
|
|
float Re = R0 + (1 - R0) * c * c * c * c*c;
|
|
float Tr = 1.f - Re;
|
|
float P = .25f + .5f * Re;
|
|
float RP = Re / P;
|
|
float TP = Tr / (1.f - P);
|
|
|
|
if (GetRandom(seed0, seed1) < P) { /* R.R. */
|
|
vsmul(throughput, RP, throughput);
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
rassign(currentRay, reflRay);
|
|
continue;
|
|
} else {
|
|
vsmul(throughput, TP, throughput);
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
rinit(currentRay, hitPoint, transDir);
|
|
continue;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
static void RadianceDirectLighting(
|
|
#ifdef GPU_KERNEL
|
|
OCL_CONSTANT_BUFFER
|
|
#endif
|
|
const Sphere *spheres,
|
|
const unsigned int sphereCount,
|
|
const Ray *startRay,
|
|
unsigned int *seed0, unsigned int *seed1,
|
|
Vec *result) {
|
|
Ray currentRay; rassign(currentRay, *startRay);
|
|
Vec rad; vinit(rad, 0.f, 0.f, 0.f);
|
|
Vec throughput; vinit(throughput, 1.f, 1.f, 1.f);
|
|
|
|
unsigned int depth = 0;
|
|
int specularBounce = 1;
|
|
for (;; ++depth) {
|
|
// Removed Russian Roulette in order to improve execution on SIMT
|
|
if (depth > 6) {
|
|
*result = rad;
|
|
return;
|
|
}
|
|
|
|
float t; /* distance to intersection */
|
|
unsigned int id = 0; /* id of intersected object */
|
|
if (!Intersect(spheres, sphereCount, ¤tRay, &t, &id)) {
|
|
*result = rad; /* if miss, return */
|
|
return;
|
|
}
|
|
|
|
#ifdef GPU_KERNEL
|
|
OCL_CONSTANT_BUFFER
|
|
#endif
|
|
const Sphere *obj = &spheres[id]; /* the hit object */
|
|
|
|
Vec hitPoint;
|
|
vsmul(hitPoint, t, currentRay.d);
|
|
vadd(hitPoint, currentRay.o, hitPoint);
|
|
|
|
Vec normal;
|
|
vsub(normal, hitPoint, obj->p);
|
|
vnorm(normal);
|
|
|
|
const float dp = vdot(normal, currentRay.d);
|
|
|
|
Vec nl;
|
|
// SIMT optimization
|
|
const float invSignDP = -1.f * sign(dp);
|
|
vsmul(nl, invSignDP, normal);
|
|
|
|
/* Add emitted light */
|
|
Vec eCol; vassign(eCol, obj->e);
|
|
if (!viszero(eCol)) {
|
|
if (specularBounce) {
|
|
vsmul(eCol, fabs(dp), eCol);
|
|
vmul(eCol, throughput, eCol);
|
|
vadd(rad, rad, eCol);
|
|
}
|
|
|
|
*result = rad;
|
|
return;
|
|
}
|
|
|
|
if (obj->refl == DIFF) { /* Ideal DIFFUSE reflection */
|
|
specularBounce = 0;
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
/* Direct lighting component */
|
|
|
|
Vec Ld;
|
|
SampleLights(spheres, sphereCount, seed0, seed1, &hitPoint, &nl, &Ld);
|
|
vmul(Ld, throughput, Ld);
|
|
vadd(rad, rad, Ld);
|
|
|
|
*result = rad;
|
|
return;
|
|
} else if (obj->refl == SPEC) { /* Ideal SPECULAR reflection */
|
|
specularBounce = 1;
|
|
|
|
Vec newDir;
|
|
vsmul(newDir, 2.f * vdot(normal, currentRay.d), normal);
|
|
vsub(newDir, currentRay.d, newDir);
|
|
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
rinit(currentRay, hitPoint, newDir);
|
|
continue;
|
|
} else {
|
|
specularBounce = 1;
|
|
|
|
Vec newDir;
|
|
vsmul(newDir, 2.f * vdot(normal, currentRay.d), normal);
|
|
vsub(newDir, currentRay.d, newDir);
|
|
|
|
Ray reflRay; rinit(reflRay, hitPoint, newDir); /* Ideal dielectric REFRACTION */
|
|
int into = (vdot(normal, nl) > 0); /* Ray from outside going in? */
|
|
|
|
float nc = 1.f;
|
|
float nt = 1.5f;
|
|
float nnt = into ? nc / nt : nt / nc;
|
|
float ddn = vdot(currentRay.d, nl);
|
|
float cos2t = 1.f - nnt * nnt * (1.f - ddn * ddn);
|
|
|
|
if (cos2t < 0.f) { /* Total internal reflection */
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
rassign(currentRay, reflRay);
|
|
continue;
|
|
}
|
|
|
|
float kk = (into ? 1 : -1) * (ddn * nnt + sqrt(cos2t));
|
|
Vec nkk;
|
|
vsmul(nkk, kk, normal);
|
|
Vec transDir;
|
|
vsmul(transDir, nnt, currentRay.d);
|
|
vsub(transDir, transDir, nkk);
|
|
vnorm(transDir);
|
|
|
|
float a = nt - nc;
|
|
float b = nt + nc;
|
|
float R0 = a * a / (b * b);
|
|
float c = 1 - (into ? -ddn : vdot(transDir, normal));
|
|
|
|
float Re = R0 + (1 - R0) * c * c * c * c*c;
|
|
float Tr = 1.f - Re;
|
|
float P = .25f + .5f * Re;
|
|
float RP = Re / P;
|
|
float TP = Tr / (1.f - P);
|
|
|
|
if (GetRandom(seed0, seed1) < P) { /* R.R. */
|
|
vsmul(throughput, RP, throughput);
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
rassign(currentRay, reflRay);
|
|
continue;
|
|
} else {
|
|
vsmul(throughput, TP, throughput);
|
|
vmul(throughput, throughput, obj->c);
|
|
|
|
rinit(currentRay, hitPoint, transDir);
|
|
continue;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
#endif
|
|
|
|
#endif /* _GEOMFUNC_H */
|
|
|