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
00008
00009
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;
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
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,
00102 double x1,double y1,double x2,double y2)
00103 {
00104 double dx,dy;
00105 double x,y;
00106 double xq,yq,ra,rb;
00107 int t;
00108
00109
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
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
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
00141 if (rb >= 0 && rb <= 1.0)
00142 return G_distance (x,y,xp,yp);
00143
00144
00145
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 }