Reducing points in polygon

jeffthompson · August 25, 2019, 3:03pm

I’ve been working with this example code over on the old forum, which reduces the number of points in a polygon using the Ramer-Douglas-Peucker algorithm. It works great for my purposes, but for one weird thing: it creates duplicates of all the points

Here’s my demo code; simplification functions below with a few notes:

int numPoints = 12;     // reduce to this many points

float[][] pointsArray = { 
  { 1687, 1667 }, { 1684, 1708 }, { 1663, 1948 }, { 1662, 1955 }, { 1642, 2030 }, { 1640, 2037 }, { 1639, 2040 }, { 1572, 2225 }, { 1570, 2230 }, { 1457, 2472 }, { 1456, 2474 }, { 1444, 2497 }, { 1378, 2623 }, { 1359, 2656 }, { 1355, 2662 }, { 1332, 2694 }, { 1307, 2717 }, { 1297, 2725 }, { 1283, 2736 }, { 1280, 2738 }, { 1265, 2747 }, { 1261, 2749 }, { 1235, 2756 }, { 1172, 2759 }, { 1110, 2759 }, { 1104, 2756 }, { 1096, 2751 }, { 1095, 2750 }, { 796, 2213 }, { 772, 2164 }, { 753, 2118 }, { 749, 2108 }, { 695, 1960 }, { 674, 1899 }, { 672, 1891 }, { 622, 1631 }, { 610, 1567 }, { 588, 1414 }, { 588, 1350 }, { 591, 1292 }, { 595, 1224 }, { 603, 1134 }, { 628, 964 }, { 643, 911 }, { 660, 863 }, { 664, 854 }, { 768, 670 }, { 1039, 207 }, { 1099, 207 }, { 1101, 209 }, { 1386, 634 }, { 1431, 711 }, { 1522, 870 }, { 1597, 1055 }, { 1600, 1063 }, { 1602, 1069 }, { 1648, 1217 }, { 1650, 1224 }, { 1651, 1228 }, { 1669, 1302 }, { 1670, 1312 }, { 1687, 1495 } 
};


void setup() {
  size(800, 600);
  background(255);

  // scale points to fit screen
  ArrayList<PVector> points = new ArrayList<PVector>();
  for (int i=0; i<pointsArray.length; i++) {
    float x = pointsArray[i][0] * 0.2;
    float y = pointsArray[i][1] * 0.2;
    points.add( new PVector(x,y) );
  }
  println("original pts:\t\t" + points.size());

  // draw original polygon
  noFill();
  stroke(0, 150, 255);
  strokeWeight(2);
  beginShape();
  for (PVector point : points) {
    vertex(point.x, point.y);
  }
  endShape(CLOSE);

  fill(0, 100);
  noStroke();
  for (PVector point : points) {
    ellipse(point.x, point.y, 15, 15);
  }
  
  // reduce the number of points in the polygon until it reaches the desired number
  int epsilon = 0;
  ArrayList<PVector> smoothedPoints = douglasPeucker(points, epsilon);
  while (true) {
    epsilon += 1;
    smoothedPoints = douglasPeucker(points, epsilon);
    println("- epsilon " + epsilon + ":\t\t" + smoothedPoints.size());
    if (smoothedPoints.size() <= numPoints) {
      break;
    }
  }
  
  // draw the result
  translate(width/2, 0);
  noFill();
  stroke(0, 150, 255);
  strokeWeight(2);
  beginShape();
  for (PVector point : smoothedPoints) {
    vertex(point.x, point.y);
  }
  endShape(CLOSE);

  fill(0, 100);
  noStroke();
  for (PVector point : smoothedPoints) {
    ellipse(point.x, point.y, 15, 15);
  }
}

And here’s the reduction code – my messy solution is noted about halfway down, where basically I just return the first half of the array list.

ArrayList<PVector> douglasPeucker(ArrayList<PVector> points, float epsilon) {
  return douglasPeucker(points, 0, points.size()-1, epsilon);
}

