|
PDL::Transform - Image transformations and N-D functions |
PDL::Transform - Image transformations and N-D functions
use PDL::Transform;
my $t = new PDL::Transform::<type>(<opt>)
$out = $t->apply($in) # Apply transform to some N-vectors (Transform method) $out = $in->apply($t) # Apply transform to some N-vectors (PDL method)
$im1 = $t->map($im); # Transform image coordinates (Transform method) $im1 = $im->map($t); # Transform image coordinates (PDL method)
$t2 = $t->compose($t1); # compose two transforms $t2 = $t x $t1; # compose two transforms (by analogy to matrix mult.)
$t3 = $t2->inverse(); # invert a transform $t3 = !$t2; # invert a transform (by analogy to logical "not")
PDL::Transform is a convenient way to represent coordinate transformations and resample images. It embodies functions mapping R^N -> R^M, both with and without inverses. Provision exists for parametrizing functions, and for composing them. You can use this part of the Transform object to keep track of arbitrary functions mapping R^N -> R^M with or without inverses.
The simplest way to use a Transform object is to transform vector data between coordinate systems. The apply method accepts a PDL whose 0th dimension is coordinate index (all other dimensions are threaded over) and transforms the vectors into the new coordinate system.
Transform also includes image resampling, via the map method, You define a coordinate transform using a Transform object, then use it to remap an image PDL. The output is a remapped, resampled image.
You can define and compose several transformations, then apply them all at once to an image. The image is interpolated only once, when all the composed transformations are applied.
In keeping with standard practice, but somewhat counterintuitively, the map engine uses the inverse transform to map coordinates FROM the destination dataspace (or image plane) TO the source dataspace; hence PDL::Transform keeps track of both the forward and inverse transform.
For terseness and convenience, most of the constructors are exported
into the current package with the name t_<transform>, so the following
(for example) are synonyms:
$t = new PDL::Transform::Radial(); # Long way
$t = t_radial(); # Short way
Several math operators are overloaded, so that you can compose and invert functions with expression syntax instead of method syntax (see below).
Coordinate transformations and mappings are a little counterintuitive at first. Here are some examples of transforms in action:
use PDL::Transform;
$a = rfits('m51.fits'); # Substitute path if necessary!
$ts = t_linear(Scale=>3); # Scaling transform
$w = pgwin(xs); $w->imag($a);
## Grow m51 by a factor of 3; origin is at lower left.
$b = $ts->map($a,{pix=>1}); # pix option uses direct pixel coord system
$w->imag($b);
## Shrink m51 by a factor of 3; origin still at lower left.
$c = $ts->unmap($a, {pix=>1});
$w->imag($c);
## Grow m51 by a factor of 3; origin is at scientific origin. $d = $ts->map($a,$a->hdr); # FITS hdr template prevents autoscaling $w->imag($d);
## Shrink m51 by a factor of 3; origin is still at sci. origin. $e = $ts->unmap($a,$a->hdr); $w->imag($e);
## A no-op: shrink m51 by a factor of 3, then autoscale back to size $f = $ts->map($a); # No template causes autoscaling of output
$f->inverse()->compose($g)->compose($f) # long way !$f x $g x $f # short way
Both of those expressions are equivalent to the mathematical expression f^-1 o g o f, or f^-1(g(f(x))).
compose($f) # long way
$f**3 # short way
Transforms are perl hashes. Here's a list of the meaning of each key:
funcs should support
inplace operation to save memory when the data are flagged
inplace. But func should always return its result even when
flagged to compute in-place.
func should treat the 0th dimension of its input as a dimensional
index (running 0..N-1 for R^N operation) and thread over all other input
dimensions.
compose() method.
Transforms should be inplace-aware where possible, to prevent excessive memory usage.
If you define a new type of transform, consider generating a new stringify method for it. Just define the sub ``stringify'' in the subclass package. It should call SUPER::stringify to generate the first line (though the PDL::Transform::Composition bends this rule by tweaking the top-level line), then output (indented) additional lines as necessary to fully describe the transformation.
Transforms have a mechanism for labeling the units and type of each coordinate, but it is just advisory. A routine to identify and, if necessary, modify units by scaling would be a good idea. Currently, it just assumes that the coordinates are correct for (e.g.) FITS scientific-to-pixel transformations.
Composition works OK but should probably be done in a more sophisticated way so that, linear transformations are combined at the matrix level instead of just strung together pixel-to-pixel.
Linear transformations should probably be handled with heterogeneous Hershey matrices, rather than matrix-multiply-and-add operations.
There are both operators and constructors. The constructors are all exported, all begin with ``t_'', and all return objects that are subclasses of PDL::Transform.
The apply, invert, map,
and unmap methods are also exported to the PDL package: they
are both Transform methods and PDL methods.
Signature: (data(); PDL::Transform t)
$out = $data->apply($t); $out = $t->apply($data);
Apply a transformation to some input coordinates.
In the example, $t is a PDL::Transform and $data is a piddle to
be interpreted as a collection of N-vectors (with index in the 0th
dimension). The output is a similar but transformed piddle.
For convenience, this is both a PDL method and a Transform method.
Signature: (data(); PDL::Transform t)
$out = $t->invert($data); $out = $data->invert($t);
Apply an inverse transformation to some input coordinates.
In the example, $t is a PDL::Transform and $data is a piddle to
be interpreted as a collection of N-vectors (with index in the 0th
dimension). The output is a similar piddle.
For convenience this is both a PDL method and a PDL::Transform method.
Signature: (data(); PDL::Transform a; template(); \%opt)
$output = $input->map($transform,[<template>],[<options>]); $output = $transform->map($input,[<template>],[<options>]);
Resample an image or N-D dataset using a coordinate transform.
The data are resampled so that the new pixel indices are proportional to the transformed coordinates rather than the original ones.
The operation uses the inverse transform: each output pixel location is inverse-transformed back to a location in the original dataset, and the value is interpolated or sampled appropriately and copied into the output domain. A variety of sampling options are available, trading off speed and mathematical correctness.
For convenience, this is both a PDL method and a PDL::Transform method.
map is FITS-aware: if there is a FITS header in the input data,
then the coordinate transform acts on the scientific coordinate system
rather than the pixel coordinate system.
By default, the output coordinates are separated from pixel coordinates
by a single layer of indirection. You can specify the mapping between
output transform (scientific) coordinates to pixel coordinates using
the orange and irange options (see below), or by supplying a
FITS header in the template.
If you don't specify an output transform, then the output is
autoscaled: map transforms a few vectors in the forward direction
to generate a mapping that will put most of the data on the image
plane, for most transformations. The calculated mapping gets stuck in the
output's FITS header.
You can operate in pixel space, and avoid autoscaling of the output,
by setting the nofits option (see below).
The output has the same data type as the input. This is a feature, but it can lead to strange-looking banding behaviors if you use interpolation on an integer input variable.
The template can be one of:
OPTIONS:
The following options are interpreted:
pdl([[-1,2],[3,4]])
limites the map to the quadrilateral in input space defined by the
four points (-1,3), (-1,4), (2,4), and (2,3).
As with plain autoscaling, the quadrilateral gets sparsely sampled by the autoranger, so pathological transformations can give you strange results.
This parameter is overridden by orange, below.
irange, except that it specifies
a quadrilateral in output space. Since the output pixel array is itself
a quadrilateral, you get pretty much exactly what you asked for.
This parameter overrides irange, if both are specified.
=item * s, sample (default for integers)
Pixel values in the output plane are sampled from the closest data value in the input plane. This is very fast but not very accurate for either magnification or decimation (shrinking). It is the default for templates of integer type.
JACOBIAN TRACKING:
This method of interpolation gives photometrically accurate resampling
of the input data for arbitrary transformations. The Jacobian of the
reverse transformation is the matrix J_ij = d X_i / d x_j, where i
and j are index variables, the X_i are the input-plane coordinates,
and the x_j are the output-plane coordinates. At each pixel, the code
generates a linear approximation to the transformation using the local
discrete Jacobian. The output pixels are treated as circles of radius
1.0, and transformed via the linear approximation to ellipses in the
input plane. The singular values of the Jacobian are padded to a
minimum of 1.4 pixels, ensuring that the transformed ellipses are fat
enough to encounter at least one sample point in the input plane.
To avoid numerical runaway, there are some limitations on the reverse-
transformed ellipses. In particular, the computational efficiency
scales inversely as the ratio between the largest and smallest
eigenvalues, so the ellipses are not allowed to get too eccentric. The
maximum eccentricity is given in the eccentricity option, and
defaults to 4.0.
A caveat about Jacobian remapping is that it assumes the transformation is continuous. Transformations that contain discontinuities will give incorrect results near the discontinuity. In particular, the 180th meridian isn't handled well in lat/lon mapping transformations (see the PDL::Transform::Cartography manpage) -- pixels along the 180th meridian get the average value of everything along the parallel occupied by the pixel. This flaw is inherent in the assumptions that underly creating a Jacobian matrix. Maybe someone will write code to work around it. Maybe that someone is you.
NOTES:
Jacobian tracking can be a memory hog, especially if the transformation includes very large regions outside of the original input plane. Some sort of memory guard needs to be put in place.
Signature: (data(); PDL::Transform a; template(); \%opt)
$out_image = $in_image->unmap($t,[<options>],[<template>]); $out_image = $t->unmap($in_image,[<options>],[<template>]);
Map an image or N-D dataset using the inverse as a coordinate transform.
This convenience function just inverts $t and calls map on the inverse; everything works the same otherwise. For convenience, it is both a PDL method and a PDL::Transform method.
$t2 = t_inverse($t); $t2 = $t->inverse; $t2 = $t ** -1; $t2 = !$t;
Return the inverse of a PDL::Transform. This just reverses the func/inv, idim/odim, itype/otype, and iunit/ounit pairs. Note that sometimes you end up with a transform that cannot be applied or mapped, because either the mathematical inverse doesn't exist or the inverse func isn't implemented.
You can invert a transform by raising it to a negative power, or by negating it with '!'.
The inverse transform remains connected to the main transform because they both point to the original parameters hash. That turns out to be useful.
$f2 = t_compose($f, $g,[...]); $f2 = $f->compose($g[,$h,$i,...]); $f2 = $f x $g x ...;
Function composition: f(g(x)), f(g(h(x))), ...
You can also compose transforms using the overloaded matrix-multiplication (nee repeat) operator 'x'.
This is accomplished by inserting a splicing code ref into the func
and inv slots. It combines multiple compositions into a single
list of transforms to be executed in order, fram last to first (in
keeping with standard mathematical notation). If one of the functions is
itself a composition, it is interpolated into the list rather than left
separate. Ultimately, linear transformations may also be combined within
the list.
No checking is done that the itype/otype and iunit/ounit fields are compatible -- that may happen later, or you can implement it yourself if you like.
$g1fg = $f->wrap($g); $g1fg = t_wrap($f,$g);
Shift a transform into a different space by 'wrapping' it with a second.
This is just a convenience function for two
compose calls. $a-wrap($b)> is the same as
(!$b) x $a x $b: the resulting transform first hits the data with
$b, then with $a, then with the inverse of $b.
For example, to shift the origin of rotation, do this:
$im = rfits('m51.fits');
$tf = t_fits($im);
$tr = t_linear({rot=>30});
$im1 = $tr->map($tr); # Rotate around pixel origin
$im2 = $tr->map($tr->wrap($tf)); # Rotate round FITS scientific origin
my $xform = t_identity my $xform = new PDL::Transform;
Generic constructor generates the identity transform.
This constructor really is trivial -- it is mainly used by the other transform constructors. It takes no parameters and returns the identity transform.
$f = t_lookup($lookup, {<options>});
Transform by lookup into an explicit table.
You specify an N+1-D PDL that is interpreted as an N-D lookup table of column vectors (vector index comes last). The last dimension has order equal to the output dimensionality of the transform.
For added flexibility in data space, You can specify pre-lookup linear scaling and offset of the data. Of course you can specify the interpolation method to be used. The linear scaling stuff is a little primitive; if you want more, try composing the linear transform with this one.
The prescribed values in the lookup table are treated as pixel-centered: that is, if your input array has N elements per row then valid data exist between the locations (-0.5) and (N-0.5) in lookup pixel space, because the pixels (which are numbered from 0 to N-1) are centered on their locations.
Lookup is done using interpND, so the boundary conditions and threading behaviour follow from that.
The indexed-over dimensions come first in the table, followed by a
single dimension containing the column vector to be output for each
set of other dimensions -- ie to output 2-vectors from 2 input
parameters, each of which can range from 0 to 49, you want an index
that has dimension list (50,50,2). For the identity lookup table
you could use cat(xvals(50,50),yvals(50,50)).
If you want to output a single value per input vector, you still need
that last index threading dimension -- if necessary, use dummy(-1,1).
The lookup index scaling is: out = lookup[ (scale * data) + offset ].
The inverse transform is calculated.
Options are listed below; there are several synonyms for each.
EXAMPLE
To scale logarithmically the Y axis of m51, try:
$a = rfits('m51.fits');
$lookup = xvals(256,256) -> cat( 10**(yvals(256,256)/100) * 256/10**2.55 );
$t = t_lookup($lookup);
$b = $t->map($a);
To do the same thing but with a smaller lookup table, try:
$lookup = 16 * xvals(17,17)->cat(10**(yvals(17,17)/(100/16)) * 16/10**2.55);
$t = t_lookup($lookup,{scale=>1/16.0});
$b = $t->map($a);
(Notice that, although the lookup table coordinates are is divided by 16, it is a 17x17 -- so linear interpolation works right to the edge of the original domain.)
NOTES
Inverses are not yet implemented -- the best way to do it might be by
judicious use of map() on the forward transformation.
the type/unit fields are ignored.
$f = t_linear({options});
Heterogeneous linear coordinate transformations.
You specify the linear transformation with pre-offset, a mixing matrix, and a post-offset. That overspecifies the transformation, so you can choose your favorite method to specify the transform you want. The inverse transform is automagically generated, provided that it actually exists (the transform matrix is invertible). Otherwise, the inverse transform just croaks.
Extra dimensions in the input vector are ignored, so if you pass a 3xN vector into a 3-D linear transformation, the final dimension is passed through unchanged.
The options you can usefully pass in are:
The rotation matrix is left-multiplied with the transformation matrix you specify, so that if you specify both rotation and a general matrix the rotation happens after the more general operation -- though that is deprecated.
Of course, you can duplicate this functionality -- and get more
general -- by generating your own rotation matrix and feeding it in
with the matrix option.
pre and post options, those define
the number of dimensions. But if you only supply scalars, there is no way
to tell and the default number of dimensions is 2. This provides a way
to do, e.g., 3-D scaling: just set {s=<scale-factor>, dims=>3}> and
you are on your way.
NOTES
the type/unit fields are currently ignored by t_linear.
$f = t_scale(<scale>)
Convenience interface to t_linear.
t_scale produces a tranform that scales around the origin by a fixed
amount. It acts exactly the same as t_linear(Scale=\<scale\>)>.
$f = t_offset(<shift>)
Convenience interface to t_linear.
t_offset produces a transform that shifts the origin to a new location.
It acts exactly the same as t_linear(Pre=\<shift\>)>.
$f = t_rot(<rotation-in-degrees>)
Convenience interface to t_linear.
t_rot produces a rotation transform in 2-D (scalar), 3-D (3-vector), or
N-D (matrix). It acts exactly the same as t_linear(Rot=\<shift\>)>.
$f = t_fits($fits,[option]);
FITS pixel-to-scientific transformation with inverse
You feed in a hash ref or a PDL with one of those as a header, and you get back a transform that converts 0-originated, pixel-centered coordinates into scientific coordinates via the transformation in the FITS header. For most FITS headers, the transform is reversible, so applying the inverse goes the other way. This is just a convenience subclass of PDL::Transform::Linear, but with unit/type support using the FITS header you supply.
For now, this transform is rather limited -- it really ought to accept units differences and stuff like that, but they are just ignored for now. Probably that would require putting units into the whole transform framework.
This transform implements the linear transform part of the WCS FITS standard outlined in Greisen & Calabata 2002 (A&A in press; find it at ``http://arxiv.org/abs/astro-ph/0207407'').
As a special case, you can pass in the boolean option ``ignore_rgb'' (default 0), and if you pass in a 3-D FITS header in which the last dimension has exactly 3 elements, it will be ignored in the output transformation. That turns out to be handy for handling rgb images.
$f = t_code(<func>,[<inv>],[options]);
Transform implementing arbitrary perl code.
This is a way of getting quick-and-dirty new transforms. You pass in anonymous (or otherwise) code refs pointing to subroutines that implement the forward and, optionally, inverse transforms. The subroutines should accept a data PDL followed by a parameter hash ref, and return the transformed data PDL. The parameter hash ref can be set via the options, if you want to.
Options that are accepted are:
The code variables are executable perl code, either as a code ref or as a string that will be eval'ed to produce code refs. If you pass in a string, it gets eval'ed at call time to get a code ref. If it compiles OK but does not return a code ref, then it gets re-evaluated with ``sub { ... }'' wrapped around it, to get a code ref.
Note that code callbacks like this can be used to do really weird things and generate equally weird results -- caveat scriptor!
$f = t_radial(<options>);
Convert Cartesian to radial/cylindrical coordinates. (2-D/3-D; with inverse)
Converts 2-D Cartesian to radial (theta,r) coordinates. You can choose direct or conformal conversion. Direct conversion preserves radial distance from the origin; conformal conversion preserves local angles, so that each small-enough part of the image only appears to be scaled and rotated, not stretched. Conformal conversion puts the radius on a logarithmic scale, so that scaling of the original image plane is equivalent to a simple offset of the transformed image plane.
If you use three or more dimensions, the higher dimensions are ignored, yielding a conversion from Cartesian to cylindrical coordinates, which is why there are two aliases for the same transform. If you use higher dimensionality than 2, you must manually specify the origin or you will get dimension mismatch errors when you apply the transform.
Theta runs clockwise instead of the more usual counterclockwise; that is to preserve the mirror sense of small structures.
OPTIONS:
ln(r/r0)) coordinates out. Theta is in radians, and the
radial coordinate varies by 1 for each e-folding of the r0-scaled
distance from the input origin. The logarithmic scaling is useful for
viewing both large and small things at the same time, and for keeping
shapes of small things preserved in the image.
EXAMPLES
These examples do transformations back into the same size image as they started from; by suitable use of the ``transform'' option to unmap you can send them to any size array you like.
Examine radial structure in M51: Here, we scale the output to stretch 2*pi radians out to the full image width in the horizontal direction, and to stretch 1 radius out to a diameter in the vertical direction.
$a = rfits('m51.fits');
$ts = t_linear(s => [250/2.0/3.14159, 2]); # Scale to fill orig. image
$tu = t_radial(o => [130,130]); # Expand around galactic core
$b = $a->map($ts x $tu);
Examine radial structure in M51 (conformal): Here, we scale the output to stretch 2*pi radians out to the full image width in the horizontal direction, and scale the vertical direction by the exact same amount to preserve conformality of the operation. Notice that each piece of the image looks ``natural'' -- only scaled and not stretched.
$a = rfits('m51.fits')
$ts = t_linear(s=> 250/2.0/3.14159); # Note scalar (heh) scale.
$tu = t_radial(o=> [130,130], r0=>5); # 5 pix. radius -> bottom of image
$b = $ts->compose($tu)->unmap($a);
$t = t_quadratic(<options>);
Quadratic scaling -- cylindrical pincushion (n-d; with inverse)
Quadratic scaling emulates pincushion in a cylindrical optical system: separate quadratic scaling is applied to each axis. You can apply separate distortion along any of the principal axes. If you want different axes, use wrap and t_linear to rotate them to the correct angle. The scaling options may be scalars or vectors; if they are scalars then the expansion is isotropic.
The formula for the expansion is:
f(a) = ( <a> + <strength> * a^2/<L_0> ) / (abs(<strength>) + 1)
where <strength> is a scaling coefficient and <L_0> is a fundamental length scale. Negative values of <strength> result in a pincushion contraction.
OPTIONS
$t = t_spherical(<options>);
Convert Cartesian to spherical coordinates. (3-D; with inverse)
Convert 3-D Cartesian to spherical (theta, phi, r) coordinates. Theta is longitude, centered on 0, and phi is latitude, also centered on 0. Unless you specify Euler angles, the pole points in the +Z direction and the prime meridian is in the +X direction. The default is for theta and phi to be in radians; you can select degrees if you want them.
Just as the t_radial 2-D transform acts like a 3-D cylindrical transform by ignoring third and higher dimensions, Spherical acts like a hypercylindrical transform in four (or higher) dimensions. Also as with t_radial, you must manually specify the origin if you want to use more dimensions than 3.
To deal with latitude & longitude on the surface of a sphere (rather than full 3-D coordinates), see t_unitsphere.
OPTIONS:
Copyright 2002 Craig DeForest. This module may be modified and distributed under the same terms as PDL itself. The module comes with NO WARRANTY.
|
PDL::Transform - Image transformations and N-D functions |