home *** CD-ROM | disk | FTP | other *** search
/ Nebula 1995 August / NEBULA.mdf / SourceCode / Classes / ContourView / sortContour.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-12-02  |  12.2 KB  |  371 lines
  1. /*  Sort contours in the order of outer most to inner most.
  2.     This is necessary for filling the inside of contours with color.
  3.     See below.
  4.  
  5.       William Grosso 92-05-27
  6. */
  7.  
  8. #include <stdio.h>
  9. #include <stdlib.h>
  10.  
  11. #include "contour.h"
  12.  
  13.  
  14. #ifndef TRUE
  15. #define TRUE 1
  16. #endif
  17. #ifndef FALSE
  18. #define FALSE 0
  19. #endif
  20.  
  21.  
  22. /* Function prototypes private to this module */
  23. int contoursub1(float xp, float yp, ContourPath *currentj, int npts);
  24. int contoursub2(double x1, double y1, double x2, double y2, double xp, double yp);
  25.  
  26.  
  27. /*     
  28.     =====================================================================================
  29.     A brief explanation: Our goal for this subroutine is this: to determine
  30.     which contours lie inside other contours (we assume no contours
  31.     intersect). The algorithm is essentially an adaptation of the pl proof
  32.     of the Jordan Curve Theorem: we draw lines out to infinity and check
  33.     parity. E.g. Given closed curves C1 and C2 we cant to check whether C1
  34.     contains C2 (C2 is "inside of" C1). To do so, we choose a point [xp,yp]
  35.     on C2, draw a line from [xp,yp] to infinity, and count the number of
  36.     intersections this line has  C1. If even, C2 is not inside C1. If odd,
  37.     C2 is inside C1. If the above was  incomprehensible, may I suggest "The
  38.     Topology of Plane Sets" by Newman. If you're curious about how this
  39.     works in higher dimensions, Moise's "Geometric Topology in Dimensions 2
  40.     & 3" is a fine book (as is Zeeman's "Seminar on Combinatorial topology").
  41.     
  42.     functions: contoursub1 draws the line and checks whether the line
  43.     intersects each of the line segments that make up C1 (by using
  44.     contoursub2--which solves the linear equation).
  45.     
  46.     optimization notes: the code is fairly incomprehensible at points
  47.     because I wanted to optimize. The matrix "matrixie" is superfluous to
  48.     the algorithm, but helps cut down the number of "containment" checks (if
  49.     C1 is  contained in C2, then C2 cannot be contained in C1). A 1 in the
  50.     i,j place indicates Ci is contained in Cj. the (i,nContours+1) place
  51.     holds the total number of contours that contain Ci
  52.     
  53.     syntax notes:
  54.     C1 above is the "i" contour below. C2 above is the "j" contour below.
  55.     =====================================================================================
  56.     Usage Notes:
  57.     The last argument must have the storage for pointers allocated as below before
  58.     calling the function:
  59.     SCntrPtr = (ContourPath **)(malloc((size_t)(sizeof(ContourPath *)*nContours)));
  60.  
  61.     cntrList: linked list of contours containing all the information about contours.
  62.     nContours: number of contours in cntrList.
  63.     SCntrPtr: array of pointers to contours in cntrList, sorted in outer to inner order.
  64.     =====================================================================================
  65. */
  66. void sort_contourList(ContourPath *cntrList, float xxmin, float xxmax, float yymin, float yymax,
  67.             int nContours, ContourPath **SCntrPtr)
  68. {
  69. int i, j, k, jj=0;
  70. float xp=0.0, yp= 0.0;
  71. int *matrixie;        /* where we'll store the containment data (temporarily).
  72.                 matrixie (i,j) indicates whether contour i is inside countour j
  73.                 Why is this here at all ? Well, we first check to
  74.                 see if contour i is contained in contour j (for i<j).
  75.                 Then, we use the fact that two contours can not each
  76.                 contain the other to cut down on our checking time. */
  77. ContourPath *currenti, *currentj;
  78.             /* pointers to the current paths for the i and j loops */
  79.     
  80.     
  81.     /*
  82.     We've declared our variables...Now we must initialize them. 
  83.     Only the matrixie must be initialized (the others are done in loops).
  84.     */
  85.     
  86.     matrixie = (int *)calloc((size_t)((nContours+1)*nContours), (size_t)sizeof(int) );
  87.  
  88.     
  89.     currenti=cntrList;
  90.     for (i=0;i<nContours;i++) {
  91.     *(matrixie+(i+1)*(nContours+1)-1)=0; /* how many guys contain contour i */
  92.     currentj=cntrList;
  93.     /* pick the point on contour-i which is inside the domain as the point from which
  94.        we draw the line out to infinity.  Start with the x[3], y[3] in an attempt to
  95.        get a point that is not on the edge.
  96.      */
  97.     for(j=3; j<(currenti->num_pts+3); j++) {
  98.         jj = j % currenti->num_pts;
  99.         xp=currenti->x[jj];
  100.         yp=currenti->y[jj];
  101.         if(xp >= xxmin && xp <= xxmax && yp >= yymin && yp <= yymax)
  102.         break;    /* we got it */
  103.     }
  104.  
  105.     for (j=0;j<nContours;j++) {
  106.         if ((j==i) || ((j<i) && (*(matrixie+j*(nContours+1)+i)==1))) {    
  107.         currentj=currentj->next;
  108.         continue;
  109.         }
  110.         *(matrixie+i*(nContours+1)+j)=contoursub1(xp,yp,currentj,currentj->num_pts);
  111.         *(matrixie +(i+1)*(nContours+1)-1)+=*(matrixie +i*(nContours+1)+j);    
  112.         currentj=currentj->next;
  113.     }
  114.     currenti=currenti->next;
  115.     }
  116.         
  117.     
  118. /*
  119.     If all went well, we now have a nContours by nContours matrix
  120.     holding all the "containment" information. Now, we create a list telling
  121.     the drawing routine which ones to draw first (array of pointers to cntrList: SCntrPtr).
  122.     
  123.     The methodology is this:  We create a list that has the contours which
  124.     aren't contained in any other contour at the head of the list. Then come
  125.     the contours that are contained in only one other contour. Then come the
  126.     contours that are contained in only two other contours... 
  127. */
  128.  
  129.     k=0;
  130.     for (i=0;i<nContours;i++) {
  131.     currentj=cntrList;
  132.     for(j=0;j<nContours;j++) {
  133.         if (*(matrixie+(j+1)*(nContours +1)-1) == i) {
  134.         SCntrPtr[k]=currentj;
  135.         k++;
  136.         }
  137.         currentj=currentj->next;    
  138.     }
  139.     if (k==nContours) break;
  140.     }
  141.     free(matrixie);
  142.  
  143. }
  144.  
  145.  
  146.  
  147. /* See commments above: this checks to see whether the line from [xp,yp] hits 
  148.    contourj an odd or an even number of times.
  149.    Return val = 0    (xp,yp) is outside ContourPath *currentj
  150.    Return val = 1    (xp,yp) is inside  ContourPath *currentj
  151. */
  152. int contoursub1 (float xp , float yp, ContourPath *currentj , int npts)
  153. {        
  154. int k, counter=0;
  155.     for (k=1;k<=npts;k++) {
  156.         counter += contoursub2(    currentj->x[k],currentj->y[k],
  157.                     currentj->x[k-1],currentj->y[k-1],xp,yp);
  158.     }
  159.     counter += contoursub2( currentj->x[npts],currentj->y[npts],
  160.                 currentj->x[0],currentj->y[0],xp,yp);
  161.     return(counter%2);
  162.     
  163. }
  164.  
  165.  
  166.  
  167. /* solves the linear equation for contoursub1
  168.  * return val = 1: lines intersect
  169.  * return val = 0: don't intersect.
  170.  * Tests whether a line from (xp,yp) intersects line (x1,y1) --- (x2,y2).
  171.  */
  172.  
  173. int contoursub2 (double x1, double y1, double x2, double y2, double xp, double yp)
  174. {
  175. double dummy;
  176.     if(y1 == y2) return (0);
  177.     dummy=(yp-y1)/(y2-y1);
  178.     if ((0.0<=dummy) && (dummy<1.0) && ((dummy*(x2-x1)+x1-xp) >= 0.0)) return(1);
  179.     else return(0);        
  180. }
  181.  
  182.  
  183.  
  184. /* --------------------------------------------------------------------------------
  185.     Routine returns (2) if (xp,yp) is inside the
  186.     polygon x[n_path], y[n_path], (0) if outside, and (1) if exactly on the
  187.     path edge. Uses non-zero winding rule in Adobe PostScript Language
  188.     reference manual, section 4.6 on Painting. Path should have been closed
  189.     first, so that x[n_path-1] = x[0], and y[n_path-1] = y[0].
  190.     
  191.     This is version 2, trying to kill a bug in above routine:  If point is
  192.     on X edge, fails to discover that it is on edge.
  193.     
  194.     We are imagining a ray extending "up" from the point (xp,yp); the ray
  195.     has equation x = xp, y >= yp. Starting with crossing_count = 0, we add 1
  196.     each time the path crosses the ray in the +x direction, and subtract 1
  197.     each time the ray crosses in the -x direction. After traversing the
  198.     entire polygon, if (crossing_count) then point is inside.  We also watch
  199.     for edge conditions.
  200.     
  201.     If two or more points on the path have x[i] == xp, then we have an
  202.     ambiguous case, and we have to find the points in the path before and
  203.     after this section, and check them. 
  204. --------------------------------------------------------------------------------- */
  205.  
  206. int non_zero_winding(float xp, float yp, float *x, float *y, int n_path)
  207. {    
  208. int    i, j, k, jend, crossing_count, above;
  209. double    y_sect;
  210.  
  211.     above = FALSE;
  212.     crossing_count = 0;
  213.  
  214.     /* First make sure first point in path is not a special case:  */
  215.     j = jend = n_path - 1;
  216.     if (x[j] == xp) {
  217.     /* Trouble already.  We might get lucky:  */
  218.     if (y[j] == yp) return(1);
  219.  
  220.     /* Go backward down the polygon until x[i] != xp:  */
  221.     if (y[j] > yp) above = TRUE;
  222.     i = j - 1;
  223.     while (x[i] == xp && i > 0) {
  224.         if (y[i] == yp) return (1);
  225.         if (!(above) && y[i] > yp) above = TRUE;
  226.         i--;
  227.     }
  228.  
  229.     /* Now if i == 0 polygon is degenerate line x=xp;
  230.        since we know xp,yp is inside bounding box,
  231.        it must be on edge:  */
  232.     if (i == 0) return(1);
  233.  
  234.     /* Now we want to mark this as the end, for later:  */
  235.     jend = i;
  236.  
  237.     /* Now if (j-i)>1 there are some segments the point could be exactly on:  */
  238.     for (k = i+1; k < j; k++) {
  239.         if((y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp))
  240.             return (1);
  241.     }
  242.  
  243.  
  244.     /* Now we have arrived where i is > 0 and < n_path-1, and x[i] != xp.
  245.        We have been using j = n_path-1.  Now we need to move j forward 
  246.        from the origin:  */
  247.     j = 1;
  248.     while (x[j] == xp) {
  249.         if (y[j] == yp) return (1);
  250.         if (!(above) && y[j] > yp) above = TRUE;
  251.         j++;
  252.     }
  253.  
  254.     /* Now at the worst, j == jstop, and we have a polygon with only 1 vertex
  255.         not at x = xp.  But now it doesn't matter, that would end us at
  256.         the main while below.  Again, if j>=2 there are some segments to check:  */
  257.     for (k = 0; k < j-1; k++) {
  258.         if((y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp))
  259.             return (1);
  260.     }
  261.  
  262.  
  263.     /* Finally, we have found an i and j with points != xp.
  264.      * If (above) we may have crossed the ray: 
  265.      */
  266.     if(above && x[i] < xp && x[j] > xp) 
  267.         crossing_count++;
  268.     else if(above && x[i] > xp && x[j] < xp) 
  269.         crossing_count--;
  270.  
  271.     /* End nightmare scenario for x[0] == xp.  */
  272.     }
  273.     else {
  274.     /* Get here when x[0] != xp:  */
  275.     i = 0;
  276.     j = 1;
  277.     while (x[j] == xp && j < jend) {
  278.         if (y[j] == yp) return (1);
  279.         if (!(above) && y[j] > yp) above = TRUE;
  280.         j++;
  281.     }
  282.     /* Again, if j==jend, (i.e., 0) then we have a polygon with only 1 vertex
  283.        not on xp and we will branch out below.  */
  284.  
  285.     /* if ((j-i)>2) the point could be on intermediate segments:  */
  286.     for (k = i+1; k < j-1; k++) {
  287.         if ( (y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp) )
  288.             return (1);
  289.     }
  290.  
  291.     /* Now we have x[i] != xp and x[j] != xp.
  292.        If (above) and x[i] and x[j] on opposite sides, we are certain to have crossed the ray.
  293.        If not (above) and (j-i)>1, then we have not crossed it.
  294.        If not (above) and j-i == 1, then we have to check the intersection point.
  295.     */
  296.     if (x[i] < xp && x[j] > xp) {
  297.         if(above) 
  298.             crossing_count++;
  299.         else if( (j-i) == 1) {
  300.         y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
  301.         if (y_sect == yp) return (1);
  302.         if (y_sect > yp) crossing_count++;
  303.         }
  304.     }
  305.     if (x[i] > xp && x[j] < xp) {
  306.         if (above) 
  307.             crossing_count--;
  308.         else if ( (j-i) == 1) {
  309.         y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
  310.         if (y_sect == yp) return (1);
  311.         if (y_sect > yp) crossing_count--;
  312.         }
  313.     }
  314.                 
  315.     /* End easier case for x[0] != xp  */
  316.     }
  317.  
  318.     /* Now MAIN WHILE LOOP begins:
  319.         Set i = j, and search for a new j, and do as before.  */
  320.  
  321.     i = j;
  322.     while (i < jend) {
  323.     above = FALSE;
  324.     j = i+1;
  325.     while (x[j] == xp) {
  326.         if (y[j] == yp) return (1);
  327.         if (!(above) && y[j] > yp) above = TRUE;
  328.         j++;
  329.     }
  330.     /* if ((j-i)>2) the point could be on intermediate segments:  */
  331.     for (k = i+1; k < j-1; k++) {
  332.         if ( (y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp) )
  333.             return (1);
  334.     }
  335.  
  336.     /* Now we have x[i] != xp and x[j] != xp.
  337.        If (above) and x[i] and x[j] on opposite sides, we are certain to have crossed the ray.
  338.        If not (above) and (j-i)>1, then we have not crossed it.
  339.        If not (above) and j-i == 1, then we have to check the intersection point.  */
  340.  
  341.     if (x[i] < xp && x[j] > xp) {
  342.         if (above) 
  343.             crossing_count++;
  344.         else if ( (j-i) == 1) {
  345.         y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
  346.         if (y_sect == yp) return (1);
  347.         if (y_sect > yp) crossing_count++;
  348.         }
  349.     }
  350.     if (x[i] > xp && x[j] < xp) {
  351.         if (above) 
  352.         crossing_count--;
  353.         else if ( (j-i) == 1) {
  354.         y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
  355.         if (y_sect == yp) return (1);
  356.         if (y_sect > yp) crossing_count--;
  357.         }
  358.     }
  359.  
  360.     /* That's it for this piece.  Advance i:  */
  361.     i = j;
  362.     }
  363.  
  364.     /* End of MAIN WHILE.  Get here when we have gone all around without landing on edge.  */
  365.  
  366.     if (crossing_count)
  367.         return(2);
  368.     else
  369.         return(0);
  370. }
  371.