Reducing points in polygon

Processing Coding Questions
1

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 :frowning:

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!

1 Like
2

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:

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

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.

3 Likes
3

I’ve made a Python Mode version too: :snake:

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:

337,333
337,342
333,390
332,391
328,406
328,407
328,408
314,445
314,446
291,494
291,495
289,499
276,525
272,531
271,532
266,539
261,543
259,545
257,547
256,548
253,549
252,550
247,551
234,552
222,552
221,551
219,550
219,550
159,443
154,433
151,424
150,422
139,392
135,380
134,378
124,326
122,313
118,283
118,270
118,258
119,245
121,227
126,193
129,182
132,173
133,171
154,134
208,41
220,41
220,42
277,127
286,142
304,174
319,211
320,213
320,214
330,243
330,245
330,246
334,260
334,262
337,299

simple.csv:

92,158
80,158
69,153
67,147
64,136
63,131
66,126
70,124
75,120
77,118
84,114
96,109
103,105
104,101
103,96
98,92
91,92
83,91
67,91
59,89
55,88
49,83
43,79
43,74
47,66
55,63
61,62
71,64
82,68
93,70
105,73
111,66
112,59
118,55
127,55
129,64
132,74
136,77
147,66
154,41
166,84
169,52
174,64
178,56
183,57
152,13
15,27
70,37
123,23
91,52
84,48
78,55
70,50
65,54
57,52
27,74

backtracking.csv:

14,111
18,114
22,110
25,106
30,106
37,108
44,112
48,115
56,118
57,195
61,118
69,119
73,120
82,120
83,117
88,113
90,47
92,112
99,115
104,114
111,110
113,101
111,91
112,84
115,80
190,77
116,75
117,67
120,64
117,56
116,53
115,46
114,37
110,31
10,27
56,22
74,25
89,22
91,13
107,22
110,26
115,21
117,14

diff.csv:

24,173
26,170
24,166
27,162
37,161
45,157
48,152
46,143
40,140
34,137
26,134
24,130
24,125
28,121
36,118
46,117
63,121
76,125
82,120
86,111
88,103
90,91
95,87
107,89
107,104
106,117
109,129
119,131
131,131
139,134
138,143
131,152
119,154
111,149
105,143
91,139
80,142
81,152
76,163
67,161
59,149
63,138
3 Likes
4

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?

5

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

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[]:

  1. shorten() / Reference / Processing.org
  2. concat() / Reference / Processing.org

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);

  1. List (Java SE 11 & JDK 11 )
  2. Collection (Java SE 11 & JDK 11 )
2 Likes
6

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


7

For an additional completeness’ sake, here’s the PyScript + p5*js version: :snake:
https://PyScript.com/@gotoloop/ramer-douglas-peucker

https://GoToloop.PyScriptApps.com/ramer-douglas-peucker
https://GoToloop.PyScriptApps.com/ramer-douglas-peucker/latest/js.html

1 Like