[mapserver-commits] r9353 - trunk/mapserver
svn at osgeo.org
svn at osgeo.org
Mon Sep 28 12:06:10 EDT 2009
Author: aboudreault
Date: 2009-09-28 12:06:09 -0400 (Mon, 28 Sep 2009)
New Revision: 9353
Modified:
trunk/mapserver/HISTORY.TXT
trunk/mapserver/mapprimitive.c
trunk/mapserver/mapsearch.c
trunk/mapserver/mapserver.h
Log:
Reduce use of sqrt() calls when determining distances (#3134)
Modified: trunk/mapserver/HISTORY.TXT
===================================================================
--- trunk/mapserver/HISTORY.TXT 2009-09-28 12:30:06 UTC (rev 9352)
+++ trunk/mapserver/HISTORY.TXT 2009-09-28 16:06:09 UTC (rev 9353)
@@ -14,12 +14,14 @@
Current Version (SVN trunk):
----------------------------
- - support axis odering for WFS 1.1 (#2899)
+- Reduce use of sqrt() calls when determining distances (#3134)
- - const changes to avoid warnings with msLoadProjectionString()
+- support axis odering for WFS 1.1 (#2899)
- - mapgd.c: removed unused drawVectorSymbolGD() function.
+- const changes to avoid warnings with msLoadProjectionString()
+- mapgd.c: removed unused drawVectorSymbolGD() function.
+
Version 5.6.0-beta1 (2009-09-23):
---------------------------------
Modified: trunk/mapserver/mapprimitive.c
===================================================================
--- trunk/mapserver/mapprimitive.c 2009-09-28 12:30:06 UTC (rev 9352)
+++ trunk/mapserver/mapprimitive.c 2009-09-28 16:06:09 UTC (rev 9353)
@@ -1109,10 +1109,11 @@
/* compute a distance to the polygon */
for(j=0;j<p->numlines;j++) {
for(i=1; i<p->line[j].numpoints; i++) {
- dist = msDistancePointToSegment(lp, &(p->line[j].point[i-1]), &(p->line[j].point[i]));
+ dist = msSquareDistancePointToSegment(lp, &(p->line[j].point[i-1]), &(p->line[j].point[i]));
if((dist < min_dist) || (min_dist < 0)) min_dist = dist;
}
}
+ min_dist = sqrt(min_dist);
if(min_dist > .1*MS_MAX(maxx-minx, maxy-miny))
return(MS_SUCCESS); /* point is not too close to the edge */
Modified: trunk/mapserver/mapsearch.c
===================================================================
--- trunk/mapserver/mapsearch.c 2009-09-28 12:30:06 UTC (rev 9352)
+++ trunk/mapserver/mapsearch.c 2009-09-28 16:06:09 UTC (rev 9353)
@@ -295,34 +295,50 @@
double msDistancePointToPoint(pointObj *a, pointObj *b)
{
double d;
+
+ d = sqrt(msSquareDistancePointToPoint(a, b));
+
+ return(d);
+}
+
+/*
+** Quickly compute the square of the distance; avoids expensive sqrt() call on each invocation
+*/
+double msSquareDistancePointToPoint(pointObj *a, pointObj *b)
+{
double dx, dy;
dx = a->x - b->x;
dy = a->y - b->y;
- d = sqrt(dx*dx + dy*dy);
- return(d);
+
+ return(dx*dx + dy*dy);
}
double msDistancePointToSegment(pointObj *p, pointObj *a, pointObj *b)
{
- double l; /* length of line ab */
+ return (sqrt(msSquareDistancePointToSegment(p, a, b)));
+}
+
+double msSquareDistancePointToSegment(pointObj *p, pointObj *a, pointObj *b)
+{
+ double l_squared; /* squared length of line ab */
double r,s;
- l = msDistancePointToPoint(a,b);
+ l_squared = msSquareDistancePointToPoint(a,b);
- if(l == 0.0) /* a = b */
- return( msDistancePointToPoint(a,p));
+ if(l_squared == 0.0) /* a = b */
+ return(msSquareDistancePointToPoint(a,p));
- r = ((a->y - p->y)*(a->y - b->y) - (a->x - p->x)*(b->x - a->x))/(l*l);
+ r = ((a->y - p->y)*(a->y - b->y) - (a->x - p->x)*(b->x - a->x))/(l_squared);
if(r > 1) /* perpendicular projection of P is on the forward extention of AB */
- return(MS_MIN(msDistancePointToPoint(p, b),msDistancePointToPoint(p, a)));
+ return(MS_MIN(msSquareDistancePointToPoint(p, b),msSquareDistancePointToPoint(p, a)));
if(r < 0) /* perpendicular projection of P is on the backward extention of AB */
- return(MS_MIN(msDistancePointToPoint(p, b),msDistancePointToPoint(p, a)));
+ return(MS_MIN(msSquareDistancePointToPoint(p, b),msSquareDistancePointToPoint(p, a)));
- s = ((a->y - p->y)*(b->x - a->x) - (a->x - p->x)*(b->y - a->y))/(l*l);
+ s = ((a->y - p->y)*(b->x - a->x) - (a->x - p->x)*(b->y - a->y))/l_squared;
- return(fabs(s*l));
+ return(fabs(s*s*l_squared));
}
#define SMALL_NUMBER 0.00000001
@@ -430,6 +446,18 @@
double msDistancePointToShape(pointObj *point, shapeObj *shape)
{
+ double d;
+
+ d = msSquareDistancePointToShape(point, shape);
+
+ return(sqrt(d));
+}
+
+/*
+** As msDistancePointToShape; avoid expensive sqrt calls
+*/
+double msSquareDistancePointToShape(pointObj *point, shapeObj *shape)
+{
int i, j;
double dist, minDist=-1;
@@ -437,7 +465,7 @@
case(MS_SHAPE_POINT):
for(j=0;j<shape->numlines;j++) {
for(i=0; i<shape->line[j].numpoints; i++) {
- dist = msDistancePointToPoint(point, &(shape->line[j].point[i]));
+ dist = msSquareDistancePointToPoint(point, &(shape->line[j].point[i]));
if((dist < minDist) || (minDist < 0)) minDist = dist;
}
}
@@ -445,7 +473,7 @@
case(MS_SHAPE_LINE):
for(j=0;j<shape->numlines;j++) {
for(i=1; i<shape->line[j].numpoints; i++) {
- dist = msDistancePointToSegment(point, &(shape->line[j].point[i-1]), &(shape->line[j].point[i]));
+ dist = msSquareDistancePointToSegment(point, &(shape->line[j].point[i-1]), &(shape->line[j].point[i]));
if((dist < minDist) || (minDist < 0)) minDist = dist;
}
}
@@ -456,7 +484,7 @@
else { /* treat shape just like a line */
for(j=0;j<shape->numlines;j++) {
for(i=1; i<shape->line[j].numpoints; i++) {
- dist = msDistancePointToSegment(point, &(shape->line[j].point[i-1]), &(shape->line[j].point[i]));
+ dist = msSquareDistancePointToSegment(point, &(shape->line[j].point[i-1]), &(shape->line[j].point[i]));
if((dist < minDist) || (minDist < 0)) minDist = dist;
}
}
@@ -478,22 +506,24 @@
case(MS_SHAPE_POINT): /* shape1 */
for(i=0;i<shape1->numlines;i++) {
for(j=0; j<shape1->line[i].numpoints; j++) {
- dist = msDistancePointToShape(&(shape1->line[i].point[j]), shape2);
+ dist = msSquareDistancePointToShape(&(shape1->line[i].point[j]), shape2);
if((dist < minDist) || (minDist < 0))
minDist = dist;
}
}
+ minDist = sqrt(minDist);
break;
case(MS_SHAPE_LINE): /* shape1 */
switch(shape2->type) {
case(MS_SHAPE_POINT):
for(i=0;i<shape2->numlines;i++) {
for(j=0; j<shape2->line[i].numpoints; j++) {
- dist = msDistancePointToShape(&(shape2->line[i].point[j]), shape1);
+ dist = msSquareDistancePointToShape(&(shape2->line[i].point[j]), shape1);
if((dist < minDist) || (minDist < 0))
minDist = dist;
}
}
+ minDist = sqrt(minDist);
break;
case(MS_SHAPE_LINE):
for(i=0;i<shape1->numlines;i++) {
@@ -545,11 +575,12 @@
case(MS_SHAPE_POINT):
for(i=0;i<shape2->numlines;i++) {
for(j=0; j<shape2->line[i].numpoints; j++) {
- dist = msDistancePointToShape(&(shape2->line[i].point[j]), shape1);
+ dist = msSquareDistancePointToShape(&(shape2->line[i].point[j]), shape1);
if((dist < minDist) || (minDist < 0))
minDist = dist;
}
}
+ minDist = sqrt(minDist);
break;
case(MS_SHAPE_LINE):
/* shape1 (the polygon) could contain shape2 or one of it's parts */
Modified: trunk/mapserver/mapserver.h
===================================================================
--- trunk/mapserver/mapserver.h 2009-09-28 12:30:06 UTC (rev 9352)
+++ trunk/mapserver/mapserver.h 2009-09-28 16:06:09 UTC (rev 9353)
@@ -1704,8 +1704,11 @@
MS_DLL_EXPORT void msMergeRect(rectObj *a, rectObj *b);
MS_DLL_EXPORT double msDistancePointToPoint(pointObj *a, pointObj *b);
+MS_DLL_EXPORT double msSquareDistancePointToPoint(pointObj *a, pointObj *b);
MS_DLL_EXPORT double msDistancePointToSegment(pointObj *p, pointObj *a, pointObj *b);
+MS_DLL_EXPORT double msSquareDistancePointToSegment(pointObj *p, pointObj *a, pointObj *b);
MS_DLL_EXPORT double msDistancePointToShape(pointObj *p, shapeObj *shape);
+MS_DLL_EXPORT double msSquareDistancePointToShape(pointObj *p, shapeObj *shape);
MS_DLL_EXPORT double msDistanceSegmentToSegment(pointObj *pa, pointObj *pb, pointObj *pc, pointObj *pd);
MS_DLL_EXPORT double msDistanceShapeToShape(shapeObj *shape1, shapeObj *shape2);
MS_DLL_EXPORT int msIntersectSegments(pointObj *a, pointObj *b, pointObj *c, pointObj *d);
More information about the mapserver-commits
mailing list