open-axiom repository from github
/***********************************************************************
contour.c
begun 19 September 1992, Jim Wen
***********************************************************************/
#include "header.h"
#include "draw.h"
#define segmentDEBUG_X
#define contourDEBUG_X
#define use_min
#define realloc_bug_fixed_NOT
/*=====================================================================*
Static variables
*=====================================================================*/
int noo=0;
poly *contour_poly_list;
poly *active_list_first, *active_list_last, *active_list_current;
segment_list_struct *tmp_segment_list;
/*=====================================================================*
macro definitions
*=====================================================================*/
#define foreach_poly(p) for ((p)=contour_poly_list; p != NIL(poly); p=p->next)
#define foreach_active(p,f,l) if (f != NIL(poly)) for (p=f; p!=l; p=(p)->next)
#define in_range(x,a,b) ( ((x>=a) && (x<b)) || ((x<a) && (x>=b)) )
#define poly_in_plane(p,z) ( (p->contour_min==z) && (p->contour_max==z) )
/*=====================================================================*
local function declarations
*=====================================================================*/
int contour_compare();
void add_segment();
segment_struct *make_the_segment();
viewTriple *mkpoint();
CONTOUR_float get_t_from_pts();
void make_active_list();
int maintain_active_list();
void contour_minMaxPolygons();
/*=====================================================================*
contour_compare()
The compare function passed to msort.
*=====================================================================*/
int contour_compare(p1,p2)
poly *p1, *p2;
{
#ifdef use_min
if (p1->contour_min < p2->contour_min) return(-1);
else if (p1->contour_min == p2->contour_min) return(0);
else return(1);
#else
if (p1->contour_max > p2->contour_max) return(-1);
else if (p1->contour_max == p2->contour_max) return(0);
else return(1);
#endif /* use_min */
} /* contour_compare() */
/*=====================================================================*
do_contour_map()
*=====================================================================*/
int do_contour_map()
{
poly *pp;
poly *ap;
CONTOUR_float z_now;
int got_more;
int *anIndex;
viewTriple *aPt, *daPt, *one_point;
int jj, segment_index;
int got_one_intersection;
int done;
segment_list_struct *sl;
#ifdef contour_object_rotate
float rotMat[4][4];
float transformed_zmax, transformed_zmin;
#endif
#ifdef conDEBUG
segment_struct *seg;
#endif
/*---------------------------------------------------*
the flag "did_contour" should be set to "no" whenever
the user modifies one of the parameters
*---------------------------------------------------*/
#ifdef contourDEBUG_x
fprintf(stderr,"Contour is %s\n",(did_contour)?"yes":"no");
#endif
if (did_contour) return 1;
did_contour = yes;
#ifdef contour_object_rotate
/*---------------------------------------------------*
transform all the points for arbitrary plane
slicing. this includes all the viewTriples
being referenced as well as the boundaries
in viewData (to get contour_min, contour_max).
*---------------------------------------------------*/
#ifdef oldie
rot_theta = 0;
rot_phi = pi/2;
#endif
sinTheta = sin(-rot_theta);
cosTheta = cos(-rot_theta);
sinPhi = sin(rot_phi);
cosPhi = cos(rot_phi);
ROTATE1(rotMat);
/*---------------------------------------------------*
transform all the points.
the zmin and zmax values need to be recalculated
now that the object has been transformed. note
that transforming the z extreme points will not
work correctly (i know, i've tried it).
*---------------------------------------------------*/
{
int i,j,k;
LLPoint *anLLPoint;
LPoint *anLPoint;
int *anIndex;
viewTriple *daPoint;
float v_in[4], v_out[4];
int first_time = yes;
anLLPoint = viewData.lllp.llp;
for (i=0; i<viewData.lllp.numOfComponents; i++,anLLPoint++) {
anLPoint = anLLPoint->lp;
for (j=0; j<anLLPoint->numOfLists; j++,anLPoint++) {
anIndex = anLPoint->indices;
for (k=0; k<anLPoint->numOfPoints; k++,anIndex++) {
daPoint = refPt3D(viewData, *anIndex);
/*---------------------------------------------------*
inefficient code so that the vector package could
be used to see if things work;
should change viewTriple's <x,y,z> to an array(?)
*---------------------------------------------------*/
v_in[0] = daPoint->x; v_in[1] = daPoint->y;
v_in[2] = daPoint->z; v_in[3] = 1.0;
vectorMatrix4(v_in, rotMat, v_out);
daPoint->contour_x = v_out[0];
daPoint->contour_y = v_out[1];
daPoint->contour_z = v_out[2];
if (first_time) {
first_time = no;
transformed_zmin = transformed_zmax = v_out[2];
} else {
if (v_out[2] < transformed_zmin)
transformed_zmin = v_out[2];
else if (v_out[2] > transformed_zmax)
transformed_zmax = v_out[2];
}
} /* for points in LPoints (k) */
} /* for LPoints in LLPoints (j) */
} /* for LLPoints in LLLPoints (i) */
}
#endif
/*---------------------------------------------------*
set up the step size - it should be user adjustable
and include all the slices possible.
max_cuts is #slices seen
z_step is #slices made
cuts_used is #slices actually displayed
*---------------------------------------------------*/
z_step = (transformed_zmax - transformed_zmin)/max_cuts;
#ifdef oldie
cuts_used = max_cuts;
#endif
/*---------------------------------------------------*
calculate the bounds of the polygons - the
contour_minMaxPolygons routine looks at the
contour points rather than the object or
projected points
*---------------------------------------------------*/
contour_poly_list = copyPolygons(viewData.polygons);
contour_minMaxPolygons(contour_poly_list);
/*---------------------------------------------------*
sort the polygons by the zmax value
(or whatever transformed value, in general)
*---------------------------------------------------*/
contour_poly_list = msort(contour_poly_list, 0,
viewData.numPolygons, contour_compare);
/*---------------------------------------------------*
having figured out how many cuts we need (should
be a one time overhead so the following stuff
should be in the initialization routine), we
allocate an array of segment lists and initialize
them (this part can stay here).
*---------------------------------------------------*/
#ifdef oldie
segment_list = saymem("contour.c: segment_list",
max_cuts, sizeof(segment_list_struct));
#else
/*---------------------------------------------------*
if the append flag is set, then we want to add
the new segment stuff onto the end of the old
stuff (to build up a model of the surface).
tmp_segment_list is use to keep track of the
head of the new list while the routine goes
through its paces.
*---------------------------------------------------*/
if (contour_append_lists_flag) {
noo = 0;
#ifdef segmentDEBUG
fprintf(stderr,"======= series %d ========\n",noo);
fprintf(stderr," ---> before: sl->num=%d [%x]\n",
segment_list->num, segment_list);
#endif
#ifdef realloc_bug_fixed
realloc(segment_list, (cuts_used + max_cuts) * sizeof(segment_list_struct));
tmp_segment_list = segment_list;
segment_list += cuts_used; /* shift to end of list */
cuts_used += max_cuts; /* size of new list */
#ifdef segmentDEBUG
fprintf(stderr,"> %d cuts => seg at %x [old=%x]\n",
cuts_used, segment_list, tmp_segment_list);
fprintf(stderr," sl->num=%d, tsl->num=%d\n",
segment_list->num, tmp_segment_list->num);
#endif
#else /* DONT_WORK */
/*---------------------------------------------------*
Because realloc doesn't seem to work properly,
we need to do this by hand
*---------------------------------------------------*/
/*---------------------------------------------------*
allocate new space
*---------------------------------------------------*/
tmp_segment_list = saymem("contour.c: segment_list, append",
cuts_used + max_cuts,
sizeof(segment_list_struct));
/*---------------------------------------------------*
copy over old data (1..cuts_used)
*---------------------------------------------------*/
{
segment_list_struct *tsl;
for (segment_index=0, sl=segment_list, tsl=tmp_segment_list;
segment_index<cuts_used;
segment_index++, sl++, tsl++) {
tsl->num = sl->num;
tsl->max_num = sl->max_num;
tsl->num_segs = sl->num_segs;
tsl->segments = sl->segments;
}
}
/*---------------------------------------------------*
free the old stuff
*---------------------------------------------------*/
/* free(segment_list); */
/*---------------------------------------------------*
now set segment_list to point to the point
where tmp_segment_list stops - there ought
to me max_cuts storage spaces left.
*---------------------------------------------------*/
segment_list = tmp_segment_list + cuts_used;
/*---------------------------------------------------*
update cuts_used to have everything for a possible
next iteration
*---------------------------------------------------*/
cuts_used += max_cuts;
#endif /* realloc_BUG */
} else {
if (contour_allocated) {
/* free(segment_list); */
} else {
contour_allocated = yes;
}
noo = 0;
segment_list = saymem("contour.c: segment_list",
max_cuts, sizeof(segment_list_struct));
cuts_used = max_cuts;
#ifdef segmentDEBUG
fprintf(stderr,"======= series %d ========\n",noo);
fprintf(stderr,"%d cuts => seg at %x\n",cuts_used, segment_list);
#endif
}
#endif
for (segment_index=0, sl=segment_list;
segment_index<max_cuts;
segment_index++, sl++) {
sl->num = ++noo;
sl->max_num = max_cuts;
#ifdef segmentDEBUG
fprintf(stderr,"Made segment list %d [%x]\n",noo,sl);
/* fprintf(stderr," ...(tsl->num=%d [%x]\n",
tmp_segment_list->num, tmp_segment_list); */
#endif
sl->num_segs = 0;
sl->segments = NIL(segment_struct);
}
/*---------------------------------------------------*
create an "active_list" of polygons such that
this slice step intersects all and only those
polygons in the active list.
*---------------------------------------------------*/
make_active_list(contour_poly_list,
#ifdef use_min
transformed_zmin, transformed_zmin+z_step,
#else
transformed_zmax, transformed_zmax-z_step,
#endif
&active_list_first, &active_list_last, &active_list_current);
/*---------------------------------------------------*
iterate from zmax down to zmin with z_step increments
*---------------------------------------------------*/
segment_index = 0;
#ifdef use_min
for (z_now=transformed_zmin; z_now<transformed_zmax; /* see below for incr*/) {
#else
for (z_now=transformed_zmax; z_now>transformed_zmin; /* see below for incr*/) {
#endif
/*---------------------------------------------------*
for each of the polygons on the active list,
intersect each of line equation for the sides
with the plane equation for the plane at z_now.
one of the following may occur:
no intersections: haha - can't happen coz active
one intersection : at the point, create point line
two intersections: create line connecting points
lies in the plane: create three segments
note that this is a fairly inefficient approach but
i'm just throwing this stuff together in an afternoon
to see how it looks
*---------------------------------------------------*/
foreach_active(ap, active_list_first, active_list_last) {
/*---------------------------------------------------*
do line-plane equation for each side of the
polygon
for now - just 3+ sided polygons (no degenerates)
*---------------------------------------------------*/
if (ap->numpts >= 3) {
if (poly_in_plane(ap, z_now)) {
/*---------------------------------------------------*
re-create all the segments of the polygon
*---------------------------------------------------*/
daPt = refPt3D(viewData, *(ap->indexPtr + (ap->numpts - 1)));
for (jj=0, anIndex=ap->indexPtr; jj<ap->numpts; jj++, anIndex++) {
aPt = refPt3D(viewData, *anIndex);
add_segment(segment_list, segment_index,
make_the_segment(aPt, daPt));
daPt = aPt;
}
} else {
/*---------------------------------------------------*
find the line that defines the intersection of
the polygon with the z-plane
*---------------------------------------------------*/
got_one_intersection = no;
done = no;
daPt = refPt3D(viewData, *(ap->indexPtr + (ap->numpts - 1)));
for (jj=0, anIndex=ap->indexPtr;
!done && jj<ap->numpts;
jj++, anIndex++) {
aPt = refPt3D(viewData, *anIndex);
if (in_range(z_now, aPt->contour_z, daPt->contour_z)) {
if (got_one_intersection) {
add_segment(segment_list, segment_index,
make_the_segment(one_point,
mkpoint(aPt, daPt, z_now)));
done = yes;
} else {
one_point = mkpoint(aPt, daPt, z_now);
got_one_intersection = yes;
}
}
daPt = aPt;
} /* for */
} /* else not lie in plane */
}
} /* foreach_active(ap) */
/*---------------------------------------------------*
maintain/update the active list, pruning off things
the fall off the top (beyond z_now - z_step) and
===> adding on things to the bottom that now fall into
the range [z_now ---> z_now-z_step].
*---------------------------------------------------*/
segment_index++;
#ifdef use_min
z_now += z_step;
got_more = maintain_active_list(&active_list_first,
&active_list_last,
&active_list_current,
z_now + z_step,
z_now);
#else
z_now -= z_step;
got_more = maintain_active_list(&active_list_first,
&active_list_last,
&active_list_current,
z_now,
z_now - z_step);
#endif
} /* for z_now from zmax to zmin */
/*---------------------------------------------------*
if the segment lists have been appended, reset
the global segment lists pointer to the top of
the lists
*---------------------------------------------------*/
if (contour_append_lists_flag) {
segment_list = tmp_segment_list;
#ifdef segmentDEBUG
#ifdef oldie
fprintf(stderr," setting seg to old=%x\n", segment_list);
fprintf(stderr," num is %d\n",segment_list->num);
#endif
{
for (segment_index=0, sl=segment_list;
segment_index<cuts_used;
segment_index++, sl++) {
fprintf(stderr," sl->num = %d\n",sl->num);
}
}
#endif
}
} /* do_contour_map() */
/*=====================================================================*
make_active_list(da_list, z_min, z_max, af, al, ac)
*=====================================================================*/
void make_active_list(da_list, z_min, z_max, af, al, ac)
poly *da_list;
CONTOUR_float z_min, z_max;
poly **af, **al, **ac;
{
poly *tmp_p;
/*---------------------------------------------------*
the first active polygon is the first one in the
given, sorted list. note that if it doesn't fall
===> inside the z_max --> z_min range, af is set to NIL.
*---------------------------------------------------*/
#ifdef use_min
if (da_list->contour_min > z_max) {
#else
if (da_list->contour_max < z_min) {
#endif
*af = NIL(poly);
return;
} else {
*af = da_list;
}
/*---------------------------------------------------*
the current active polygon is set to "af" at this
point but it could be the case that it is set to
"al" if, for example, the current span of z-values
has no polygons - af is set to NIL, and the next
time we look we start at ac=af.
*---------------------------------------------------*/
*ac = *af;
/*---------------------------------------------------*
the last active polygon is the polygon right before
===> the first one whose zmax is too small to make the
list
*---------------------------------------------------*/
*al = da_list;
tmp_p = da_list->next;
for (; tmp_p != NIL(poly); tmp_p = tmp_p->next) {
#ifdef use_min
if (tmp_p->contour_min > z_max) return;
#else
if (tmp_p->contour_max < z_min) return;
#endif
*al = tmp_p;
}
} /* make_active_list() */
/*=====================================================================*
maintain_active_list(af, al, ac, z_max, z_min)
*=====================================================================*/
int maintain_active_list(af, al, ac, z_max, z_min)
poly **al, **af, **ac;
CONTOUR_float z_max, z_min;
{
poly *tmp_p;
/*---------------------------------------------------*
first, get the lower boundary to be within range,
pruning elements from the head of the list
*---------------------------------------------------*/
*af = *ac;
#ifdef use_min
while ((*af) && (*af)->contour_max < z_min) {
#else
while ((*af) && (*af)->contour_min > z_max) {
#endif
*af = (*af)->next;
}
/*---------------------------------------------------*
check to see if the upper boundary is still in
range
*---------------------------------------------------*/
#ifdef use_min
if ((*af) == NIL(poly) || ((*af)->contour_min > z_max)) {
#else
if ((*af) == NIL(poly) || ((*af)->contour_max < z_min)) {
#endif
/* --- nope, it wasn't --- */
*ac = *af;
*af = NIL(poly);
return(0);
}
/*---------------------------------------------------*
upper boundary is okay...see if we need to add to
the list on the lower bound side
*---------------------------------------------------*/
tmp_p = (*al)->next;
for (; tmp_p != NIL(poly); tmp_p = tmp_p->next) {
#ifdef use_min
if (tmp_p->contour_min > z_max) return;
#else
if (tmp_p->contour_max < z_min) return;
#endif
*al = tmp_p;
}
return(1);
} /* maintain_active_list() */
/*=====================================================================*
add_segment(seg_list, index, seg)
*=====================================================================*/
void add_segment(seg_list, index, seg)
segment_list_struct *seg_list;
int index;
segment_struct *seg;
{
segment_list_struct *sl;
sl = seg_list + index;
seg->next = sl->segments;
sl->segments = seg;
sl->num_segs++;
} /* add_segment() */
/*=====================================================================*
make_the_segment(pt1, pt2)
*=====================================================================*/
segment_struct *make_the_segment(pt1, pt2)
viewTriple *pt1, *pt2;
{
segment_struct *seg;
seg = (segment_struct *)saymem("contour.c: segment",1,sizeof(segment_struct));
seg->point1 = pt1;
seg->point2 = pt2;
return(seg);
} /* make_the_segment() */
/*=====================================================================*
viewTriple *mkpoint(vt1, vt2, z_val)
*=====================================================================*/
viewTriple *mkpoint(vt1, vt2, z_val)
viewTriple *vt1, *vt2;
CONTOUR_float z_val;
{
viewTriple *vt;
CONTOUR_float t;
vt = (viewTriple *)saymem("contour.c: viewTriple",1,sizeof(viewTriple));
t = get_t_from_pts(vt1->contour_z, vt2->contour_z, z_val);
#ifdef waitaminute
vt->x = vt1->contour_x + (vt2->contour_x - vt1->contour_x) * t;
vt->y = vt1->contour_y + (vt2->contour_y - vt1->contour_y) * t;
vt->z = z_val;
#else
vt->x = vt1->x + (vt2->x - vt1->x) * t;
vt->y = vt1->y + (vt2->y - vt1->y) * t;
vt->z = vt1->z + (vt2->z - vt1->z) * t;
#endif
vt->contour_x = vt1->contour_x + (vt2->contour_x - vt1->contour_x) * t;
vt->contour_y = vt1->contour_y + (vt2->contour_y - vt1->contour_y) * t;
vt->contour_z = z_val;
return(vt);
} /* mkpoint() */
/*=====================================================================*
get_t_from_pts(z_min, z_max, z_val)
*=====================================================================*/
CONTOUR_float get_t_from_pts(z_min, z_max, z_val)
CONTOUR_float z_min, z_max, z_val;
{
CONTOUR_float t;
if (z_min == z_max) return 0;
t = (z_val - z_min)/(z_max - z_min);
return(t);
} /* get_t_from_pts() */
void contour_minMaxPolygons(aPoly)
poly *aPoly;
{
int *anIndex;
int i;
for (; aPoly != NIL(poly); aPoly = aPoly->next) {
anIndex = aPoly->indexPtr;
aPoly->contour_min = aPoly->contour_max =
refPt3D(viewData,*anIndex)->contour_z;
for (i=1,anIndex++; i<aPoly->numpts; i++,anIndex++) {
if (refPt3D(viewData,*anIndex)->contour_z < aPoly->contour_min)
aPoly->contour_min = refPt3D(viewData,*anIndex)->contour_z;
else if (refPt3D(viewData,*anIndex)->contour_z > aPoly->contour_max)
aPoly->contour_max = refPt3D(viewData,*anIndex)->contour_z;
}
}
} /* contour_minMaxPolygons */