Skip to content

Instantly share code, notes, and snippets.

@wassname
Created July 12, 2015 04:20
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save wassname/39eb361450b0cf4ce8d9 to your computer and use it in GitHub Desktop.
Simplify a polygon using perpendicular distance
# -*- coding: utf-8 -*-
"""
Simplify a polygon using perpendicular distanc, see help(simplify_perp_dist)
"""
import numpy as np
import pylab
def simplify_perp_dist(M, trgt_points = np.inf):
"""
Simply a polygon based on perpendicular distance. This is not as accurate as other algorithms
but it will let us get an exact number of points and will handle large datasets with out exponential
increases in speed
http://psimpl.sourceforge.net/perpendicular-distance.html
Note that I found RDP to be better and faster than this method for large polygons.
Link: https://github.com/fhirschmann/rdp
Inputs:
:param M: an array
:type M: Nx2 numpy array
:param trgt_points: number of points to keep
:type trgt_points: int
Output:
:param Keep: a boolean array describing which points to keep
"""
# initialize
keep=np.ones(M.shape[0],dtype=np.bool)
dists=np.ones(M.shape[0])*np.inf
indices=np.arange(M.shape[0])
#evaluate initial distances
for i in xrange(1, M.shape[0]-1):
dists[i] = pylab.dist_point_to_segment(M[i], M[i-1], M[i+1])
def update_dist(l,dists,M=M,keep=keep):
"""
Update the distance matrix for point i
"""
if l==0 or l==M.shape[0]-1:
return dists
else:
kept_points=indices[keep]
a=np.searchsorted(kept_points,l)
h=kept_points[a+1]
j=kept_points[a-1]
dists[l] = pylab.dist_point_to_segment(M[l], M[h], M[j])
return dists
def update_neighbouring_dists(k,dists,M=M,keep=keep):
"""
Update the distance matrix for the neighbours around i
"""
kept_points=indices[keep]
a=np.searchsorted(kept_points,k)
h=kept_points[a]
j=kept_points[a-1]
if k != 0:
dists=update_dist(h,dists,M=M,keep=keep)
if k!=M.shape[0]-1:
dists=update_dist(j,dists,M=M,keep=keep)
return dists
# now remove points, and update neighbouring distance
# until we reach target number of points
for l in xrange(M.shape[0]-trgt_points):
k=np.argmin(dists)
keep[k]=False
dists[k]=np.inf
# now update the distances to the next and previous
# kept points
update_neighbouring_dists(k,dists,M=M,keep=keep)
return keep
import unittest
class TestHelperMethods(unittest.TestCase):
def setUp(self):
self.downhill=np.array([ 0. , 226.98202787, 239.34162503, 210.28424477,
244.00490321, 236.71432386, 80.58639256, 168.83050901,
87.30252713, 224.55499072, 70.1743241 , 137.15605142,
132.9436124 , 109.60138946, 334.73187408, 86.52304 ,
54.82528 ])
self.line=np.array([
[-1932992.26880429, -3270025.96855725],
[-1931930.09345217, -3272920.84696409],
[-1930344.97925264, -3275496.74141852],
[-1929133.93518365, -3278317.75904836],
[-1927396.82809749, -3280552.00678365],
[-1926418.42355616, -3283372.4500318 ],
[-1925121.10039399, -3281980.65474085],
[-1925313.25647357, -3279323.35443727],
[-1923839.21092072, -3278520.38315222],
[-1922721.63008425, -3279125.67676915],
[-1921923.44219151, -3278300.85260667],
[-1919433.38215566, -3278051.01533315],
[-1918060.75992851, -3275817.24249265],
[-1917607.9779197 , -3272946.305138 ],
[-1919676.42925548, -3272875.92827008],
[-1920216.44875936, -3271654.86850963],
[-1921745.32379784, -3270425.70908153]
])
def test_simplify_perp_dist(self):
keep=simplify_perp_dist(self.line,trgt_points=10)
sline=self.line[keep]
self.assertEqual(len(sline),10) # line and rec should be equal
if __name__=="__main__":
unittest.main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment