blob: 6fccd2040e31aa673f5bf11e125a044e59b08360 [file] [log] [blame]
/*************************************************************************/
/* */
/* Copyright (c) 1994 Stanford University */
/* */
/* All rights reserved. */
/* */
/* Permission is given to use, copy, and modify this software for any */
/* non-commercial purpose as long as this copyright notice is not */
/* removed. All other uses, including redistribution in whole or in */
/* part, are forbidden without prior written permission. */
/* */
/* This software is provided with absolutely no warranty and no */
/* support. */
/* */
/*************************************************************************/
/************************************************************************
*
* Patch management.
*
* This module has the following functions:
* (1) Create/initialize a new instance of the patch object.
* (2) Management of BSP tree (insertion,traversal)
*
*************************************************************************/
#include <stdio.h>
EXTERN_ENV;
include(radiosity.h)
static void _foreach_patch(Patch *node, void (*func)(), long arg1, long process_id);
static void _foreach_d_s_patch(Vertex *svec, Patch *node, void (*func)(), long arg1, long process_id);
static void split_into_3(Patch *patch, ElemVertex *ev1, ElemVertex *ev2, ElemVertex *ev3, Edge *e12, Edge *e23, Edge *e31, Patch *parent, long process_id);
static void split_into_2(Patch *patch, ElemVertex *ev1, ElemVertex *ev2, ElemVertex *ev3, Edge *e12, Edge *e23, Edge *e31, Patch *parent, long process_id);
/***************************************************************************
****************************************************************************
*
* Methods of Patch object
*
****************************************************************************
****************************************************************************
****************************************************************************
*
* foreach_patch_in_bsp()
*
* General purpose driver. For each patch in the BSP tree, apply specified
* function. The function is passed a pointer to the patch.
* Traversal is in-order, that is, subtree(-) -> node -> subtree(+).
*
****************************************************************************/
void foreach_patch_in_bsp(void (*func)(), long arg1, long process_id)
{
_foreach_patch( global->bsp_root, func, arg1, process_id ) ;
}
static void _foreach_patch(Patch *node, void (*func)(), long arg1, long process_id)
{
if( node == 0 )
return ;
/* Process subtree(-) */
if( node->bsp_negative )
_foreach_patch( node->bsp_negative, func, arg1, process_id ) ;
/* Apply function to this node */
func( node, arg1, process_id ) ;
/* Process subtree(+) */
if( node->bsp_positive )
_foreach_patch( node->bsp_positive, func, arg1, process_id ) ;
}
/***************************************************************************
*
* foreach_depth_sorted_patch()
*
* For each patch in the BSP tree, apply specified function. In the depth
* sorted order along the given vector (from tail to arrow head of the
* vector).
* The function is passed a pointer to the patch.
*
****************************************************************************/
void foreach_depth_sorted_patch(Vertex *sort_vec, void (*func)(), long arg1, long process_id)
{
_foreach_d_s_patch( sort_vec, global->bsp_root, func, arg1, process_id ) ;
}
static void _foreach_d_s_patch(Vertex *svec, Patch *node, void (*func)(), long arg1, long process_id)
{
float sign ;
if( node == 0 )
return ;
/* Compute inner product */
sign = inner_product( svec, &node->plane_equ.n ) ;
if( sign >= 0.0 )
{
/* The vector is approaching from the negative side of the patch */
/* Process subtree(-) */
if( node->bsp_negative )
_foreach_d_s_patch( svec, node->bsp_negative, func, arg1, process_id ) ;
/* Apply function to this node */
func( node, arg1, process_id ) ;
/* Process subtree(+) */
if( node->bsp_positive )
_foreach_d_s_patch( svec, node->bsp_positive, func, arg1, process_id ) ;
}
else
{
/* Process subtree(+) */
if( node->bsp_positive )
_foreach_d_s_patch( svec, node->bsp_positive, func, arg1, process_id ) ;
/* Apply function to this node */
func( node, arg1, process_id ) ;
/* Process subtree(-) */
if( node->bsp_negative )
_foreach_d_s_patch( svec, node->bsp_negative, func, arg1, process_id ) ;
}
}
/***************************************************************************
*
* define_patch()
*
* Insert a new patch in the BSP tree and put refinement task in the queue.
*
****************************************************************************/
void define_patch(Patch *patch, Patch *root, long process_id)
{
Patch *parent = root ;
long xing_code ;
/* Lock the BSP tree */
LOCK(global->bsp_tree_lock);
/* If this is the first patch, link directly */
if( parent == 0 )
{
if( global->bsp_root == 0 )
{
/* This is really the first patch */
global->bsp_root = patch ;
patch->bsp_positive = 0 ;
patch->bsp_negative = 0 ;
patch->bsp_parent = 0 ;
attach_element( patch, process_id ) ;
UNLOCK(global->bsp_tree_lock);
return ;
}
else
/* Race condition. The root was NULL when the task was
created */
parent = global->bsp_root ;
}
/* Traverse the BSP tree and get to the leaf node */
while( 1 )
{
/* Check the sign */
xing_code = patch_intersection( &parent->plane_equ, &patch->p1,
&patch->p2, &patch->p3, process_id ) ;
/* Traverse down the tree according to the sign */
if( POSITIVE_SIDE( xing_code ) )
{
if( parent->bsp_positive == 0 )
{
/* Insert the patch */
parent->bsp_positive = patch ;
patch->bsp_parent = parent ;
attach_element( patch, process_id ) ;
UNLOCK(global->bsp_tree_lock);
foreach_patch_in_bsp( refine_newpatch, (long)patch, process_id ) ;
return ;
}
else
/* Traverse down to the subtree(+) */
parent = parent->bsp_positive ;
}
else if( NEGATIVE_SIDE( xing_code ) )
{
if( parent->bsp_negative == 0 )
{
/* Insert the patch */
parent->bsp_negative = patch ;
patch->bsp_parent = parent ;
attach_element( patch, process_id ) ;
UNLOCK(global->bsp_tree_lock);
foreach_patch_in_bsp( refine_newpatch, (long)patch, process_id ) ;
return ;
}
else
/* Traverse down to the subtree(-) */
parent = parent->bsp_negative ;
}
else
{
/* The patch must be split. Insertion is taken care of by
split_patch(). */
UNLOCK(global->bsp_tree_lock);
split_patch( patch, parent, xing_code, process_id ) ;
return ;
}
}
}
/***************************************************************************
*
* split_patch()
* split_into_3()
* split_into_2()
*
* Split a patch and insert in the BSP tree.
* split_into_3() Split a patch into 3 patches. The routine assumes both
* P1-P2 and P1-P3 intersect the plane.
* split_into_2() Split a patch into 2 patches. The routine assuems P1 is
* on the plane and P2-P3 intersects the plane.
* split_patch() Classify intersection type, rename vertices, and
* call split_into_X().
*
****************************************************************************/
void split_patch(Patch *patch, Patch *node, long xing_code, long process_id)
{
long c1, c2, c3 ;
c1 = P1_CODE( xing_code ) ;
c2 = P2_CODE( xing_code ) ;
c3 = P3_CODE( xing_code ) ;
/* Classify intersection type */
if( c1 == c2 )
/* P3 is on the oposite side */
split_into_3( patch, patch->ev3, patch->ev1, patch->ev2,
patch->e31, patch->e12, patch->e23, node, process_id) ;
else if( c1 == c3 )
/* P2 is on the oposite side */
split_into_3( patch, patch->ev2, patch->ev3, patch->ev1,
patch->e23, patch->e31, patch->e12, node, process_id ) ;
else if( c2 == c3 )
/* P1 is on the oposite side */
split_into_3( patch, patch->ev1, patch->ev2, patch->ev3,
patch->e12, patch->e23, patch->e31, node, process_id ) ;
else if( c1 == POINT_ON_PLANE )
/* P1 is on the plane. P2 and P3 are on the oposite side */
split_into_2( patch, patch->ev1, patch->ev2, patch->ev3,
patch->e12, patch->e23, patch->e31, node, process_id ) ;
else if( c2 == POINT_ON_PLANE )
/* P2 is on the plane. P3 and P1 are on the oposite side */
split_into_2( patch, patch->ev2, patch->ev3, patch->ev1,
patch->e23, patch->e31, patch->e12, node, process_id ) ;
else
/* P3 is on the plane. P1 and P2 are on the oposite side */
split_into_2( patch, patch->ev3, patch->ev1, patch->ev2,
patch->e31, patch->e12, patch->e23, node, process_id ) ;
}
static void split_into_3(Patch *patch, ElemVertex *ev1, ElemVertex *ev2, ElemVertex *ev3, Edge *e12, Edge *e23, Edge *e31, Patch *parent, long process_id)
{
ElemVertex *ev_a ; /* Intersection of P1-P2 & the patch */
ElemVertex *ev_b ; /* Intersection of P1-P3 & the patch */
float h1, h2, h3 ;
float u2, u3 ;
Patch *new ;
Edge *e_ab, *e_3a ;
long rev_e12, rev_e31 ;
/* Compute intersection in terms of parametarized distance from P1 */
h1 = plane_equ( &parent->plane_equ, &ev1->p, process_id ) ;
h2 = plane_equ( &parent->plane_equ, &ev2->p, process_id ) ;
h3 = plane_equ( &parent->plane_equ, &ev3->p, process_id ) ;
/* NOTE: Length of P1-P2 and P1-P3 are at least 2*F_COPLANAR.
So, no check is necessary before division */
u2 = h1 / (h1 - h2) ;
if( (rev_e12 = EDGE_REVERSE( e12, ev1, ev2 )) )
subdivide_edge( e12, u2, process_id ) ;
else
subdivide_edge( e12, (float)1.0 - u2, process_id ) ;
ev_a = e12->ea->pb ;
u3 = h1 / (h1 - h3) ;
if( (rev_e31 = EDGE_REVERSE( e31, ev3, ev1 )) )
subdivide_edge( e31, (float)1.0 - u3, process_id ) ;
else
subdivide_edge( e31, u3, process_id ) ;
ev_b = e31->ea->pb ;
/* Now insert patches in the tree */
/* (1) Put P1-Pa-Pb */
new = get_patch(process_id) ;
new->p1 = ev1->p ;
new->p2 = ev_a->p ;
new->p3 = ev_b->p ;
new->ev1 = ev1 ;
new->ev2 = e12->ea->pb ;
new->ev3 = e31->ea->pb ;
new->e12 = (!rev_e12)? e12->ea : e12->eb ;
new->e23 = e_ab = create_edge(ev_a, ev_b, process_id ) ;
new->e31 = (!rev_e31)? e31->eb : e31->ea ;
new->plane_equ = patch->plane_equ ;
new->area = u2 * u3 * patch->area ;
new->color = patch->color ;
new->emittance = patch->emittance ;
define_patch( new, parent, process_id ) ;
/* (2) Put Pa-P2-P3 */
new = get_patch(process_id) ;
new->p1 = ev_a->p ;
new->p2 = ev2->p ;
new->p3 = ev3->p ;
new->ev1 = ev_a ;
new->ev2 = ev2 ;
new->ev3 = ev3 ;
new->e12 = (!rev_e12)? e12->eb : e12->ea ;
new->e23 = e23 ;
new->e31 = e_3a = create_edge( ev3, ev_a, process_id ) ;
new->plane_equ = patch->plane_equ ;
new->area = (1.0 - u2) * patch->area ;
new->color = patch->color ;
new->emittance = patch->emittance ;
define_patch( new, parent, process_id ) ;
/* (3) Put Pa-P3-Pb. Reuse the original patch */
patch->p1 = ev_a->p ;
patch->p2 = ev3->p ;
patch->p3 = ev_b->p ;
patch->ev1 = ev_a ;
patch->ev2 = ev3 ;
patch->ev3 = ev_b ;
patch->e12 = e_3a ;
patch->e23 = (!rev_e31)? e31->ea : e31->eb ;
patch->e31 = e_ab ;
patch->area = u2 * (1.0 - u3) * patch->area ;
define_patch( patch, parent, process_id ) ;
}
static void split_into_2(Patch *patch, ElemVertex *ev1, ElemVertex *ev2, ElemVertex *ev3, Edge *e12, Edge *e23, Edge *e31, Patch *parent, long process_id)
{
ElemVertex *ev_a ;
Edge *e_a1 ;
float h2, h3 ;
float u2 ;
Patch *new ;
long rev_e23 ;
/* Compute intersection in terms of parameterized distance from P2 */
h2 = plane_equ( &parent->plane_equ, &ev2->p, process_id ) ;
h3 = plane_equ( &parent->plane_equ, &ev3->p, process_id ) ;
/* NOTE: Length of P2-P3 is at least 2*F_COPLANAR.
So, no check is necessary before division */
u2 = h2 / (h2 - h3) ;
if( (rev_e23 = EDGE_REVERSE( e23, ev2, ev3 )) )
subdivide_edge( e23, u2, process_id ) ;
else
subdivide_edge( e23, (float)1.0 - u2, process_id ) ;
ev_a = e23->ea->pb ;
/* Now put patches in the tree */
/* (1) Put P1-P2-Pa */
new = get_patch(process_id) ;
new->p1 = ev1->p ;
new->p2 = ev2->p ;
new->p3 = ev_a->p ;
new->ev1 = ev1 ;
new->ev2 = ev2 ;
new->ev3 = ev_a ;
new->e12 = e12 ;
new->e23 = (!rev_e23)? e23->ea : e23->eb ;
new->e31 = e_a1 = create_edge( ev_a, ev1, process_id ) ;
new->plane_equ = patch->plane_equ ;
new->area = u2 * patch->area ;
new->color = patch->color ;
new->emittance = patch->emittance ;
define_patch( new, parent, process_id ) ;
/* (2) Put P1-Pa-P3. Reuse the original patch */
patch->p1 = ev1->p ;
patch->p2 = ev_a->p ;
patch->p3 = ev3->p ;
patch->ev1 = ev1 ;
patch->ev2 = ev_a ;
patch->ev3 = ev3 ;
patch->e12 = e_a1 ;
patch->e23 = (!rev_e23)? e23->eb : e23->ea ;
patch->e31 = e31 ;
patch->area = (1.0 - u2) * patch->area ;
define_patch( patch, parent, process_id ) ;
}
/***************************************************************************
*
* attach_element()
*
* Attach an element to the patch. This element becomes the
* root of the quad tree.
*
****************************************************************************/
void attach_element(Patch *patch, long process_id)
{
Element *pelem ;
/* Create and link an element to the patch */
pelem = get_element(process_id) ;
patch->el_root = pelem ;
/* Initialization of the element */
pelem->patch = patch ;
pelem->ev1 = patch->ev1 ;
pelem->ev2 = patch->ev2 ;
pelem->ev3 = patch->ev3 ;
pelem->e12 = patch->e12 ;
pelem->e23 = patch->e23 ;
pelem->e31 = patch->e31 ;
pelem->area = patch->area ;
pelem->rad = patch->emittance ;
}
/***************************************************************************
*
* refine_newpatch()
*
* Recursively subdivide
*
****************************************************************************/
void refine_newpatch(Patch *patch, long newpatch, long process_id)
{
long cc ;
Patch *new_patch = (Patch *)newpatch ;
/* Check sequence number */
if( patch->seq_no >= new_patch->seq_no )
/* Racing condition due to multiprocessing */
return ;
/* Check visibility */
cc = patch_intersection( &patch->plane_equ,
&new_patch->p1, &new_patch->p2, &new_patch->p3, process_id ) ;
if( NEGATIVE_SIDE(cc) )
/* If negative or on the plane, then do nothing */
return ;
cc = patch_intersection( &new_patch->plane_equ,
&patch->p1, &patch->p2, &patch->p3, process_id ) ;
if( NEGATIVE_SIDE(cc) )
/* If negative or on the plane, then do nothing */
return ;
/* Create a new task or do it by itself */
create_ff_refine_task( patch->el_root, new_patch->el_root, 0, process_id ) ;
}
/***************************************************************************
*
* get_patch()
*
* Returns a new instance of the Patch object.
*
****************************************************************************/
Patch *get_patch(long process_id)
{
Patch *p ;
/* LOCK the free list */
LOCK(global->free_patch_lock);
/* Test pointer */
if( global->free_patch == 0 )
{
printf( "Fatal: Ran out of patch buffer\n" ) ;
UNLOCK(global->free_patch_lock);
exit( 1 ) ;
}
/* Get a patch data structure */
p = global->free_patch ;
global->free_patch = p->bsp_positive ;
global->n_total_patches++ ;
global->n_free_patches-- ;
/* Unlock the list */
UNLOCK(global->free_patch_lock);
/* Clear pointers just in case.. */
p->el_root = 0 ;
p->bsp_positive = 0 ;
p->bsp_negative = 0 ;
p->bsp_parent = 0 ;
return( p ) ;
}
/***************************************************************************
*
* init_patchlist()
*
* Initialize patch free list.
*
****************************************************************************/
void init_patchlist(long process_id)
{
long i ;
/* Initialize Patch free list */
for( i = 0 ; i < MAX_PATCHES-1 ; i++ )
{
global->patch_buf[i].bsp_positive = &global->patch_buf[i+1] ;
global->patch_buf[i].seq_no = i ;
}
global->patch_buf[ MAX_PATCHES-1 ].bsp_positive = 0 ;
global->patch_buf[ MAX_PATCHES-1 ].seq_no = MAX_PATCHES - 1 ;
global->free_patch = global->patch_buf ;
global->n_total_patches = 0 ;
global->n_free_patches = MAX_PATCHES ;
LOCKINIT(global->free_patch_lock) ;
#if PATCH_ASSIGNMENT == PATCH_ASSIGNMENT_COSTBASED
/* Initialize Patch_Cost structure */
for( i = 0 ; i < MAX_PATCHES ; i++ )
{
global->patch_cost[i].patch = &global->patch_buf[i] ;
global->patch_cost[i].cost_lock
= get_sharedlock( SHARED_LOCK_SEGANY, process_id ) ;
global->patch_cost[i].n_bsp_node = 0 ;
global->patch_cost[i].n_total_inter = 0 ;
global->patch_cost[i].cost_estimate = 0 ;
global->patch_cost[i].cost_history[0] = 0 ;
global->patch_cost[i].cost_history[1] = 0 ;
global->patch_cost[i].cost_history[2] = 0 ;
global->patch_cost[i].cost_history[3] = 0 ;
}
#endif
}
/***************************************************************************
*
* print_patch()
*
* Print patch information.
*
****************************************************************************/
void print_patch(Patch *patch, long process_id)
{
printf( "Patch (%ld)\n", (long)patch ) ;
print_point( &patch->p1 ) ;
print_point( &patch->p2 ) ;
print_point( &patch->p3 ) ;
print_plane_equ( &patch->plane_equ, process_id ) ;
printf( "\tArea %f\n", patch->area ) ;
}
/***************************************************************************
*
* print_bsp_tree()
*
* Print BSP tree
*
****************************************************************************/
void print_bsp_tree(long process_id)
{
printf( "**** BSP TREE ***\n" ) ;
foreach_patch_in_bsp( _pr_patch, 0, process_id ) ;
printf( "\n\n" ) ;
}
void _pr_patch(Patch *patch, long dummy, long process_id)
{
print_patch( patch, process_id ) ;
}
/***************************************************************************
****************************************************************************
*
* Methods for PlaneEqu object
*
****************************************************************************
****************************************************************************
****************************************************************************
*
*
* plane_equ()
*
* Returns the value H(Px, Py, Pz) where H is the equation of the plane
*
****************************************************************************/
float plane_equ(PlaneEqu *plane, Vertex *point, long process_id)
{
float h ;
h = plane->c + point->x * plane->n.x
+ point->y * plane->n.y
+ point->z * plane->n.z ;
return( h ) ;
}
/***************************************************************************
*
* comp_plane_equ()
*
* Compute plane equation from the three vertices on the plane.
*
****************************************************************************/
float comp_plane_equ(PlaneEqu *pln, Vertex *p1, Vertex *p2, Vertex *p3, long process_id)
{
float length ;
/* Compute normal vector */
length = plane_normal( &pln->n, p1, p2, p3 ) ;
/* Calculate constant factor */
pln->c = -inner_product( &pln->n, p1 ) ;
return( length ) ;
}
/***************************************************************************
*
* point_intersection()
* patch_intersection()
*
* Returns the intersection code according to the relationship between the
* point/patch and the plane.
* point_intersection() returns 2 bits code that represents:
* 01: Point is on the positive side (H(x,y,z) > 0)
* 10: Point is on the negative side (H(x,y,z) < 0)
* 00: Point is on the plane (H(x,y,z) = 0)
*
* patch_intersection() returns 3 sets of 2 bits code each represents the
* relationship of each vertex of the triangle patch.
*
****************************************************************************/
long point_intersection(PlaneEqu *plane, Vertex *point, long process_id)
{
float h ;
long result_code = 0 ;
/* Compare H(x,y,z) against allowance */
if( (h = plane_equ( plane, point, process_id )) < -F_COPLANAR )
result_code |= POINT_NEGATIVE_SIDE ;
if( h > F_COPLANAR )
result_code |= POINT_POSITIVE_SIDE ;
return( result_code ) ;
}
long patch_intersection(PlaneEqu *plane, Vertex *p1, Vertex *p2, Vertex *p3, long process_id)
{
long c1, c2, c3 ;
c1 = point_intersection( plane, p1, process_id ) ;
c2 = point_intersection( plane, p2, process_id ) ;
c3 = point_intersection( plane, p3, process_id ) ;
return( (c3 << 4) | (c2 << 2) | c1 ) ;
}
/***************************************************************************
*
* print_plane_equ()
*
* Print plane equation
*
****************************************************************************/
void print_plane_equ(PlaneEqu *peq, long process_id)
{
printf( "\tPLN: %.3f x + %.3f y + %.3f z + %.3f\n",
peq->n.x, peq->n.y, peq->n.z, peq->c ) ;
}