http://qs1969.pair.com?node_id=1219239

dideod.yang has asked for the wisdom of the Perl Monks concerning the following question:

Hi Monks:). I know monks are professional perler:) My question is where is easy and nice module to combine polygon.. I searched all of perl module site. and I got Math::polygon. This module gave me many functions. There are calculating area, perimeter, searching point in polygon etc.. But I want to achieve to combine polygons. Unfortunately that module has not that function.. Below examples, I show only rectangular but it is not always rectangular, sometimes polygons.. Also that polygons have not hole and always have minimum 5 points. Is there any nice idea? or modules you know to combine polygons? I couldn't sleep because of that problems.. Please help me
#case 1 => There is common area. merge case my @poly1 = ([0,0],[0,2],[2,2],[2,0],[0,0]); my @poly2 = ([0,0],[3,0],[3,1],[0,1],[0,0]); my @merge_poly12; print ("@merge_poly12\n"); #then I want to print [0,0],[0,2],[2,2],[2, +1],[3,1],[3,0],[0,0]; #case 2 => There is not common area, No merge my @poly1 = ([0,0],[0,2],[2,2],[2,0],[0,0]); my @poly3 = ([3,0],[3,2],[4,2],[4,0],[3,0]); my @merge_poly13; print ("@merge_poly13\n"); # then I want to print nothing or only @pol +y1. #case 3 => There is not common area, But have line, Merge case my @poly1 = ([0,0],[0,2],[2,2],[2,0],[0,0]); my @poly4 = ([0,2],[0,3],[1,3],[1,2],[0,2]); my @merge_poly13; print ("@merge_poly14\n"); # then I want to print [0,0],[0,3],[1,3],[ +1,2],[2,2],[2,0],[0,0]

