/* RAY_RAD - A Very Simple Ray-Traced Radiosity Renderer */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define FLOOR 0 /* X-Y plane, Z = 0.0 */
#define SOUTH_WALL 1 /* X-Z plane, Y = 0.0 */
#define EAST_WALL 2 /* Y-Z plane, X = 8.0 */
#define NORTH_WALL 3 /* X-Z plane, Y = 8.0 */
#define WEST_WALL 4 /* X-Z plane, X = 0.0 */
#define MAX_RAYS 25000 /* Maximum number of rays */
#define MAX_VAL 1.0e+6 /* Maximum floating point value */
#define MIN_VAL 1.0e-6 /* Minimum floating point value */
#define NUM_LOOP 20 /* Number of iterations */
#define NUM_SURF 10 /* Number of surfaces */
#define PI 3.141592654
#define DRAND() ((double) rand() / (double) RAND_MAX)
typedef struct /* 3-D point coordinates */
{
double x; /* X-axis coordinate */
double y; /* Y-axis coordinate */
double z; /* Z-axis coordinate */
} PT_3D;
typedef struct /* 3-D ray */
{
PT_3D org; /* Origin */
PT_3D dir; /* Direction */
} RAY;
typedef struct /* Element */
{
double total; /* Total flux */
double unsent; /* Unsent flux */
} ELEM;
typedef struct /* Surface */
{
PT_3D vertex[4]; /* Vertex array */
int num_row; /* Number of element rows */
int num_col; /* Number of element columns */
double col_dim; /* Element column dimension */
double row_dim; /* Element row dimension */
double emittance; /* Emittance */
double reflectance; /* Reflectance */
ELEM *elemp; /* Element array pointer */
} SURF;
static int send_flux(int, int, int);
static void display_surface(void);
static void find_element(RAY *, int *, int *, int *);
static void global_to_local(int, PT_3D *, PT_3D *);
static void local_to_global(int, RAY *, RAY *);
static void select_ray(int, int, int, RAY *);
static double MaxEmittance = 0.0; /* Maximum emittance */
static SURF RoomSurf[NUM_SURF] = { /* Room surfaces */
{ { { 0.0, 0.0, 0.0 }, { 8.0, 0.0, 0.0 }, /* Floor */
{ 8.0, 8.0, 0.0 }, { 0.0, 8.0, 0.0 } }, 8, 8, 1.0, 1.0,
0.0, 0.2, NULL },
{ { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 8.0 }, /* South wall */
{ 8.0, 0.0, 8.0 }, { 8.0, 0.0, 0.0 } }, 8, 8, 1.0, 1.0,
0.0, 0.5, NULL },
{ { { 8.0, 0.0, 0.0 }, { 8.0, 0.0, 8.0 }, /* East wall */
{ 8.0, 8.0, 8.0 }, { 8.0, 8.0, 0.0 } }, 8, 8, 1.0, 1.0,
0.0, 0.5, NULL },
{ { { 8.0, 8.0, 0.0 }, { 8.0, 8.0, 8.0 }, /* North wall */
{0.0, 8.0, 8.0 }, { 0.0, 8.0, 0.0 } }, 8, 8, 1.0, 1.0,
0.0, 0.5, NULL },
{ { { 0.0, 8.0, 0.0 }, { 0.0, 8.0, 8.0 }, /* West wall */
{ 0.0, 0.0, 8.0 }, { 0.0, 0.0, 0.0 } }, 8, 8, 1.0, 1.0,
0.0, 0.5, NULL },
{ { { 0.0, 0.0, 8.0 }, { 0.0, 3.0, 8.0 }, /* South ceil. */
{ 8.0, 3.0, 8.0 }, { 8.0, 0.0, 8.0 } }, 8, 3, 1.0, 1.0,
0.0, 0.8, NULL },
{ { { 5.0, 3.0, 8.0 }, { 5.0, 5.0, 8.0 }, /* East ceiling */
{8.0, 5.0, 8.0 }, { 8.0, 3.0, 8.0 } }, 3, 2, 1.0, 1.0,
0.0, 0.8, NULL },
{ { { 0.0, 5.0, 8.0 }, { 0.0, 8.0, 8.0 }, /* North ceil. */
{ 8.0, 8.0, 8.0 }, { 8.0, 5.0, 8.0 } }, 8, 3, 1.0, 1.0,
0.0, 0.8, NULL },
{ { { 0.0, 3.0, 8.0 }, { 0.0, 5.0, 8.0 }, /* West ceiling */
{ 3.0, 5.0, 8.0 }, { 3.0, 3.0, 8.0 } }, 3, 2, 1.0, 1.0,
0.0, 0.8, NULL },
{ { { 3.0, 3.0, 8.0 }, { 3.0, 5.0, 8.0 }, /* Ceil. light */
{ 5.0, 5.0, 8.0 }, { 5.0, 3.0, 8.0 } }, 2, 2, 1.0, 1.0,
5000.0, 0.8, NULL } };
int main(void)
{
int col; /* Column counter */
int elem; /* Element counter */
int maximum; /* Maximum number of rays */
int num_elem; /* Number of elements */
int num_rays; /* Number of rays */
int row; /* Row counter */
int surf; /* Surface counter */
ELEM *elemp; /* Element pointer */
SURF *surfp; /* Surface pointer */
for (surf = 0; surf < NUM_SURF; surf++) {
/* Instantiate elements */
surfp = &(RoomSurf[surf]);
num_elem = surfp->num_row * surfp->num_col;
if ((surfp->elemp = (ELEM *) malloc(num_elem * sizeof(ELEM)))
== NULL) {
fputs("Out of memory!\n", stderr);
return (2);
}
elemp = surfp->elemp;
for (elem = 0; elem < num_elem; elem++, elemp++)
elemp->total = elemp->unsent = surfp->emittance;
if (surfp->emittance > MaxEmittance)
MaxEmittance = surfp->emittance;
}
do {
/* Distribute flux between elements */
num_rays = 0;
for (surf = 0; surf < NUM_SURF; surf++) {
surfp = &(RoomSurf[surf]);
for (row = 0; row < surfp->num_row; row++)
for (col = 0; col < surfp->num_col; col++} {
if ((maximum = send_flux(surf, row, col)) > num_rays)
num_rays = maximum;
}
}
display_surface(); /* Display intermediate results */
} while (num_rays > 0); /* Repeat until convergence */
for (surf = 0; surf < NUM_SURF; surf++) /* Free memory */
free(RoomSurf[surf].elemp);
return (0);
}
static int send_flux ( /* Send flux to other elements */
int s_surf, /* Sending surface index */
int s_row, /* Sending element row */
int s_col ) /* Sending element column */
{
int h_col; /* Hit element column */
int h_row; /* Hit element row */
int h_surf; /* Hit surface index */
int num_rays; /* Number of rays */
int ray; /* Ray counter */
double inc_flux; /* Incident (ray) flux */
double ref_flux; /* Reflected flux */
RAY shoot; /* Shooting ray */
ELEM *h_elemp; /* Hit element pointer */
ELEM *s_elemp; /* Sending element pointer */
SURF *h_surfp; /* Hit surface pointer */
SURF *s_surfp; /* Sending surface pointer */
s_surfp = &(RoomSurf[s_surf]);
s_elemp = s_surfp->elemp + s_row * s_surfp->num_col + s_col;
/* Number of rays to shoot proportional to unsent flux */
if ((num_rays = (int) (MAX_RAYS * s_elemp->unsent /
MaxEmittance)) == 0)
return (0);
inc_flux = s_elemp->unsent / num_rays; /* Get ray flux */
For (ray = 0; ray < num_rays; ray++) { /* Distribute flux */
select_ray(s_surf, s_row, s_col, &shoot); /* Select ray */
/* Find hit surface and element */
find_element(&shoot, &h_surf, &h_row, &h_col);
/* Calculate flux reflected from hit element */
h_surfp = &(RoomSurf[h_surf]);
h_elemp = h_surfp->elemp + h_row * h_surfp->num_col + h_col;
ref_flux = h_surfp->reflectance * inc_flux;
h_elemp->total += ref_flux; /* Update total flux */
h_elemp->unsent += ref_flux; /* Update unsent flux */
}
s_elemp->unsent = 0.0; /* Reset sender's unsent flux */
return (num_rays);
}
static void select_ray ( /* Randomly select ray */
int surf, /* Surface index */
int row, /* Element row */
int col, /* Element column */
RAY *rayp ) /* Ray pointer */
{
double horz; /* Ray direction local horizontal angle */
double vert; /* Ray direction local vertical angle */
RAY local; /* Ray local coordinates */
SURF *surfp; /* Surface pointer */
surfp = &(RoomSurf[surf]);
/* Randomly select local (surface) ray origin coordinates */
/* within sending element boundary */
local.org.x = ((double) col + DRAND()) / surfp->col_dim;
local.org.y = ((double) row + DRAND()) / surfp->row_dim;
local.org.z = 0.0;
/* Randomly select local (surface) ray direction coordinates */
/* with cosine probability distribution in vertical plane */
horz = 2.0 * PI * DRAND(); /* Horizontal angle */
vert = acos(sqrt(1.0 - DRAND())); /* Vertical angle */
/* Convert from spherical to rectilinear coordinates */
local.dir.x = sin(vert) * cos(horz);
local.dir.y = sin(vert) * sin(horz);
local.dir.z = cos(vert);
/* Convert ray coordinates from local to global */
local_to_global(surf, &local, rayp);
}
static void find_element ( /* Find intersected element */
RAY *rayp, /* Ray pointer */
int *h_surfp, /* Hit surface index pointer */
int *h_rowp, /* Hit element row pointer */
int *h_colp ) /* Hit element column pointer */
{
int surf; /* Surface counter */
double s = MAX_VAL; /* Smallest parametric value */
double t; /* Current parametric value */
PT_3D hit; /* Hit point coordinates */
PT_3D temp; /* Temporary point coordinates */
SURF *surfp; /* Surface pointer */
for (surf = 0; surf < NUM_SURF; surf++) {
surfp = &(RoomSurf[surf]);
/* Determine if and where ray intersects surface plane */
switch (surf) {
case EAST_WALL:
case WEST_WALL:
if (fabs(rayp->dir.x) > MIN_VAL) {
t = (surfp->vertex[0].x - rayp->org.x) / rayp->dir.x;
if (t > MIN_VAL && t < s) {
s = t;
*h_surfp = surf;
}
}
break;
case SOUTH_WALL:
case NORTH_WALL:
if (fabs(rayp->dir.y) > MIN_VAL) {
t = (surfp->vertex[0].y - rayp->org.y) / rayp->dir.y;
if (t > MIN_VAL && t < s) {
s = t;
*h_surfp = surf;
}
}
break;
case FLOOR:
if (fabs(rayp->dir.z) > MIN_VAL) {
t = (surfp->vertex[0].z - rayp->org.z) / rayp->dir.z;
if (t > MIN_VAL && t < s) {
S = t;
*h_surfp = surf;
}
}
break;
default:
if (fabs(rayp->dir.z) > MIN_VAL) {
t = (surfp->vertex[0].z - rayp->org.z) / rayp->dir.z;
if (t = MIN_VAL && t < s) {
/* Check for intersection inside surface boundary */
temp.x = rayp->org.x + t * rayp->dir.x;
temp.y = rayp->org.y + t * rayp->dir.y;
if ((temp.x >= surfp->vertex[0].x && temp.x < surfp->vertex[3].x) &&
(temp.y >= surfp->vertex[0].y && temp.y < surfp->vertex[1].y)) {
s = t;
*h_surfp = surf;
}
}
}
break;
}
}
/* Calculate hit point global coordinates */
temp.x = rayp->org.x + s * rayp->dir.x;
temp.y = rayp->org.y + s * rayp->dir.y;
temp.z = rayp->org.z + s * rayp->dir.z;
/* Convert global to local (surface) coordinates */
global_to_local(*h_surfp, &temp, &hit);
/* Calculate hit element row and column indices */
*h_rowp = (int) (hit.y / surfp->row_dim);
*h_colp = (int) (hit.x / surfp->col_dim);
}
static void global_to_local ( /* Convert global to local */
int surf, /* Surface index */
PT_3D *globalp, /* Global coordinates pointer */
PT_3D *localp ) /* Local coordinates pointer */
{
PT_3D *vp; /* Vertex pointer */
vp = &(RoomSurf[surf].vertex[0]);
switch (surf) {
case FLOOR:
localp->x = globalp->x - vp->x;
localp->y = globalp->y - vp->y;
break;
case SOUTH_WALL:
localp->x = globalp->z - vp->z;
localp->y = globalp->x - vp->x;
break;
case EAST_WALL:
localp->x = globalp->z - vp->z;
localp->y = globalp->y - vp->y;
break;
case NORTH_WALL:
localp->x = globalp->z - vp->z;
localp->y = vp->x - globalp->x;
break;
case WEST_WALL:
localp->x = globalp->z - vp->z;
localp->y = vp->y - globalp->y;
break;
default:
localp->x = globalp->y - vp->y;
localp->y = globalp->x - vp->x;
break;
}
}
static void local_to_global ( /* Convert ray coordinates */
int surf, /* Surface index */
RAY *localp, /* Local coordinates pointer */
RAY *globalp ) /* Global coordinates pointer */
{
PT_3D *vp; /* Vertex pointer */
vp = &(RoomSurf[surf].vertex[0]);
switch (surf) {
case FLOOR:
globalp->org.x = vp->x + localp->org.x;
globalp->org.y = vp->y + localp->org.y;
globalp->org.z = vp->z + localp->org.z;
globalp->dir.x = localp->dir.x;
globalp->dir.y = localp->dir.y;
globalp->dir.z = localp->dir.z;
break;
case SOUTH_WALL:
globalp->org.x = vp->x + localp->org.y;
globalp->org.y = vp->y + localp->org.z;
globalp->org.z = vp->z + localp->org.x;
globalp->dir.x = localp->dir.y;
globalp->dir.y = localp->dir.z;
globalp->dir.z = localp->dir.x;
break;
case EAST_WALL:
globalp->org.x = vp->x - localp->org.z;
globalp->org.y = vp->y + localp->org.y;
globalp->org.z = vp->z + localp->org.x;
globalp->dir.x = -localp->dir.z;
globalp->dir.y = localp->dir.y;
globalp->dir.z = localp->dir.x;
break;
case NORTH_WALL:
globalp->org.x = vp->x - localp->org.y;
globalp->org.y = vp->y - localp->org.z;
globalp->org.z = vp->z + localp->org.x;
globalp->dir.x = -localp->dir.y;
globalp->dir.y = -localp->dir.z;
globalp->dir.z = localp->dir.x;
break;
case WEST_WALL:
globalp->org.x = vp->x + localp->org.z;
globalp->org.y = vp->y - localp->org.y;
globalp->org.z = vp->z + localp->org.x;
globalp->dir.x = localp->dir.z;
globalp->dir.y = -localp->dir.y;
globalp->dir.z = localp->dir.x;
break;
default:
globalp->org.x = vp->x + localp->org.y;
globalp->org.y = vp->y + localp->org.x;
globalp->org.z = vp->z - localp->org.z;
globalp->dir.x = localp->dir.y;
globalp->dir.y = localp->dir.x;
globalp->dir.z = -localp->dir.z;
break;
}
}
static void display_surface(void) /* Display surfaces */
{
int col; /* Element column index */
int row; /* Element row index */
int surf; /* Surface counter */
ELEM *elemp; /* Element pointer */
SURF *surfp; /* Surface pointer */
static int loop = 0; /* Iteration counter */
printf("Loop %d\n", loop++);
for (surf = 0; surf < NUM_SURF; surf++) {
printf("Surface %d\n", surf);
surfp = &(RoomSurf[surf]);
for (row = 0; row < surfp->num_row; row++) {
for (col = 0; col < surfp->num_col; col++) {
elemp = surfp->elemp + row * surfp->num_col + col;
printf("% 8.2f ", elemp->total / (PI * surfp->row_dim *
surfp->col_dim));
}
putchar('\n');
}
}
}
/* End of File */