1Transform(3) User Contributed Perl Documentation Transform(3)
2
6 PDL::Transform - Coordinate transforms, image warping, and N-D
7 functions
8
10 use PDL::Transform;
11 my $t = new PDL::Transform::<type>(<opt>)
13 $out = $t->apply($in) # Apply transform to some N-vectors (Transform method)
15 $out = $in->apply($t) # Apply transform to some N-vectors (PDL method)
16 $im1 = $t->map($im); # Transform image coordinates (Transform method)
18 $im1 = $im->map($t); # Transform image coordinates (PDL method)
19 $t2 = $t->compose($t1); # compose two transforms
21 $t2 = $t x $t1; # compose two transforms (by analogy to matrix mult.)
22 $t3 = $t2->inverse(); # invert a transform
24 $t3 = !$t2; # invert a transform (by analogy to logical "not")
25
27 PDL::Transform is a convenient way to represent coordinate
28 transformations and resample images. It embodies functions mapping R^N
29 -> R^M, both with and without inverses. Provision exists for
30 parametrizing functions, and for composing them. You can use this part
31 of the Transform object to keep track of arbitrary functions mapping
32 R^N -> R^M with or without inverses.
33 The simplest way to use a Transform object is to transform vector data
35 between coordinate systems. The "apply" method accepts a PDL whose 0th
36 dimension is coordinate index (all other dimensions are threaded over)
37 and transforms the vectors into the new coordinate system.
38 Transform also includes image resampling, via the "map" method. You
40 define a coordinate transform using a Transform object, then use it to
41 remap an image PDL. The output is a remapped, resampled image.
42 You can define and compose several transformations, then apply them all
44 at once to an image. The image is interpolated only once, when all the
45 composed transformations are applied.
46 In keeping with standard practice, but somewhat counterintuitively, the
48 "map" engine uses the inverse transform to map coordinates FROM the
49 destination dataspace (or image plane) TO the source dataspace; hence
50 PDL::Transform keeps track of both the forward and inverse transform.
51 For terseness and convenience, most of the constructors are exported
53 into the current package with the name "t_<transform>", so the
54 following (for example) are synonyms:
55 $t = new PDL::Transform::Radial(); # Long way
57 $t = t_radial(); # Short way
59 Several math operators are overloaded, so that you can compose and
61 invert functions with expression syntax instead of method syntax (see
62 below).
63
65 Coordinate transformations and mappings are a little counterintuitive
66 at first. Here are some examples of transforms in action:
67 use PDL::Transform;
69 $x = rfits('m51.fits'); # Substitute path if necessary!
70 $ts = t_linear(Scale=>3); # Scaling transform
71 $w = pgwin(xs);
73 $w->imag($x);
74 ## Grow m51 by a factor of 3; origin is at lower left.
76 $y = $ts->map($x,{pix=>1}); # pix option uses direct pixel coord system
77 $w->imag($y);
78 ## Shrink m51 by a factor of 3; origin still at lower left.
80 $c = $ts->unmap($x, {pix=>1});
81 $w->imag($c);
82 ## Grow m51 by a factor of 3; origin is at scientific origin.
84 $d = $ts->map($x,$x->hdr); # FITS hdr template prevents autoscaling
85 $w->imag($d);
86 ## Shrink m51 by a factor of 3; origin is still at sci. origin.
88 $e = $ts->unmap($x,$x->hdr);
89 $w->imag($e);
90 ## A no-op: shrink m51 by a factor of 3, then autoscale back to size
92 $f = $ts->map($x); # No template causes autoscaling of output
93
95 '!'
96 The bang is a unary inversion operator. It binds exactly as tightly
97 as the normal bang operator.
98 'x'
100 By analogy to matrix multiplication, 'x' is the compose operator, so
101 these two expressions are equivalent:
102 $f->inverse()->compose($g)->compose($f) # long way
104 !$f x $g x $f # short way
105 Both of those expressions are equivalent to the mathematical
107 expression f^-1 o g o f, or f^-1(g(f(x))).
108 '**'
110 By analogy to numeric powers, you can apply an operator a positive
111 integer number of times with the ** operator:
112 $f->compose($f)->compose($f) # long way
114 $f**3 # short way
115
117 Transforms are perl hashes. Here's a list of the meaning of each key:
118 func
120 Ref to a subroutine that evaluates the transformed coordinates.
121 It's called with the input coordinate, and the "params" hash. This
122 springboarding is done via explicit ref rather than by subclassing,
123 for convenience both in coding new transforms (just add the
124 appropriate sub to the module) and in adding custom transforms at
125 run-time. Note that, if possible, new "func"s should support inplace
126 operation to save memory when the data are flagged inplace. But
127 "func" should always return its result even when flagged to compute
128 in-place.
129 "func" should treat the 0th dimension of its input as a dimensional
131 index (running 0..N-1 for R^N operation) and thread over all other
132 input dimensions.
133 inv
135 Ref to an inverse method that reverses the transformation. It must
136 accept the same "params" hash that the forward method accepts. This
137 key can be left undefined in cases where there is no inverse.
138 idim, odim
140 Number of useful dimensions for indexing on the input and output
141 sides (ie the order of the 0th dimension of the coordinates to be
142 fed in or that come out). If this is set to 0, then as many are
143 allocated as needed.
144 name
146 A shorthand name for the transformation (convenient for debugging).
147 You should plan on using UNIVERAL::isa to identify classes of
148 transformation, e.g. all linear transformations should be subclasses
149 of PDL::Transform::Linear. That makes it easier to add smarts to,
150 e.g., the compose() method.
151 itype
153 An array containing the name of the quantity that is expected from
154 the input ndarray for the transform, for each dimension. This field
155 is advisory, and can be left blank if there's no obvious quantity
156 associated with the transform. This is analogous to the CTYPEn
157 field used in FITS headers.
158 oname
160 Same as itype, but reporting what quantity is delivered for each
161 dimension.
162 iunit
164 The units expected on input, if a specific unit (e.g. degrees) is
165 expected. This field is advisory, and can be left blank if there's
166 no obvious unit associated with the transform.
167 ounit
169 Same as iunit, but reporting what quantity is delivered for each
170 dimension.
171 params
173 Hash ref containing relevant parameters or anything else the func
174 needs to work right.
175 is_inverse
177 Bit indicating whether the transform has been inverted. That is
178 useful for some stringifications (see the PDL::Transform::Linear
179 stringifier), and may be useful for other things.
180 Transforms should be inplace-aware where possible, to prevent excessive
182 memory usage.
183 If you define a new type of transform, consider generating a new
185 stringify method for it. Just define the sub "stringify" in the
186 subclass package. It should call SUPER::stringify to generate the
187 first line (though the PDL::Transform::Composition bends this rule by
188 tweaking the top-level line), then output (indented) additional lines
189 as necessary to fully describe the transformation.
190
192 Transforms have a mechanism for labeling the units and type of each
193 coordinate, but it is just advisory. A routine to identify and, if
194 necessary, modify units by scaling would be a good idea. Currently, it
195 just assumes that the coordinates are correct for (e.g.) FITS
196 scientific-to-pixel transformations.
197 Composition works OK but should probably be done in a more
199 sophisticated way so that, for example, linear transformations are
200 combined at the matrix level instead of just strung together pixel-to-
201 pixel.
202
204 There are both operators and constructors. The constructors are all
205 exported, all begin with "t_", and all return objects that are
206 subclasses of PDL::Transform.
207 The "apply", "invert", "map", and "unmap" methods are also exported to
209 the "PDL" package: they are both Transform methods and PDL methods.
210
212 apply
213 Signature: (data(); PDL::Transform t)
214 $out = $data->apply($t);
216 $out = $t->apply($data);
217 Apply a transformation to some input coordinates.
219 In the example, $t is a PDL::Transform and $data is a PDL to be
221 interpreted as a collection of N-vectors (with index in the 0th
222 dimension). The output is a similar but transformed PDL.
223 For convenience, this is both a PDL method and a Transform method.
225 invert
227 Signature: (data(); PDL::Transform t)
228 $out = $t->invert($data);
230 $out = $data->invert($t);
231 Apply an inverse transformation to some input coordinates.
233 In the example, $t is a PDL::Transform and $data is an ndarray to be
235 interpreted as a collection of N-vectors (with index in the 0th
236 dimension). The output is a similar ndarray.
237 For convenience this is both a PDL method and a PDL::Transform method.
239 map
241 Signature: (k0(); SV *in; SV *out; SV *map; SV *boundary; SV *method;
242 SV *big; SV *blur; SV *sv_min; SV *flux; SV *bv)
243 match
245 $y = $x->match($c); # Match $c's header and size
246 $y = $x->match([100,200]); # Rescale to 100x200 pixels
247 $y = $x->match([100,200],{rect=>1}); # Rescale and remove rotation/skew.
248 Resample a scientific image to the same coordinate system as another.
250 The example above is syntactic sugar for
252 $y = $x->map(t_identity, $c, ...);
254 it resamples the input PDL with the identity transformation in
256 scientific coordinates, and matches the pixel coordinate system to $c's
257 FITS header.
258 There is one difference between match and map: match makes the
260 "rectify" option to "map" default to 0, not 1. This only affects
261 matching where autoscaling is required (i.e. the array ref example
262 above). By default, that example simply scales $x to the new size and
263 maintains any rotation or skew in its scientific-to-pixel coordinate
264 transform.
265 map
267 $y = $x->map($xform,[<template>],[\%opt]); # Distort $x with transform $xform
268 $y = $x->map(t_identity,[$pdl],[\%opt]); # rescale $x to match $pdl's dims.
269 Resample an image or N-D dataset using a coordinate transform.
271 The data are resampled so that the new pixel indices are proportional
273 to the transformed coordinates rather than the original ones.
274 The operation uses the inverse transform: each output pixel location is
276 inverse-transformed back to a location in the original dataset, and the
277 value is interpolated or sampled appropriately and copied into the
278 output domain. A variety of sampling options are available, trading
279 off speed and mathematical correctness.
280 For convenience, this is both a PDL method and a PDL::Transform method.
282 "map" is FITS-aware: if there is a FITS header in the input data, then
284 the coordinate transform acts on the scientific coordinate system
285 rather than the pixel coordinate system.
286 By default, the output coordinates are separated from pixel coordinates
288 by a single layer of indirection. You can specify the mapping between
289 output transform (scientific) coordinates to pixel coordinates using
290 the "orange" and "irange" options (see below), or by supplying a FITS
291 header in the template.
292 If you don't specify an output transform, then the output is
294 autoscaled: "map" transforms a few vectors in the forward direction to
295 generate a mapping that will put most of the data on the image plane,
296 for most transformations. The calculated mapping gets stuck in the
297 output's FITS header.
298 Autoscaling is especially useful for rescaling images -- if you specify
300 the identity transform and allow autoscaling, you duplicate the
301 functionality of rescale2d, but with more options for interpolation.
302 You can operate in pixel space, and avoid autoscaling of the output, by
304 setting the "nofits" option (see below).
305 The output has the same data type as the input. This is a feature, but
307 it can lead to strange-looking banding behaviors if you use
308 interpolation on an integer input variable.
309 The "template" can be one of:
311 • a PDL
313 The PDL and its header are copied to the output array, which is then
315 populated with data. If the PDL has a FITS header, then the FITS
316 transform is automatically applied so that $t applies to the output
317 scientific coordinates and not to the output pixel coordinates. In
318 this case the NAXIS fields of the FITS header are ignored.
319 • a FITS header stored as a hash ref
321 The FITS NAXIS fields are used to define the output array, and the
323 FITS transformation is applied to the coordinates so that $t applies
324 to the output scientific coordinates.
325 • a list ref
327 This is a list of dimensions for the output array. The code
329 estimates appropriate pixel scaling factors to fill the available
330 space. The scaling factors are placed in the output FITS header.
331 • nothing
333 In this case, the input image size is used as a template, and
335 scaling is done as with the list ref case (above).
336 OPTIONS:
338 The following options are interpreted:
340 b, bound, boundary, Boundary (default = 'truncate')
342 This is the boundary condition to be applied to the input image; it
343 is passed verbatim to range or interpND in the sampling or
344 interpolating stage. Other values are 'forbid','extend', and
345 'periodic'. You can abbreviate this to a single letter. The
346 default 'truncate' causes the entire notional space outside the
347 original image to be filled with 0.
348 pix, Pixel, nf, nofits, NoFITS (default = 0)
350 If you set this to a true value, then FITS headers and
351 interpretation are ignored; the transformation is treated as being
352 in raw pixel coordinates.
353 j, J, just, justify, Justify (default = 0)
355 If you set this to 1, then output pixels are autoscaled to have unit
356 aspect ratio in the output coordinates. If you set it to a non-1
357 value, then it is the aspect ratio between the first dimension and
358 all subsequent dimensions -- or, for a 2-D transformation, the
359 scientific pixel aspect ratio. Values less than 1 shrink the scale
360 in the first dimension compared to the other dimensions; values
361 greater than 1 enlarge it compared to the other dimensions. (This
362 is the same sense as in the PGPLOT interface.)
363 ir, irange, input_range, Input_Range
365 This is a way to modify the autoscaling. It specifies the range of
366 input scientific (not necessarily pixel) coordinates that you want
367 to be mapped to the output image. It can be either a nested array
368 ref or an ndarray. The 0th dim (outside coordinate in the array
369 ref) is dimension index in the data; the 1st dim should have order
370 2. For example, passing in either [[-1,2],[3,4]] or
371 pdl([[-1,2],[3,4]]) limites the map to the quadrilateral in input
372 space defined by the four points (-1,3), (-1,4), (2,4), and (2,3).
373 As with plain autoscaling, the quadrilateral gets sparsely sampled
375 by the autoranger, so pathological transformations can give you
376 strange results.
377 This parameter is overridden by "orange", below.
379 or, orange, output_range, Output_Range
381 This sets the window of output space that is to be sampled onto the
382 output array. It works exactly like "irange", except that it
383 specifies a quadrilateral in output space. Since the output pixel
384 array is itself a quadrilateral, you get pretty much exactly what
385 you asked for.
386 This parameter overrides "irange", if both are specified. It forces
388 rectification of the output (so that scientific axes align with the
389 pixel grid).
390 r, rect, rectify
392 This option defaults TRUE and controls how autoscaling is performed.
393 If it is true or undefined, then autoscaling adjusts so that pixel
394 coordinates in the output image are proportional to individual
395 scientific coordinates. If it is false, then autoscaling
396 automatically applies the inverse of any input FITS transformation
397 *before* autoscaling the pixels. In the special case of linear
398 transformations, this preserves the rectangular shape of the
399 original pixel grid and makes output pixel coordinate proportional
400 to input coordinate.
401 m, method, Method
403 This option controls the interpolation method to be used.
404 Interpolation greatly affects both speed and quality of output. For
405 most cases the option is directly passed to interpND for
406 interpolation. Possible options, in order from fastest to slowest,
407 are:
408 • s, sample (default for ints)
410 Pixel values in the output plane are sampled from the closest
412 data value in the input plane. This is very fast but not very
413 accurate for either magnification or decimation (shrinking). It
414 is the default for templates of integer type.
415 • l, linear (default for floats)
417 Pixel values are linearly interpolated from the closest data
419 value in the input plane. This is reasonably fast but only
420 accurate for magnification. Decimation (shrinking) of the image
421 causes aliasing and loss of photometry as features fall between
422 the samples. It is the default for floating-point templates.
423 • c, cubic
425 Pixel values are interpolated using an N-cubic scheme from a
427 4-pixel N-cube around each coordinate value. As with linear
428 interpolation, this is only accurate for magnification.
429 • f, fft
431 Pixel values are interpolated using the term coefficients of the
433 Fourier transform of the original data. This is the most
434 appropriate technique for some kinds of data, but can yield
435 undesired "ringing" for expansion of normal images. Best suited
436 to studying images with repetitive or wavelike features.
437 • h, hanning
439 Pixel values are filtered through a spatially-variable filter
441 tuned to the computed Jacobian of the transformation, with
442 hanning-window (cosine) pixel rolloff in each dimension. This
443 prevents aliasing in the case where the image is distorted or
444 shrunk, but allows small amounts of aliasing at pixel edges
445 wherever the image is enlarged.
446 • g, gaussian, j, jacobian
448 Pixel values are filtered through a spatially-variable filter
450 tuned to the computed Jacobian of the transformation, with radial
451 Gaussian rolloff. This is the most accurate resampling method,
452 in the sense of introducing the fewest artifacts into a properly
453 sampled data set. This method uses a lookup table to speed up
454 calculation of the Gaussian weighting.
455 • G
457 This method works exactly like 'g' (above), except that the
459 Gaussian values are computed explicitly for every sample -- which
460 is considerably slower.
461 blur, Blur (default = 1.0)
463 This value scales the input-space footprint of each output pixel in
464 the gaussian and hanning methods. It's retained for historical
465 reasons. Larger values yield blurrier images; values significantly
466 smaller than unity cause aliasing. The parameter has slightly
467 different meanings for Gaussian and Hanning interpolation. For
468 Hanning interpolation, numbers smaller than unity control the
469 sharpness of the border at the edge of each pixel (so that blur=>0
470 is equivalent to sampling for non-decimating transforms). For
471 Gaussian interpolation, the blur factor parameter controls the width
472 parameter of the Gaussian around the center of each pixel.
473 sv, SV (default = 1.0)
475 This value lets you set the lower limit of the transformation's
476 singular values in the hanning and gaussian methods, limiting the
477 minimum radius of influence associated with each output pixel.
478 Large numbers yield smoother interpolation in magnified parts of the
479 image but don't affect reduced parts of the image.
480 big, Big (default = 0.5)
482 This is the largest allowable input spot size which may be mapped to
483 a single output pixel by the hanning and gaussian methods, in units
484 of the largest non-thread input dimension. (i.e. the default won't
485 let you reduce the original image to less than 5 pixels across).
486 This places a limit on how long the processing can take for
487 pathological transformations. Smaller numbers keep the code from
488 hanging for a long time; larger numbers provide for photometric
489 accuracy in more pathological cases. Numbers larer than 1.0 are
490 silly, because they allow the entire input array to be compressed
491 into a region smaller than a single pixel.
492 Wherever an output pixel would require averaging over an area that
494 is too big in input space, it instead gets NaN or the bad value.
495 phot, photometry, Photometry
497 This lets you set the style of photometric conversion to be used in
498 the hanning or gaussian methods. You may choose:
499 • 0, s, surf, surface, Surface (default)
501 (this is the default): surface brightness is preserved over the
503 transformation, so features maintain their original intensity.
504 This is what the sampling and interpolation methods do.
505 • 1, f, flux, Flux
507 Total flux is preserved over the transformation, so that the
509 brightness integral over image regions is preserved. Parts of
510 the image that are shrunk wind up brighter; parts that are
511 enlarged end up fainter.
512 VARIABLE FILTERING:
514 The 'hanning' and 'gaussian' methods of interpolation give
516 photometrically accurate resampling of the input data for arbitrary
517 transformations. At each pixel, the code generates a linear
518 approximation to the input transformation, and uses that linearization
519 to estimate the "footprint" of the output pixel in the input space.
520 The output value is a weighted average of the appropriate input spaces.
521 A caveat about these methods is that they assume the transformation is
523 continuous. Transformations that contain discontinuities will give
524 incorrect results near the discontinuity. In particular, the 180th
525 meridian isn't handled well in lat/lon mapping transformations (see
526 PDL::Transform::Cartography) -- pixels along the 180th meridian get the
527 average value of everything along the parallel occupied by the pixel.
528 This flaw is inherent in the assumptions that underly creating a
529 Jacobian matrix. Maybe someone will write code to work around it.
530 Maybe that someone is you.
531 BAD VALUES:
533 "map()" supports bad values in the data array. Bad values in the input
535 array are propagated to the output array. The 'g' and 'h' methods will
536 do some smoothing over bad values: if more than 1/3 of the weighted
537 input-array footprint of a given output pixel is bad, then the output
538 pixel gets marked bad.
539 map does not process bad values. It will set the bad-value flag of all
541 output ndarrays if the flag is set for any of the input ndarrays.
542 unmap
544 Signature: (data(); PDL::Transform a; template(); \%opt)
545 $out_image = $in_image->unmap($t,[<options>],[<template>]);
547 $out_image = $t->unmap($in_image,[<options>],[<template>]);
548 Map an image or N-D dataset using the inverse as a coordinate
550 transform.
551 This convenience function just inverts $t and calls "map" on the
553 inverse; everything works the same otherwise. For convenience, it is
554 both a PDL method and a PDL::Transform method.
555 t_inverse
557 $t2 = t_inverse($t);
558 $t2 = $t->inverse;
559 $t2 = $t ** -1;
560 $t2 = !$t;
561 Return the inverse of a PDL::Transform. This just reverses the
563 func/inv, idim/odim, itype/otype, and iunit/ounit pairs. Note that
564 sometimes you end up with a transform that cannot be applied or mapped,
565 because either the mathematical inverse doesn't exist or the inverse
566 func isn't implemented.
567 You can invert a transform by raising it to a negative power, or by
569 negating it with '!'.
570 The inverse transform remains connected to the main transform because
572 they both point to the original parameters hash. That turns out to be
573 useful.
574 t_compose
576 $f2 = t_compose($f, $g,[...]);
577 $f2 = $f->compose($g[,$h,$i,...]);
578 $f2 = $f x $g x ...;
579 Function composition: f(g(x)), f(g(h(x))), ...
581 You can also compose transforms using the overloaded matrix-
583 multiplication (nee repeat) operator 'x'.
584 This is accomplished by inserting a splicing code ref into the "func"
586 and "inv" slots. It combines multiple compositions into a single list
587 of transforms to be executed in order, fram last to first (in keeping
588 with standard mathematical notation). If one of the functions is
589 itself a composition, it is interpolated into the list rather than left
590 separate. Ultimately, linear transformations may also be combined
591 within the list.
592 No checking is done that the itype/otype and iunit/ounit fields are
594 compatible -- that may happen later, or you can implement it yourself
595 if you like.
596 t_wrap
598 $g1fg = $f->wrap($g);
599 $g1fg = t_wrap($f,$g);
600 Shift a transform into a different space by 'wrapping' it with a
602 second.
603 This is just a convenience function for two "t_compose" calls.
605 "$x->wrap($y)" is the same as "(!$y) x $x x $y": the resulting
606 transform first hits the data with $y, then with $x, then with the
607 inverse of $y.
608 For example, to shift the origin of rotation, do this:
610 $im = rfits('m51.fits');
612 $tf = t_fits($im);
613 $tr = t_linear({rot=>30});
614 $im1 = $tr->map($tr); # Rotate around pixel origin
615 $im2 = $tr->map($tr->wrap($tf)); # Rotate round FITS scientific origin
616 t_identity
618 my $xform = t_identity
619 my $xform = new PDL::Transform;
620 Generic constructor generates the identity transform.
622 This constructor really is trivial -- it is mainly used by the other
624 transform constructors. It takes no parameters and returns the
625 identity transform.
626 t_lookup
628 $f = t_lookup($lookup, {<options>});
629 Transform by lookup into an explicit table.
631 You specify an N+1-D PDL that is interpreted as an N-D lookup table of
633 column vectors (vector index comes last). The last dimension has order
634 equal to the output dimensionality of the transform.
635 For added flexibility in data space, You can specify pre-lookup linear
637 scaling and offset of the data. Of course you can specify the
638 interpolation method to be used. The linear scaling stuff is a little
639 primitive; if you want more, try composing the linear transform with
640 this one.
641 The prescribed values in the lookup table are treated as pixel-
643 centered: that is, if your input array has N elements per row then
644 valid data exist between the locations (-0.5) and (N-0.5) in lookup
645 pixel space, because the pixels (which are numbered from 0 to N-1) are
646 centered on their locations.
647 Lookup is done using interpND, so the boundary conditions and threading
649 behaviour follow from that.
650 The indexed-over dimensions come first in the table, followed by a
652 single dimension containing the column vector to be output for each set
653 of other dimensions -- ie to output 2-vectors from 2 input parameters,
654 each of which can range from 0 to 49, you want an index that has
655 dimension list (50,50,2). For the identity lookup table you could use
656 "cat(xvals(50,50),yvals(50,50))".
657 If you want to output a single value per input vector, you still need
659 that last index threading dimension -- if necessary, use "dummy(-1,1)".
660 The lookup index scaling is: out = lookup[ (scale * data) + offset ].
662 A simplistic table inversion routine is included. This means that you
664 can (for example) use the "map" method with "t_lookup" transformations.
665 But the table inversion is exceedingly slow, and not practical for
666 tables larger than about 100x100. The inversion table is calculated in
667 its entirety the first time it is needed, and then cached until the
668 object is destroyed.
669 Options are listed below; there are several synonyms for each.
671 s, scale, Scale
673 (default 1.0) Specifies the linear amount of scaling to be done
674 before lookup. You can feed in a scalar or an N-vector; other
675 values may cause trouble. If you want to save space in your table,
676 then specify smaller scale numbers.
677 o, offset, Offset
679 (default 0.0) Specifies the linear amount of offset before lookup.
680 This is only a scalar, because it is intended to let you switch to
681 corner-centered coordinates if you want to (just feed in o=-0.25).
682 b, bound, boundary, Boundary
684 Boundary condition to be fed to interpND
685 m, method, Method
687 Interpolation method to be fed to interpND
688 EXAMPLE
690 To scale logarithmically the Y axis of m51, try:
692 $x = float rfits('m51.fits');
694 $lookup = xvals(128,128) -> cat( 10**(yvals(128,128)/50) * 256/10**2.55 );
695 $t = t_lookup($lookup);
696 $y = $t->map($x);
697 To do the same thing but with a smaller lookup table, try:
699 $lookup = 16 * xvals(17,17)->cat(10**(yvals(17,17)/(100/16)) * 16/10**2.55);
701 $t = t_lookup($lookup,{scale=>1/16.0});
702 $y = $t->map($x);
703 (Notice that, although the lookup table coordinates are is divided by
705 16, it is a 17x17 -- so linear interpolation works right to the edge of
706 the original domain.)
707 NOTES
709 Inverses are not yet implemented -- the best way to do it might be by
711 judicious use of map() on the forward transformation.
712 the type/unit fields are ignored.
714 t_linear
716 $f = t_linear({options});
717 Linear (affine) transformations with optional offset
719 t_linear implements simple matrix multiplication with offset, also
721 known as the affine transformations.
722 You specify the linear transformation with pre-offset, a mixing matrix,
724 and a post-offset. That overspecifies the transformation, so you can
725 choose your favorite method to specify the transform you want. The
726 inverse transform is automagically generated, provided that it actually
727 exists (the transform matrix is invertible). Otherwise, the inverse
728 transform just croaks.
729 Extra dimensions in the input vector are ignored, so if you pass a 3xN
731 vector into a 3-D linear transformation, the final dimension is passed
732 through unchanged.
733 The options you can usefully pass in are:
735 s, scale, Scale
737 A scaling scalar (heh), vector, or matrix. If you specify a vector
738 it is treated as a diagonal matrix (for convenience). It gets left-
739 multiplied with the transformation matrix you specify (or the
740 identity), so that if you specify both a scale and a matrix the
741 scaling is done after the rotation or skewing or whatever.
742 r, rot, rota, rotation, Rotation
744 A rotation angle in degrees -- useful for 2-D and 3-D data only. If
745 you pass in a scalar, it specifies a rotation from the 0th axis
746 toward the 1st axis. If you pass in a 3-vector as either a PDL or
747 an array ref (as in "rot=>[3,4,5]"), then it is treated as a set of
748 Euler angles in three dimensions, and a rotation matrix is generated
749 that does the following, in order:
750 • Rotate by rot->(2) degrees from 0th to 1st axis
752 • Rotate by rot->(1) degrees from the 2nd to the 0th axis
754 • Rotate by rot->(0) degrees from the 1st to the 2nd axis
756 The rotation matrix is left-multiplied with the transformation
758 matrix you specify, so that if you specify both rotation and a
759 general matrix the rotation happens after the more general operation
760 -- though that is deprecated.
761 Of course, you can duplicate this functionality -- and get more
763 general -- by generating your own rotation matrix and feeding it in
764 with the "matrix" option.
765 m, matrix, Matrix
767 The transformation matrix. It does not even have to be square, if
768 you want to change the dimensionality of your input. If it is
769 invertible (note: must be square for that), then you automagically
770 get an inverse transform too.
771 pre, preoffset, offset, Offset
773 The vector to be added to the data before they get multiplied by the
774 matrix (equivalent of CRVAL in FITS, if you are converting from
775 scientific to pixel units).
776 post, postoffset, shift, Shift
778 The vector to be added to the data after it gets multiplied by the
779 matrix (equivalent of CRPIX-1 in FITS, if youre converting from
780 scientific to pixel units).
781 d, dim, dims, Dims
783 Most of the time it is obvious how many dimensions you want to deal
784 with: if you supply a matrix, it defines the transformation; if you
785 input offset vectors in the "pre" and "post" options, those define
786 the number of dimensions. But if you only supply scalars, there is
787 no way to tell and the default number of dimensions is 2. This
788 provides a way to do, e.g., 3-D scaling: just set
789 "{s="<scale-factor>, dims=>3}> and you are on your way.
790 NOTES
792 the type/unit fields are currently ignored by t_linear.
794 t_scale
796 $f = t_scale(<scale>)
797 Convenience interface to "t_linear".
799 t_scale produces a transform that scales around the origin by a fixed
801 amount. It acts exactly the same as "t_linear(Scale="\<scale\>)>.
802 t_offset
804 $f = t_offset(<shift>)
805 Convenience interface to "t_linear".
807 t_offset produces a transform that shifts the origin to a new location.
809 It acts exactly the same as "t_linear(Pre="\<shift\>)>.
810 t_rot
812 $f = t_rot(<rotation-in-degrees>)
813 Convenience interface to "t_linear".
815 t_rot produces a rotation transform in 2-D (scalar), 3-D (3-vector), or
817 N-D (matrix). It acts exactly the same as "t_linear(Rot="\<shift\>)>.
818 t_fits
820 $f = t_fits($fits,[option]);
821 FITS pixel-to-scientific transformation with inverse
823 You feed in a hash ref or a PDL with one of those as a header, and you
825 get back a transform that converts 0-originated, pixel-centered
826 coordinates into scientific coordinates via the transformation in the
827 FITS header. For most FITS headers, the transform is reversible, so
828 applying the inverse goes the other way. This is just a convenience
829 subclass of PDL::Transform::Linear, but with unit/type support using
830 the FITS header you supply.
831 For now, this transform is rather limited -- it really ought to accept
833 units differences and stuff like that, but they are just ignored for
834 now. Probably that would require putting units into the whole
835 transform framework.
836 This transform implements the linear transform part of the WCS FITS
838 standard outlined in Greisen & Calabata 2002 (A&A in press; find it at
839 "http://arxiv.org/abs/astro-ph/0207407").
840 As a special case, you can pass in the boolean option "ignore_rgb"
842 (default 0), and if you pass in a 3-D FITS header in which the last
843 dimension has exactly 3 elements, it will be ignored in the output
844 transformation. That turns out to be handy for handling rgb images.
845 t_code
847 $f = t_code(<func>,[<inv>],[options]);
848 Transform implementing arbitrary perl code.
850 This is a way of getting quick-and-dirty new transforms. You pass in
852 anonymous (or otherwise) code refs pointing to subroutines that
853 implement the forward and, optionally, inverse transforms. The
854 subroutines should accept a data PDL followed by a parameter hash ref,
855 and return the transformed data PDL. The parameter hash ref can be set
856 via the options, if you want to.
857 Options that are accepted are:
859 p,params
861 The parameter hash that will be passed back to your code (defaults
862 to the empty hash).
863 n,name
865 The name of the transform (defaults to "code").
866 i, idim (default 2)
868 The number of input dimensions (additional ones should be passed
869 through unchanged)
870 o, odim (default 2)
872 The number of output dimensions
873 itype
875 The type of the input dimensions, in an array ref (optional and
876 advisiory)
877 otype
879 The type of the output dimension, in an array ref (optional and
880 advisory)
881 iunit
883 The units that are expected for the input dimensions (optional and
884 advisory)
885 ounit
887 The units that are returned in the output (optional and advisory).
888 The code variables are executable perl code, either as a code ref or as
890 a string that will be eval'ed to produce code refs. If you pass in a
891 string, it gets eval'ed at call time to get a code ref. If it compiles
892 OK but does not return a code ref, then it gets re-evaluated with "sub
893 { ... }" wrapped around it, to get a code ref.
894 Note that code callbacks like this can be used to do really weird
896 things and generate equally weird results -- caveat scriptor!
897 t_cylindrical
899 "t_cylindrical" is an alias for "t_radial"
900 t_radial
902 $f = t_radial(<options>);
903 Convert Cartesian to radial/cylindrical coordinates. (2-D/3-D; with
905 inverse)
906 Converts 2-D Cartesian to radial (theta,r) coordinates. You can choose
908 direct or conformal conversion. Direct conversion preserves radial
909 distance from the origin; conformal conversion preserves local angles,
910 so that each small-enough part of the image only appears to be scaled
911 and rotated, not stretched. Conformal conversion puts the radius on a
912 logarithmic scale, so that scaling of the original image plane is
913 equivalent to a simple offset of the transformed image plane.
914 If you use three or more dimensions, the higher dimensions are ignored,
916 yielding a conversion from Cartesian to cylindrical coordinates, which
917 is why there are two aliases for the same transform. If you use higher
918 dimensionality than 2, you must manually specify the origin or you will
919 get dimension mismatch errors when you apply the transform.
920 Theta runs clockwise instead of the more usual counterclockwise; that
922 is to preserve the mirror sense of small structures.
923 OPTIONS:
925 d, direct, Direct
927 Generate (theta,r) coordinates out (this is the default);
928 incompatible with Conformal. Theta is in radians, and the radial
929 coordinate is in the units of distance in the input plane.
930 r0, c, conformal, Conformal
932 If defined, this floating-point value causes t_radial to generate
933 (theta, ln(r/r0)) coordinates out. Theta is in radians, and the
934 radial coordinate varies by 1 for each e-folding of the r0-scaled
935 distance from the input origin. The logarithmic scaling is useful
936 for viewing both large and small things at the same time, and for
937 keeping shapes of small things preserved in the image.
938 o, origin, Origin [default (0,0,0)]
940 This is the origin of the expansion. Pass in a PDL or an array ref.
941 u, unit, Unit [default 'radians']
943 This is the angular unit to be used for the azimuth.
944 EXAMPLES
946 These examples do transformations back into the same size image as they
948 started from; by suitable use of the "transform" option to "unmap" you
949 can send them to any size array you like.
950 Examine radial structure in M51: Here, we scale the output to stretch
952 2*pi radians out to the full image width in the horizontal direction,
953 and to stretch 1 radius out to a diameter in the vertical direction.
954 $x = rfits('m51.fits');
956 $ts = t_linear(s => [250/2.0/3.14159, 2]); # Scale to fill orig. image
957 $tu = t_radial(o => [130,130]); # Expand around galactic core
958 $y = $x->map($ts x $tu);
959 Examine radial structure in M51 (conformal): Here, we scale the output
961 to stretch 2*pi radians out to the full image width in the horizontal
962 direction, and scale the vertical direction by the exact same amount to
963 preserve conformality of the operation. Notice that each piece of the
964 image looks "natural" -- only scaled and not stretched.
965 $x = rfits('m51.fits')
967 $ts = t_linear(s=> 250/2.0/3.14159); # Note scalar (heh) scale.
968 $tu = t_radial(o=> [130,130], r0=>5); # 5 pix. radius -> bottom of image
969 $y = $ts->compose($tu)->unmap($x);
970 t_quadratic
972 $t = t_quadratic(<options>);
973 Quadratic scaling -- cylindrical pincushion (n-d; with inverse)
975 Quadratic scaling emulates pincushion in a cylindrical optical system:
977 separate quadratic scaling is applied to each axis. You can apply
978 separate distortion along any of the principal axes. If you want
979 different axes, use "t_wrap" and "t_linear" to rotate them to the
980 correct angle. The scaling options may be scalars or vectors; if they
981 are scalars then the expansion is isotropic.
982 The formula for the expansion is:
984 f(a) = ( <a> + <strength> * a^2/<L_0> ) / (abs(<strength>) + 1)
986 where <strength> is a scaling coefficient and <L_0> is a fundamental
988 length scale. Negative values of <strength> result in a pincushion
989 contraction.
990 Note that, because quadratic scaling does not have a strict inverse for
992 coordinate systems that cross the origin, we cheat slightly and use $x
993 * abs($x) rather than $x**2. This does the Right thing for pincushion
994 and barrel distortion, but means that t_quadratic does not behave
995 exactly like t_cubic with a null cubic strength coefficient.
996 OPTIONS
998 o,origin,Origin
1000 The origin of the pincushion. (default is the, er, origin).
1001 l,l0,length,Length,r0
1003 The fundamental scale of the transformation -- the radius that
1004 remains unchanged. (default=1)
1005 s,str,strength,Strength
1007 The relative strength of the pincushion. (default = 0.1)
1008 honest (default=0)
1010 Sets whether this is a true quadratic coordinate transform. The
1011 more common form is pincushion or cylindrical distortion, which
1012 switches branches of the square root at the origin (for symmetric
1013 expansion). Setting honest to a non-false value forces true
1014 quadratic behavior, which is not mirror-symmetric about the origin.
1015 d, dim, dims, Dims
1017 The number of dimensions to quadratically scale (default is the
1018 dimensionality of your input vectors)
1019 t_cubic
1021 $t = t_cubic(<options>);
1022 Cubic scaling - cubic pincushion (n-d; with inverse)
1024 Cubic scaling is a generalization of t_quadratic to a purely cubic
1026 expansion.
1027 The formula for the expansion is:
1029 f(a) = ( a' + st * a'^3/L_0^2 ) / (1 + abs(st)) + origin
1031 where a'=(a-origin). That is a simple pincushion expansion/contraction
1033 that is fixed at a distance of L_0 from the origin.
1034 Because there is no quadratic term the result is always invertible with
1036 one real root, and there is no mucking about with complex numbers or
1037 multivalued solutions.
1038 OPTIONS
1040 o,origin,Origin
1042 The origin of the pincushion. (default is the, er, origin).
1043 l,l0,length,Length,r0
1045 The fundamental scale of the transformation -- the radius that
1046 remains unchanged. (default=1)
1047 d, dim, dims, Dims
1049 The number of dimensions to treat (default is the dimensionality of
1050 your input vectors)
1051 t_quartic
1053 $t = t_quartic(<options>);
1054 Quartic scaling -- cylindrical pincushion (n-d; with inverse)
1056 Quartic scaling is a generalization of t_quadratic to a quartic
1058 expansion. Only even powers of the input coordinates are retained, and
1059 (as with t_quadratic) sign is preserved, making it an odd function
1060 although a true quartic transformation would be an even function.
1061 You can apply separate distortion along any of the principal axes. If
1063 you want different axes, use "t_wrap" and "t_linear" to rotate them to
1064 the correct angle. The scaling options may be scalars or vectors; if
1065 they are scalars then the expansion is isotropic.
1066 The formula for the expansion is:
1068 f(a) = ( <a> + <strength> * a^2/<L_0> ) / (abs(<strength>) + 1)
1070 where <strength> is a scaling coefficient and <L_0> is a fundamental
1072 length scale. Negative values of <strength> result in a pincushion
1073 contraction.
1074 Note that, because quadratic scaling does not have a strict inverse for
1076 coordinate systems that cross the origin, we cheat slightly and use $x
1077 * abs($x) rather than $x**2. This does the Right thing for pincushion
1078 and barrel distortion, but means that t_quadratic does not behave
1079 exactly like t_cubic with a null cubic strength coefficient.
1080 OPTIONS
1082 o,origin,Origin
1084 The origin of the pincushion. (default is the, er, origin).
1085 l,l0,length,Length,r0
1087 The fundamental scale of the transformation -- the radius that
1088 remains unchanged. (default=1)
1089 s,str,strength,Strength
1091 The relative strength of the pincushion. (default = 0.1)
1092 honest (default=0)
1094 Sets whether this is a true quadratic coordinate transform. The
1095 more common form is pincushion or cylindrical distortion, which
1096 switches branches of the square root at the origin (for symmetric
1097 expansion). Setting honest to a non-false value forces true
1098 quadratic behavior, which is not mirror-symmetric about the origin.
1099 d, dim, dims, Dims
1101 The number of dimensions to quadratically scale (default is the
1102 dimensionality of your input vectors)
1103 t_spherical
1105 $t = t_spherical(<options>);
1106 Convert Cartesian to spherical coordinates. (3-D; with inverse)
1108 Convert 3-D Cartesian to spherical (theta, phi, r) coordinates. Theta
1110 is longitude, centered on 0, and phi is latitude, also centered on 0.
1111 Unless you specify Euler angles, the pole points in the +Z direction
1112 and the prime meridian is in the +X direction. The default is for
1113 theta and phi to be in radians; you can select degrees if you want
1114 them.
1115 Just as the "t_radial" 2-D transform acts like a 3-D cylindrical
1117 transform by ignoring third and higher dimensions, Spherical acts like
1118 a hypercylindrical transform in four (or higher) dimensions. Also as
1119 with "t_radial", you must manually specify the origin if you want to
1120 use more dimensions than 3.
1121 To deal with latitude & longitude on the surface of a sphere (rather
1123 than full 3-D coordinates), see t_unit_sphere.
1124 OPTIONS:
1126 o, origin, Origin [default (0,0,0)]
1128 This is the Cartesian origin of the spherical expansion. Pass in a
1129 PDL or an array ref.
1130 e, euler, Euler [default (0,0,0)]
1132 This is a 3-vector containing Euler angles to change the angle of
1133 the pole and ordinate. The first two numbers are the (theta, phi)
1134 angles of the pole in a (+Z,+X) spherical expansion, and the last is
1135 the angle that the new prime meridian makes with the meridian of a
1136 simply tilted sphere. This is implemented by composing the output
1137 transform with a PDL::Transform::Linear object.
1138 u, unit, Unit (default radians)
1140 This option sets the angular unit to be used. Acceptable values are
1141 "degrees","radians", or reasonable substrings thereof (e.g. "deg",
1142 and "rad", but "d" and "r" are deprecated). Once genuine unit
1143 processing comes online (a la Math::Units) any angular unit should
1144 be OK.
1145 t_projective
1147 $t = t_projective(<options>);
1148 Projective transformation
1150 Projective transforms are simple quadratic, quasi-linear
1152 transformations. They are the simplest transformation that can
1153 continuously warp an image plane so that four arbitrarily chosen points
1154 exactly map to four other arbitrarily chosen points. They have the
1155 property that straight lines remain straight after transformation.
1156 You can specify your projective transformation directly in homogeneous
1158 coordinates, or (in 2 dimensions only) as a set of four unique points
1159 that are mapped one to the other by the transformation.
1160 Projective transforms are quasi-linear because they are most easily
1162 described as a linear transformation in homogeneous coordinates (e.g.
1163 (x',y',w) where w is a normalization factor: x = x'/w, etc.). In those
1164 coordinates, an N-D projective transformation is represented as simple
1165 multiplication of an N+1-vector by an N+1 x N+1 matrix, whose lower-
1166 right corner value is 1.
1167 If the bottom row of the matrix consists of all zeroes, then the
1169 transformation reduces to a linear affine transformation (as in
1170 "t_linear").
1171 If the bottom row of the matrix contains nonzero elements, then the
1173 transformed x,y,z,etc. coordinates are related to the original
1174 coordinates by a quadratic polynomial, because the normalization factor
1175 'w' allows a second factor of x,y, and/or z to enter the equations.
1176 OPTIONS:
1178 m, mat, matrix, Matrix
1180 If specified, this is the homogeneous-coordinate matrix to use. It
1181 must be N+1 x N+1, for an N-dimensional transformation.
1182 p, point, points, Points
1184 If specified, this is the set of four points that should be mapped
1185 one to the other. The homogeneous-coordinate matrix is calculated
1186 from them. You should feed in a 2x2x4 PDL, where the 0 dimension
1187 runs over coordinate, the 1 dimension runs between input and output,
1188 and the 2 dimension runs over point. For example, specifying
1189 p=> pdl([ [[0,1],[0,1]], [[5,9],[5,8]], [[9,4],[9,3]], [[0,0],[0,0]] ])
1191 maps the origin and the point (0,1) to themselves, the point (5,9)
1193 to (5,8), and the point (9,4) to (9,3).
1194 This is similar to the behavior of fitwarp2d with a quadratic
1196 polynomial.
1197
1199 Copyright 2002, 2003 Craig DeForest. There is no warranty. You are
1200 allowed to redistribute this software under certain conditions. For
1201 details, see the file COPYING in the PDL distribution. If this file is
1202 separated from the PDL distribution, the copyright notice should be
1203 included in the file.
1204perl v5.34.0 2021-08-16 Transform(3)