[mapserver-dev] reduce sqrt() usage in mapserver distance calculations

Bret S. Lambert bret.lambert at freecodeint.com
Fri Sep 25 06:09:20 EDT 2009


Looking through mapsearch.c, I noticed a pretty textbook case of
sqrt() overuse in distance calculation functions. Given that these
functions are called many times in tight loops, deferring the actual
calculation of the square root until necessary should reduce the
computational overhead needed.

The attached diff creates a pair of 'Quick' function versions,
which require the caller to call sqrt() on the returned value in
order to get the actual distance; however, when finding a minimum
distance, comparing the squares of the distances gives the same
result as comparing the distances themselves, so calling the
Quick version in a loop and then doing a single sqrt() on the
value found should result in a measurable improvement.

Caveats:

Tested by compilation and inspection only; I have no access to
a testbed until after the weekend, and am relying on the kindness
of strangers ;) (however, I will be available via email for the
foreseeable future, and have 6 paid days in the next 2 weeks that
I can use to hack on this for a client).

If this diff is lost/mangled by the mailer daemon, I'll find an
alternate way of getting it to the list.
-------------- next part --------------
--- ../mapserver-5.6.0-beta1/mapsearch.c	Thu Oct 23 05:19:03 2008
+++ mapsearch.c	Fri Sep 25 09:18:31 2009
@@ -295,34 +295,47 @@ int msIntersectPolygons(shapeObj *p1, shapeObj *p2) {
 double msDistancePointToPoint(pointObj *a, pointObj *b)
 {
   double d;
+
+  d = sqrt(msDistancePointToPointQuick(a, b));
+
+  return(d);
+}
+
+/*
+** Quickly compute the square of the distance; avoids expensive sqrt() call on each invocation
+** NOTE: requires caller to compute sqrt() of value
+*/
+double msDistancePointToPointQuick(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 */
+  double l_squared; /* squared length of line ab */
   double r,s;
 
-  l = msDistancePointToPoint(a,b);
+  l_squared = msDistancePointToPointQuick(a,b);
 
-  if(l == 0.0) /* a = b */
+  if(l_squared == 0.0) /* a = b */
     return( msDistancePointToPoint(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(sqrt(MS_MIN(msDistancePointToPointQuick(p, b),msDistancePointToPointQuick(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(sqrt(MS_MIN(msDistancePointToPointQuick(p, b),msDistancePointToPointQuick(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*sqrt(l_squared)));
 }
 
 #define SMALL_NUMBER 0.00000001
@@ -430,6 +443,19 @@ double msDistanceSegmentToSegment(pointObj *pa, pointO
 
 double msDistancePointToShape(pointObj *point, shapeObj *shape)
 {
+  double d;
+
+  d = msDistancePointToShapeQuick(point, shape);
+
+  return(sqrt(d));
+}
+
+/*
+** As msDistancePointToShape; avoid expensive sqrt calls
+** NOTE: requires caller to compute sqrt() of returned value
+*/
+double msDistancePointToShapeQuick(pointObj *point, shapeObj *shape)
+{
   int i, j;
   double dist, minDist=-1;
 
@@ -437,7 +463,7 @@ double msDistancePointToShape(pointObj *point, shapeOb
   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 = msDistancePointToPointQuick(point, &(shape->line[j].point[i]));
         if((dist < minDist) || (minDist < 0)) minDist = dist;
       }
     }
@@ -445,7 +471,7 @@ double msDistancePointToShape(pointObj *point, shapeOb
   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 = msDistancePointToSegmentQuick(point, &(shape->line[j].point[i-1]), &(shape->line[j].point[i]));
         if((dist < minDist) || (minDist < 0)) minDist = dist;
       }
     }
@@ -456,7 +482,7 @@ double msDistancePointToShape(pointObj *point, shapeOb
     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 = msDistancePointToSegmentQuick(point, &(shape->line[j].point[i-1]), &(shape->line[j].point[i]));
           if((dist < minDist) || (minDist < 0)) minDist = dist;
         }
       }
@@ -478,22 +504,24 @@ double msDistanceShapeToShape(shapeObj *shape1, shapeO
   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 = msDistancePointToShapeQuick(&(shape1->line[i].point[j]), shape2);
         if((dist < minDist) || (minDist < 0)) 
     minDist = dist;
       }
     }
+    minDist = sqrt(minDist); /* msDistancePointToShapeQuick() returns the square of the distance */
     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 = msDistancePointToShapeQuick(&(shape2->line[i].point[j]), shape1);
           if((dist < minDist) || (minDist < 0)) 
       minDist = dist;
         }
       }
+      minDist = sqrt(minDist); /* msDistancePointToShapeQuick() returns the square of the distance */
       break;
     case(MS_SHAPE_LINE):
       for(i=0;i<shape1->numlines;i++) {
@@ -545,11 +573,12 @@ double msDistanceShapeToShape(shapeObj *shape1, shapeO
     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 = msDistancePointToShapeQuick(&(shape2->line[i].point[j]), shape1);
           if((dist < minDist) || (minDist < 0)) 
       minDist = dist;
         }
       }
+      minDist = sqrt(minDist); /* msDistancePointToShapeQuick() returns the square of the distance */
       break;
     case(MS_SHAPE_LINE):
       /* shape1 (the polygon) could contain shape2 or one of it's parts       */
--- ../mapserver-5.6.0-beta1/mapserver.h	Thu Sep 24 01:46:24 2009
+++ mapserver.h	Fri Sep 25 11:47:23 2009
@@ -1707,8 +1707,10 @@ MS_DLL_EXPORT void msPointToFormattedString(pointObj *
 
 MS_DLL_EXPORT void msMergeRect(rectObj *a, rectObj *b);
 MS_DLL_EXPORT double msDistancePointToPoint(pointObj *a, pointObj *b);
+MS_DLL_EXPORT double msDistancePointToPointQuick(pointObj *a, pointObj *b);
 MS_DLL_EXPORT double msDistancePointToSegment(pointObj *p, pointObj *a, pointObj *b);
 MS_DLL_EXPORT double msDistancePointToShape(pointObj *p, shapeObj *shape);
+MS_DLL_EXPORT double msDistancePointToShapeQuick(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-dev mailing list