plucky (3) PDL::ImageND.3pm.gz
NAME
PDL::ImageND - useful image processing in N dimensions
DESCRIPTION
These routines act on PDLs as N-dimensional objects, not as broadcasted sets of 0-D or 1-D
objects. The file is sort of a catch-all for broadly functional routines, most of which
could legitimately be filed elsewhere (and probably will, one day).
ImageND is not a part of the PDL core (v2.4) and hence must be explicitly loaded.
SYNOPSIS
use PDL::ImageND;
$y = $x->convolveND($kernel,{bound=>'periodic'});
$y = $x->rebin(50,30,10);
FUNCTIONS
convolve
Signature: (a(m); b(n); indx adims(p); indx bdims(q); [o]c(m))
N-dimensional convolution (Deprecated; use convolveND)
$new = convolve $x, $kernel
Convolve an array with a kernel, both of which are N-dimensional. This routine does
direct convolution (by copying) but uses quasi-periodic boundary conditions: each dim
"wraps around" to the next higher row in the next dim.
This routine is kept for backwards compatibility with earlier scripts; for most purposes
you want convolveND instead: it runs faster and handles a variety of boundary conditions.
convolve does not process bad values. It will set the bad-value flag of all output
ndarrays if the flag is set for any of the input ndarrays.
ninterpol()
N-dimensional interpolation routine
Signature: ninterpol(point(),data(n),[o]value())
$value = ninterpol($point, $data);
"ninterpol" uses "interpol" to find a linearly interpolated value in N dimensions,
assuming the data is spread on a uniform grid. To use an arbitrary grid distribution,
need to find the grid-space point from the indexing scheme, then call "ninterpol" -- this
is far from trivial (and ill-defined in general).
See also interpND, which includes boundary conditions and allows you to switch the method
of interpolation, but which runs somewhat slower.
rebin
Signature: (a(m); [o]b(n); int ns => n)
N-dimensional rebinning algorithm
$new = rebin $x, $dim1, $dim2,..;.
$new = rebin $x, $template;
$new = rebin $x, $template, {Norm => 1};
Rebin an N-dimensional array to newly specified dimensions. Specifying `Norm' keeps the
sum constant, otherwise the intensities are kept constant. If more template dimensions
are given than for the input pdl, these dimensions are created; if less, the final
dimensions are maintained as they were.
So if $x is a 10 x 10 pdl, then "rebin($x,15)" is a 15 x 10 pdl, while
"rebin($x,15,16,17)" is a 15 x 16 x 17 pdl (where the values along the final dimension are
all identical).
Expansion is performed by sampling; reduction is performed by averaging. If you want
different behavior, use PDL::Transform::map instead. PDL::Transform::map runs slower but
is more flexible.
rebin does not process bad values. It will set the bad-value flag of all output ndarrays
if the flag is set for any of the input ndarrays.
circ_mean_p
Calculates the circular mean of an n-dim image and returns the projection. Optionally
takes the center to be used.
$cmean=circ_mean_p($im);
$cmean=circ_mean_p($im,{Center => [10,10]});
circ_mean
Smooths an image by applying circular mean. Optionally takes the center to be used.
circ_mean($im);
circ_mean($im,{Center => [10,10]});
kernctr
`centre' a kernel (auxiliary routine to fftconvolve)
$kernel = kernctr($image,$smallk);
fftconvolve($image,$kernel);
kernctr centres a small kernel to emulate the behaviour of the direct convolution
routines.
convolveND
Signature: (k0(); pdl *k; pdl *aa; pdl *a)
Speed-optimized convolution with selectable boundary conditions
$new = convolveND($x, $kernel, [ {options} ]);
Convolve an array with a kernel, both of which are N-dimensional.
If the kernel has fewer dimensions than the array, then the extra array dimensions are
broadcasted over. There are options that control the boundary conditions and method used.
The kernel's origin is taken to be at the kernel's center. If your kernel has a dimension
of even order then the origin's coordinates get rounded up to the next higher pixel (e.g.
(1,2) for a 3x4 kernel). This mimics the behavior of the earlier "convolve" and
fftconvolve routines, so convolveND is a drop-in replacement for them.
The kernel may be any size compared to the image, in any dimension.
The kernel and the array are not quite interchangeable (as in mathematical convolution):
the code is inplace-aware only for the array itself, and the only allowed boundary
condition on the kernel is truncation.
convolveND is inplace-aware: say "convolveND(inplace $x ,$k)" to modify a variable in-
place. You don't reduce the working memory that way -- only the final memory.
OPTIONS
Options are parsed by PDL::Options, so unique abbreviations are accepted.
boundary (default: 'truncate')
The boundary condition on the array, which affects any pixel closer to the edge than
the half-width of the kernel.
The boundary conditions are the same as those accepted by range, because this option is
passed directly into range. Useful options are 'truncate' (the default), 'extend', and
'periodic'. You can select different boundary conditions for different axes -- see
range for more detail.
The (default) truncate option marks all the near-boundary pixels as BAD if you have bad
values compiled into your PDL and the array's badflag is set.
method (default: 'auto')
The method to use for the convolution. Acceptable alternatives are 'direct', 'fft', or
'auto'. The direct method is an explicit copy-and-multiply operation; the fft method
takes the Fourier transform of the input and output kernels. The two methods give the
same answer to within double-precision numerical roundoff. The fft method is much
faster for large kernels; the direct method is faster for tiny kernels. The tradeoff
occurs when the array has about 400x more pixels than the kernel.
The default method is 'auto', which chooses direct or fft convolution based on the size
of the input arrays.
NOTES
At the moment there's no way to broadcast over kernels. That could/should be fixed.
The broadcasting over input is cheesy and should probably be fixed: currently the kernel
just gets dummy dimensions added to it to match the input dims. That does the right thing
tersely but probably runs slower than a dedicated broadcastloop.
The direct copying code uses PP primarily for the generic typing: it includes its own
broadcastloops.
convolveND does not process bad values. It will set the bad-value flag of all output
ndarrays if the flag is set for any of the input ndarrays.
contour_segments
Signature: (c(); data(m,n); points(d,m,n);
[o] segs(d,q=CALC(($SIZE(m)-1)*($SIZE(n)-1)*4)); indx [o] cnt();)
Finds a contour in given data. Takes 3 ndarrays as input:
$c is the contour value (broadcast with this)
$data is an [m,n] array of values at each point
$points is a list of [d,m,n] points. It should be a grid monotonically increasing with m
and n.
Returns $segs, and $cnt which is the highest 2nd-dim index in $segs that's defined. The
contours are a collection of disconnected line segments rather than a set of closed
polygons.
The data array represents samples of some field observed on the surface described by
points. This uses a variant of the Marching Squares algorithm, though without being data-
driven.
contour_segments does not process bad values. It will set the bad-value flag of all
output ndarrays if the flag is set for any of the input ndarrays.
contour_polylines
Signature: (c(); data(m,n); points(d,m,n);
indx [o] pathendindex(q=CALC(($SIZE(m)-1)*($SIZE(n)-1)*5)); [o] paths(d,q);
byte [t] seenmap(m,n))
Finds polylines describing contours in given data. Takes 3 ndarrays as input:
$c is the contour value (broadcast with this)
$data is an [m,n] array of values at each point
$points is a list of [d,m,n] points. It should be a grid monotonically increasing with m
and n.
Returns $pathendindex, and $paths. Any $pathendindex entries after the pointers to the
ends of polylines are negative.
Algorithm
Has two modes: scanning, and line-walking. Scanning is done from the top left, along each
row. Each point can be considered as, at "a":
a|b
+-+-
c|d|e
Every potential boundary above, or to the left of (including the bottom boundaries), "a"
has been cleared (marked with a space above).
Boundary detection
This is done by first checking both points' coordinates are within bounds, then checking
if the boundary is marked seen, then detecting whether the two cells' values cross the
contour threshold.
Scanning
If detect boundary between "a"-"b", and also "a"-"d", "d"-"e", or "b"-"e", line-walking
starts "a"-"b" facing south.
If not, mark "a"-"b" seen.
If detect boundary "a"-"d" and "c"-"d", line-walking starts "a"-"d" facing west.
If detect boundary "a"-"d" and also "d"-"e" or "b"-"e", line-walking starts "a"-"d" facing
east.
If not, mark "a"-"d" seen, and continue scanning.
Line-walking
The conditions above guarantee that any line started will have at least two points, since
two connected "points" (boundaries between two cells) have been detected. The coordinates
of the back end of the starting "point" (boundary with direction) are recorded.
At each, a line-point is emitted and that "point" is marked seen. The coordinates emitted
are linearly interpolated between the coordinates of the two cells similarly to the
Marching Squares algorithm.
The next "point" is sought, looking in order right, straight ahead, then left. Each one
not detected is marked seen. That order means the walked boundary will always turn as much
right (go clockwise) as available, thereby guaranteeing enclosing the area, which deals
with saddle points.
If a next "point" is found, move to that and repeat.
If not, then if the front of the ending "point" (boundary plus direction) is identical to
the back of the starting point, a final point is emitted to close the shape. Then the
polyline is closed by emitting the current point-counter into "polyendindex".
use PDL;
use PDL::ImageND;
use PDL::Graphics::Simple;
$SIZE = 500;
$vals = rvals($SIZE,$SIZE)->divide($SIZE/12.5)->sin;
@cntr_threshes = zeroes(9)->xlinvals($vals->minmax)->list;
$win = pgswin();
$xrange = [0,$vals->dim(0)-1]; $yrange = [0,$vals->dim(1)-1];
$win->plot(with=>'image', $vals, {xrange=>$xrange,yrange=>$yrange,j=>1},);
for $thresh (@cntr_threshes) {
($pi, $p) = contour_polylines($thresh, $vals, $vals->ndcoords);
$pi_max = $pi->max;
next if $pi_max < 0;
$pi = $pi->where($pi > -1);
$p = $p->slice(',0:'.$pi_max);
@paths = path_segs($pi, $p->mv(0,-1));
$win->oplot(
(map +(with=>'lines', $_->dog), @paths),
{xrange=>$xrange,yrange=>$yrange,j=>1},
);
}
print "ret> "; <>;
contour_polylines does not process bad values. It will set the bad-value flag of all
output ndarrays if the flag is set for any of the input ndarrays.
path_join
Signature: (e(v=2,n);
indx [o] pathendindex(n); indx [o] paths(nout=CALC($SIZE(n)*2));
indx [t] highestoutedge(d); indx [t] outedges(d,d); byte [t] hasinward(d);
indx [t] sourceids(d);
; PDL_Indx d => d; int directed)
Links a (by default directed) graph's edges into paths.
The outputs are the indices into "paths" ending each path. Past the last path, the indices
are set to -1.
path_join does not process bad values. It will set the bad-value flag of all output
ndarrays if the flag is set for any of the input ndarrays.
path_segs
Divide a path into segments.
@segments = path_segs($pathindices, $paths);
Returns a series of slices of the "paths", such as those created by "path_join", stopping
at the first negative index. Currently does not broadcast.
use PDL;
use PDL::ImageND;
use PDL::Graphics::Simple;
$SIZE = 500;
$vals = rvals($SIZE,$SIZE)->divide($SIZE/12.5)->sin;
@cntr_threshes = zeroes(9)->xlinvals($vals->minmax)->list;
$win = pgswin();
$xrange = [0,$vals->dim(0)-1]; $yrange = [0,$vals->dim(1)-1];
$win->plot(with=>'image', $vals, {xrange=>$xrange,yrange=>$yrange,j=>1},);
for $thresh (@cntr_threshes) {
my ($segs, $cnt) = contour_segments($thresh, $vals, $vals->ndcoords);
my $segscoords = $segs->slice(',0:'.$cnt->max)->clump(-1)->splitdim(0,4);
$linesegs = $segscoords->splitdim(0,2);
$uniqcoords = $linesegs->uniqvec;
next if $uniqcoords->dim(1) < 2;
$indexed = vsearchvec($linesegs, $uniqcoords)->uniqvec;
@paths = path_segs(path_join($indexed, $uniqcoords->dim(1), 0));
@paths = map $uniqcoords->dice_axis(1, $_)->mv(0,-1), @paths;
$win->oplot(
(map +(with=>'lines', $_->dog), @paths),
{xrange=>$xrange,yrange=>$yrange,j=>1},
);
}
print "ret> "; <>;
AUTHORS
Copyright (C) Karl Glazebrook and Craig DeForest, 1997, 2003 All rights reserved. There is
no warranty. You are allowed to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL distribution. If this file is
separated from the PDL distribution, the copyright notice should be included in the file.