/*
* raytracer.cc: Raytracer class
*
* This file is part of Pyrit Ray Tracer.
*
* Copyright 2006, 2007, 2008 Radek Brich
*
* 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.
*/
#include <pthread.h>
#include <stdio.h>
#include <malloc.h>
#include "raytracer.h"
int pyrit_verbosity = 2;
// Hammersley spherical point distribution
// http://www.cse.cuhk.edu.hk/~ttwong/papers/udpoint/udpoints.html
Vector3 Raytracer::SphereDistribute(int i, int n, Float extent, Vector3 &normal)
{
Float p, t, st, phi, phirad;
int kk;
t = 0;
for (p=0.5, kk=i; kk; p*=0.5, kk>>=1)
if (kk & 1)
t += p;
t = 1.0 + (t - 1.0)*extent;
phi = (i + 0.5) / n;
phirad = phi * 2.0 * M_PI;
st = sqrt(1.0 - t*t);
Float x, y, z, xx, yy, zz, q;
x = st * cos(phirad);
y = st * sin(phirad);
z = t;
// rotate against Y axis
q = acos(normal.z);
zz = z*cos(q) - x*sin(q);
xx = z*sin(q) + x*cos(q);
yy = y;
// rotate against Z axis
q = atan2f(normal.y, normal.x);
x = xx*cos(q) - yy*sin(q);
y = xx*sin(q) + yy*cos(q);
z = zz;
return Vector3(x, y, z);
}
// ---- tyto dve funkce budou v budouci verzi metody objektu PhongShader
// calculate shader function
// P is point of intersection, N normal in this point
Colour PhongShader_ambient(Material &mat, Vector3 &P)
{
Colour col;
if (mat.texture)
col = mat.texture->evaluate(P);
else
col = mat.colour;
// ambient
return mat.ambient * col;
}
/*
P is point of intersection,
N normal in this point,
R direction of reflected ray,
V direction to the viewer
*/
Colour PhongShader_calculate(Material &mat, Vector3 &P, Vector3 &N, Vector3 &R, Vector3 &V,
Light &light)
{
Colour I = Colour();
Vector3 L = light.pos - P;
L.normalize();
Float L_dot_N = dot(L, N);
Float R_dot_V = dot(R, V);
Colour col;
if (mat.texture)
col = mat.texture->evaluate(P);
else
col = mat.colour;
// diffuse
I = mat.diffuse * col * light.colour * L_dot_N;
// specular
if (R_dot_V > 0)
I += mat.specular * light.colour * powf(R_dot_V, mat.shininess);
return I;
}
Colour Raytracer::raytrace(Ray &ray, int depth, Shape *origin_shape)
{
Float nearest_distance = Inf;
Shape *nearest_shape = top->nearest_intersection(origin_shape, ray, nearest_distance);
if (nearest_shape == NULL) {
return bg_colour;
} else {
Colour col = Colour();
Vector3 P = ray.o + ray.dir * nearest_distance; // point of intersection
Vector3 normal = nearest_shape->normal(P);
bool from_inside = false;
// make shapes double sided
if (dot(normal, ray.dir) > 0.0)
{
normal = - normal;
from_inside = true;
}
col = PhongShader_ambient(*nearest_shape->material, P);
vector<Light*>::iterator light;
for (light = lights.begin(); light != lights.end(); light++) {
Vector3 jo, L = (*light)->pos - P; // direction vector to light
L.normalize();
Float L_dot_N = dot(L, normal);
if (L_dot_N > 0) {
// test if this light is occluded (sharp shadows)
if ((*light)->cast_shadows) {
Ray shadow_ray = Ray(P, L);
Float dist = FLT_MAX;
if (top->nearest_intersection(nearest_shape, shadow_ray, dist))
continue;
}
// shading function
Vector3 R = L - 2.0 * L_dot_N * normal;
col += PhongShader_calculate(*nearest_shape->material,
P, normal, R, ray.dir, **light);
}
}
if (depth < max_depth)
{
Colour trans_col, refl_col;
Float trans = nearest_shape->material->transmissivity;
Float refl = nearest_shape->material->reflectivity;
const Float cos_i = - dot(normal, ray.dir);
// reflection
if (refl > 0.01)
{
Vector3 newdir = ray.dir + 2.0 * cos_i * normal;
Ray newray = Ray(P, newdir);
refl_col = raytrace(newray, depth + 1, nearest_shape);
}
// refraction
if (trans > 0.01)
{
Float n, n1, n2;
if (from_inside)
{
n1 = nearest_shape->material->refract_index;
n2 = 1.0;
n = n1;
}
else
{
n1 = 1.0;
n2 = nearest_shape->material->refract_index;
n = 1.0 / n2;
}
const Float sin2_t = n*n * (1 - cos_i*cos_i);
if (sin2_t >= 1.0)
{
// totally reflected
refl += trans;
trans = 0;
}
else
{
const Float cos_t = sqrtf(1 - sin2_t);
const Float Rdiv = 1.0/(n1*cos_i + n2*cos_t);
const Float Rper = (n1*cos_i - n2*cos_t)*Rdiv;
const Float Rpar = (n2*cos_i - n1*cos_t)*Rdiv;
const Float R = (Rper*Rper + Rpar*Rpar)/2;
refl += R*trans;
trans = (1-R)*trans;
Vector3 newdir = n * ray.dir + (n*cos_i - cos_t) * normal;
Ray newray = Ray(P, newdir);
trans_col = raytrace(newray, depth + 1, nearest_shape);
}
}
col = (1-refl-trans)*col + refl*refl_col + trans*trans_col;
}
// ambient occlusion
if (!from_inside && ao_samples)
{
Float miss = 0;
for (int i = 0; i < ao_samples; i++) {
Vector3 dir = SphereDistribute(i, ao_samples, ao_angle, normal);
Ray ao_ray = Ray(P, dir);
Float dist = ao_distance;
Shape *shape_in_way = top->nearest_intersection(nearest_shape, ao_ray, dist);
if (shape_in_way == NULL)
miss += 1.0;
else
miss += dist / ao_distance;
}
Float ao_intensity = miss / ao_samples;
col = col * ao_intensity;
}
return col;
}
}
#if 0
static void *renderrow(void *data)
{
RenderrowData *d = (RenderrowData*) data;
const int subsample = d->rt->getSubsample();
const Float subsample2 = 1.0/(subsample*subsample);
const int oversample = d->rt->getOversample();
const int ww = d->w*3;
Vector3 dir = d->dfix;
for (int x = 0; x < d->w; x += subsample) {
// generate a ray from eye passing through this pixel
if (subsample > 1)
{
Colour ic;
// top-left
dir.normalize();
Ray ray(d->eye, dir);
Colour c1 = d->rt->raytrace(ray, 0, NULL);
// top-right
Vector3 tmpdir = dir + (subsample-1)*d->dx;
tmpdir.normalize();
ray.dir = tmpdir;
Colour c2 = d->rt->raytrace(ray, 0, NULL);
// bottom right
tmpdir += (subsample-1)*d->dy;
tmpdir.normalize();
ray.dir = tmpdir;
Colour c4 = d->rt->raytrace(ray, 0, NULL);
// bottom left
tmpdir = dir + (subsample-1)*d->dy;
tmpdir.normalize();
ray.dir = tmpdir;
Colour c3 = d->rt->raytrace(ray, 0, NULL);
// are the colors similar?
Float m = (c1-c2).mag2();
m = max(m, (c2-c3).mag2());
m = max(m, (c3-c4).mag2());
m = max(m, (c4-c1).mag2());
if (m < 0.001)
{
// interpolate
Float *i = d->iter;
for (int x = 0; x < subsample; x++)
{
for (int y = 0; y < subsample; y++)
{
ic = c1*(subsample-x)*(subsample-y)*subsample2
+ c2*(x)*(subsample-y)*subsample2
+ c3*(subsample-x)*(y)*subsample2
+ c4*(x)*(y)*subsample2;
*(i + ww*y) = ic.r;
*(i + ww*y + 1) = ic.g;
*(i + ww*y + 2) = ic.b;
}
i += 3;
}
d->iter = i;
}
else
{
// render all pixels
Vector3 tmpdir = dir;
if (oversample)
{
for (int x = 0; x < subsample; x++)
{
for (int y = 0; y < subsample; y++)
{
Vector3 tmp2dir = tmpdir + y*d->dy;
samplepixel(ic, tmp2dir, d, oversample);
*(d->iter + ww*y) = ic.r;
*(d->iter + ww*y + 1) = ic.g;
*(d->iter + ww*y + 2) = ic.b;
}
d->iter += 3;
tmpdir += d->dx;
}
}
else
{
/* this is so complex because it tries to reuse
already computed corner pixels
though, above code will also work for non-oversampling... */
// first column
*(d->iter) = c1.r;
*(d->iter + 1) = c1.g;
*(d->iter + 2) = c1.b;
for (int y = 1; y < subsample-1; y++)
{
Vector3 tmp2dir = tmpdir + y*d->dy;
tmp2dir.normalize();
ray.dir = tmp2dir;
ic = d->rt->raytrace(ray, 0, NULL);
*(d->iter + ww*y) = ic.r;
*(d->iter + ww*y + 1) = ic.g;
*(d->iter + ww*y + 2) = ic.b;
}
*(d->iter + ww*(subsample-1)) = c3.r;
*(d->iter + ww*(subsample-1) + 1) = c3.g;
*(d->iter + ww*(subsample-1) + 2) = c3.b;
d->iter += 3;
tmpdir += d->dx;
// middle
for (int x = 1; x < subsample-1; x++)
{
for (int y = 0; y < subsample; y++)
{
Vector3 tmp2dir = tmpdir + y*d->dy;
tmp2dir.normalize();
ray.dir = tmp2dir;
ic = d->rt->raytrace(ray, 0, NULL);
*(d->iter + ww*y) = ic.r;
*(d->iter + ww*y + 1) = ic.g;
*(d->iter + ww*y + 2) = ic.b;
}
d->iter += 3;
tmpdir += d->dx;
}
// last column
*(d->iter) = c2.r;
*(d->iter + 1) = c2.g;
*(d->iter + 2) = c2.b;
for (int y = 1; y < subsample-1; y++)
{
Vector3 tmp2dir = tmpdir + y*d->dy;
tmp2dir.normalize();
ray.dir = tmp2dir;
ic = d->rt->raytrace(ray, 0, NULL);
*(d->iter + ww*y) = ic.r;
*(d->iter + ww*y + 1) = ic.g;
*(d->iter + ww*y + 2) = ic.b;
}
*(d->iter + ww*(subsample-1)) = c4.r;
*(d->iter + ww*(subsample-1) + 1) = c4.g;
*(d->iter + ww*(subsample-1) + 2) = c4.b;
d->iter += 3;
}
}
}
else // subsample <= 1
{
Colour c;
samplepixel(c, dir, d, oversample);
// write color to buffer
*d->iter++ = c.r;
*d->iter++ = c.g;
*d->iter++ = c.b;
}
dir += d->dx*subsample;
}
#ifdef PTHREADS
pthread_exit((void *)d);
#endif
return (void *)d;
}
#endif
void *Raytracer::raytrace_worker(void *d)
{
static const int my_queue_size = 256;
Raytracer *rt = (Raytracer*)d;
Sample my_queue[my_queue_size];
Colour my_colours[my_queue_size];
int my_count;
Ray ray;
for (;;)
{
pthread_mutex_lock(&rt->sample_queue_mutex);
while (rt->sample_queue_count == 0)
{
if (rt->end_of_samples)
{
dbgmsg(4, "T thread [%d] exiting\n", pthread_self());
pthread_mutex_unlock(&rt->sample_queue_mutex);
pthread_exit(NULL);
}
pthread_cond_wait(&rt->sample_queue_cond, &rt->sample_queue_mutex);
}
if (rt->sample_queue_count >= my_queue_size)
my_count = my_queue_size;
else
my_count = rt->sample_queue_count;
rt->sample_queue_count -= my_count;
// copy samples to local queue
if (rt->sample_queue_pos + my_count >= rt->sample_queue_size)
{
register int c = rt->sample_queue_size - rt->sample_queue_pos;
memcpy(my_queue, rt->sample_queue + rt->sample_queue_pos, c*sizeof(Sample));
memcpy(my_queue + c, rt->sample_queue, (my_count - c)*sizeof(Sample));
rt->sample_queue_pos = my_count - c;
}
else
{
memcpy(my_queue, rt->sample_queue + rt->sample_queue_pos,
my_count*sizeof(Sample));
rt->sample_queue_pos += my_count;
}
if (rt->sample_queue_count <= my_queue_size*2)
pthread_cond_signal(&rt->worker_ready_cond);
pthread_mutex_unlock(&rt->sample_queue_mutex);
// do the work
for (int i = 0; i < my_count; i++)
{
ray = rt->camera->makeRay(my_queue[i]);
my_colours[i] = rt->raytrace(ray, 0, NULL);
}
// save the results
pthread_mutex_lock(&rt->sampler_mutex);
for (int i = 0; i < my_count; i++)
rt->sampler->saveSample(my_queue[i], my_colours[i]);
pthread_mutex_unlock(&rt->sampler_mutex);
}
}
void Raytracer::render()
{
if (!sampler || !camera || !top)
return;
static const int my_count_max = 256;
sample_queue_size = my_count_max*2*num_threads;
sample_queue = new Sample [sample_queue_size];
sample_queue_pos = 0;
sample_queue_count = 0;
end_of_samples = false;
int my_count, my_pos = 0;
int sampnum = 0, sampdone;
bool more_samples;
sampler->init();
// create workers
dbgmsg(1, "* running %d threads\n", num_threads);
pthread_t threads[num_threads];
for (int t = 0; t < num_threads; t++)
{
int rc = pthread_create(&threads[t], NULL, raytrace_worker, (void*)this);
if (rc) {
dbgmsg(0, "\nE pthread_create unsuccessful, return code was %d\n", rc);
exit(1);
}
}
dbgmsg(1, "* raytracing...\n");
dbgmsg(2, "- 0%% done");
pthread_mutex_lock(&sampler_mutex);
while ( (sampnum = sampler->initSampleSet()) > 0 )
{
sampdone = 0;
for (;;)
{
my_count = 0;
while ( more_samples = sampler->nextSample(&sample_queue[my_pos++]) )
{
my_count++;
if (my_pos >= sample_queue_size)
my_pos = 0;
if (my_count >= my_count_max)
break;
}
if (!more_samples && !my_count)
break;
pthread_mutex_unlock(&sampler_mutex);
pthread_mutex_lock(&sample_queue_mutex);
sample_queue_count += my_count;
// wait for workers if there is enough samples ready on queue
while (sample_queue_count > (2*num_threads-1)*my_count_max)
pthread_cond_wait(&worker_ready_cond, &sample_queue_mutex);
pthread_cond_signal(&sample_queue_cond);
pthread_mutex_unlock(&sample_queue_mutex);
sampdone += my_count;
dbgmsg(2, "\b\b\b\b\b\b\b\b%2d%% done", (sampdone - sample_queue_count)*100/sampnum);
pthread_mutex_lock(&sampler_mutex);
}
dbgmsg(2, "\b\b\b\b\b\b\b\b100%% done\n");
}
pthread_mutex_unlock(&sampler_mutex);
// wait for workers
dbgmsg(2, "- waiting for threads to finish\n");
pthread_mutex_lock(&sample_queue_mutex);
end_of_samples = true;
while (sample_queue_count)
{
pthread_cond_broadcast(&sample_queue_cond);
pthread_mutex_unlock(&sample_queue_mutex);
pthread_mutex_lock(&sample_queue_mutex);
}
pthread_cond_broadcast(&sample_queue_cond);
pthread_mutex_unlock(&sample_queue_mutex);
for (int t = 0; t < num_threads; t++)
pthread_join(threads[t], NULL);
delete[] sample_queue;
}
void Raytracer::addlight(Light *light)
{
lights.push_back(light);
}
void Raytracer::ambientocclusion(int samples, Float distance, Float angle)
{
ao_samples = samples;
ao_distance = distance;
ao_angle = angle;
if (ao_distance == 0)
/* 0 ==> Inf */
ao_distance = Inf;
}