Beefy Boxes and Bandwidth Generously Provided by pair Networks
The stupid question is the question not asked
 
PerlMonks  

Marching Squares (for contouring) with a PDL convolution

by etj (Curate)
on May 05, 2024 at 16:12 UTC ( [id://11159294]=CUFP: print w/replies, xml ) Need Help??

This implements a partial Marching Squares algorithm (see https://en.wikipedia.org/wiki/Marching_squares). Limitations:
  • It doesn't do the linear interpolation bit, because I couldn't figure a lazy way of getting it to do that. Probably doubling the coordinate offsets and using those as index offsets would work.
  • Making a bunch of individual line-segments and drawing each one is very slow in PGS. Joining them into polylines is possible with the not-yet-released next version of PDL (there's a path_join which allows this), which goes much quicker.
If you change the if (0) to 1, it shows you its lookup table instead of drawing contours.
use strict; use warnings; use PDL; use PDL::Image2D; use PDL::Graphics::Simple; my $LOOKUP = pdl( # relative to cell, x1,y1,x2,y2 for each line; 0 is invalid: lines s +tart edge [[ 0, 0, 0, 0],[ 0, 0, 0, 0]], # 0 [[-0.5, 0, 0,-0.5],[ 0, 0, 0, 0]], [[ 0,-0.5, 0.5, 0],[ 0, 0, 0, 0]], # 2 [[-0.5, 0, 0.5, 0],[ 0, 0, 0, 0]], [[ 0, 0.5, 0.5, 0],[ 0, 0, 0, 0]], # 4 [[ 0,-0.5, 0.5, 0],[-0.5, 0, 0, 0.5]], [[ 0,-0.5, 0, 0.5],[ 0, 0, 0, 0]], # 6 [[-0.5, 0, 0, 0.5],[ 0, 0, 0, 0]], [[-0.5, 0, 0, 0.5],[ 0, 0, 0, 0]], # 8 [[ 0,-0.5, 0, 0.5],[ 0, 0, 0, 0]], [[-0.5, 0, 0,-0.5],[ 0, 0.5, 0.5, 0]], # 10 [[ 0, 0.5, 0.5, 0],[ 0, 0, 0, 0]], [[-0.5, 0, 0.5, 0],[ 0, 0, 0, 0]], # 12 [[ 0,-0.5, 0.5, 0],[ 0, 0, 0, 0]], [[-0.5, 0, 0,-0.5],[ 0, 0, 0, 0]], # 14 [[ 0, 0, 0, 0],[ 0, 0, 0, 0]], ); sub marching_squares { my ($c, $data, $points) = @_; my $kernel = pdl q[4 8; 2 1]; my $contcells = conv2d($data < $c, $kernel)->slice(':-2,:-2'); my $segs = $LOOKUP->slice([],[],$contcells->flat)->clump(1..2); my $segsinds = $segs->orover; my $segsmask = $segsinds->dummy(0,4); my $contcoords = +($contcells->ndcoords->inflateN(1,2)->dupN(2) + 0. +5)->clump(1,2); my $segscoords = ($segs + $contcoords)->whereND($segsmask); $segscoords->splitdim(0,4)->clump(1,2); } if (0) { my $win = pgswin(multi=>[4,4]); my $xrange = [-0.5,0.5]; my $yrange = [-0.5,0.5]; my $i = 0; for my $lines ($LOOKUP->dog) { $win->plot( (map +(with=>'lines', $_->splitdim(0,2)->mv(0,-1)->dog), $lines->d +og), {xrange=>$xrange,yrange=>$yrange,j=>1,title=>$i++}, ); } print "ret> "; <>; exit; } my $SIZE = 50; my $vals = rvals($SIZE,$SIZE)->divide($SIZE/12.5)->sin; my $cntr_cnt = 9; my @cntr_threshes = zeroes($cntr_cnt+2)->xlinvals($vals->minmax)->list +; @cntr_threshes = @cntr_threshes[1..$#cntr_threshes-1]; my $win = pgswin(); my $xrange = [0,$vals->dim(0)-1]; my $yrange = [0,$vals->dim(1)-1]; $win->plot(with=>'image', $vals, {xrange=>$xrange,yrange=>$yrange,j=>1 +},); for my $thresh (@cntr_threshes) { my $segscoords = marching_squares($thresh, $vals); $win->oplot( (map +(with=>'lines', $_->splitdim(0,2)->mv(0,-1)->dog), $segscoor +ds->splitdim(0,4)->clump(1,2)->dog), {xrange=>$xrange,yrange=>$yrange,j=>1}, ); } print "ret> "; <>;

Log In?
Username:
Password:

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

How do I use this?Last hourOther CB clients
Other Users?
Others admiring the Monastery: (3)
As of 2024-06-16 21:14 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found

    Notices?
    erzuuli‥ 🛈The London Perl and Raku Workshop takes place on 26th Oct 2024. If your company depends on Perl, please consider sponsoring and/or attending.