Replies are listed 'Best First'.
Re: combine(merge) polygons
by haj (Vicar) on Jul 25, 2018 at 11:47 UTC

    This problem looks so easy if you tackle it with a pencil and paper, and yet it is rather complex to do by a computer program. There's a lot of tedious geometry to do. On the other hand, geometry isn't Perl specific. There's an Article on stackoverflow which explains the steps.

    While writing this, I saw that choroba recommends Math::Geometry::Planar which helps a lot.

    There are some more issues to consider:

    • Be aware of rounding errors, for example if vertices are "almost parallel".
    • Be aware that the union of the two shapes can contain a hole unless both polygons are convex. A merge of two polygons isn't necessarily convex, as your examples show. Therefore, the library returns a structure which contain more than just a list of points.
    • Be aware of self-intersecting polygons, like these:
      my @poly1 = ([0,0], [1,0], [0,1], [1,1], [0,0]); my @poly2 = ([0,0], [0,1], [1,0], [1,1], [0,0]); # A possible result could be: my @merge = ([0,0], [1,0], [1,1], [0,1], [0,0]);

    Fun fact: Math::Polygon::Calc will return 0 as the area of @poly1 and @poly2, while the merged polygon has an area of 1.

      ... the union of the two shapes can contain a hole unless both polygons are concave. [emphasis added]

      Shouldn't "concave" be "convex" in the quoted passage?


      Give a man a fish:  <%-{-{-{-<

        Ooops... Thanks, fixed!
      Thanks a lot:) I search Math::Geometry::Planar in CPAN. But I don't know which function make merge polygons. Did I miss function??? I am not used to work perl modules yet.. Can you tell me which function is merge polygon. Thanks you:)
Re: combine(merge) polygons
by choroba (Cardinal) on Jul 25, 2018 at 11:27 UTC
Re: combine(merge) polygons
by hda (Chaplain) on Jul 25, 2018 at 13:02 UTC

    You have already received good Perl answers. I nevertheless want to flag about the potential complexity of this issue. Think about what do you mean exactly by "combine polygons", since there are several different ways to "combine" polygons (see, for example, this link). As you already noted, these operations can soon get very complicated for non regular shapes.

    Reinventing the wheel is surely educative and I believe it is always worth the effort. However, if you are in a hurry and need to get job done, I would recommend exploring any of the solutions that are extensively available in Geographic Information Systems.

    There are a lot of good free and open source GIS packages out there, some more advanced than others. See, for example: GRASS, QGIS, and, if you like operating on databases: Spatialite or, on a PostgreSQL database, PostGIS.

    All these apply the underlying engine for spatial operations implemented in GEOS, which in turn implements specifications by Open Geospatial Consortium

    Good luck!

Re: combine(merge) polygons
by vr (Curate) on Jul 25, 2018 at 13:16 UTC

    Looks like everything (and everyone) points into direction of GPC and its Perl bindings. In case its license is OK (if not, whoever brave enough follow links on that page and bind to e.g. Fortran lib), but here's for a start. Your requirements are a bit odd, in that you want union if there is intersection or touch. It's solved with check for only single sub-polyline in union (assuming you neither use nor create "holes" or self-intersecting shapes). The "cleanup" is to "remove colinear points", but note that start point is not duplicated as end point.

    use strict; use warnings; use feature 'say'; use Data::Dump 'dd'; my @sets = ( [[0,0],[0,2],[2,2],[2,0],[0,0]], [[0,0],[3,0],[3,1],[0,1],[0,0]], [[3,0],[3,2],[4,2],[4,0],[3,0]], [[0,2],[0,3],[1,3],[1,2],[0,2]], ); ######################################## #use Math::Geometry::Planar::GPC::PolygonXS; # ... skipped, don't need it ######################################## use Math::Geometry::Planar; my $poly = Math::Geometry::Planar-> new; $poly-> points( $sets[ 0 ]); for ( 1 .. 3 ) { say "********* case #$_:"; my $clip = Math::Geometry::Planar-> new; $clip-> points( $sets[ $_ ]); my $gpc = GpcClip( 'UNION', $poly-> convert2gpc, $clip-> convert2gpc, ); my @pgons = Gpc2Polygons( $gpc ); dd $pgons[ 0 ]-> cleanup unless $#pgons } __END__ ********* case #1: [[[2, 2], [0, 2], [0, 0], [3, 0], [3, 1], [2, 1]]] ********* case #2: ********* case #3: [[[1, 3], [0, 3], [0, 0], [2, 0], [2, 2], [1, 2]]]

      TL;DR: Bound to functions in compiled FORTRAN source. FORTRAN still rocks (speed!), vs C :) Though, seriously, I didn't investigate the reason for such difference. I also hope I didn't make stupid mistakes, it's not that I accomplish such things often, and some are here for the 1st time. Sorry in advance, if anything. Please report any bugs :)

      bind to e.g. Fortran lib

      I had some spare time yesterday, so why not... Initially thought to write my first one in "Meditations" or rather "CUFP" here, but then who am I to presume my sentimental reminiscences of punching FORTRAN code onto cards some 30+ years ago are "cool" :). Therefore, writing this as follow-up comment. And yes, wanted to do something FORTRAN-related for a couple years anyway, so I'm quite pleased that (looks like) it, actually, works, however clumsily and incompletely implemented.

      OK, to business. This experiment was done on Windows (Strawberry Perl + MinGW) and later today confirmed in Linux, with only -fPIC option additionally required (see below).

      Looks like the source for the POLYPACK library is only available (publicly, now) as part of somewhat large distribution, but just a few (4, to be precise) small files are required from it. The discussion how to compile this FORTRAN source can be found here, a few remarks contain valuable information, one link leads further to here.

      So, not only this library was written 20+ years ago, but the stars were aligned so that the author had chosen to use non-standard, 35+ 45+ years old extension (IFTRAN) i.e. source file filter, to translate to "proper" FORTRAN. Looks like compiling this stand-alone filter is most difficult part.

      1. As per information in links above, this header (re-named to c.h) has to be edited -- "SED_NGCALLF" replaced with "reg##_" -- and made available for next step.

      2. C library, required for Iftran pre-processor, is compiled next:

      gcc -I. -c -o rwchinfl.o rwchinfl.c
      

      3. Same for Iftran itself:

      gfortran -c -o Iftran.o Iftran.f
      gfortran -o Iftran.exe Iftran.o rwchinfl.o
      

      4. Then, at last, pre-process the library file:

      Iftran.exe < CodeIftran > polypack.f
      *** IFTRAN SUMMARY - 10067 CARDS,   0 ERRORS ***
      

      (Nice output... 10000+ cards...)

      5. Then just compile it ("-fPIC" was required in my Linux):

      gfortran -c -o polypack.o polypack.f
      

      6. The next step is per instructions for Inline::C::Cookbook:

      ar cru libpolypack.a polypack.o
      

      And we are done with compiling POLYPACK library. Perl package that uses this library is as follows:

      package MyUnion; use strict; use warnings; use FindBin; our $VERSION = 0.001; our @EXPORT_OK = qw/ wrap_ppunpo /; use Inline C => Config => LIBS => "-L$FindBin::Bin/lib -lgfortran -lpolypack"; use Inline C => <<'EOC'; extern void ppunpo_( float* xclip, float* yclip, int* nclip, float* xpoly, float* ypoly, int* npoly, float* rwork, void* iwork, int* nwork, void* cback, int* error ); AV* ary; // AoA; collect output here, array (i.e. one call // to call-back function, below) per sub-path static void urpp( float* xcra, float* ycra, int* ncra ) { int i; AV* pgone = newAV(); for ( i = 0; i < *ncra; i ++ ) { AV* point = newAV(); av_push( point, newSVnv( xcra[ i ])); av_push( point, newSVnv( ycra[ i ])); av_push( pgone, newRV_noinc( point )); } av_push( ary, newRV_noinc( pgone )); } SV* wrap_ppunpo( SV* ccp, SV* csp ) { AV* ccp_ary = SvRV( ccp ); AV* csp_ary = SvRV( csp ); int nccp = av_top_index( ccp_ary ) + 1; int ncsp = av_top_index( csp_ary ) + 1; int nwrk = 20 * ( nccp + ncsp ); float* xccp = malloc( nccp * sizeof( float )); float* yccp = malloc( nccp * sizeof( float )); float* xcsp = malloc( ncsp * sizeof( float )); float* ycsp = malloc( ncsp * sizeof( float )); float* work = malloc( nwrk * sizeof( float )); int i; AV* point; for ( i = 0; i < nccp; i ++ ) { point = SvRV( *av_fetch( ccp_ary, i, 0 )); xccp[ i ] = SvNV( *av_fetch( point, 0, 0 )); yccp[ i ] = SvNV( *av_fetch( point, 1, 0 )); } for ( i = 0; i < ncsp; i ++ ) { point = SvRV( *av_fetch( csp_ary, i, 0 )); xcsp[ i ] = SvNV( *av_fetch( point, 0, 0 )); ycsp[ i ] = SvNV( *av_fetch( point, 1, 0 )); } ary = newAV(); int ierr; ppunpo_( &xccp[0], &yccp[0], &nccp, &xcsp[0], &ycsp[0], &ncsp, &work[0], &work[0], &nwrk, &urpp, &ierr ); SV* ret; if ( ierr ) { av_undef( ary ); ret = newSV( 0 ); } else { ret = newRV_noinc( ary ); } free( xccp ); free( yccp ); free( xcsp ); free( ycsp ); free( work ); // Note, this also frees memory // used by data passed to call-back return( ret ); } EOC 1;

      It's purely experimental module with just single wrapper function, with conveniences and many other things missing, e.g. error "reporting" (rather, "taking note of") quite rudimental. Unfortunately, I have already "purged" it from quite a few typecasts, which, I now see, are required for "clean" i.e. warning-free compilation, but if everything goes OK, no-one will see these warnings. So let it be "as is"...

      I also see, in retrospect, I've changed formal parameters names in "extern" declaration to something perhaps more meaningful, but elsewhere they are as in description, and don't they add a nice retro touch :).

      This function expects 2 "normal" Perl arrayrefs (one per polygon) of "usual" 2D X-Y points (refs to 2-element Perl arrays) and translates them to somewhat exotic representation that POLYPACK requires. Beware, also, that, if polygons are "compound" (more than single sub-path per path), this library wants them in unconventional form (X-Y pairs have to be arranged accordingly, see documentation), but, if such compound path happens in output, they are returned in "standard" form of, well, multiple polygons.

      The test data for benchmarking was found in pages about another project. I have extracted 2 data files from 69839 folder.

      Directory structure for test script is expected as follows:

      180731_polypack
          |___lib
          |    |___libpolypack.a
          |    |___MyUnion.pm
          |
          |___sample
          |    |___c.wlr
          |    |___s.wlr
          |
          |___180731.pl
          
      

      The "wlr" files (simple kind, a polygon per file) are as

      1
      1
      31084
      26381,10485
      27120,10807
      27120,15466
      26837,15466
      26835,15465
      

      ...and so on.

      use strict; use warnings; use feature 'say'; use FindBin; use Data::Dump 'dd'; use autodie; use Benchmark 'cmpthese'; use lib "$FindBin::Bin/lib"; use MyUnion 'wrap_ppunpo'; use Math::Geometry::Planar; { ###### Demo data output block my @sets = ( [[0,0],[0,2],[2,2],[2,0],[0,0]], [[0,0],[3,0],[3,1],[0,1],[0,0]], [[3,0],[3,2],[4,2],[4,0],[3,0]], [[0,2],[0,3],[1,3],[1,2],[0,2]], ); say "\n**** Demo output for GPC:\n"; my $clip_gpc_obj = Math::Geometry::Planar-> new; my $poly_gpc_obj = Math::Geometry::Planar-> new; $clip_gpc_obj-> points( $sets[ 0 ]); for ( 1.. 3 ) { $poly_gpc_obj-> points( $sets[ $_ ]); dd [ map @{ $_-> polygons }, Gpc2Polygons( GpcClip( 'UNION', $poly_gpc_obj-> convert2gpc, $clip_gpc_obj-> convert2gpc, ) ) ] } say "\n\n*** Demo output for POLYPACK:\n"; dd wrap_ppunpo( @sets[ 0, $_ ]) for 1 .. 3; } ####################### Benchmark my $clip_points = _get_data( 'c' ); # clipmask polygon points my $poly_points = _get_data( 's' ); # sample polygon points my $clip_gpc_obj = Math::Geometry::Planar-> new; my $poly_gpc_obj = Math::Geometry::Planar-> new; $clip_gpc_obj-> points( $clip_points ); $poly_gpc_obj-> points( $poly_points ); my $clip_gpc_struct = $clip_gpc_obj-> convert2gpc; my $poly_gpc_struct = $poly_gpc_obj-> convert2gpc; say "\n\n*** Benchmark (union of polygons with " . @$clip_points . ' and ' . @$poly_points . " points):\n"; cmpthese( -5, { GPC_lib => sub { my $union_gpc_struct = GpcClip( 'UNION', $poly_gpc_struct, $clip_gpc_struct, ); my ( $contour ) = Gpc2Polygons( $union_gpc_struct ); return $contour-> polygons; }, POLYPACK => sub { return wrap_ppunpo( $clip_points, $poly_points ); }, }); sub _get_data { # only simple wlr files expected! my $name = shift; open my $fh, '<', "sample/$name.wlr"; <$fh>; <$fh>; <$fh>; my @a; push @a, [ map { 0 + $_ } split ',', <$fh> ] until eof $fh; return \@a } __END__

      Output:

      **** Demo output for GPC: [[[2, 2], [0, 2], [0, 1], [0, 0], [3, 0], [3, 1], [2, 1]]] [ [[4, 2], [3, 2], [3, 0], [4, 0]], [[2, 2], [0, 2], [0, 0], [2, 0]], ] [[[1, 3], [0, 3], [0, 2], [0, 0], [2, 0], [2, 2], [1, 2]]] *** Demo output for POLYPACK: [ [[2, 2], [2, 1], [3, 1], [3, 0], [2, 0], [0, 0], [0, 2], [2, 2]], ] [ [[2, 2], [2, 0], [0, 0], [0, 2], [2, 2]], [[4, 2], [4, 0], [3, 0], [3, 2], [4, 2]], ] [ [[1, 3], [1, 2], [2, 2], [2, 0], [0, 0], [0, 2], [0, 3], [1, 3]], ] *** Benchmark (union of polygons with 31084 and 38755 points): s/iter GPC_lib POLYPACK GPC_lib 1.26 -- -98% POLYPACK 2.74e-002 4491% --

      The "demo" block can be deleted, it serves as rudimentary test suite, and also to show a few peculiarities: (1) POLYPACK polygons are always "closed" (as per documentation); (2) direction, and presence (or absence) of co-linear points may differ. (3) I had also to do some acrobatics with Math::Geometry::Planar, in that, e.g. Gpc2Polygons should have been named "Gpc2Contours", as it returns list of "contours", per its parlance.

      In benchmarking section, I gave quite a handicap to Math::Geometry::Planar, as you see. "Points" are converted to internal representation only once (as opposed to my module). Perhaps it's fair, because, unlike call to POLYPACK routine, GPC object data can be retained for subsequent calls (e.g., same mask to clip a list of polygons).

      Even so, POLYPACK seems to be a lot faster. OTOH, GPC (as wrapped into Math::Geometry::Planar) offers plenty of other useful functions, while POLYPACK does only a few boolean operations. And, I don't know if "usual" for 25+ years ago use of low precision "float" can lead to nasty errors in practical application. And of course, single test case is not enough for any serious judgement.

      Plus, like I said, licence for GPC may or may not be important, TWIMC.