ArrayList<PVector> douglasPeucker(ArrayList<PVector> points, int startIndex, int endIndexInc, float epsilon) {
  
  // find point with max dist
  float dMax = 0;
  int index = -1;
  for (int i = startIndex+1; i <= endIndexInc; i++) {
    PVector current = points.get(i);
    PVector start =   points.get(startIndex);
    PVector end =     points.get(endIndexInc);
    float d = distToSegmentSquared(current.x, current.y, start.x, start.y, end.x, end.y);
    if (d > dMax) {
      index = i;
      dMax = d;
    }
  }
  dMax = sqrt(dMax);

  // if it's greater we simplify
  if (dMax > epsilon) {
    ArrayList<PVector> resultFront = douglasPeucker(points, startIndex, index, epsilon);
    ArrayList<PVector> resultBack =  douglasPeucker(points, index, endIndexInc, epsilon);

    // combine
    ArrayList<PVector> result = new ArrayList<PVector>(resultFront);
    result.addAll(resultBack);
    return result;
  } 
  else {
    ArrayList<PVector> result = new ArrayList<PVector>();
    result.add( points.get(startIndex) );
    result.add( points.get(endIndexInc) );

    // ***** PROBLEM HERE!!!! ******
    // for some reason gives us duplicates of all the points,
    // so we only return every other one
    ArrayList<PVector> out = new ArrayList<PVector>();
    for (int i=0; i<result.size()/2; i++) {
      out.add(result.get(i));
    }
    return out;
  }
}

float dist2(float x1, float y1, float x2, float y2) { 
  return sq(x1 - x2) + sq(y1 - y2);
}

float distToSegment(float px, float  py, float lx1, float  ly1, float lx2, float ly2) { 
  return sqrt(distToSegmentSquared(px, py, lx1, ly1, lx2, ly2));
}

// inspired by stackoverflow.com/a/1501725/1022707
float distToSegmentSquared(float px, float py, float lx1, float ly1, float lx2, float ly2) {
  float lineDist = dist2(lx1, ly1, lx2, ly2);
  if (lineDist == 0) return dist2(px, py, lx1, ly1);
  float t = ((px - lx1) * (lx2 - lx1) + (py - ly1) * (ly2 - ly1)) / lineDist;
  if (t < 0) return dist2(px, py, lx1, ly1);
  if (t > 1) return dist2(px, py, lx2, ly2);
  return dist2(px, py, lx1 + t * (lx2 - lx1), ly1 + t * (ly2 - ly1));
}

I’m pretty lost as to how this is working, so any help would be much appreciated!

GoToLoop · August 28, 2019, 8:45am

Here’s my take on converting the code from the Wikipedia article you’ve posted:

And its JS implementation: https://Karthaus.nl/rdp/js/rdp2.js

Here’s the online code:

And its corresponding gist link:

gist.github.com

https://gist.github.com/GoToLoop/a5db257be4d7756a00220a3e97066dd5

.block

height: 600
scrolling: no
border: yes
license: cc-by-4.0

Calc.pde

static final String[] ERRORS = {
  "Polygon's newSize parameter can't be less than 2!",
  "Parameter dots[] must contain at least 2 PVector objects!",
  "Parameter epsilon can't be a negative value!"
};

