GIS编程

开发平台：
Visual C++

intersect.h：源码内容
							/*
 * FUNCTION:
 * This file contains a number of utilities useful to 3D graphics in
 * general, and to the generation of tubing and extrusions in particular
 * 
 * HISTORY:
 * Written by Linas Vepstas, August 1991
 * Updated to correctly handle degenerate cases, Linas,  February 1993 
 */
#include <math.h>
#include "port.h"
#include "vvector.h"
#define BACKWARDS_INTERSECT (2)
/* ========================================================== */
/*
 * the Degenerate_Tolerance token represents the greatest amount by
 * which different scales in a graphics environment can differ before
 * they should be considered "degenerate".   That is, when one vector is
 * a million times longer than another, changces are that the second will
 * be less than a pixel long, and therefore was probably meant to be
 * degenerate (by the CAD package, etc.)  But what should this tolerance
 * be?  At least 1 in onethousand (since screen sizes are 1K pixels), but
 * les than 1 in 4 million (since this is the limit of single-precision
 * floating point accuracy).  Of course, if double precision were used,
 * then the tolerance could be increased.
 * 
 * Potentially, this naive assumption could cause problems if the CAD
 * package attempts to zoom in on small details, and turns out, certain
 * points should not hvae been degenerate.  The problem presented here
 * is that the tolerance could run out before single-precision ran
 * out, and so the CAD packages would perceive this as a "bug".
 * One alternative is to fiddle around & try to tighten the tolerance.
 * However, the right alternative is to code the graphics pipeline in
 * double-precision (and tighten the tolerance).
 *
 * By the way, note that Degernate Tolerance is a "dimensionless"
 * quantitiy -- it has no units -- it does not measure feet, inches,
 * millimeters or pixels.  It is used only in the computations of ratios
 * and relative lengths.
 */
/* 
 * Right now, the tolerance is set to 2 parts in a million, which
 * corresponds to a 19-bit distinction of mantissas. Note that
 * single-precsion numbers have 24 bit mantissas.
 */
#define DEGENERATE_TOLERANCE   (0.000002)
/* ========================================================== */
/* 
 * The macro and subroutine INTERSECT are designed to compute the
 * intersection of a line (defined by the points v1 and v2) and a plane
 * (defined as plane which is normal to the vector n, and contains the
 * point p).  Both return the point sect, which is the point of
 * interesection.
 *
 * This MACRO attemps to be fairly robust by checking for a divide by
 * zero.
 */
/* ========================================================== */
/*
 * HACK ALERT
 * The intersection parameter t has the nice property that if t>1,
 * then the intersection is "in front of" p1, and if t<0, then the
 * intersection is "behind" p2. Unfortunately, as the intersecting plane
 * and the line become parallel, t wraps through infinity -- i.e. t can
 * become so large that t becomes "greater than infinity" and comes back 
 * as a negative number (i.e. winding number hopped by one unit).  We 
 * have no way of detecting this situation without adding gazzillions 
 * of lines of code of topological algebra to detect the winding number;
 * and this would be incredibly difficult, and ruin performance.
 * 
 * Thus, we've installed a cheap hack for use by the "cut style" drawing
 * routines. If t proves to be a large negative number (more negative
 * than -5), then we assume that t was positive and wound through
 * infinity.  This makes most cuts look good, without introducing bogus
 * cuts at infinity.
 */
/* ========================================================== */
#define INTERSECT(sect,p,n,v1,v2)			
{							
   gleDouble deno, numer, t, omt;			
							
   deno = (v1[0] - v2[0]) * n[0];			
   deno += (v1[1] - v2[1]) * n[1];			
   deno += (v1[2] - v2[2]) * n[2];			
   							
   if (deno == 0.0) {					
      VEC_COPY (n, v1);					
      /* printf ("Intersect: Warning: line is coplanar with plane n"); */ 
   } else {						
							
      numer = (p[0] - v2[0]) * n[0];			
      numer += (p[1] - v2[1]) * n[1];			
      numer += (p[2] - v2[2]) * n[2];			
							
      t = numer / deno;					
      omt = 1.0 - t;					
							
      sect[0] = t * v1[0] + omt * v2[0];		
      sect[1] = t * v1[1] + omt * v2[1];		
      sect[2] = t * v1[2] + omt * v2[2];		
   }							
}
/* ========================================================== */
/* 
 * The macro and subroutine BISECTING_PLANE compute a normal vector that
 * describes the bisecting plane between three points (v1, v2 and v3).  
 * This bisecting plane has the following properties:
 * 1) it contains the point v2
 * 2) the angle it makes with v21 == v2 - v1 is equal to the angle it 
 *    makes with v32 == v3 - v2 
 * 3) it is perpendicular to the plane defined by v1, v2, v3.
 *
 * Having input v1, v2, and v3, it returns a unit vector n.
 *
 * In some cases, the user may specify degenerate points, and still
 * expect "reasonable" or "obvious" behaviour.  The "expected"
 * behaviour for these degenerate cases is:
 *
 * 1) if v1 == v2 == v3, then return n=0
 * 2) if v1 == v2, then return v32 (normalized).
 * 3) if v2 == v3, then return v21 (normalized).
 * 4) if v1, v2 and v3 co-linear, then return v21 (normalized).
 *
 * Mathematically, these special cases "make sense" -- we just have to
 * code around potential divide-by-zero's in the code below.
 */
