Implement deviation metric callback functionality.

Author: mikebolt

DeviationMetrics.c

@@ -0,0 +1,97 @@
+/*
+   DistanceMetrics.c
+   08/29/2015
+   Authors: Michael Casebolt, Brett Casebolt
+*/
+
+/*
+The MIT License (MIT)
+
+Copyright (c) 2015 Michael Casebolt and Brett Casebolt
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+*/
+
+// This file uses lines up to 100 characters long. If this 100 character line fits then you're good.
+
+#include "PathCompacter.h"
+
+// This one just returns the shortest distance to the infinite extension of the line segment.
+static double perpendicularDistance(DVector2D start, DVector2D end, DVector2D mid,
+                                    double dSquareSegmentLength)
+   {
+   double dArea;
+
+   dArea = 0.5 * (start.dX * (mid.dY - end.dY) +
+                  mid.dX * (end.dY - start.dY) +
+                  end.dX * (start.dY - mid.dY));
+
+   return dArea * dArea / dSquareSegmentLength;
+   }
+DeviationMetric perpendicularDistanceDeviationMetric = &perpendicularDistance;
+
+static double shortestDistanceToSegment(DVector2D start, DVector2D end, DVector2D mid,
+                                        double dSquareSegmentLength)
+   {
+   double dAX, dAY, dBX, dBY, dCX, dCY, dAdotB, dBdotC, dArea;
+
+   // Start->End forms vector A.
+   dAX = end.dX - start.dX;
+   dAY = end.dY - start.dY;
+
+   // Start->Mid forms vector B.
+   dBX = mid.dX - start.dX;
+   dBY = mid.dY - start.dY;
+
+   // End->Mid forms vector C;
+   dCX = mid.dX - end.dX;
+   dCY = mid.dY - end.dY;
+
+   // Find the dot product of A and B.
+   dAdotB = dAX * dBX + dAY * dBY;
+
+   // Find the dot product of B and C;
+   dBdotC = dBX * dCX + dBY * dCY;
+
+   // If the signs are different, the closest point on the segment is not an endpoint.
+   if (dAdotB > 0.0 && dBdotC < 0.0)
+      {
+      dArea = 0.5 * (start.dX * (mid.dY - end.dY) +
+                     mid.dX * (end.dY - start.dY) +
+                     end.dX * (start.dY - mid.dY));
+
+      return dArea * dArea / dSquareSegmentLength;
+      }
+
+   // Otherwise, figure out which endpoint it is closer to.
+   else
+      {
+      if (dAdotB < 0.0 && dBdotC < 0.0)
+         {
+         // It is closer to the start point.
+         return dAX * dAX + dAY * dAY;
+         }
+      else
+         {
+         // It is closer to the end point.
+         return dCX * dCX + dCY * dCY;
+         }
+      }
+   }
+DeviationMetric shortestDistanceToSegmentDeviationMetric = &shortestDistanceToSegment;

PathCompacter.c

@@ -55,16 +55,17 @@ typedef enum CompactPathResultCode
 #define FAILURE 0
 #define SUCCESS 1
 
-static CompactPathResultCode CompactPathSubproblemSolver(DVector2D *pPointArray,
+static CompactPathResultCode compactPathSubproblemSolver(DVector2D *pPointArray,
       int iPointsInCurrentPath, DVector2D *pResultPointArray,
-      int *piPointsInResultPath, int *piDivisionIndex, double dEpsilon);
+      int *piPointsInResultPath, int *piDivisionIndex, double dEpsilon,
+      DeviationMetric deviationMetric);
 
 // The call stack will start able to hold this many calls and grow by this amount whenever it needs
 // to grow in size.
 #define COMPACT_PATH_CALL_STACK_UNIT 2048
 
 // callStackBase is a double pointer because realloc might move the base pointer.
