distance.c

Go to the documentation of this file.
00001 #include <math.h>
00002 #include <grass/gis.h>
00003 static double min4(double, double, double, double);
00004 static double min2(double, double);
00005 static double dabs (double);
00006 
00007 /* WARNING: this code is preliminary and may be changed,
00008  * including calling sequences to any of the functions
00009  * defined here
00010  */
00011 
00012 static int projection = 0;
00013 static double factor = 1.0;
00014 
00015 
00028 int G_begin_distance_calculations()
00029 {
00030     double a, e2;
00031 
00032     factor = 1.0;
00033     switch (projection = G_projection())
00034     {
00035     case PROJECTION_LL:
00036         G_get_ellipsoid_parameters (&a, &e2);
00037         G_begin_geodesic_distance (a, e2);
00038         return 2;
00039     default:
00040         factor = G_database_units_to_meters_factor();
00041         if (factor <= 0.0)
00042         {
00043             factor = 1.0;          /* assume meter grid */
00044             return 0;
00045         }
00046         return 1;
00047     }
00048 }
00049 
00050 
00066 double G_distance (double e1,double n1,double e2,double n2)
00067 {
00068     double hypot();
00069 
00070     if (projection == PROJECTION_LL)
00071         return G_geodesic_distance (e1, n1, e2, n2);
00072     else
00073         return factor * hypot (e1-e2,n1-n2);
00074 }
00075 
00076 double
00077 G_distance_between_line_segments (
00078     double ax1,double ay1,
00079     double ax2,double ay2,
00080     double bx1,double by1,
00081     double bx2,double by2)
00082 {
00083     double ra, rb;
00084     double x, y;
00085 
00086 /* if the segments intersect, then the distance is zero */
00087     if(G_intersect_line_segments(ax1,ay1,ax2,ay2,
00088                                bx1,by1,bx2,by2,
00089                                &ra,&rb,&x,&y) > 0)
00090             return 0.0;
00091     return min4 (
00092         G_distance_point_to_line_segment (ax1,ay1,bx1,by1,bx2,by2) ,
00093         G_distance_point_to_line_segment (ax2,ay2,bx1,by1,bx2,by2) ,
00094         G_distance_point_to_line_segment (bx1,by1,ax1,ay1,ax2,ay2) ,
00095         G_distance_point_to_line_segment (bx2,by2,ax1,ay1,ax2,ay2)
00096                 ) ;
00097 }
00098 
00099 double
00100 G_distance_point_to_line_segment (
00101     double xp,double yp,         /* the point */
00102     double x1,double y1,double x2,double y2)   /* the line segment */
00103 {
00104     double dx,dy;
00105     double x,y;
00106     double xq,yq,ra,rb;
00107     int t;
00108 
00109 /* define the perpendicular to the segment thru the point */
00110     dx = x1 - x2;
00111     dy = y1 - y2;
00112 
00113     if (dx == 0.0 && dy == 0.0)
00114         return G_distance (x1,y1,xp,yp);
00115 
00116     if (dabs(dy) > dabs(dx))
00117     {
00118         xq = xp + dy;
00119         yq = (dx/dy) * (xp - xq) + yp;
00120     }
00121     else
00122     {
00123         yq = yp + dx;
00124         xq = (dy/dx) * (yp - yq) + xp;
00125     }
00126 
00127 /* find the intersection of the perpendicular with the segment */
00128     switch(t=G_intersect_line_segments(xp,yp,xq,yq, x1,y1,x2,y2, &ra,&rb,&x,&y))
00129     {
00130     case 0:
00131     case 1:
00132         break;
00133     default:
00134         /* parallel/colinear cases shouldn't occur with perpendicular lines */
00135         fprintf (stderr,"G_distance_point_to_line_segment: shouldn't happen\n");
00136 fprintf(stderr," code=%d P=(%f,%f) S=(%f,%f)(%f,%f)\n",t,xp,yp,x1,y1,x2,y2);
00137         return -1.0;
00138     }
00139 
00140 /* if x,y lies on the segment, then the distance is from to x,y */
00141     if (rb >= 0 && rb <= 1.0)
00142         return G_distance (x,y,xp,yp);
00143 
00144 /* otherwise the distance is the short of the distances to the endpoints
00145  * of the segment
00146  */
00147     return min2(G_distance (x1,y1,xp,yp), G_distance (x2,y2,xp,yp));
00148 }
00149 
00150 static double dabs ( double x)
00151 {
00152     return x < 0 ? -x : x;
00153 }
00154 
00155 static double min4(double a,double b,double c,double d)
00156 {
00157     return min2 (min2(a,b), min2(c,d));
00158 }
00159 
00160 static double min2(double a,double b)
00161 {
00162     return a < b ? a : b;
00163 }

Generated on Wed Dec 19 14:59:05 2007 for GRASS by  doxygen 1.5.4