NAME¶
PDL::Transform - Coordinate transforms, image warping, and N-D functions
SYNOPSIS¶
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")
DESCRIPTION¶
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).
EXAMPLE¶
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
OPERATOR OVERLOADS¶
- '!'
- The bang is a unary inversion operator. It binds exactly as tightly as the
normal bang operator.
- 'x'
- By analogy to matrix multiplication, 'x' is the compose operator, so these
two expressions are equivalent:
$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))).
- '**'
- By analogy to numeric powers, you can apply an operator a positive integer
number of times with the ** operator:
$f->compose($f)->compose($f) # long way
$f**3 # short way
INTERNALS¶
Transforms are perl hashes. Here's a list of the meaning of each key:
- func
- Ref to a subroutine that evaluates the transformed coordinates. It's
called with the input coordinate, and the "params" hash. This
springboarding is done via explicit ref rather than by subclassing, for
convenience both in coding new transforms (just add the appropriate sub to
the module) and in adding custom transforms at run-time. Note that, if
possible, new "func"s 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.
- inv
- Ref to an inverse method that reverses the transformation. It must accept
the same "params" hash that the forward method accepts. This key
can be left undefined in cases where there is no inverse.
- idim, odim
- Number of useful dimensions for indexing on the input and output sides (ie
the order of the 0th dimension of the coordinates to be fed in or that
come out). If this is set to 0, then as many are allocated as needed.
- name
- A shorthand name for the transformation (convenient for debugging). You
should plan on using UNIVERAL::isa to identify classes of transformation,
e.g. all linear transformations should be subclasses of
PDL::Transform::Linear. That makes it easier to add smarts to, e.g., the
compose() method.
- itype
- An array containing the name of the quantity that is expected from the
input piddle for the transform, for each dimension. This field is
advisory, and can be left blank if there's no obvious quantity associated
with the transform. This is analogous to the CTYPEn field used in FITS
headers.
- oname
- Same as itype, but reporting what quantity is delivered for each
dimension.
- iunit
- The units expected on input, if a specific unit (e.g. degrees) is
expected. This field is advisory, and can be left blank if there's no
obvious unit associated with the transform.
- ounit
- Same as iunit, but reporting what quantity is delivered for each
dimension.
- params
- Hash ref containing relevant parameters or anything else the func needs to
work right.
- is_inverse
- Bit indicating whether the transform has been inverted. That is useful for
some stringifications (see the PDL::Transform::Linear stringifier), and
may be useful for other things.
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.
NOTES¶
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, for example, linear transformations are combined at the matrix level
instead of just strung together pixel-to-pixel.
FUNCTIONS¶
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.
FUNCTIONS¶
apply¶
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 PDL to be interpreted as a
collection of N-vectors (with index in the 0th dimension). The output is a
similar but transformed PDL.
For convenience, this is both a PDL method and a Transform method.
invert¶
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.
map¶
Signature: (k0(); SV *in; SV *out; SV *map; SV *boundary; SV *method;
SV *big; SV *blur; SV *sv_min; SV *flux; SV *bv)
PDL::match¶
$b = $a->match($c); # Match $c's header and size
$b = $a->match([100,200]); # Rescale to 100x200 pixels
$b = $a->match([100,200],{rect=>1}); # Rescale and remove rotation/skew.
Resample a scientific image to the same coordinate system as another.
The example above is syntactic sugar for
$b = $a->map(t_identity, $c, ...);
it resamples the input PDL with the identity transformation in scientific
coordinates, and matches the pixel coordinate system to $c's FITS header.
There is one difference between match and map: match makes the
"rectify" option to "map" default to 0, not 1. This only
affects matching where autoscaling is required (i.e. the array ref example
above). By default, that example simply scales $a to the new size and
maintains any rotation or skew in its scientiic-to-pixel coordinate transform.
map¶
$b = $a->map($xform,[<template>],[\%opt]); # Distort $a with tranform $xform
$b = $a->map(t_identity,[$pdl],[\%opt]); # rescale $a to match $pdl's dims.
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.
Autoscaling is especially useful for rescaling images -- if you specify the
identity transform and allow autoscaling, you duplicate the functionality of
rescale2d, but with more options for interpolation.
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:
- •
- a PDL
The PDL and its header are copied to the output array, which is then
populated with data. If the PDL has a FITS header, then the FITS transform
is automatically applied so that $t applies to the output scientific
coordinates and not to the output pixel coordinates. In this case the
NAXIS fields of the FITS header are ignored.
- •
- a FITS header stored as a hash ref
The FITS NAXIS fields are used to define the output array, and the FITS
transformation is applied to the coordinates so that $t applies to the
output scientific coordinates.
- •
- a list ref
This is a list of dimensions for the output array. The code estimates
appropriate pixel scaling factors to fill the available space. The scaling
factors are placed in the output FITS header.
- •
- nothing
In this case, the input image size is used as a template, and scaling is
done as with the list ref case (above).
OPTIONS:
The following options are interpreted:
- b, bound, boundary, Boundary (default = 'truncate')
- This is the boundary condition to be applied to the input image; it is
passed verbatim to range or interpND in the sampling or interpolating
stage. Other values are 'forbid','extend', and 'periodic'. You can
abbreviate this to a single letter. The default 'truncate' causes the
entire notional space outside the original image to be filled with 0.
- pix, Pixel, nf, nofits, NoFITS (default = 0)
- If you set this to a true value, then FITS headers and interpretation are
ignored; the transformation is treated as being in raw pixel
coordinates.
- j, J, just, justify, Justify (default = 0)
- If you set this to 1, then output pixels are autoscaled to have unit
aspect ratio in the output coordinates. If you set it to a non-1 value,
then it is the aspect ratio between the first dimension and all subsequent
dimensions -- or, for a 2-D transformation, the scientific pixel aspect
ratio. Values less than 1 shrink the scale in the first dimension compared
to the other dimensions; values greater than 1 enlarge it compared to the
other dimensions. (This is the same sense as in the PGPLOTinterface.)
- ir, irange, input_range, Input_Range
- This is a way to modify the autoscaling. It specifies the range of input
scientific (not necessarily pixel) coordinates that you want to be mapped
to the output image. It can be either a nested array ref or a piddle. The
0th dim (outside coordinate in the array ref) is dimension index in the
data; the 1st dim should have order 2. For example, passing in either
[[-1,2],[3,4]] or 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.
- or, orange, output_range, Output_Range
- This sets the window of output space that is to be sampled onto the output
array. It works exactly like "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. It
forces rectification of the output (so that scientific axes align with the
pixel grid).
- r, rect, rectify
- This option defaults TRUE and controls how autoscaling is performed. If it
is true or undefined, then autoscaling adjusts so that pixel coordinates
in the output image are proportional to individual scientific coordinates.
If it is false, then autoscaling automatically applies the inverse of any
input FITS transformation *before* autoscaling the pixels. In the special
case of linear transformations, this preserves the rectangular shape of
the original pixel grid and makes output pixel coordinate proportional to
input coordinate.
- m, method, Method
- This option controls the interpolation method to be used. Interpolation
greatly affects both speed and quality of output. For most cases the
option is directly passed to interpND for interpolation. Possible options,
in order from fastest to slowest, are:
- •
- s, sample (default for ints)
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.
- •
- l, linear (default for floats)
Pixel values are linearly interpolated from the closest data value in the
input plane. This is reasonably fast but only accurate for magnification.
Decimation (shrinking) of the image causes aliasing and loss of photometry
as features fall between the samples. It is the default for floating-point
templates.
- •
- c, cubic
Pixel values are interpolated using an N-cubic scheme from a 4-pixel N-cube
around each coordinate value. As with linear interpolation, this is only
accurate for magnification.
- •
- f, fft
Pixel values are interpolated using the term coefficients of the Fourier
transform of the original data. This is the most appropriate technique for
some kinds of data, but can yield undesired "ringing" for
expansion of normal images. Best suited to studying images with repetitive
or wavelike features.
- •
- h, hanning
Pixel values are filtered through a spatially-variable filter tuned to the
computed Jacobian of the transformation, with hanning-window (cosine)
pixel rolloff in each dimension. This prevents aliasing in the case where
the image is distorted or shrunk, but allows small amounts of aliasing at
pixel edges wherever the image is enlarged.
- •
- g, gaussian, j, jacobian
Pixel values are filtered through a spatially-variable filter tuned to the
computed Jacobian of the transformation, with radial Gaussian rolloff.
This is the most accurate resampling method, in the sense of introducing
the fewest artifacts into a properly sampled data set.
- blur, Blur (default = 1.0)
- This value scales the input-space footprint of each output pixel in the
gaussian and hanning methods. It's retained for historical reasons. Larger
values yield blurrier images; values significantly smaller than unity
cause aliasing.
- sv, SV (default = 1.0)
- This value lets you set the lower limit of the transformation's singular
values in the hanning and gaussian methods, limiting the minimum radius of
influence associated with each output pixel. Large numbers yield smoother
interpolation in magnified parts of the image but don't affect reduced
parts of the image.
- big, Big (default = 0.2)
- This is the largest allowable input spot size which may be mapped to a
single output pixel by the hanning and gaussian methods, in units of the
largest non-thread input dimension. (i.e. the default won't let you reduce
the original image to less than 5 pixels across). This places a limit on
how long the processing can take for pathological transformations. Smaller
numbers keep the code from hanging for a long time; larger numbers provide
for photometric accuracy in more pathological cases. Numbers larer than
1.0 are silly, because they allow the entire input array to be compressed
into a region smaller than a single pixel.
Wherever an output pixel would require averaging over an area that is too
big in input space, it instead gets NaN or the equivalent (bad values are
not yet supported).
- phot, photometry, Photometry
- This lets you set the style of photometric conversion to be used in the
hanning or gaussian methods. You may choose:
- •
- 0, s, surf, surface, Surface (default)
(this is the default): surface brightness is preserved over the
transformation, so features maintain their original intensity. This is
what the sampling and interpolation methods do.
- •
- 1, f, flux, Flux
Total flux is preserved over the transformation, so that the brightness
integral over image regions is preserved. Parts of the image that are
shrunk wind up brighter; parts that are enlarged end up fainter.
VARIABLE FILTERING:
The 'hanning' and 'gaussian' methods of interpolation give photometrically
accurate resampling of the input data for arbitrary transformations. At each
pixel, the code generates a linear approximation to the input transformation,
and uses that linearization to estimate the "footprint" of the
output pixel in the input space. The output value is a weighted average of the
appropriate input spaces.
A caveat about these methods is that they assume 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
PDL::Transform::Cartography) -- 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.
map does not process bad values. It will set the bad-value flag of all output
piddles if the flag is set for any of the input piddles.
unmap¶
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.
t_inverse¶
$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.
t_compose¶
$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.
t_wrap¶
$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 t_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
t_identity¶
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.
t_lookup¶
$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 ].
A simplistic table inversion routine is included. This means that you can (for
example) use the "map" method with "t_lookup"
transformations. But the table inversion is exceedingly slow, and not
practical for tables larger than about 100x100. The inversion table is
calculated in its entirety the first time it is needed, and then cached until
the object is destroyed.
Options are listed below; there are several synonyms for each.
- s, scale, Scale
- (default 1.0) Specifies the linear amount of scaling to be done before
lookup. You can feed in a scalar or an N-vector; other values may cause
trouble. If you want to save space in your table, then specify smaller
scale numbers.
- o, offset, Offset
- (default 0.0) Specifies the linear amount of offset before lookup. This is
only a scalar, because it is intended to let you switch to corner-centered
coordinates if you want to (just feed in o=-0.25).
- b, bound, boundary, Boundary
- Boundary condition to be fed to interpND
- m, method, Method
- Interpolation method to be fed to interpND
EXAMPLE
To scale logarithmically the Y axis of m51, try:
$a = float rfits('m51.fits');
$lookup = xvals(128,128) -> cat( 10**(yvals(128,128)/50) * 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.
t_linear¶
$f = t_linear({options});
Linear (affine) transformations with optional offset
t_linear implements simple matrix multiplication with offset, also known as the
affine 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:
- s, scale, Scale
- A scaling scalar (heh), vector, or matrix. If you specify a vector it is
treated as a diagonal matrix (for convenience). It gets left-multiplied
with the transformation matrix you specify (or the identity), so that if
you specify both a scale and a matrix the scaling is done after the
rotation or skewing or whatever.
- r, rot, rota, rotation, Rotation
- A rotation angle in degrees -- useful for 2-D and 3-D data only. If you
pass in a scalar, it specifies a rotation from the 0th axis toward the 1st
axis. If you pass in a 3-vector as either a PDL or an array ref (as in
"rot=>[3,4,5]"), then it is treated as a set of Euler angles
in three dimensions, and a rotation matrix is generated that does the
following, in order:
- •
- Rotate by rot->(2) degrees from 0th to 1st axis
- •
- Rotate by rot->(1) degrees from the 2nd to the 0th axis
- •
- Rotate by rot->(0) degrees from the 1st to the 2nd axis
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.
- m, matrix, Matrix
- The transformation matrix. It does not even have to be square, if you want
to change the dimensionality of your input. If it is invertible (note:
must be square for that), then you automagically get an inverse transform
too.
- pre, preoffset, offset, Offset
- The vector to be added to the data before they get multiplied by the
matrix (equivalent of CRVAL in FITS, if you are converting from scientific
to pixel units).
- post, postoffset, shift, Shift
- The vector to be added to the data after it gets multiplied by the matrix
(equivalent of CRPIX-1 in FITS, if youre converting from scientific to
pixel units).
- d, dim, dims, Dims
- Most of the time it is obvious how many dimensions you want to deal with:
if you supply a matrix, it defines the transformation; if you input offset
vectors in the "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.
t_scale¶
$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\>)>.
t_offset¶
$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\>)>.
t_rot¶
$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\>)>.
t_fits¶
$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.
t_code¶
$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:
- p,params
- The parameter hash that will be passed back to your code (defaults to the
empty hash).
- n,name
- The name of the transform (defaults to "code").
- i, idim (default 2)
- The number of input dimensions (additional ones should be passed through
unchanged)
- o, odim (default 2)
- The number of output dimensions
- itype
- The type of the input dimensions, in an array ref (optional and
advisiory)
- otype
- The type of the output dimension, in an array ref (optional and
advisory)
- iunit
- The units that are expected for the input dimensions (optional and
advisory)
- ounit
- The units that are returned in the output (optional and advisory).
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!
t_cylindrical¶
t_radial¶
$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:
- d, direct, Direct
- Generate (theta,r) coordinates out (this is the default); incompatible
with Conformal. Theta is in radians, and the radial coordinate is in the
units of distance in the input plane.
- r0, c, conformal, Conformal
- If defined, this floating-point value causes t_radial to generate (theta,
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.
- o, origin, Origin [default (0,0,0)]
- This is the origin of the expansion. Pass in a PDL or an array ref.
- u, unit, Unit [default 'radians']
- This is the angular unit to be used for the azimuth.
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_quadratic¶
$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 t_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.
Note that, because quadratic scaling does not have a strict inverse for
coordinate systems that cross the origin, we cheat slightly and use $x *
abs($x) rather than $x**2. This does the Right thing for pincushion and barrel
distortion, but means that t_quadratic does not behave exactly like t_cubic
with a null cubic strength coefficient.
OPTIONS
- o,origin,Origin
- The origin of the pincushion. (default is the, er, origin).
- l,l0,length,Length,r0
- The fundamental scale of the transformation -- the radius that remains
unchanged. (default=1)
- s,str,strength,Strength
- The relative strength of the pincushion. (default = 0.1)
- honest (default=0)
- Sets whether this is a true quadratic coordinate transform. The more
common form is pincushion or cylindrical distortion, which switches
branches of the square root at the origin (for symmetric expansion).
Setting honest to a non-false value forces true quadratic behavior, which
is not mirror-symmetric about the origin.
- d, dim, dims, Dims
- The number of dimensions to quadratically scale (default is the
dimensionality of your input vectors)
t_cubic¶
$t = t_cubic(<options>);
Cubic scaling - cubic pincushion (n-d; with inverse)
Cubic scaling is a generalization of t_quadratic to a purely cubic expansion.
The formula for the expansion is:
f(a) = ( a' + st * a'^3/L_0^2 ) / (1 + abs(st)) + origin
where a'=(a-origin). That is a simple pincushion expansion/contraction that is
fixed at a distance of L_0 from the origin.
Because there is no quadratic term the result is always invertible with one real
root, and there is no mucking about with complex numbers or multivalued
solutions.
OPTIONS
- o,origin,Origin
- The origin of the pincushion. (default is the, er, origin).
- l,l0,length,Length,r0
- The fundamental scale of the transformation -- the radius that remains
unchanged. (default=1)
- d, dim, dims, Dims
- The number of dimensions to treat (default is the dimensionality of your
input vectors)
t_quartic¶
$t = t_quartic(<options>);
Quartic scaling -- cylindrical pincushion (n-d; with inverse)
Quartic scaling is a generalization of t_quadratic to a quartic expansion. Only
even powers of the input coordinates are retained, and (as with t_quadratic)
sign is preserved, making it an odd function although a true quartic
transformation would be an even function.
You can apply separate distortion along any of the principal axes. If you want
different axes, use t_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.
Note that, because quadratic scaling does not have a strict inverse for
coordinate systems that cross the origin, we cheat slightly and use $x *
abs($x) rather than $x**2. This does the Right thing for pincushion and barrel
distortion, but means that t_quadratic does not behave exactly like t_cubic
with a null cubic strength coefficient.
OPTIONS
- o,origin,Origin
- The origin of the pincushion. (default is the, er, origin).
- l,l0,length,Length,r0
- The fundamental scale of the transformation -- the radius that remains
unchanged. (default=1)
- s,str,strength,Strength
- The relative strength of the pincushion. (default = 0.1)
- honest (default=0)
- Sets whether this is a true quadratic coordinate transform. The more
common form is pincushion or cylindrical distortion, which switches
branches of the square root at the origin (for symmetric expansion).
Setting honest to a non-false value forces true quadratic behavior, which
is not mirror-symmetric about the origin.
- d, dim, dims, Dims
- The number of dimensions to quadratically scale (default is the
dimensionality of your input vectors)
t_spherical¶
$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_unit_sphere.
OPTIONS:
- o, origin, Origin [default (0,0,0)]
- This is the Cartesian origin of the spherical expansion. Pass in a PDL or
an array ref.
- e, euler, Euler [default (0,0,0)]
- This is a 3-vector containing Euler angles to change the angle of the pole
and ordinate. The first two numbers are the (theta, phi) angles of the
pole in a (+Z,+X) spherical expansion, and the last is the angle that the
new prime meridian makes with the meridian of a simply tilted sphere. This
is implemented by composing the output transform with a
PDL::Transform::Linear object.
- u, unit, Unit (default radians)
- This option sets the angular unit to be used. Acceptable values are
"degrees","radians", or reasonable substrings thereof
(e.g. "deg", and "rad", but "d" and
"r" are deprecated). Once genuine unit processing comes online
(a la Math::Units) any angular unit should be OK.
t_projective¶
$t = t_projective(<options>);
Projective transformation
Projective transforms are simple quadratic, quasi-linear transformations. They
are the simplest transformation that can continuously warp an image plane so
that four arbitrarily chosen points exactly map to four other arbitrarily
chosen points. They have the property that straight lines remain straight
after transformation.
You can specify your projective transformation directly in homogeneous
coordinates, or (in 2 dimensions only) as a set of four unique points that are
mapped one to the other by the transformation.
Projective transforms are quasi-linear because they are most easily described as
a linear transformation in homogeneous coordinates (e.g. (x',y',w) where w is
a normalization factor: x = x'/w, etc.). In those coordinates, an N-D
projective transformation is represented as simple multiplication of an
N+1-vector by an N+1 x N+1 matrix, whose lower-right corner value is 1.
If the bottom row of the matrix consists of all zeroes, then the transformation
reduces to a linear affine transformation (as in t_linear).
If the bottom row of the matrix contains nonzero elements, then the transformed
x,y,z,etc. coordinates are related to the original coordinates by a quadratic
polynomial, because the normalization factor 'w' allows a second factor of
x,y, and/or z to enter the equations.
OPTIONS:
- m, mat, matrix, Matrix
- If specified, this is the homogeneous-coordinate matrix to use. It must be
N+1 x N+1, for an N-dimensional transformation.
- p, point, points, Points
- If specified, this is the set of four points that should be mapped one to
the other. The homogeneous-coordinate matrix is calculated from them. You
should feed in a 2x2x4 PDL, where the 0 dimension runs over coordinate,
the 1 dimension runs between input and output, and the 2 dimension runs
over point. For example, specifying
p=> pdl([ [[0,1],[0,1]], [[5,9],[5,8]], [[9,4],[9,3]], [[0,0],[0,0]] ])
maps the origin and the point (0,1) to themselves, the point (5,9) to (5,8),
and the point (9,4) to (9,3).
This is similar to the behavior of fitwarp2d with a quadratic
polynomial.
AUTHOR¶
Copyright 2002, 2003 Craig DeForest. There is no warranty. You are allowed to
redistribute this software 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.