/* ========================================================== */
#define BISECTING_PLANE(valid,n,v1,v2,v3)			
{								
   double v21[3], v32[3];					
   double len21, len32;						
   double dot;							
								
   VEC_DIFF (v21, v2, v1);					
   VEC_DIFF (v32, v3, v2);					
								
   VEC_LENGTH (len21, v21);					
   VEC_LENGTH (len32, v32);					
								
   if (len21 <= DEGENERATE_TOLERANCE * len32) {			
								
      if (len32 == 0.0) {					
         /* all three points lie ontop of one-another */	
         VEC_ZERO (n);						
         valid = FALSE;						
      } else {							
         /* return a normalized copy of v32 as bisector */	
         len32 = 1.0 / len32;					
         VEC_SCALE (n, len32, v32);				
         valid = TRUE;						
      }								
								
   } else {							
								
      valid = TRUE;						
								
      if (len32 <= DEGENERATE_TOLERANCE * len21) {		
         /* return a normalized copy of v21 as bisector */	
         len21 = 1.0 / len21;					
         VEC_SCALE (n, len21, v21);				
								
      } else {							
								
         /* normalize v21 to be of unit length */		
         len21 = 1.0 / len21;					
         VEC_SCALE (v21, len21, v21);				
								
         /* normalize v32 to be of unit length */		
         len32 = 1.0 / len32;					
         VEC_SCALE (v32, len32, v32);				
								
         VEC_DOT_PRODUCT (dot, v32, v21);			
								
         /* if dot == 1 or -1, then points are colinear */	
         if ((dot >= (1.0-DEGENERATE_TOLERANCE)) || 		
             (dot <= (-1.0+DEGENERATE_TOLERANCE))) {		
            VEC_COPY (n, v21);					
         } else {						
   								
            /* go do the full computation */ 			
            n[0] = dot * (v32[0] + v21[0]) - v32[0] - v21[0];	
            n[1] = dot * (v32[1] + v21[1]) - v32[1] - v21[1];	
            n[2] = dot * (v32[2] + v21[2]) - v32[2] - v21[2];	
								
            /* if above if-test's passed, 			
             * n should NEVER be of zero length */		
            VEC_NORMALIZE (n);					
         } 							
      } 							
   } 								
}
/* ========================================================== */
/*
 * The block of code below is ifdef'd out, and is here for reference
 * purposes only.  It performs the "mathematically right thing" for
 * computing a bisecting plane, but is, unfortunately, subject ot noise
 * in the presence of near degenerate points.  Since computer graphics,
 * due to sloppy coding, laziness, or correctness, is filled with
 * degenerate points, we can't really use this version.  The code above
 * is far more appropriate for graphics.
 */
#ifdef MATHEMATICALLY_EXACT_GRAPHICALLY_A_KILLER
#define BISECTING_PLANE(n,v1,v2,v3)				
{								
   double v21[3], v32[3];					
   double len21, len32;						
   double dot;							
								
   VEC_DIFF (v21, v2, v1);					
   VEC_DIFF (v32, v3, v2);					
								
   VEC_LENGTH (len21, v21);					
   VEC_LENGTH (len32, v32);					
								
   if (len21 == 0.0) {						
								
      if (len32 == 0.0) {					
         /* all three points lie ontop of one-another */	
         VEC_ZERO (n);						
         valid = FALSE;						
      } else {							
         /* return a normalized copy of v32 as bisector */	
         len32 = 1.0 / len32;					
         VEC_SCALE (n, len32, v32);				
      }								
								
   } else {							
								
      /* normalize v21 to be of unit length */			
      len21 = 1.0 / len21;					
      VEC_SCALE (v21, len21, v21);				
								
      if (len32 == 0.0) {					
         /* return a normalized copy of v21 as bisector */	
         VEC_COPY (n, v21);					
      } else {							
								
         /* normalize v32 to be of unit length */		
         len32 = 1.0 / len32;					
         VEC_SCALE (v32, len32, v32);				
								
         VEC_DOT_PRODUCT (dot, v32, v21);			
								
         /* if dot == 1 or -1, then points are colinear */	
         if ((dot == 1.0) || (dot == -1.0)) {			
            VEC_COPY (n, v21);					
         } else {						
   								
            /* go do the full computation */ 			
            n[0] = dot * (v32[0] + v21[0]) - v32[0] - v21[0];	
            n[1] = dot * (v32[1] + v21[1]) - v32[1] - v21[1];	
            n[2] = dot * (v32[2] + v21[2]) - v32[2] - v21[2];	
								
            /* if above if-test's passed, 			
             * n should NEVER be of zero length */		
            VEC_NORMALIZE (n);					
         } 							
      } 							
   } 								
}
#endif
/* ========================================================== */
/*
 * This macro computes the plane perpendicular to the the plane
 * defined by three points, and whose normal vector is givven as the
 * difference between the two vectors ...
 * 
 * (See way below for the "math" model if you want to understand this.
 * The comments about relative errors above apply here.)
 */
