This is sortContour.c in view mode; [Download] [Up]
/* Sort contours in the order of outer most to inner most.
This is necessary for filling the inside of contours with color.
See below.
William Grosso 92-05-27
*/
#include <stdio.h>
#include <stdlib.h>
#include "contour.h"
#ifndef TRUE
#define TRUE 1
#endif
#ifndef FALSE
#define FALSE 0
#endif
/* Function prototypes private to this module */
int contoursub1(float xp, float yp, ContourPath *currentj, int npts);
int contoursub2(double x1, double y1, double x2, double y2, double xp, double yp);
/*
=====================================================================================
A brief explanation: Our goal for this subroutine is this: to determine
which contours lie inside other contours (we assume no contours
intersect). The algorithm is essentially an adaptation of the pl proof
of the Jordan Curve Theorem: we draw lines out to infinity and check
parity. E.g. Given closed curves C1 and C2 we cant to check whether C1
contains C2 (C2 is "inside of" C1). To do so, we choose a point [xp,yp]
on C2, draw a line from [xp,yp] to infinity, and count the number of
intersections this line has C1. If even, C2 is not inside C1. If odd,
C2 is inside C1. If the above was incomprehensible, may I suggest "The
Topology of Plane Sets" by Newman. If you're curious about how this
works in higher dimensions, Moise's "Geometric Topology in Dimensions 2
& 3" is a fine book (as is Zeeman's "Seminar on Combinatorial topology").
functions: contoursub1 draws the line and checks whether the line
intersects each of the line segments that make up C1 (by using
contoursub2--which solves the linear equation).
optimization notes: the code is fairly incomprehensible at points
because I wanted to optimize. The matrix "matrixie" is superfluous to
the algorithm, but helps cut down the number of "containment" checks (if
C1 is contained in C2, then C2 cannot be contained in C1). A 1 in the
i,j place indicates Ci is contained in Cj. the (i,nContours+1) place
holds the total number of contours that contain Ci
syntax notes:
C1 above is the "i" contour below. C2 above is the "j" contour below.
=====================================================================================
Usage Notes:
The last argument must have the storage for pointers allocated as below before
calling the function:
SCntrPtr = (ContourPath **)(malloc((size_t)(sizeof(ContourPath *)*nContours)));
cntrList: linked list of contours containing all the information about contours.
nContours: number of contours in cntrList.
SCntrPtr: array of pointers to contours in cntrList, sorted in outer to inner order.
=====================================================================================
*/
void sort_contourList(ContourPath *cntrList, float xxmin, float xxmax, float yymin, float yymax,
int nContours, ContourPath **SCntrPtr)
{
int i, j, k, jj=0;
float xp=0.0, yp= 0.0;
int *matrixie; /* where we'll store the containment data (temporarily).
matrixie (i,j) indicates whether contour i is inside countour j
Why is this here at all ? Well, we first check to
see if contour i is contained in contour j (for i<j).
Then, we use the fact that two contours can not each
contain the other to cut down on our checking time. */
ContourPath *currenti, *currentj;
/* pointers to the current paths for the i and j loops */
/*
We've declared our variables...Now we must initialize them.
Only the matrixie must be initialized (the others are done in loops).
*/
matrixie = (int *)calloc((size_t)((nContours+1)*nContours), (size_t)sizeof(int) );
currenti=cntrList;
for (i=0;i<nContours;i++) {
*(matrixie+(i+1)*(nContours+1)-1)=0; /* how many guys contain contour i */
currentj=cntrList;
/* pick the point on contour-i which is inside the domain as the point from which
we draw the line out to infinity. Start with the x[3], y[3] in an attempt to
get a point that is not on the edge.
*/
for(j=3; j<(currenti->num_pts+3); j++) {
jj = j % currenti->num_pts;
xp=currenti->x[jj];
yp=currenti->y[jj];
if(xp >= xxmin && xp <= xxmax && yp >= yymin && yp <= yymax)
break; /* we got it */
}
for (j=0;j<nContours;j++) {
if ((j==i) || ((j<i) && (*(matrixie+j*(nContours+1)+i)==1))) {
currentj=currentj->next;
continue;
}
*(matrixie+i*(nContours+1)+j)=contoursub1(xp,yp,currentj,currentj->num_pts);
*(matrixie +(i+1)*(nContours+1)-1)+=*(matrixie +i*(nContours+1)+j);
currentj=currentj->next;
}
currenti=currenti->next;
}
/*
If all went well, we now have a nContours by nContours matrix
holding all the "containment" information. Now, we create a list telling
the drawing routine which ones to draw first (array of pointers to cntrList: SCntrPtr).
The methodology is this: We create a list that has the contours which
aren't contained in any other contour at the head of the list. Then come
the contours that are contained in only one other contour. Then come the
contours that are contained in only two other contours...
*/
k=0;
for (i=0;i<nContours;i++) {
currentj=cntrList;
for(j=0;j<nContours;j++) {
if (*(matrixie+(j+1)*(nContours +1)-1) == i) {
SCntrPtr[k]=currentj;
k++;
}
currentj=currentj->next;
}
if (k==nContours) break;
}
free(matrixie);
}
/* See commments above: this checks to see whether the line from [xp,yp] hits
contourj an odd or an even number of times.
Return val = 0 (xp,yp) is outside ContourPath *currentj
Return val = 1 (xp,yp) is inside ContourPath *currentj
*/
int contoursub1 (float xp , float yp, ContourPath *currentj , int npts)
{
int k, counter=0;
for (k=1;k<=npts;k++) {
counter += contoursub2( currentj->x[k],currentj->y[k],
currentj->x[k-1],currentj->y[k-1],xp,yp);
}
counter += contoursub2( currentj->x[npts],currentj->y[npts],
currentj->x[0],currentj->y[0],xp,yp);
return(counter%2);
}
/* solves the linear equation for contoursub1
* return val = 1: lines intersect
* return val = 0: don't intersect.
* Tests whether a line from (xp,yp) intersects line (x1,y1) --- (x2,y2).
*/
int contoursub2 (double x1, double y1, double x2, double y2, double xp, double yp)
{
double dummy;
if(y1 == y2) return (0);
dummy=(yp-y1)/(y2-y1);
if ((0.0<=dummy) && (dummy<1.0) && ((dummy*(x2-x1)+x1-xp) >= 0.0)) return(1);
else return(0);
}
/* --------------------------------------------------------------------------------
Routine returns (2) if (xp,yp) is inside the
polygon x[n_path], y[n_path], (0) if outside, and (1) if exactly on the
path edge. Uses non-zero winding rule in Adobe PostScript Language
reference manual, section 4.6 on Painting. Path should have been closed
first, so that x[n_path-1] = x[0], and y[n_path-1] = y[0].
This is version 2, trying to kill a bug in above routine: If point is
on X edge, fails to discover that it is on edge.
We are imagining a ray extending "up" from the point (xp,yp); the ray
has equation x = xp, y >= yp. Starting with crossing_count = 0, we add 1
each time the path crosses the ray in the +x direction, and subtract 1
each time the ray crosses in the -x direction. After traversing the
entire polygon, if (crossing_count) then point is inside. We also watch
for edge conditions.
If two or more points on the path have x[i] == xp, then we have an
ambiguous case, and we have to find the points in the path before and
after this section, and check them.
--------------------------------------------------------------------------------- */
int non_zero_winding(float xp, float yp, float *x, float *y, int n_path)
{
int i, j, k, jend, crossing_count, above;
double y_sect;
above = FALSE;
crossing_count = 0;
/* First make sure first point in path is not a special case: */
j = jend = n_path - 1;
if (x[j] == xp) {
/* Trouble already. We might get lucky: */
if (y[j] == yp) return(1);
/* Go backward down the polygon until x[i] != xp: */
if (y[j] > yp) above = TRUE;
i = j - 1;
while (x[i] == xp && i > 0) {
if (y[i] == yp) return (1);
if (!(above) && y[i] > yp) above = TRUE;
i--;
}
/* Now if i == 0 polygon is degenerate line x=xp;
since we know xp,yp is inside bounding box,
it must be on edge: */
if (i == 0) return(1);
/* Now we want to mark this as the end, for later: */
jend = i;
/* Now if (j-i)>1 there are some segments the point could be exactly on: */
for (k = i+1; k < j; k++) {
if((y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp))
return (1);
}
/* Now we have arrived where i is > 0 and < n_path-1, and x[i] != xp.
We have been using j = n_path-1. Now we need to move j forward
from the origin: */
j = 1;
while (x[j] == xp) {
if (y[j] == yp) return (1);
if (!(above) && y[j] > yp) above = TRUE;
j++;
}
/* Now at the worst, j == jstop, and we have a polygon with only 1 vertex
not at x = xp. But now it doesn't matter, that would end us at
the main while below. Again, if j>=2 there are some segments to check: */
for (k = 0; k < j-1; k++) {
if((y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp))
return (1);
}
/* Finally, we have found an i and j with points != xp.
* If (above) we may have crossed the ray:
*/
if(above && x[i] < xp && x[j] > xp)
crossing_count++;
else if(above && x[i] > xp && x[j] < xp)
crossing_count--;
/* End nightmare scenario for x[0] == xp. */
}
else {
/* Get here when x[0] != xp: */
i = 0;
j = 1;
while (x[j] == xp && j < jend) {
if (y[j] == yp) return (1);
if (!(above) && y[j] > yp) above = TRUE;
j++;
}
/* Again, if j==jend, (i.e., 0) then we have a polygon with only 1 vertex
not on xp and we will branch out below. */
/* if ((j-i)>2) the point could be on intermediate segments: */
for (k = i+1; k < j-1; k++) {
if ( (y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp) )
return (1);
}
/* Now we have x[i] != xp and x[j] != xp.
If (above) and x[i] and x[j] on opposite sides, we are certain to have crossed the ray.
If not (above) and (j-i)>1, then we have not crossed it.
If not (above) and j-i == 1, then we have to check the intersection point.
*/
if (x[i] < xp && x[j] > xp) {
if(above)
crossing_count++;
else if( (j-i) == 1) {
y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
if (y_sect == yp) return (1);
if (y_sect > yp) crossing_count++;
}
}
if (x[i] > xp && x[j] < xp) {
if (above)
crossing_count--;
else if ( (j-i) == 1) {
y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
if (y_sect == yp) return (1);
if (y_sect > yp) crossing_count--;
}
}
/* End easier case for x[0] != xp */
}
/* Now MAIN WHILE LOOP begins:
Set i = j, and search for a new j, and do as before. */
i = j;
while (i < jend) {
above = FALSE;
j = i+1;
while (x[j] == xp) {
if (y[j] == yp) return (1);
if (!(above) && y[j] > yp) above = TRUE;
j++;
}
/* if ((j-i)>2) the point could be on intermediate segments: */
for (k = i+1; k < j-1; k++) {
if ( (y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp) )
return (1);
}
/* Now we have x[i] != xp and x[j] != xp.
If (above) and x[i] and x[j] on opposite sides, we are certain to have crossed the ray.
If not (above) and (j-i)>1, then we have not crossed it.
If not (above) and j-i == 1, then we have to check the intersection point. */
if (x[i] < xp && x[j] > xp) {
if (above)
crossing_count++;
else if ( (j-i) == 1) {
y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
if (y_sect == yp) return (1);
if (y_sect > yp) crossing_count++;
}
}
if (x[i] > xp && x[j] < xp) {
if (above)
crossing_count--;
else if ( (j-i) == 1) {
y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
if (y_sect == yp) return (1);
if (y_sect > yp) crossing_count--;
}
}
/* That's it for this piece. Advance i: */
i = j;
}
/* End of MAIN WHILE. Get here when we have gone all around without landing on edge. */
if (crossing_count)
return(2);
else
return(0);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.