-static int CompactPathCallStackPush(CompactPathSubproblemCall **ppCallStackBase,
+static int compactPathCallStackPush(CompactPathSubproblemCall **ppCallStackBase,
                                     int *piCallStackCapacity, int *piNumCallsInStack,
                                     CompactPathSubproblemCall *pCall)
    {
@@ -89,7 +90,7 @@ static int CompactPathCallStackPush(CompactPathSubproblemCall **ppCallStackBase,
       return SUCCESS;
    }
 
-static int CompactPathCallStackPop(CompactPathSubproblemCall *pCallStackBase,
+static int compactPathCallStackPop(CompactPathSubproblemCall *pCallStackBase,
                                    int *piNumCallsInStack, CompactPathSubproblemCall *pPoppedCall)
    {
    // Check if there is a call to pop.
@@ -124,8 +125,9 @@ static int CompactPathCallStackPop(CompactPathSubproblemCall *pCallStackBase,
 // and resultPointArray, keeping in mind that doing so will likely alter pointArray.
 // Returns a true value (1) on successful completion, and returns a false value (0) otherwise.
 // On failure, pointsInResultPath and the contents of resultPointArray are undefined.
-int CompactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
-                DVector2D *pResultPointArray, int *piPointsInResultPath, double dEpsilon)
+int compactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
+                DVector2D *pResultPointArray, int *piPointsInResultPath, double dEpsilon,
+                DeviationMetric deviationMetric)
    {
    CompactPathSubproblemCall current, firstSubproblem, secondSubproblem;
    int iDivisionIndex; // Where should we split the problem into subproblems?
@@ -164,22 +166,22 @@ int CompactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
    current.iPointsInCurrentPath = iPointsInCurrentPath;
 
    // Add the first instance to the stack
-   if (!CompactPathCallStackPush(&pCallStackBase, &iCallStackCapacity, &iNumCallsInStack, &current))
+   if (!compactPathCallStackPush(&pCallStackBase, &iCallStackCapacity, &iNumCallsInStack, &current))
       {
       COMPACT_PATH_RETURN(FAILURE);
       }
 
    // As long as there are calls on the stack, pop one and process it.
    while (iNumCallsInStack > 0)
       {
-      if (!CompactPathCallStackPop(pCallStackBase, &iNumCallsInStack, &current))
+      if (!compactPathCallStackPop(pCallStackBase, &iNumCallsInStack, &current))
          {
          COMPACT_PATH_RETURN(FAILURE);
          }
 
-      subproblemResultCode = CompactPathSubproblemSolver(current.pPointArray,
+      subproblemResultCode = compactPathSubproblemSolver(current.pPointArray,
          current.iPointsInCurrentPath, current.pResultPointArray, &iPointsInResultPath,
-         &iDivisionIndex, dEpsilon);
+         &iDivisionIndex, dEpsilon, deviationMetric);
 
       if (subproblemResultCode == COMPACT_PATH_RESULT_CODE_DIVIDE)
          {
@@ -199,7 +201,7 @@ int CompactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
             secondSubproblem.pPointArray = current.pPointArray + iDivisionIndex;
             secondSubproblem.pResultPointArray = current.pResultPointArray + iDivisionIndex;
             secondSubproblem.iPointsInCurrentPath = current.iPointsInCurrentPath - iDivisionIndex;
-            if (!CompactPathCallStackPush(&pCallStackBase, &iCallStackCapacity,
+            if (!compactPathCallStackPush(&pCallStackBase, &iCallStackCapacity,
                                           &iNumCallsInStack, &secondSubproblem))
                {
                COMPACT_PATH_RETURN(FAILURE);
@@ -208,7 +210,7 @@ int CompactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
             firstSubproblem.pPointArray = current.pPointArray;
             firstSubproblem.pResultPointArray = current.pResultPointArray;
             firstSubproblem.iPointsInCurrentPath = iDivisionIndex + 1;
-            if (!CompactPathCallStackPush(&pCallStackBase, &iCallStackCapacity,
+            if (!compactPathCallStackPush(&pCallStackBase, &iCallStackCapacity,
                                           &iNumCallsInStack, &firstSubproblem))
                {
                COMPACT_PATH_RETURN(FAILURE);
@@ -251,12 +253,12 @@ int CompactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
 // If the result is COMPACT_PATH_RESULT_CODE_SOLVED, it means that the algorithm does not need to
 // do any further work on the subproblem, because it is already solved. In this case, divisionIndex
 // is not set.
-static CompactPathResultCode CompactPathSubproblemSolver(DVector2D *pPointArray,
+static CompactPathResultCode compactPathSubproblemSolver(DVector2D *pPointArray,
       int iPointsInCurrentPath, DVector2D *pResultPointArray,
-      int *piPointsInResultPath, int *piDivisionIndex, double dEpsilon)
+      int *piPointsInResultPath, int *piDivisionIndex, double dEpsilon,
+      DeviationMetric deviationMetric)
    {
-   double dSquareSegLen, dDX, dDY, dArea, dSquareDeviation, dMaxSquareDeviationInThisSegment,
-      ax, ay, bx, by, cx, cy;
+   double dSquareSegLen, dSquareDeviation, dMaxSquareDeviationInThisSegment, dDX, dDY;
    int i, iMaxPointIndex;
 
    // If there are fewer than three points provided, the problem is solved already.
@@ -272,22 +274,16 @@ static CompactPathResultCode CompactPathSubproblemSolver(DVector2D *pPointArray,
       }
 
    dMaxSquareDeviationInThisSegment = 0.0;
-   ax = pPointArray[0].dX;
-   ay = pPointArray[0].dY;
-   cx = pPointArray[iPointsInCurrentPath-1].dX;
-   cy = pPointArray[iPointsInCurrentPath-1].dY;
-   dDX = cx - ax;
-   dDY = cy - ay;
+
+   dDX = pPointArray[iPointsInCurrentPath - 1].dX - pPointArray[0].dX;
+   dDY = pPointArray[iPointsInCurrentPath - 1].dY - pPointArray[0].dY;
    dSquareSegLen = dDX * dDX + dDY * dDY;
 
    for (i = 1; i < iPointsInCurrentPath - 1; ++i)
       {
-      bx = pPointArray[i].dX;
-      by = pPointArray[i].dY;
-      // TODO: optionally use the endpoint distance if the triangle is skew.
-      dArea = ax * (by - cy) + bx * (cy - ay) + cx * (ay - by);
-      
-      dSquareDeviation = dArea * dArea / dSquareSegLen;
+      dSquareDeviation = deviationMetric(pPointArray[0],
+         pPointArray[iPointsInCurrentPath - 1], pPointArray[i], dSquareSegLen);
+
       if (dSquareDeviation > dMaxSquareDeviationInThisSegment)
          {
          iMaxPointIndex = i;
@@ -311,4 +307,3 @@ static CompactPathResultCode CompactPathSubproblemSolver(DVector2D *pPointArray,
       return COMPACT_PATH_RESULT_CODE_DIVIDE;
       }
    }
-

PathCompacter.h

@@ -44,8 +44,8 @@ typedef struct DVector2D
 // callback does not have to recalculate it every time it gets called.
 // These callbacks must return the square of the actual value that they intend to be compared
 // against epsilon. This is because squaring often is better than sqrt'ing sometimes.
-typedef double (*DeviationMetric)(DVector2D startOfSegment, DVector2D endOfSegment,
-                                  DVector2D point, double dSquareSegmentLength);
+typedef double (*DeviationMetric)(DVector2D /*startOfSegment*/, DVector2D /*endOfSegment*/,
+                                  DVector2D /*point*/, double /*dSquareSegmentLength*/);
 
 // This function iteratively simulates the recursive Ramer-Douglas-Peucker algorithm.
 // https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm
@@ -56,13 +56,13 @@ typedef double (*DeviationMetric)(DVector2D startOfSegment, DVector2D endOfSegme
 // and resultPointArray, keeping in mind that doing so will likely alter pointArray.
 // Returns a true value (1) on successful completion, and returns a false value (0) otherwise.
 // On failure, pointsInResultPath and the contents of resultPointArray are undefined.
-extern int compactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
-                       DVector2D *pResultPointArray, int *piPointsInResultPath,
-                       double dEpsilon, DeviationMetric deviationMetric);
+int compactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
+                DVector2D *pResultPointArray, int *piPointsInResultPath,
+                double dEpsilon, DeviationMetric deviationMetric);
 
 // Declare several metric function implementations.
 
-extern DeviationMetric perpendicularDistanceDeviationMetric;
-extern DeviationMetric shortestDistanceToSegmentDeviationMetric;
+DeviationMetric perpendicularDistanceDeviationMetric;
+DeviationMetric shortestDistanceToSegmentDeviationMetric;
 
 #endif

PathCompacterRecursive.c

@@ -41,8 +41,12 @@ THE SOFTWARE.
 // array where the next points to be solved should be placed.
 // Only CompactPath and CompactPathRecursive should touch this variable.
 static DVector2D *pCompacterLocation;
+
+// Save the deviation metric function pointer so that it doesn't have to be copied on the stack
+// a bunch of times.
+static DeviationMetric currentDeviationMetric;
 
-static void CompactPathRecursive(DVector2D *pPointArray, int iPointsInCurrentPath, double dEpsilon);
+static void compactPathRecursive(DVector2D *pPointArray, int iPointsInCurrentPath, double dEpsilon);
 
 // This function iteratively simulates the recursive Ramer-Douglas-Peucker algorithm.
 // https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm
@@ -53,8 +57,9 @@ static void CompactPathRecursive(DVector2D *pPointArray, int iPointsInCurrentPat
 // and resultPointArray, keeping in mind that doing so will likely alter pointArray.
 // Returns a true value (1) on successful completion, and returns a false value (0) otherwise.
 // This function might blow up the stack and crash your program.
-int CompactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
-                DVector2D *pResultPointArray, int *piPointsInResultPath, double dEpsilon)
+extern int compactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
+                       DVector2D *pResultPointArray, int *piPointsInResultPath,
+                       double dEpsilon, DeviationMetric deviationMetric)
    {
    int iFinalSuccessValue;
 
@@ -74,8 +79,10 @@ int CompactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
    *pResultPointArray = *pPointArray;
    pCompacterLocation = pResultPointArray + 1;
 
+   currentDeviationMetric = deviationMetric;
+
    // Launch the first call of CompactPathRecursive.
-   CompactPathRecursive(pPointArray, iPointsInCurrentPath, dEpsilon);
+   compactPathRecursive(pPointArray, iPointsInCurrentPath, dEpsilon);
 
    // Calculate and return through pointer the number of points in the resulting path.
    *piPointsInResultPath = pCompacterLocation - pResultPointArray;
@@ -92,10 +99,9 @@ int CompactPath(DVector2D *pPointArray, int iPointsInCurrentPath,
 // If the result is COMPACT_PATH_RESULT_CODE_SOLVED, it means that the algorithm does not need to
 // do any further work on the subproblem, because it is already solved. In this case, divisionIndex
 // is not set.
-static void CompactPathRecursive(DVector2D *pPointArray, int iPointsInCurrentPath, double dEpsilon)
+static void compactPathRecursive(DVector2D *pPointArray, int iPointsInCurrentPath, double dEpsilon)
    {
-   double dSquareSegLen, dDX, dDY, dArea, dSquareDeviation, dMaxSquareDeviationInThisSegment,
-      ax, ay, bx, by, cx, cy;
+   double dSquareSegLen, dDX, dDY, dArea, dSquareDeviation, dMaxSquareDeviationInThisSegment;
    int i, iMaxPointIndex;
 
    // If there are fewer than three points provided, the problem is solved already.
@@ -113,22 +119,16 @@ static void CompactPathRecursive(DVector2D *pPointArray, int iPointsInCurrentPat
       }
 
    dMaxSquareDeviationInThisSegment = 0.0;
-   ax = pPointArray[0].dX;
-   ay = pPointArray[0].dY;
-   cx = pPointArray[iPointsInCurrentPath-1].dX;
-   cy = pPointArray[iPointsInCurrentPath-1].dY;
-   dDX = cx - ax;
-   dDY = cy - ay;
+
+   dDX = pPointArray[iPointsInCurrentPath-1].dX - pPointArray[0].dX;
+   dDY = pPointArray[iPointsInCurrentPath-1].dY - pPointArray[0].dY;
    dSquareSegLen = dDX * dDX + dDY * dDY;
 
    for (i = 1; i < iPointsInCurrentPath - 1; ++i)
       {
-      bx = pPointArray[i].dX;
-      by = pPointArray[i].dY;
-      // TODO: optionally use the endpoint distance if the triangle is skew.
-      dArea = ax * (by - cy) + bx * (cy - ay) + cx * (ay - by);
+      dSquareDeviation = currentDeviationMetric(pPointArray[0],
+         pPointArray[iPointsInCurrentPath - 1], pPointArray[i], dSquareSegLen);
 
-      dSquareDeviation = dArea * dArea / dSquareSegLen;
       if (dSquareDeviation > dMaxSquareDeviationInThisSegment)
          {
          iMaxPointIndex = i;
@@ -148,9 +148,9 @@ static void CompactPathRecursive(DVector2D *pPointArray, int iPointsInCurrentPat
       {
       // Split the subproblem. Make sure to keep this left-recursive.
 
-      CompactPathRecursive(pPointArray, iMaxPointIndex + 1, dEpsilon);
+      compactPathRecursive(pPointArray, iMaxPointIndex + 1, dEpsilon);
 
-      CompactPathRecursive(pPointArray + iMaxPointIndex, iPointsInCurrentPath - iMaxPointIndex,
+      compactPathRecursive(pPointArray + iMaxPointIndex, iPointsInCurrentPath - iMaxPointIndex,
                            dEpsilon);
       }