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 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 \$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> "; <>;

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?