static final PVector[] reducePolygon(final PVector[] dots, final int newSize) {
  if (newSize < 2)  throw new Error(ERRORS[0]);

  PVector[] reduced;

This file has been truncated. show original

Data.pde

static final String FILENAME = "xy.csv", DELIM = ",";

PVector[] loadCoordsXYAsPVectorsArray() {
  final String[] lines = loadStrings(FILENAME);
  final PVector[] vecs = new PVector[lines.length];

  for (int i = 0; i < lines.length; ++i) {
    final PVector v = vecs[i] = new PVector();
    v.set(float(lines[i].split(DELIM)));
  }

This file has been truncated. show original

There are more than three files. show original

It’s also hosted on OpenProcessing too. But it’s glitched there, probly due to the do/while loop syntax:
OpenProcessing.org/sketch/747172

In order to run it, go to the right panel and switch on & off the “Infinite Loop Protection” option.

In my attempt, it takes 9 calls to douglasPeucker() in order to reduce the array’s length to 12

Your version takes 1 less: 8 calls.

Also, my array’s length slightly varies on each epsilon value compared to yours.

GoToLoop · August 28, 2019, 8:16pm

I’ve made a Python Mode version too:

RamerDouglasPeucker.pyde:

"""
 * Ramer-Douglas-Peucker Algorithm (v1.2.3)
 * GoToLoop & JeffThompson (2019/Aug/27)
 *
 * https://Discourse.Processing.org/t/reducing-points-in-polygon/13590/3
 *
 * https://en.Wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
 * https://Karthaus.nl/rdp/js/rdp2.js
 *
 * https://Bl.ocks.org/GoToLoop/a5db257be4d7756a00220a3e97066dd5
 * https://www.OpenProcessing.org/sketch/747172
"""

from Calc import reducePolygon
from Data import loadCoordsXYAsPVectorsTuple

QTY, DIAM, BOLD = 12, 15, 2;
BG, FG, STROKE = -1, 0x64000000, 0xff0096FF

DEBUG = True

def setup():
    size(800, 600)
    smooth(), noLoop(), strokeWeight(BOLD)

    global coords, reduced

    coords = loadCoordsXYAsPVectorsTuple()
    if DEBUG: print "Polygon's original size: %d" % len(coords)

    reduced = reducePolygon(coords, QTY, DEBUG = DEBUG)
    if DEBUG: print "Polygon's final size: %d\n" % len(reduced)

    if DEBUG:
        txt, digits = '%s: %s', len(`max(0, len(reduced) - 1)`)
        for i, v in enumerate(reduced): print txt % (`i`.zfill(digits), v)


def draw():
    background(BG)
    drawPoints(coords)

    translate(width >> 1, 0)
    drawPoints(reduced)


def drawPoints(dots):
    if not dots: return

    noFill()
    stroke(STROKE)

    with beginClosedShape():
        for v in dots: vertex(v.x, v.y)

    fill(FG)
    noStroke()

    for v in dots: circle(v.x, v.y, DIAM)

Calc.py:

ERRORS = (
    "Polygon's newSize parameter can't be less than 2!",
    "Parameter dots[] must contain at least 2 PVector objects!",
    "Parameter epsilon can't be a negative value!"
)

MSG = 'Step: %d\t\tSize: %d'

def reducePolygon(dots, newSize, ep=0, DEBUG=True):
    if newSize < 2: raise ValueError(ERRORS[0])

    while True:
        reduced = shortDistRDP(dots, ep)
        currLen = len(reduced)

        if DEBUG: print MSG % (ep, currLen)

        if currLen <= newSize: return reduced
        ep += 1


def shortDistRDP(dots, epsilon):
    if not dots or len(dots) < 2: raise ValueError(ERRORS[1])
    if epsilon < 0: raise ValueError(ERRORS[2])

    qty = len(dots)
    if (qty == 2): return dots[:]

    head = dots[0]
    tail = dots[-1]
    last = qty - 1

    maxDistSq = 0
    idx = -1

    for i in xrange(1, last):
        currDistSq = dotToLineDistSq(dots[i], head, tail)

        if (currDistSq > maxDistSq):
            maxDistSq = currDistSq
            idx = i

    if sqrt(maxDistSq) > epsilon:
        recurResL = shortDistRDP(dots[:idx+1], epsilon)
        recurResR = shortDistRDP(dots[idx:], epsilon)

        return recurResL[:-1] + recurResR

    return dots[::last]


def dotToLineDistSq(p, v, w, z=PVector()):
    lineLen = vec2dDistSq(v, w)
    if not lineLen: return vec2dDistSq(p, v)

    lineX = w.x - v.x
    lineY = w.y - v.y
    t = ((p.x - v.x)*lineX + (p.y - v.y)*lineY) / lineLen

    if t < 0: return vec2dDistSq(p, v)
    if t > 1: return vec2dDistSq(p, w)

    return vec2dDistSq(p, z.set(v).add(t * lineX, t * lineY, 0))


def vec2dDistSq(a, b): return sq(a.x - b.x) + sq(a.y - b.y)

Data.py:

# 'xy.csv', 'simple.csv', 'backtracking.csv', 'diff.csv'

def loadCoordsXYAsPVectorsTuple(FILE='xy.csv', SEP=',', f=PApplet.parseFloat):
    return tuple(PVector().set(f(r.split(SEP))) for r in loadStrings(FILE))

xy.csv:

simple.csv:

backtracking.csv:

diff.csv:

jeffthompson · August 27, 2019, 12:26pm

Thanks @GoToLoop! I need to dig into your version, but with a quick look it seems that you have the two parts building the final point list, like in the version I posted:

recurRes1 = douglasPeucker((PVector[]) subset(dots, 0, idx + 1), epsilon), 
recurRes2 = douglasPeucker((PVector[]) subset(dots, idx), epsilon);

That was the part that seemed to give me duplicates – could you explain what’s going on here?

GoToLoop · August 28, 2019, 8:48am

This is your version for the split recursion + container concatenation block:

jeffthompson:

if (dMax > epsilon) {
  ArrayList<PVector> resultFront = douglasPeucker(points, startIndex, index, epsilon);
  ArrayList<PVector> resultBack = douglasPeucker(points, index, endIndexInc, epsilon);

  ArrayList<PVector> result = new ArrayList<PVector>(resultFront);
  result.addAll(resultBack);

  return result;
}

Below’s the corresponding original pseudocode block excerpt:

if ( dmax > epsilon ) {
  recResults1[] = DouglasPeucker(PointList[1...index], epsilon)
  recResults2[] = DouglasPeucker(PointList[index...end], epsilon)

  ResultList[] = {recResults1[1...length(recResults1)-1], recResults2[1...length(recResults2)]
}

And here’s my attempt conversion:

  if (sqrt(dMax) > epsilon) {
    final PVector[]
      recurResL = shortDistRDP((PVector[]) subset(dots, 0, idx + 1), epsilon), 
      recurResR = shortDistRDP((PVector[]) subset(dots, idx), epsilon);

    return (PVector[]) concat(shorten(recurResL), recurResR);
  }

As you can notice, I call shorten() in order to remove the tail element from container recurResL[] before concat() it w/ recurResR[]:

I did so b/c it seems to me the original does that: recResults1[1...length(recResults1)-1]

The pseudocode’s wiki states it’s using containers w/ their starting index = 1.
So index 1 is actually index 0.

The ... denotes a range, where both its start & end indices are inclusive, apparently.

B/c it subtracts 1 from the recResults1[] container’s length, its tail element is cut off.

At your solution here: ArrayList<PVector> result = new ArrayList<PVector>(resultFront);

You initialize another ArrayList w/ the full contents of resultFront.

I believe you need to remove() its tail before invoking addAll() for the 2nd container:
result.remove(result.size() - 1);

GoToLoop · August 28, 2019, 10:48pm

For completeness’ sake, I’ve ended up making a p5js version too:

gist.github.com

https://gist.github.com/GoSubRoutine/da4939559d8b786df13f5694ea2edd30

.block

height: 600
scrolling: no
border: yes
license: cc-by-4.0

calc.mjs

const ERRORS = [
  "Polygon's newSize parameter can't be less than 2!",
  "Parameter dots[] must contain at least 2 { x, y } objects!",
  "Parameter epsilon can't be a negative value!"
];

export default function reducePolygon(dots, newSize, debug=true) {
  if (newSize < 2)  throw RangeError(ERRORS[0]);

  let reduced, ep = 0;

This file has been truncated. show original

data.mjs

export default function mapXYTableToVectorArray(t) {
  return t.getArray().map(([ x, y ]) => createVector(+x, +y));
}

There are more than three files. show original