#define CUTTING_PLANE(valid,n,v1,v2,v3)				
{								
   double v21[3], v32[3];					
   double len21, len32;						
   double lendiff;						
								
   VEC_DIFF (v21, v2, v1);					
   VEC_DIFF (v32, v3, v2);					
								
   VEC_LENGTH (len21, v21);					
   VEC_LENGTH (len32, v32);					
								
   if (len21 <= DEGENERATE_TOLERANCE * len32) {			
								
      if (len32 == 0.0) {					
         /* all three points lie ontop of one-another */	
         VEC_ZERO (n);						
         valid = FALSE;						
      } else {							
         /* return a normalized copy of v32 as cut-vector */	
         len32 = 1.0 / len32;					
         VEC_SCALE (n, len32, v32);				
         valid = TRUE;						
      }								
								
   } else {							
								
      valid = TRUE;						
								
      if (len32 <= DEGENERATE_TOLERANCE * len21) {		
         /* return a normalized copy of v21 as cut vector */	
         len21 = 1.0 / len21;					
         VEC_SCALE (n, len21, v21);				
      } else {							
								
         /* normalize v21 to be of unit length */		
         len21 = 1.0 / len21;					
         VEC_SCALE (v21, len21, v21);				
								
         /* normalize v32 to be of unit length */		
         len32 = 1.0 / len32;					
         VEC_SCALE (v32, len32, v32);				
								
         VEC_DIFF (n, v21, v32);				
         VEC_LENGTH (lendiff, n);				
								
         /* if the perp vector is very small, then the two 	
          * vectors are darn near collinear, and the cut 	
          * vector is probably poorly defined. */		
         if (lendiff < DEGENERATE_TOLERANCE) {			
            VEC_ZERO (n);					
            valid = FALSE;					
         } else {						
            lendiff = 1.0 / lendiff;				
            VEC_SCALE (n, lendiff, n);				
         } 							
      } 							
   } 								
}
/* ========================================================== */
#ifdef MATHEMATICALLY_EXACT_GRAPHICALLY_A_KILLER
#define CUTTING_PLANE(n,v1,v2,v3)		
{						
   double v21[3], v32[3];			
						
   VEC_DIFF (v21, v2, v1);			
   VEC_DIFF (v32, v3, v2);			
						
   VEC_NORMALIZE (v21);				
   VEC_NORMALIZE (v32);				
						
   VEC_DIFF (n, v21, v32);			
   VEC_NORMALIZE (n);				
}
#endif
/* ============================================================ */
/* This macro is used in several places to cycle through a series of
 * points to find the next non-degenerate point in a series */
#define FIND_NON_DEGENERATE_POINT(inext,npoints,len,diff,point_array)   
{                                                                       
   gleDouble slen;							
   gleDouble summa[3];							
   									
   do {                                                                 
      /* get distance to next point */                                  
      VEC_DIFF (diff, point_array[inext+1], point_array[inext]);        
      VEC_LENGTH (len, diff);                                           
      VEC_SUM (summa, point_array[inext+1], point_array[inext]);        
      VEC_LENGTH (slen, summa);                                         
      slen *= DEGENERATE_TOLERANCE;					
      inext ++;                                                         
   } while ((len <= slen) && (inext < npoints-1));                      
}
/* ========================================================== */
extern int bisecting_plane (gleDouble n[3],    /* returned */
                      gleDouble v1[3],  /* input */
                      gleDouble v2[3],  /* input */
                      gleDouble v3[3]);  /* input */

						