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