Select to view content in your preferred language

Sorting points and sweepline

84
0
Thursday
DanPatterson
MVP Esteemed Contributor
2 0 84

Take two polygons, an area of interest rectangle/square and an overlay polygon with the same extent as the first polygon, but rotated within it's bounds.  Now this could be the basis of clipping or more generally boolean operations on polygons.  From this data we can determine all of the associated overlay cases, like symmetrical difference, intersection, union, etcetera.

aoi_aoi0_sweep_example.pngPoints 0, 1, 2, 3 belong to the overlay polygon.  A 4th point is the closing duplicate of point 0.

Points 5, 6, 7, 8 belong to the polygon on the bottom.  Again, a 9th point is a duplicate of the first.

 

As arrays, you can slice the first and last values using [0, -1] where -1 means count backwards to the last.

 

 

 

aoi[[0, -1]], aoi0[[0, -1]]  # slice the first and last points

(array([[  0.00,   0.00],
        [  0.00,   0.00]]),
 array([[  0.00,   4.00],
        [  0.00,   4.00]]))

When these shapes are intersected you end up with a combined geometry with each part consisting of from-to points as shown below (extra information on clockwise/counterclockwise and part and/or bit IDs are included)

prn_tbl(geom2)
...    OID_    Fr_pnt    To_pnt    CW/CCW    Part_ID    Bit_ID  
----------------------------------------------------------------
 000      0         0         4         1          1         0
 001      1         4         8         1          1         0
 002      2         8        12         1          1         0
 003      3        12        16         1          1         0
 004      4        16        21         1          1         0

The points are arranged in clockwise order and are numbered as shown in the first figure.

 pnt shape  part  X       Y     
--------------------------------
 000     0          0.00    4.00
 001     0          0.00   10.00
 002     0          4.00   10.00
 003     0  ___     0.00    4.00
 004     1  -o      4.00   10.00
 005     1         10.00   10.00
 006     1         10.00    6.00
 007     1  ___     4.00   10.00
 008     2  -o     10.00    6.00
 009     2         10.00    0.00
 010     2          6.00    0.00
 011     2  ___    10.00    6.00
 012     3  -o      6.00    0.00
 013     3          0.00    0.00
 014     3          0.00    4.00
 015     3  ___     6.00    0.00
 016     4  -o      0.00    4.00
 017     4          4.00   10.00
 018     4         10.00    6.00
 019     4          6.00    0.00
 020     4          0.00    4.00

The points, and hence their segments, can be classed as either on, out, or in.  Duplicate ids are kept in this example so you can follow the construction of the sub-polygons.

on    :    [0, 1, 2, 3]
out   :   [5, 6, 7, 8]
in_    :   []

So an `out` point would be used to identify polygon parts that are exterior to the common overlay/clipped area.

segs_out :
  [[0, 5, 1, 1, 0],
   [1, 6, 2, 2, 1],
   [2, 7, 3, 3, 2],
   [3, 8, 0, 0, 3]]

segs_in :
  [0, 1, 1, 2, 2, 3, 3, 0]

In array format

# sequence ID values
[array([0, 1]),
 array([0, 5, 1]),
 array([1, 2]),
 array([1, 6, 2]),
 array([2, 3]),
 array([2, 7, 3]),
 array([3, 0]),
 array([3, 8, 0])]

# lexicographic sort.  Increasing x and decreasing y

[array([0, 1]),
 array([0, 5, 1]),
 array([0, 3]),
 array([0, 8, 3]),
 array([1, 2]),
 array([1, 6, 2]),
 array([3, 2]),
 array([3, 7, 2])]

In visual form

aoi_aoi0_sweep_example_0.pngThis forms the basis of "sweep-line" or "plane-sweep" algorithms.

In the figure to the right, there are 4 sweep lines which are identified where the two polygons share a common point.  The common points are 0, 1, 2, 3).

Note that the original order differs from the lexicographic order since point id 3 is encounted before number 2.

 

 

 

 

Lexicographic sorting uses "keys" to specify which values to use in considering the sort.  Consider these examples of how the key value (in this case, x or y column) and their value (positive or negative)

affects point sort order.

kys0 = np.lexsort((aoi[:, -1], aoi[:, 0]))   # increasing x, decreasing y
kys1 = np.lexsort((-aoi[:, -1], aoi[:, 0]))  # increasing x, increasing y
kys2 = np.lexsort((aoi[:, 0], aoi[:, -1]))   # increasing y, increasing x
kys3 = np.lexsort((-aoi[:, 0], aoi[:, -1]))  # increasing y, increasing y
kys0, kys1, kys2, kys3
 
(array([0, 4, 1, 3, 2]),  # -- note 0 and 4 are duplicate
 array([1, 0, 4, 2, 3]),
 array([0, 4, 3, 1, 2]),
 array([3, 0, 4, 2, 1]))

aoi[kys0]
array([[  0.00,   0.00],
       [  0.00,   0.00],
       [  0.00,  10.00],
       [ 10.00,   0.00],
       [ 10.00,  10.00]])

aoi[kys1]
array([[  0.00,  10.00],
       [  0.00,   0.00],
       [  0.00,   0.00],
       [ 10.00,  10.00],
       [ 10.00,   0.00]])

aoi[kys2]
array([[  0.00,   0.00],
       [  0.00,   0.00],
       [ 10.00,   0.00],
       [  0.00,  10.00],
       [ 10.00,  10.00]])

aoi[kys3]
array([[ 10.00,   0.00],
       [  0.00,   0.00],
       [  0.00,   0.00],
       [ 10.00,  10.00],
       [  0.00,  10.00]])

Normally in a lexicograph sort, the last key is the primary sort key and preceding keys are used as the secondary/tertiary/etc sort keys.

So don't just use "plain" sorting... consider "plane" sort or its fancy term "lexicographic" sorting.

Things can get rather complex when multiple points share common x and/or y values... have a look.

E_d1_base_image_1.pngE_d1_base_image_2.png

 

 

 

 

 

 

 

 

Or more simply.  Four lines and sorted with various orders.

a  # -- common start point (0, 0), increasing end points
array([[ 0,  0, 10,  0],
       [ 0,  0, 10,  3],
       [ 0,  0, 10,  6],
       [ 0,  0, 10, 10]])

# sort by decreasing y value of the end point
wh1_ = np.lexsort((a[:, 2], -a[:, 3], a[:, 1], a[:, 0]))  # the keys
a[wh1_]
array([[ 0,  0, 10, 10],
       [ 0,  0, 10,  6],
       [ 0,  0, 10,  3],
       [ 0,  0, 10,  0]])
Tags (3)
Contributors
About the Author
Retired Geomatics Instructor (also DanPatterson_Retired). Currently working on geometry projects (various) as they relate to GIS and spatial analysis. I use NumPy, python and kin and interface with ArcGIS Pro.
Labels