From: David Blasby Date: Wed, 4 Feb 2004 02:53:20 +0000 (+0000) Subject: applied patricia tozer's patch (distance function was taking acos of something X-Git-Tag: pgis_0_8_2~120 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=640943a3b05156c85ba3cd2d8606a870f95d487d;p=postgis applied patricia tozer's patch (distance function was taking acos of something just slightly outside [-1,1]). git-svn-id: http://svn.osgeo.org/postgis/trunk@439 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/postgis_proj.c b/postgis_proj.c index d76f3db8c..e3cf68045 100644 --- a/postgis_proj.c +++ b/postgis_proj.c @@ -8,9 +8,13 @@ * * This is free software; you can redistribute and/or modify it under * the terms of hte GNU General Public Licence. See the COPYING file. - * + * ********************************************************************** * $Log$ + * Revision 1.9 2004/02/04 02:53:20 dblasby + * applied patricia tozer's patch (distance function was taking acos of something + * just slightly outside [-1,1]). + * * Revision 1.8 2003/12/04 18:58:35 dblasby * changed david skae to skea * @@ -50,9 +54,9 @@ -//use the WKT definition of an ellipsoid -// ie. SPHEROID["name",A,rf] or SPHEROID("name",A,rf) -// SPHEROID["GRS_1980",6378137,298.257222101] +//use the WKT definition of an ellipsoid +// ie. SPHEROID["name",A,rf] or SPHEROID("name",A,rf) +// SPHEROID["GRS_1980",6378137,298.257222101] // wkt says you can use "(" or "[" PG_FUNCTION_INFO_V1(ellipsoid_in); @@ -66,11 +70,11 @@ Datum ellipsoid_in(PG_FUNCTION_ARGS) memset(sphere,0, sizeof(SPHEROID)); - if (strstr(str,"SPHEROID") != str ) + if (strstr(str,"SPHEROID") != str ) { elog(ERROR,"SPHEROID parser - doesnt start with SPHEROID"); pfree(sphere); - PG_RETURN_NULL(); + PG_RETURN_NULL(); } nitems = sscanf(str,"SPHEROID[\"%19[^\"]\",%lf,%lf]",sphere->name,&sphere->a,&rf); @@ -91,7 +95,7 @@ Datum ellipsoid_in(PG_FUNCTION_ARGS) sphere->e = sqrt(sphere->e_sq); PG_RETURN_POINTER(sphere); - + } @@ -99,7 +103,7 @@ Datum ellipsoid_in(PG_FUNCTION_ARGS) PG_FUNCTION_INFO_V1(ellipsoid_out); Datum ellipsoid_out(PG_FUNCTION_ARGS) { - SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(0); + SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(0); char *result; result = palloc(MAX_DIGS_DOUBLE + MAX_DIGS_DOUBLE + 20 + 9 + 2); @@ -123,12 +127,12 @@ double deltaLongitude(double azimuth, double sigma, double tsm,SPHEROID *sphere) das = cos(azimuth)*cos(azimuth); C = sphere->f/16.0 * das * (4.0 + sphere->f * (4.0 - 3.0 * das)); // compute the difference in longitude - + ctsm = cos(tsm); DL = ctsm + C * cos(sigma) * (-1.0 + 2.0 * ctsm*ctsm); DL = sigma + C * sin(sigma) * DL; return (1.0 - C) * sphere->f * sin(azimuth) * DL; -} +} //support function for distance calc @@ -142,7 +146,7 @@ double mu2(double azimuth,SPHEROID *sphere) e2 = sqrt(sphere->a*sphere->a-sphere->b*sphere->b)/sphere->b; return cos(azimuth)*cos(azimuth) * e2*e2; -} +} //support function for distance calc @@ -153,7 +157,7 @@ double mu2(double azimuth,SPHEROID *sphere) double bigA(double u2) { return 1.0 + u2/256.0 * (64.0 + u2 * (-12.0 + 5.0 * u2)); -} +} //support function for distance calc @@ -164,7 +168,7 @@ double bigA(double u2) double bigB(double u2) { return u2/512.0 * (128.0 + u2 * (-64.0 + 37.0 * u2)); -} +} //given 2 lat/longs and ellipse, find the distance @@ -184,7 +188,9 @@ double distance_ellipse(double lat1, double long1, double u2,A,B; double dsigma; - int iterations; + double TEMP; + + int iterations; L1 = atan((1.0 - sphere->f ) * tan( lat1) ); L2 = atan((1.0 - sphere->f ) * tan( lat2) ); @@ -202,7 +208,21 @@ double distance_ellipse(double lat1, double long1, cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1; sigma = acos(cosSigma); azimuthEQ = asin((cosU1 * cosU2 * sindl1)/sin(sigma)); - tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ))); + + // patch from patrica tozer to handle minor mathematical stability problem + TEMP = cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ)); + if(TEMP > 1) + { + TEMP = 1; + } + else if(TEMP < -1) + { + TEMP = -1; + } + tsm = acos(TEMP); + + + //tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ))); dl2 = deltaLongitude(azimuthEQ, sigma, tsm,sphere); dl3 = dl1 - (dl + dl2); dl1 = dl + dl2; @@ -211,12 +231,12 @@ double distance_ellipse(double lat1, double long1, iterations++; } while ( (iterations<999) && (fabs(dl3) > 1.0e-032)); - // compute expansions A and B + // compute expansions A and B u2 = mu2(azimuthEQ,sphere); A = bigA(u2); B = bigB(u2); - // compute length of geodesic + // compute length of geodesic dsigma = B * sin(sigma) * (cos(tsm) + (B*cosSigma*(-1.0 + 2.0 * (cos(tsm)*cos(tsm))))/4.0); return sphere->b * (A * (sigma - dsigma)); } @@ -237,7 +257,7 @@ double length2d_ellipse_linestring(LINE3D *line, SPHEROID *sphere) for (i=1; inpoints;i++) { to = &line->points[i]; - + dist += distance_ellipse(frm->y*M_PI/180.0 , frm->x*M_PI/180.0, to->y*M_PI/180.0 , to->x*M_PI/180.0, sphere); @@ -262,7 +282,7 @@ double length3d_ellipse_linestring(LINE3D *line, SPHEROID *sphere) for (i=1; inpoints;i++) { to = &line->points[i]; - + dist_ellipse = distance_ellipse(frm->y*M_PI/180.0 , frm->x*M_PI/180.0, to->y*M_PI/180.0 , to->x*M_PI/180.0, sphere); @@ -281,7 +301,7 @@ double length3d_ellipse_linestring(LINE3D *line, SPHEROID *sphere) // find the "length of a geometry" // length2d(point) = 0 // length2d(line) = length of line -// length2d(polygon) = 0 +// length2d(polygon) = 0 // uses ellipsoidal math to find the distance //// x's are longitude, and y's are latitude - both in decimal degrees @@ -289,8 +309,8 @@ PG_FUNCTION_INFO_V1(length_ellipsoid); Datum length_ellipsoid(PG_FUNCTION_ARGS) { GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1); - + SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1); + int32 *offsets1; char *o1; int32 type1,j; @@ -304,7 +324,7 @@ Datum length_ellipsoid(PG_FUNCTION_ARGS) for (j=0; j< geom->nobjs; j++) //for each object in geom1 { - o1 = (char *) geom +offsets1[j] ; + o1 = (char *) geom +offsets1[j] ; type1= geom->objType[j]; if (type1 == LINETYPE) //LINESTRING { @@ -322,7 +342,7 @@ Datum length_ellipsoid(PG_FUNCTION_ARGS) // find the "length of a geometry" // length3d(point) = 0 // length3d(line) = length of line -// length3d(polygon) = 0 +// length3d(polygon) = 0 // uses ellipsoidal math to find the distance on the XY plane, then // uses simple Pythagoras theorm to find the 3d distance on each segment // x's are longitude, and y's are latitude - both in decimal degrees @@ -331,8 +351,8 @@ PG_FUNCTION_INFO_V1(length3d_ellipsoid); Datum length3d_ellipsoid(PG_FUNCTION_ARGS) { GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1); - + SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1); + int32 *offsets1; char *o1; int32 type1,j; @@ -346,7 +366,7 @@ Datum length3d_ellipsoid(PG_FUNCTION_ARGS) for (j=0; j< geom->nobjs; j++) //for each object in geom1 { - o1 = (char *) geom +offsets1[j] ; + o1 = (char *) geom +offsets1[j] ; type1= geom->objType[j]; if (type1 == LINETYPE) //LINESTRING { @@ -362,7 +382,7 @@ Datum length3d_ellipsoid(PG_FUNCTION_ARGS) PG_FUNCTION_INFO_V1(distance_ellipsoid); Datum distance_ellipsoid(PG_FUNCTION_ARGS) { - SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(2); + SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(2); GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); @@ -377,7 +397,7 @@ Datum distance_ellipsoid(PG_FUNCTION_ARGS) elog(ERROR,"optimistic_overlap:Operation on two GEOMETRIES with different SRIDs\n"); PG_RETURN_NULL(); } - + if (geom1->type != POINTTYPE) { elog(ERROR,"optimistic_overlap: first arg isnt a point\n");