Menu

Merging polygons

2022-02-22
2022-03-21
  • Paul Higgins

    Paul Higgins - 2022-02-22

    Using Jan_Eriks Polyshape.cls I am able to am able to merge polygons (see the attachment) by following the intersections of the two polygons but it is quite complicated logic .

    But I am wondering if there is some algorithm that can be used to merge polygons. I did find a convex hull algorithm that I converted to rex but a convex hull is not what is needed.

    So my question is does anyone know of some logic to merge the points of two polygons?

     
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-02

    Found it in the link above and it's available again and may not stay so forever. Do write/port the code to ooRexx!
    https://github.com/rickbrew/GeneralPolygonClipper

    Conceptual:
    PolyShape use an Array of Points while Polygon (and subclasses) use an Array of Arrays of Points, the latter to be able to handle holes in the shape.
    gpc have code for union, intersection etc.

     

    Last edit: Jan-Erik Lärka 2022-03-02
  • Paul Higgins

    Paul Higgins - 2022-03-02

    I already look at it and found it's too hard too follow. I have downloaded many others and some converted but they do not work. For now I am sticking with the original one I had that follows the intersection points. Thanks for looking.

     
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-02

    Polygon A union with polygon B
    1. Find line section that intersect
    2. Find intersection point and insert it to A at position returned above
    3. Find side, if next point on B is outside of shape A
    4a.If true follow points in shape B forward until next intersection
    4b.if false follow points in shape B backward until next intersection
    5. Remove points in shape A from first intersection to the second intersection
    6. Add points from B until intersection and found intersection point to A
    7. Repeat until done

    Current polyshape class return only one intersection per line segment to see if shape B intersect A on the same line segment.
    use previous segments returned e.g. 2,3
    polyB~intersect( polyA, 2 )
    polyA~intersect( polyB, 3 )
    to find out.

     

    Last edit: Jan-Erik Lärka 2022-03-02
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-02
    p1 = .WGS84~new( '0.6819797,40.8142401,0 0.7656636,40.7931155,0 0.8275485,40.7320805,0 0.9003329,40.7227141,0 1.0540547,40.9405862,0 0.8796467,41.0090171,0 0.6819797,40.8142401,0' )
    
    p2 = .WGS84~new( '0.8304525,41.047812,0 0.9046103,40.9358631,0 1.1160971,40.9182243,0 1.1174704,40.994971,0 1.0172201,41.0571325,0 0.8304525,41.047812,0' )
    
    say p1~toString
    
    p1~merge( p2 )
    
    say  p1~toString
    
    ::requires 'WGS84.cls'
    

    Output:

    0.68198,40.81424,0 0.765664,40.793116,0 0.827549,40.732081,0 0.900333,40.722714,0 1.054055,40.940586,0 0.879647,41.009017,0 0.68198,40.81424,0

    0.68198,40.81424,0 0.765664,40.793116,0 0.827549,40.732081,0 0.900333,40.722714,0 1.042602,40.924354,0 1.116097,40.918224,0 1.11747,40.994971,0 1.01722,41.057133,0 0.830453,41.047812,0 0.865431,40.995009,0 0.68198,40.81424,0

    Limitations:
    1. Only capable of merge of the outer shape, no holes that is.
    2. Only polygons drawn in the same direction, clock wise or counter clock wise.
    3. and possibly not two intersections for the same line segment

    Polygon.cls

    ::method merge
    NUMERIC DIGITS self~precision
    Use Arg shape
    Do i = 1 To self~shapes~Items
       If shape~isInstanceOf( .Polygon ) Then
       Do j = 1 To shape~shapes~Items
          self~shapes[ i ]~merge( shape~shapes[ j ] )
       End
       Else If shape~isInstanceOf( .PolyShape ) Then
          self~shapes[ i ]~merge( shape )
    End
    Return self
    

    PolyShape.cls

    ::METHOD merge
    Use Arg shape
    PARSE VALUE self~intersect( shape ) WITH this_section '.' add_section
    If \DATATYPE( this_section, 'W' ) Then
       Return self
    Do i = this_section + 1 To self~point~Items
       PARSE VALUE self~intersect( shape, i ) WITH next_section '.' fin_section
       If DATATYPE( fin_section, 'W' ) Then
          LEAVE i
    End
    If \DATATYPE( next_section, 'W' ) Then Return self
    ipnt = self~intersectionPoint( shape, this_section )
    opnt = self~intersectionPoint( shape, next_section )
    If ipnt = .nil | opnt = .nil Then
       Return self
    Do next_section - this_section
       self~point~delete( this_section + 1 )
    End
    self~point~insert( ipnt, this_section )
    Do i = 1 To shape~point~Items + 1 - fin_section - add_section
       self~point~insert( shape~point[ i + add_section ], this_section + i )
    End
    self~point~insert( opnt, this_section + i )
    Return self
    
     

    Last edit: Jan-Erik Lärka 2022-03-02
  • Paul Higgins

    Paul Higgins - 2022-03-03

    Interesting. Will have a look at this and see if it can replace my existing code. Thanks

     
  • Paul Higgins

    Paul Higgins - 2022-03-03

    JanErik, I tried two polygons and get a weird result:

    p1 = .WGS84~new('117.494063,8.487999,0 117.529286,8.567740,0 117.570462,8.549952,0 117.535230,8.470211,0  117.494063,8.487999,0')
    
    p2 = .WGS84~new('117.528199,8.564692,0 117.538691,8.602737,0 117.582045,8.591045,0 117.571549,8.553000,0  117.528199,8.564692,0')
    
    p1~merge( p2 )
    
    say  p1~toString
    

    result:
    117.494063,8.487999,0 117.52863,8.566254,0 117.538691,8.602737,0 117.549874,8.558846,0 117.570462,8.549952,0 117.53523,8.470211,0 117.494063,8.487999,0

    the first attachment is the two polygons and the second is the resulting polygon

     

    Last edit: Paul Higgins 2022-03-03
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-04

    Change the calculation for this loop as it stop to soon or to late

    Do i = 1 To shape~point~Items + 1 - fin_section - add_section
       self~point~insert( shape~point[ i + add_section ], this_section + i )
    End
    

    while you're at it add another loop to take care of clock wise vs. counter clock wise direction.

    If add_section < fin_section Then
    Do i = 1 To ...
    ...
    End
    Else
    Do i = 1 To ... By -1
    ...
    End
    

    or better yet, merge into one.

     

    Last edit: Jan-Erik Lärka 2022-03-04
  • Paul Higgins

    Paul Higgins - 2022-03-04

    Are you asking me to figure out the new calculation?

     
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-04

    yes,
    insert a

    call trace '?i'
    

    right before the line you want to begin watching... then press enter to step through.
    type

    trace off

    to end before it end as the method end.

     

    Last edit: Jan-Erik Lärka 2022-03-04
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-05

    How about:

    Do i = 1 To ABS( fin_section - add_section )
       If add_section < fin_section Then
          self~point~insert( shape~point[ i + add_section ], this_section + i )
       Else
          self~point~insert( shape~point[ add_section - i + 1 ], this_section + i )
    End
    
     
  • Paul Higgins

    Paul Higgins - 2022-03-06

    Will give a try when I get a chance. Right now have other problems.

     
  • Paul Higgins

    Paul Higgins - 2022-03-06

    Its better (see attachment) but that little jog to point 6 should not exist. Point 5 should connect to point 7.

     
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-06

    no, i can't accept that, the upper polygon doesn't cover that area, thus shouldn't stretch into it.
    union/merge doesn't mean stretch, it mean use the given structure that overlap.

    is there a little pointy thing on the left side? that one should just be a smooth transition though.

     
  • Paul Higgins

    Paul Higgins - 2022-03-06

    You are right.

    The problem I get when I replace my logic with yours is on the fourth merge where it tries to merge the existing large polygon with the next one (see attachment "Two polygons to merge"). The resulting polygon (second attachment "Result Merged") does not follow the next one, probably because the polygons share a point.

    I put the two polygons into the test program:

    p1 = .WGS84~new('117.494063,8.487999,0  117.52863,8.566254,0    117.537862,8.599729,0   117.524659,8.624429,0   117.552887,8.658872,0   117.578478,8.683697,0   117.609608,8.678317,0   117.656904,8.668307,0   117.706455,8.621621,0   117.679609,8.586361,0   117.612843,8.544386,0   117.574668,8.564307,0   117.571549,8.553,0  117.549874,8.558846,0   117.570462,8.549952,0   117.53523,8.470211,0    117.494063,8.487999,0 ')  
    p2 = .WGS84~new('117.635281,8.673601,0 117.685860,8.687754,0   117.698101,8.645000,0   117.647521,8.630847  117.635281,8.673601,0')
    
    p1~merge( p2 )
    
    str =  p1~toString
    say str
    -- make a lat lon list of the result 117.494063,8.487999,0   117.52863,8.566254,0
    do n = 1 to words(str)
       parse value word(str,n) with lon ',' lat ',' .
       say lat lon
    end
    

    This is the result:
    117.494063,8.487999,0 117.52863,8.566254,0 117.537862,8.599729,0 117.524659,8.624429,0 117.552887,8.658872,0 117.578478,8.683697,0 117.609608,8.678317,0 117.6355,8.672837,0 117.647521,8.630847,0 117.68541,8.641449,0 117.706455,8.621621,0 117.679609,8.586361,0 117.612843,8.544386,0 117.574668,8.564307,0 117.571549,8.553,0 117.549874,8.558846,0 117.570462,8.549952,0 117.53523,8.470211,0 117.494063,8.487999,0
    8.487999 117.494063
    8.566254 117.52863
    8.599729 117.537862
    8.624429 117.524659
    8.658872 117.552887
    8.683697 117.578478
    8.678317 117.609608
    8.672837 117.6355
    8.630847 117.647521
    8.641449 117.68541
    8.621621 117.706455
    8.586361 117.679609
    8.544386 117.612843
    8.564307 117.574668
    8.553 117.571549
    8.558846 117.549874
    8.549952 117.570462
    8.470211 117.53523
    8.487999 117.494063

    I had the wrong data for p2 . but the result is the same

     

    Last edit: Paul Higgins 2022-03-09
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-08

    Hmm, that is a subtract, the opposite of union.
    Try to trace it and see when things go wrong.

     
  • Paul Higgins

    Paul Higgins - 2022-03-09

    I am slowly digesting your code. I have never been a fan of the trace. I prefer to simply put in say statements.
    Now that I corrected the data it makes more sense. Will let you know what I find.

     
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-09

    ok, though you can test variables and adjust them to try things out.
    What you find out is important to improve the code. :-)

     
  • Paul Higgins

    Paul Higgins - 2022-03-09

    I made the following change and it seems to work:
    if add_section > fin_section then add_section = 0
    I put that right before your new code see the attached new polyshape.cls

    I will do more testing . I think the code will fail if there are more than two intersection points.

     
  • Paul Higgins

    Paul Higgins - 2022-03-11

    I found the code would not support two cases: * When the two intersection points were on the same section (see attachment "IP on Same Section") * When there were more than two intersection points (see attachment "More than two Intersections")
    I've attached the new Polyshape.cls it works great.

     
  • Jan-Erik Lärka

    Jan-Erik Lärka - 2022-03-13

    Yes, just as I wrote about two intersections on the same section.
    Really neat that you solved it and the two or more intersection points problem.

    "merge" grew, does your improvement of the method contain one (or more) section(s) that would be useful on its own or perhaps with "subtract", "overlap" or something else?
    If so, extract it/them and write "subtract" and "overlap" as well while you're at it and we'll have a more complete and general polyshape class.

     
  • Paul Higgins

    Paul Higgins - 2022-03-14

    JanErik,

    When I put the code into my program all mt testing of it worked. But when I asked my tester too test he found it failed. I am working on a fix,

    Paul

     
  • Paul Higgins

    Paul Higgins - 2022-03-21

    Well it took awhile but finally I have the merge code that has passed all my testers.
    It follows the intersection points of the polygons to merge them .
    I had to add a new method and some routines:
    - method distance - it finds the distance between two sets of lat/lons. Used to sort the intersections.
    - routine ring_array - orders the sections of the merging polygon so they can be processed sequentially.
    - routines Create_WP and Create_Route which create gpx files of the intersection points and the polygons and the result, used for diagnostics.

    The new Polyshape.cls is attached.

     

    Last edit: Paul Higgins 2022-03-21

Log in to post a comment.