oracular (3) PDL::Image2D.3pm.gz
NAME
PDL::Image2D - Miscellaneous 2D image processing functions
DESCRIPTION
Miscellaneous 2D image processing functions - for want of anywhere else to put them.
SYNOPSIS
use PDL::Image2D;
FUNCTIONS
conv2d
Signature: (a(m,n); kern(p,q); [o]b(m,n); indx [t]mapi(isize=CALC($SIZE(p) + $SIZE(m))); indx [t]mapj(jsize=CALC($SIZE(q) + $SIZE(n))); int opt)
2D convolution of an array with a kernel (smoothing)
For large kernels, using a FFT routine, such as "fftconvolve" in PDL::FFT, will be quicker.
$new = conv2d $old, $kernel, {OPTIONS}
$smoothed = conv2d $image, ones(3,3), {Boundary => Reflect}
Boundary - controls what values are assumed for the image when kernel
crosses its edge:
=> Default - periodic boundary conditions
(i.e. wrap around axis)
=> Reflect - reflect at boundary
=> Truncate - truncate at boundary
=> Replicate - repeat boundary pixel values
Unlike the FFT routines, conv2d is able to process bad values.
med2d
Signature: (a(m,n); kern(p,q); [o]b(m,n); double+ [t]tmp(pq=CALC($SIZE(p)*$SIZE(q))); indx [t]mapi(isize=CALC($SIZE(p) + $SIZE(m))); indx [t]mapj(jsize=CALC($SIZE(q) + $SIZE(n))); int opt)
2D median-convolution of an array with a kernel (smoothing)
Note: only points in the kernel >0 are included in the median, other points are weighted by the kernel
value (medianing lots of zeroes is rather pointless)
$new = med2d $old, $kernel, {OPTIONS}
$smoothed = med2d $image, ones(3,3), {Boundary => Reflect}
Boundary - controls what values are assumed for the image when kernel
crosses its edge:
=> Default - periodic boundary conditions (i.e. wrap around axis)
=> Reflect - reflect at boundary
=> Truncate - truncate at boundary
=> Replicate - repeat boundary pixel values
Bad values are ignored in the calculation. If all elements within the kernel are bad, the output is set
bad.
med2df
Signature: (a(m,n); [o]b(m,n); indx [t]mapi(isize=CALC($SIZE(p) + $SIZE(m))); indx [t]mapj(jsize=CALC($SIZE(q) + $SIZE(n))); int p_size=>p; int q_size=>q; int opt)
2D median-convolution of an array in a pxq window (smoothing)
Note: this routine does the median over all points in a rectangular
window and is not quite as flexible as "med2d" in this regard
but slightly faster instead
$new = med2df $old, $xwidth, $ywidth, {OPTIONS}
$smoothed = med2df $image, 3, 3, {Boundary => Reflect}
Boundary - controls what values are assumed for the image when kernel
crosses its edge:
=> Default - periodic boundary conditions (i.e. wrap around axis)
=> Reflect - reflect at boundary
=> Truncate - truncate at boundary
=> Replicate - repeat boundary pixel values
med2df 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.
box2d
Signature: (a(n,m); [o] b(n,m); int wx; int wy; int edgezero)
fast 2D boxcar average
$smoothim = $im->box2d($wx,$wy,$edgezero=1);
The edgezero argument controls if edge is set to zero (edgezero=1) or just keeps the original
(unfiltered) values.
"box2d" should be updated to support similar edge options as "conv2d" and "med2d" etc.
Boxcar averaging is a pretty crude way of filtering. For serious stuff better filters are around (e.g.,
use "conv2d" with the appropriate kernel). On the other hand it is fast and computational cost grows only
approximately linearly with window size.
box2d 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.
patch2d
Signature: (a(m,n); int bad(m,n); [o]b(m,n))
patch bad pixels out of 2D images using a mask
$patched = patch2d $data, $bad;
$bad is a 2D mask array where 1=bad pixel 0=good pixel. Pixels are replaced by the average of their non-
bad neighbours; if all neighbours are bad, the original data value is copied across.
This routine does not handle bad values - use "patchbad2d" instead
patchbad2d
Signature: (a(m,n); [o]b(m,n))
patch bad pixels out of 2D images containing bad values
$patched = patchbad2d $data;
Pixels are replaced by the average of their non-bad neighbours; if all neighbours are bad, the output is
set bad. If the input ndarray contains no bad values, then a straight copy is performed (see "patch2d").
patchbad2d handles bad values. The output ndarray may contain bad values, depending on the pattern of bad
values in the input ndarray.
max2d_ind
Signature: (a(m,n); [o]val(); int [o]x(); int[o]y())
Return value/position of maximum value in 2D image
Contributed by Tim Jenness
Bad values are excluded from the search. If all pixels are bad then the output is set bad.
centroid2d
Signature: (im(m,n); x(); y(); box(); [o]xcen(); [o]ycen())
Refine a list of object positions in 2D image by centroiding in a box
$box is the full-width of the box, i.e. the window is "+/- $box/2".
Bad pixels are excluded from the centroid calculation. If all elements are bad (or the pixel sum is 0 -
but why would you be centroiding something with negatives in...) then the output values are set bad.
crop
Return bounding box of given mask in an "indx" ndarray, so it can broadcast. Use other operations (such
as "isgood" in PDL::Bad, or "eqvec" in PDL::Primitive with a colour vector) to create a mask suitable for
your application.
$x1x2y1y2 = crop($image);
cc8compt
Connected 8-component labeling of a binary image.
Connected 8-component labeling of 0,1 image - i.e. find separate segmented objects and fill object pixels
with object number. 8-component labeling includes all neighboring pixels. This is just a front-end to
ccNcompt. See also "cc4compt".
$segmented = cc8compt( $image > $threshold );
cc4compt
Connected 4-component labeling of a binary image.
Connected 4-component labeling of 0,1 image - i.e. find separate segmented objects and fill object pixels
with object number. 4-component labling does not include the diagonal neighbors. This is just a front-
end to ccNcompt. See also "cc8compt".
$segmented = cc4compt( $image > $threshold );
ccNcompt
Signature: (a(m,n); int+ [o]b(m,n); int con)
Connected component labeling of a binary image.
Connected component labeling of 0,1 image - i.e. find separate segmented objects and fill object pixels
with object number. See also "cc4compt" and "cc8compt".
The connectivity parameter must be 4 or 8.
$segmented = ccNcompt( $image > $threshold, 4);
$segmented2 = ccNcompt( $image > $threshold, 8);
where the second parameter specifies the connectivity (4 or 8) of the labeling.
ccNcompt ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output
ndarrays if the flag is set for any of the input ndarrays.
polyfill
fill the area of the given polygon with the given colour.
This function works inplace, i.e. modifies "im".
polyfill($im,$ps,$colour,[\%options]);
The default method of determining which points lie inside of the polygon used is not as strict as the
method used in "pnpoly". Often, it includes vertices and edge points. Set the "Method" option to change
this behaviour.
Method - Set the method used to determine which points lie in the polygon.
=> Default - internal PDL algorithm
=> pnpoly - use the "pnpoly" algorithm
# Make a convex 3x3 square of 1s in an image using the pnpoly algorithm
$ps = pdl([3,3],[3,6],[6,6],[6,3]);
polyfill($im,$ps,1,{'Method' =>'pnpoly'});
pnpoly
'points in a polygon' selection from a 2-D ndarray
$mask = $img->pnpoly($ps);
# Old style, do not use
$mask = pnpoly($x, $y, $px, $py);
For a closed polygon determined by the sequence of points in {$px,$py} the output of pnpoly is a mask
corresponding to whether or not each coordinate (x,y) in the set of test points, {$x,$y}, is in the
interior of the polygon. This is the 'points in a polygon' algorithm from
<http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html> and vectorized for PDL by Karl
Glazebrook.
# define a 3-sided polygon (a triangle)
$ps = pdl([3, 3], [20, 20], [34, 3]);
# $tri is 0 everywhere except for points in polygon interior
$tri = $img->pnpoly($ps);
With the second form, the x and y coordinates must also be specified. THIS IS MAINTAINED FOR BACKWARD
COMPATIBILITY ONLY .
$px = pdl( 3, 20, 34 );
$py = pdl( 3, 20, 3 );
$x = $img->xvals; # get x pixel coords
$y = $img->yvals; # get y pixel coords
# $tri is 0 everywhere except for points in polygon interior
$tri = pnpoly($x,$y,$px,$py);
polyfillv
return the (dataflowed) area of an image described by a polygon
polyfillv($im,$ps,[\%options]);
The default method of determining which points lie inside of the polygon used is not as strict as the
method used in "pnpoly". Often, it includes vertices and edge points. Set the "Method" option to change
this behaviour.
Method - Set the method used to determine which points lie in the polygon.
=> Default - internal PDL algorithm
=> pnpoly - use the "pnpoly" algorithm
# increment intensity in area bounded by $poly using the pnpoly algorithm
$im->polyfillv($poly,{'Method'=>'pnpoly'})++; # legal in perl >= 5.6
# compute average intensity within area bounded by $poly using the default algorithm
$av = $im->polyfillv($poly)->avg;
rot2d
Signature: (im(m,n); float angle(); bg(); int aa(); [o] om(p,q))
rotate an image by given "angle"
# rotate by 10.5 degrees with antialiasing, set missing values to 7
$rot = $im->rot2d(10.5,7,1);
This function rotates an image through an "angle" between -90 and + 90 degrees. Uses/doesn't use
antialiasing depending on the "aa" flag. Pixels outside the rotated image are set to "bg".
Code modified from pnmrotate (Copyright Jef Poskanzer) with an algorithm based on "A Fast Algorithm for
General Raster Rotation" by Alan Paeth, Graphics Interface '86, pp. 77-81.
Use the "rotnewsz" function to find out about the dimension of the newly created image
($newcols,$newrows) = rotnewsz $oldn, $oldm, $angle;
PDL::Transform offers a more general interface to distortions, including rotation, with various types of
sampling; but rot2d is faster.
rot2d ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output
ndarrays if the flag is set for any of the input ndarrays.
bilin2d
Signature: (Int(n,m); [io] O(q,p))
Bilinearly maps the first ndarray in the second. The interpolated values are actually added to the second
ndarray which is supposed to be larger than the first one.
bilin2d ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output
ndarrays if the flag is set for any of the input ndarrays.
rescale2d
Signature: (Int(m,n); [io] O(p,q))
The first ndarray is rescaled to the dimensions of the second (expanding or meaning values as needed) and
then added to it in place. Nothing useful is returned.
If you want photometric accuracy or automatic FITS header metadata tracking, consider using
PDL::Transform::map instead: it does these things, at some speed penalty compared to rescale2d.
rescale2d ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output
ndarrays if the flag is set for any of the input ndarrays.
fitwarp2d
Find the best-fit 2D polynomial to describe a coordinate transformation.
( $px, $py ) = fitwarp2d( $x, $y, $u, $v, $nf, { options } )
Given a set of points in the output plane ("$u,$v"), find the best-fit (using singular-value
decomposition) 2D polynomial to describe the mapping back to the image plane ("$x,$y"). The order of the
fit is controlled by the $nf parameter (the maximum power of the polynomial is "$nf - 1"), and you can
restrict the terms to fit using the "FIT" option.
$px and $py are "np" by "np" element ndarrays which describe a polynomial mapping (of order "np-1") from
the output "(u,v)" image to the input "(x,y)" image:
x = sum(j=0,np-1) sum(i=0,np-1) px(i,j) * u^i * v^j
y = sum(j=0,np-1) sum(i=0,np-1) py(i,j) * u^i * v^j
The transformation is returned for the reverse direction (ie output to input image) since that is what is
required by the "warp2d" routine. The "applywarp2d" routine can be used to convert a set of "$u,$v"
points given $px and $py.
Options:
FIT - which terms to fit? default ones(byte,$nf,$nf)
FIT "FIT" allows you to restrict which terms of the polynomial to fit: only those terms for which the FIT
ndarray evaluates to true will be evaluated. If a 2D ndarray is sent in, then it is used for the x
and y polynomials; otherwise "$fit->slice(":,:,(0)")" will be used for $px and
"$fit->slice(":,:,(1)")" will be used for $py.
The number of points must be at least equal to the number of terms to fit ("$nf*$nf" points for the
default value of "FIT").
# points in original image
$x = pdl( 0, 0, 100, 100 );
$y = pdl( 0, 100, 100, 0 );
# get warped to these positions
$u = pdl( 10, 10, 90, 90 );
$v = pdl( 10, 90, 90, 10 );
#
# shift of origin + scale x/y axis only
$fit = byte( [ [1,1], [0,0] ], [ [1,0], [1,0] ] );
( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2, { FIT => $fit } );
print "px = ${px}py = $py";
px =
[
[-12.5 1.25]
[ 0 0]
]
py =
[
[-12.5 0]
[ 1.25 0]
]
#
# Compared to allowing all 4 terms
( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2 );
print "px = ${px}py = $py";
px =
[
[ -12.5 1.25]
[ 1.110223e-16 -1.1275703e-17]
]
py =
[
[ -12.5 1.6653345e-16]
[ 1.25 -5.8546917e-18]
]
# A higher-degree polynomial should not affect the answer much, but
# will require more control points
$x = $x->glue(0,pdl(50,12.5, 37.5, 12.5, 37.5));
$y = $y->glue(0,pdl(50,12.5, 37.5, 37.5, 12.5));
$u = $u->glue(0,pdl(73,20,40,20,40));
$v = $v->glue(0,pdl(29,20,40,40,20));
( $px3, $py3 ) = fitwarp2d( $x, $y, $u, $v, 3 );
print "px3 =${px3}py3 =$py3";
px3 =
[
[-6.4981162e+08 71034917 -726498.95]
[ 49902244 -5415096.7 55945.388]
[ -807778.46 88457.191 -902.51612]
]
py3 =
[
[-6.2732159e+08 68576392 -701354.77]
[ 48175125 -5227679.8 54009.114]
[ -779821.18 85395.681 -871.27997]
]
#This illustrates an important point about singular value
#decompositions that are used in fitwarp2d: like all SVDs, the
#rotation matrices are not unique, and so the $px and $py returned
#by fitwarp2d are not guaranteed to be the "simplest" solution.
#They do still work, though:
($x3,$y3) = applywarp2d($px3,$py3,$u,$v);
print approx $x3,$x,1e-4;
[1 1 1 1 1 1 1 1 1]
print approx $y3,$y;
[1 1 1 1 1 1 1 1 1]
applywarp2d
Transform a set of points using a 2-D polynomial mapping
( $x, $y ) = applywarp2d( $px, $py, $u, $v )
Convert a set of points (stored in 1D ndarrays "$u,$v") to "$x,$y" using the 2-D polynomial with
coefficients stored in $px and $py. See "fitwarp2d" for more information on the format of $px and $py.
warp2d
Signature: (img(m,n); ldouble px(np,np); ldouble py(np,np); [o] warp(m,n); ldouble [t] poly(np); ldouble [t] kernel(ns); char *kernel_type; double noval; int nsamples => ns)
Warp a 2D image given a polynomial describing the reverse mapping.
$out = warp2d( $img, $px, $py, { options } );
Apply the polynomial transformation encoded in the $px and $py ndarrays to warp the input image $img into
the output image $out.
The format for the polynomial transformation is described in the documentation for the "fitwarp2d"
routine.
At each point "x,y", the closest 16 pixel values are combined with an interpolation kernel to calculate
the value at "u,v". The interpolation is therefore done in the image, rather than Fourier, domain. By
default, a "tanh" kernel is used, but this can be changed using the "KERNEL" option discussed below (the
choice of kernel depends on the frequency content of the input image).
The routine is based on the "warping" command from the Eclipse data-reduction package - see
http://www.eso.org/eclipse/ - and for further details on image resampling see Wolberg, G., "Digital Image
Warping", 1990, IEEE Computer Society Press ISBN 0-8186-8944-7).
Currently the output image is the same size as the input one, which means data will be lost if the
transformation reduces the pixel scale. This will (hopefully) be changed soon.
$img = rvals(byte,501,501);
imag $img, { JUSTIFY => 1 };
#
# use a not-particularly-obvious transformation:
# x = -10 + 0.5 * $u - 0.1 * $v
# y = -20 + $v - 0.002 * $u * $v
#
$px = pdl( [ -10, 0.5 ], [ -0.1, 0 ] );
$py = pdl( [ -20, 0 ], [ 1, 0.002 ] );
$wrp = warp2d( $img, $px, $py );
#
# see the warped image
imag $warp, { JUSTIFY => 1 };
The options are:
KERNEL - default value is tanh
NOVAL - default value is 0
"KERNEL" is used to specify which interpolation kernel to use (to see what these kernels look like, use
the "warp2d_kernel" routine). The options are:
tanh
Hyperbolic tangent: the approximation of an ideal box filter by the product of symmetric tanh
functions.
sinc
For a correctly sampled signal, the ideal filter in the fourier domain is a rectangle, which produces
a "sinc" interpolation kernel in the spatial domain:
sinc(x) = sin(pi * x) / (pi * x)
However, it is not ideal for the "4x4" pixel region used here.
sinc2
This is the square of the sinc function.
lanczos
Although defined differently to the "tanh" kernel, the result is very similar in the spatial domain.
The Lanczos function is defined as
L(x) = sinc(x) * sinc(x/2) if abs(x) < 2
= 0 otherwise
hann
This kernel is derived from the following function:
H(x) = a + (1-a) * cos(2*pi*x/(N-1)) if abs(x) < 0.5*(N-1)
= 0 otherwise
with "a = 0.5" and N currently equal to 2001.
hamming
This kernel uses the same H(x) as the Hann filter, but with "a = 0.54".
"NOVAL" gives the value used to indicate that a pixel in the output image does not map onto one in the
input image.
warp2d ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all output
ndarrays if the flag is set for any of the input ndarrays.
warp2d_kernel
Signature: ([o] x(n); [o] k(n); ldouble [t] kernel(n); char *name; PDL_Indx nsize => n)
Return the specified kernel, as used by "warp2d"
( $x, $k ) = warp2d_kernel( $name )
The valid values for $name are the same as the "KERNEL" option of "warp2d".
line warp2d_kernel( "hamming" );
warp2d_kernel ignores the bad-value flag of the input ndarrays. It will set the bad-value flag of all
output ndarrays if the flag is set for any of the input ndarrays.
AUTHORS
Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams (rjrw@ast.leeds.ac.uk), Tim Jenness (timj@jach.hawaii.edu), and Doug Burke (burke@ifa.hawaii.edu). 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.