Beefy Boxes and Bandwidth Generously Provided by pair Networks
Keep It Simple, Stupid
 
PerlMonks  

Re: combine(merge) polygons

by vr (Curate)
on Jul 25, 2018 at 13:16 UTC ( [id://1219256]=note: print w/replies, xml ) Need Help??


in reply to combine(merge) polygons

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]]]

Replies are listed 'Best First'.
Re^2: combine(merge) polygons (Let's do it with FORTRAN)
by vr (Curate) on Jul 31, 2018 at 20:53 UTC

    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.

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: note [id://1219256]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others scrutinizing the Monastery: (5)
As of 2024-03-28 20:04 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found