Simplify a polygon using perpendicular distance
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# -*- 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