1Transform(3) User Contributed Perl Documentation Transform(3)
2
3
4
6 PDL::Transform - Coordinate transforms, image warping, and N-D
7 functions
8
10 use PDL::Transform;
11
12 my $t = new PDL::Transform::<type>(<opt>)
13
14 $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
17 $im1 = $t->map($im); # Transform image coordinates (Transform method)
18 $im1 = $im->map($t); # Transform image coordinates (PDL method)
19
20 $t2 = $t->compose($t1); # compose two transforms
21 $t2 = $t x $t1; # compose two transforms (by analogy to matrix mult.)
22
23 $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
34 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
39 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
43 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
47 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
52 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
56 $t = new PDL::Transform::Radial(); # Long way
57
58 $t = t_radial(); # Short way
59
60 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
68 use PDL::Transform;
69 $a = rfits('m51.fits'); # Substitute path if necessary!
70 $ts = t_linear(Scale=>3); # Scaling transform
71
72 $w = pgwin(xs);
73 $w->imag($a);
74
75 ## Grow m51 by a factor of 3; origin is at lower left.
76 $b = $ts->map($a,{pix=>1}); # pix option uses direct pixel coord system
77 $w->imag($b);
78
79 ## Shrink m51 by a factor of 3; origin still at lower left.
80 $c = $ts->unmap($a, {pix=>1});
81 $w->imag($c);
82
83 ## Grow m51 by a factor of 3; origin is at scientific origin.
84 $d = $ts->map($a,$a->hdr); # FITS hdr template prevents autoscaling
85 $w->imag($d);
86
87 ## Shrink m51 by a factor of 3; origin is still at sci. origin.
88 $e = $ts->unmap($a,$a->hdr);
89 $w->imag($e);
90
91 ## A no-op: shrink m51 by a factor of 3, then autoscale back to size
92 $f = $ts->map($a); # 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
99 'x'
100 By analogy to matrix multiplication, 'x' is the compose operator, so
101 these two expressions are equivalent:
102
103 $f->inverse()->compose($g)->compose($f) # long way
104 !$f x $g x $f # short way
105
106 Both of those expressions are equivalent to the mathematical
107 expression f^-1 o g o f, or f^-1(g(f(x))).
108
109 '**'
110 By analogy to numeric powers, you can apply an operator a positive
111 integer number of times with the ** operator:
112
113 $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
119 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
130 "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
134 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
139 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
145 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
152 itype
153 An array containing the name of the quantity that is expected from
154 the input piddle 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
159 oname
160 Same as itype, but reporting what quantity is delivered for each
161 dimension.
162
163 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
168 ounit
169 Same as iunit, but reporting what quantity is delivered for each
170 dimension.
171
172 params
173 Hash ref containing relevant parameters or anything else the func
174 needs to work right.
175
176 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
181 Transforms should be inplace-aware where possible, to prevent excessive
182 memory usage.
183
184 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
198 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
208 The apply, invert, map, and unmap methods are also exported to the
209 "PDL" package: they are both Transform methods and PDL methods.
210
212 apply
213 Signature: (data(); PDL::Transform t)
214
215 $out = $data->apply($t);
216 $out = $t->apply($data);
217
218 Apply a transformation to some input coordinates.
219
220 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
224 For convenience, this is both a PDL method and a Transform method.
225
226 invert
227 Signature: (data(); PDL::Transform t)
228
229 $out = $t->invert($data);
230 $out = $data->invert($t);
231
232 Apply an inverse transformation to some input coordinates.
233
234 In the example, $t is a PDL::Transform and $data is a piddle to be
235 interpreted as a collection of N-vectors (with index in the 0th
236 dimension). The output is a similar piddle.
237
238 For convenience this is both a PDL method and a PDL::Transform method.
239
240 map
241 Signature: (k0(); SV *in; SV *out; SV *map; SV *boundary; SV *method;
242 SV *big; SV *blur; SV *sv_min; SV *flux)
243
244 PDL::match
245 $b = $a->match($c);
246
247 Resample a scientific image to the same coordinate system as another.
248
249 The example above is syntactic sugar for
250
251 $b = $a->map(t_identity, $c, ...);
252
253 it resamples the input PDL with the identity transformation in
254 scientific coordinates, and matches the pixel coordinate system to $c's
255 FITS header.
256
257 map
258 $b = $a->map($xform,[<template>],[\%opt]); # Distort $a with tranform $xform
259 $b = $a->map(t_identity,[$pdl],[\%opt]); # rescale $a to match $pdl's dims.
260
261 Resample an image or N-D dataset using a coordinate transform.
262
263 The data are resampled so that the new pixel indices are proportional
264 to the transformed coordinates rather than the original ones.
265
266 The operation uses the inverse transform: each output pixel location is
267 inverse-transformed back to a location in the original dataset, and the
268 value is interpolated or sampled appropriately and copied into the
269 output domain. A variety of sampling options are available, trading
270 off speed and mathematical correctness.
271
272 For convenience, this is both a PDL method and a PDL::Transform method.
273
274 "map" is FITS-aware: if there is a FITS header in the input data, then
275 the coordinate transform acts on the scientific coordinate system
276 rather than the pixel coordinate system.
277
278 By default, the output coordinates are separated from pixel coordinates
279 by a single layer of indirection. You can specify the mapping between
280 output transform (scientific) coordinates to pixel coordinates using
281 the "orange" and "irange" options (see below), or by supplying a FITS
282 header in the template.
283
284 If you don't specify an output transform, then the output is
285 autoscaled: "map" transforms a few vectors in the forward direction to
286 generate a mapping that will put most of the data on the image plane,
287 for most transformations. The calculated mapping gets stuck in the
288 output's FITS header.
289
290 Autoscaling is especially useful for rescaling images -- if you specify
291 the identity transform and allow autoscaling, you duplicate the
292 functionality of rescale2d, but with more options for interpolation.
293
294 You can operate in pixel space, and avoid autoscaling of the output, by
295 setting the "nofits" option (see below).
296
297 The output has the same data type as the input. This is a feature, but
298 it can lead to strange-looking banding behaviors if you use
299 interpolation on an integer input variable.
300
301 The "template" can be one of:
302
303 · a PDL
304
305 The PDL and its header are copied to the output array, which is then
306 populated with data. If the PDL has a FITS header, then the FITS
307 transform is automatically applied so that $t applies to the output
308 scientific coordinates and not to the output pixel coordinates. In
309 this case the NAXIS fields of the FITS header are ignored.
310
311 · a FITS header stored as a hash ref
312
313 The FITS NAXIS fields are used to define the output array, and the
314 FITS transformation is applied to the coordinates so that $t applies
315 to the output scientific coordinates.
316
317 · a list ref
318
319 This is a list of dimensions for the output array. The code
320 estimates appropriate pixel scaling factors to fill the available
321 space. The scaling factors are placed in the output FITS header.
322
323 · nothing
324
325 In this case, the input image size is used as a template, and
326 scaling is done as with the list ref case (above).
327
328 OPTIONS:
329
330 The following options are interpreted:
331
332 b, bound, boundary, Boundary (default = 'truncate')
333 This is the boundary condition to be applied to the input image; it
334 is passed verbatim to range or interpND in the sampling or
335 interpolating stage. Other values are 'forbid','extend', and
336 'periodic'. You can abbreviate this to a single letter. The
337 default 'truncate' causes the entire notional space outside the
338 original image to be filled with 0.
339
340 p, pix, Pixel, nf, nofits, NoFITS (default = 0)
341 If you set this to a true value, then FITS headers and
342 interpretation are ignored; the transformation is treated as being
343 in raw pixel coordinates.
344
345 j, J, just, justify, Justify (default = 0)
346 If you set this to 1, then output pixels are autoscaled to have unit
347 aspect ratio in the output coordinates. If you set it to a non-1
348 value, then it is the aspect ratio between the first dimension and
349 all subsequent dimensions -- or, for a 2-D transformation, the
350 scientific pixel aspect ratio. Values less than 1 shrink the scale
351 in the first dimension compared to the other dimensions; values
352 greater than 1 enlarge it compared to the other dimensions. (This
353 is the same sense as in the PGPLOTinterface.)
354
355 ir, irange, input_range, Input_Range
356 This is a way to modify the autoscaling. It specifies the range of
357 input scientific (not necessarily pixel) coordinates that you want
358 to be mapped to the output image. It can be either a nested array
359 ref or a piddle. The 0th dim (outside coordinate in the array ref)
360 is dimension index in the data; the 1st dim should have order 2.
361 For example, passing in either [[-1,2],[3,4]] or pdl([[-1,2],[3,4]])
362 limites the map to the quadrilateral in input space defined by the
363 four points (-1,3), (-1,4), (2,4), and (2,3).
364
365 As with plain autoscaling, the quadrilateral gets sparsely sampled
366 by the autoranger, so pathological transformations can give you
367 strange results.
368
369 This parameter is overridden by "orange", below.
370
371 or, orange, output_range, Output_Range
372 This sets the window of output space that is to be sampled onto the
373 output array. It works exactly like "irange", except that it
374 specifies a quadrilateral in output space. Since the output pixel
375 array is itself a quadrilateral, you get pretty much exactly what
376 you asked for.
377
378 This parameter overrides "irange", if both are specified.
379
380 m, method, Method
381 This option controls the interpolation method to be used.
382 Interpolation greatly affects both speed and quality of output. For
383 most cases the option is directly passed to interpND for
384 interpolation. Possible options, in order from fastest to slowest,
385 are:
386
387 · s, sample (default for ints)
388
389 Pixel values in the output plane are sampled from the closest
390 data value in the input plane. This is very fast but not very
391 accurate for either magnification or decimation (shrinking). It
392 is the default for templates of integer type.
393
394 · l, linear (default for floats)
395
396 Pixel values are linearly interpolated from the closest data
397 value in the input plane. This is reasonably fast but only
398 accurate for magnification. Decimation (shrinking) of the image
399 causes aliasing and loss of photometry as features fall between
400 the samples. It is the default for floating-point templates.
401
402 · c, cubic
403
404 Pixel values are interpolated using an N-cubic scheme from a
405 4-pixel N-cube around each coordinate value. As with linear
406 interpolation, this is only accurate for magnification.
407
408 · f, fft
409
410 Pixel values are interpolated using the term coefficients of the
411 Fourier transform of the original data. This is the most
412 appropriate technique for some kinds of data, but can yield
413 undesired "ringing" for expansion of normal images. Best suited
414 to studying images with repetitive or wavelike features.
415
416 · h, hanning
417
418 Pixel values are filtered through a spatially-variable filter
419 tuned to the computed Jacobian of the transformation, with
420 hanning-window (cosine) pixel rolloff in each dimension. This
421 prevents aliasing in the case where the image is distorted or
422 shrunk, but allows small amounts of aliasing at pixel edges
423 wherever the image is enlarged.
424
425 · g, gaussian, j, jacobian
426
427 Pixel values are filtered through a spatially-variable filter
428 tuned to the computed Jacobian of the transformation, with radial
429 Gaussian rolloff. This is the most accurate resampling method,
430 in the sense of introducing the fewest artifacts into a properly
431 sampled data set.
432
433 blur, Blur (default = 1.0)
434 This value scales the input-space footprint of each output pixel in
435 the gaussian and hanning methods. It's retained for historical
436 reasons. Larger values yield blurrier images; values significantly
437 smaller than unity cause aliasing.
438
439 sv, SV (default = 1.0)
440 This value lets you set the lower limit of the transformation's
441 singular values in the hanning and gaussian methods, limiting the
442 minimum radius of influence associated with each output pixel.
443 Large numbers yield smoother interpolation in magnified parts of the
444 image but don't affect reduced parts of the image.
445
446 big, Big (default = 0.2)
447 This is the largest allowable input spot size which may be mapped to
448 a single output pixel by the hanning and gaussian methods, in units
449 of the largest non-thread input dimension. (i.e. the default won't
450 let you reduce the original image to less than 5 pixels across).
451 This places a limit on how long the processing can take for
452 pathological transformations. Smaller numbers keep the code from
453 hanging for a long time; larger numbers provide for photometric
454 accuracy in more pathological cases. Numbers larer than 1.0 are
455 silly, because they allow the entire input array to be compressed
456 into a region smaller than a single pixel.
457
458 Wherever an output pixel would require averaging over an area that
459 is too big in input space, it instead gets NaN or the equivalent
460 (bad values are not yet supported).
461
462 p, phot, photometry, Photometry
463 This lets you set the style of photometric conversion to be used in
464 the hanning or gaussian methods. You may choose:
465
466 · 0, s, surf, surface, Surface (default)
467
468 (this is the default): surface brightness is preserved over the
469 transformation, so features maintain their original intensity.
470 This is what the sampling and interpolation methods do.
471
472 · 1, f, flux, Flux
473
474 Total flux is preserved over the transformation, so that the
475 brightness integral over image regions is preserved. Parts of
476 the image that are shrunk wind up brighter; parts that are
477 enlarged end up fainter.
478
479 VARIABLE FILTERING:
480
481 The 'hanning' and 'gaussian' methods of interpolation give
482 photometrically accurate resampling of the input data for arbitrary
483 transformations. At each pixel, the code generates a linear
484 approximation to the input transformation, and uses that linearization
485 to estimate the "footprint" of the output pixel in the input space.
486 The output value is a weighted average of the appropriate input spaces.
487
488 A caveat about these methods is that they assume the transformation is
489 continuous. Transformations that contain discontinuities will give
490 incorrect results near the discontinuity. In particular, the 180th
491 meridian isn't handled well in lat/lon mapping transformations (see
492 PDL::Transform::Cartography) -- pixels along the 180th meridian get the
493 average value of everything along the parallel occupied by the pixel.
494 This flaw is inherent in the assumptions that underly creating a
495 Jacobian matrix. Maybe someone will write code to work around it.
496 Maybe that someone is you.
497
498 map does not process bad values. It will set the bad-value flag of all
499 output piddles if the flag is set for any of the input piddles.
500
501 PDL::Transform::unmap
502 Signature: (data(); PDL::Transform a; template(); \%opt)
503
504 $out_image = $in_image->unmap($t,[<options>],[<template>]);
505 $out_image = $t->unmap($in_image,[<options>],[<template>]);
506
507 Map an image or N-D dataset using the inverse as a coordinate
508 transform.
509
510 This convenience function just inverts $t and calls map on the inverse;
511 everything works the same otherwise. For convenience, it is both a PDL
512 method and a PDL::Transform method.
513
514 t_inverse
515 $t2 = t_inverse($t);
516 $t2 = $t->inverse;
517 $t2 = $t ** -1;
518 $t2 = !$t;
519
520 Return the inverse of a PDL::Transform. This just reverses the
521 func/inv, idim/odim, itype/otype, and iunit/ounit pairs. Note that
522 sometimes you end up with a transform that cannot be applied or mapped,
523 because either the mathematical inverse doesn't exist or the inverse
524 func isn't implemented.
525
526 You can invert a transform by raising it to a negative power, or by
527 negating it with '!'.
528
529 The inverse transform remains connected to the main transform because
530 they both point to the original parameters hash. That turns out to be
531 useful.
532
533 t_compose
534 $f2 = t_compose($f, $g,[...]);
535 $f2 = $f->compose($g[,$h,$i,...]);
536 $f2 = $f x $g x ...;
537
538 Function composition: f(g(x)), f(g(h(x))), ...
539
540 You can also compose transforms using the overloaded matrix-
541 multiplication (nee repeat) operator 'x'.
542
543 This is accomplished by inserting a splicing code ref into the "func"
544 and "inv" slots. It combines multiple compositions into a single list
545 of transforms to be executed in order, fram last to first (in keeping
546 with standard mathematical notation). If one of the functions is
547 itself a composition, it is interpolated into the list rather than left
548 separate. Ultimately, linear transformations may also be combined
549 within the list.
550
551 No checking is done that the itype/otype and iunit/ounit fields are
552 compatible -- that may happen later, or you can implement it yourself
553 if you like.
554
555 t_wrap
556 $g1fg = $f->wrap($g);
557 $g1fg = t_wrap($f,$g);
558
559 Shift a transform into a different space by 'wrapping' it with a
560 second.
561
562 This is just a convenience function for two compose calls.
563 "$a-"wrap($b)> is the same as "(!$b) x $a x $b": the resulting
564 transform first hits the data with $b, then with $a, then with the
565 inverse of $b.
566
567 For example, to shift the origin of rotation, do this:
568
569 $im = rfits('m51.fits');
570 $tf = t_fits($im);
571 $tr = t_linear({rot=>30});
572 $im1 = $tr->map($tr); # Rotate around pixel origin
573 $im2 = $tr->map($tr->wrap($tf)); # Rotate round FITS scientific origin
574
575 t_identity
576 my $xform = t_identity
577 my $xform = new PDL::Transform;
578
579 Generic constructor generates the identity transform.
580
581 This constructor really is trivial -- it is mainly used by the other
582 transform constructors. It takes no parameters and returns the
583 identity transform.
584
585 t_lookup
586 $f = t_lookup($lookup, {<options>});
587
588 Transform by lookup into an explicit table.
589
590 You specify an N+1-D PDL that is interpreted as an N-D lookup table of
591 column vectors (vector index comes last). The last dimension has order
592 equal to the output dimensionality of the transform.
593
594 For added flexibility in data space, You can specify pre-lookup linear
595 scaling and offset of the data. Of course you can specify the
596 interpolation method to be used. The linear scaling stuff is a little
597 primitive; if you want more, try composing the linear transform with
598 this one.
599
600 The prescribed values in the lookup table are treated as pixel-
601 centered: that is, if your input array has N elements per row then
602 valid data exist between the locations (-0.5) and (N-0.5) in lookup
603 pixel space, because the pixels (which are numbered from 0 to N-1) are
604 centered on their locations.
605
606 Lookup is done using interpND, so the boundary conditions and threading
607 behaviour follow from that.
608
609 The indexed-over dimensions come first in the table, followed by a
610 single dimension containing the column vector to be output for each set
611 of other dimensions -- ie to output 2-vectors from 2 input parameters,
612 each of which can range from 0 to 49, you want an index that has
613 dimension list (50,50,2). For the identity lookup table you could use
614 "cat(xvals(50,50),yvals(50,50))".
615
616 If you want to output a single value per input vector, you still need
617 that last index threading dimension -- if necessary, use "dummy(-1,1)".
618
619 The lookup index scaling is: out = lookup[ (scale * data) + offset ].
620
621 A simplistic table inversion routine is included. This means that you
622 can (for example) use the "map" method with "t_lookup" transformations.
623 But the table inversion is exceedingly slow, and not practical for
624 tables larger than about 100x100. The inversion table is calculated in
625 its entirety the first time it is needed, and then cached until the
626 object is destroyed.
627
628 Options are listed below; there are several synonyms for each.
629
630 s, scale, Scale
631 (default 1.0) Specifies the linear amount of scaling to be done
632 before lookup. You can feed in a scalar or an N-vector; other
633 values may cause trouble. If you want to save space in your table,
634 then specify smaller scale numbers.
635
636 o, offset, Offset
637 (default 0.0) Specifies the linear amount of offset before lookup.
638 This is only a scalar, because it is intended to let you switch to
639 corner-centered coordinates if you want to (just feed in o=-0.25).
640
641 b, bound, boundary, Boundary
642 Boundary condition to be fed to interpND
643
644 m, method, Method
645 Interpolation method to be fed to interpND
646
647 EXAMPLE
648
649 To scale logarithmically the Y axis of m51, try:
650
651 $a = float rfits('m51.fits');
652 $lookup = xvals(128,128) -> cat( 10**(yvals(128,128)/50) * 256/10**2.55 );
653 $t = t_lookup($lookup);
654 $b = $t->map($a);
655
656 To do the same thing but with a smaller lookup table, try:
657
658 $lookup = 16 * xvals(17,17)->cat(10**(yvals(17,17)/(100/16)) * 16/10**2.55);
659 $t = t_lookup($lookup,{scale=>1/16.0});
660 $b = $t->map($a);
661
662 (Notice that, although the lookup table coordinates are is divided by
663 16, it is a 17x17 -- so linear interpolation works right to the edge of
664 the original domain.)
665
666 NOTES
667
668 Inverses are not yet implemented -- the best way to do it might be by
669 judicious use of map() on the forward transformation.
670
671 the type/unit fields are ignored.
672
673 t_linear
674 $f = t_linear({options});
675
676 Linear (affine) transformations with optional offset
677
678 t_linear implements simple matrix multiplication with offset, also
679 known as the affine transformations.
680
681 You specify the linear transformation with pre-offset, a mixing matrix,
682 and a post-offset. That overspecifies the transformation, so you can
683 choose your favorite method to specify the transform you want. The
684 inverse transform is automagically generated, provided that it actually
685 exists (the transform matrix is invertible). Otherwise, the inverse
686 transform just croaks.
687
688 Extra dimensions in the input vector are ignored, so if you pass a 3xN
689 vector into a 3-D linear transformation, the final dimension is passed
690 through unchanged.
691
692 The options you can usefully pass in are:
693
694 s, scale, Scale
695 A scaling scalar (heh), vector, or matrix. If you specify a vector
696 it is treated as a diagonal matrix (for convenience). It gets left-
697 multiplied with the transformation matrix you specify (or the
698 identity), so that if you specify both a scale and a matrix the
699 scaling is done after the rotation or skewing or whatever.
700
701 r, rot, rota, rotation, Rotation
702 A rotation angle in degrees -- useful for 2-D and 3-D data only. If
703 you pass in a scalar, it specifies a rotation from the 0th axis
704 toward the 1st axis. If you pass in a 3-vector as either a PDL or
705 an array ref (as in "rot=>[3,4,5]"), then it is treated as a set of
706 Euler angles in three dimensions, and a rotation matrix is generated
707 that does the following, in order:
708
709 · Rotate by rot->(2) degrees from 0th to 1st axis
710
711 · Rotate by rot->(1) degrees from the 2nd to the 0th axis
712
713 · Rotate by rot->(0) degrees from the 1st to the 2nd axis
714
715 The rotation matrix is left-multiplied with the transformation
716 matrix you specify, so that if you specify both rotation and a
717 general matrix the rotation happens after the more general operation
718 -- though that is deprecated.
719
720 Of course, you can duplicate this functionality -- and get more
721 general -- by generating your own rotation matrix and feeding it in
722 with the "matrix" option.
723
724 m, matrix, Matrix
725 The transformation matrix. It does not even have to be square, if
726 you want to change the dimensionality of your input. If it is
727 invertible (note: must be square for that), then you automagically
728 get an inverse transform too.
729
730 pre, preoffset, offset, Offset
731 The vector to be added to the data before they get multiplied by the
732 matrix (equivalent of CRVAL in FITS, if you are converting from
733 scientific to pixel units).
734
735 post, postoffset, shift, Shift
736 The vector to be added to the data after it gets multiplied by the
737 matrix (equivalent of CRPIX-1 in FITS, if youre converting from
738 scientific to pixel units).
739
740 d, dim, dims, Dims
741 Most of the time it is obvious how many dimensions you want to deal
742 with: if you supply a matrix, it defines the transformation; if you
743 input offset vectors in the "pre" and "post" options, those define
744 the number of dimensions. But if you only supply scalars, there is
745 no way to tell and the default number of dimensions is 2. This
746 provides a way to do, e.g., 3-D scaling: just set
747 "{s="<scale-factor>, dims=>3}> and you are on your way.
748
749 NOTES
750
751 the type/unit fields are currently ignored by t_linear.
752
753 t_scale
754 $f = t_scale(<scale>)
755
756 Convenience interface to t_linear.
757
758 t_scale produces a tranform that scales around the origin by a fixed
759 amount. It acts exactly the same as "t_linear(Scale="\<scale\>)>.
760
761 t_offset
762 $f = t_offset(<shift>)
763
764 Convenience interface to t_linear.
765
766 t_offset produces a transform that shifts the origin to a new location.
767 It acts exactly the same as "t_linear(Pre="\<shift\>)>.
768
769 t_rot
770 $f = t_rot(<rotation-in-degrees>)
771
772 Convenience interface to t_linear.
773
774 t_rot produces a rotation transform in 2-D (scalar), 3-D (3-vector), or
775 N-D (matrix). It acts exactly the same as "t_linear(Rot="\<shift\>)>.
776
777 t_fits
778 $f = t_fits($fits,[option]);
779
780 FITS pixel-to-scientific transformation with inverse
781
782 You feed in a hash ref or a PDL with one of those as a header, and you
783 get back a transform that converts 0-originated, pixel-centered
784 coordinates into scientific coordinates via the transformation in the
785 FITS header. For most FITS headers, the transform is reversible, so
786 applying the inverse goes the other way. This is just a convenience
787 subclass of PDL::Transform::Linear, but with unit/type support using
788 the FITS header you supply.
789
790 For now, this transform is rather limited -- it really ought to accept
791 units differences and stuff like that, but they are just ignored for
792 now. Probably that would require putting units into the whole
793 transform framework.
794
795 This transform implements the linear transform part of the WCS FITS
796 standard outlined in Greisen & Calabata 2002 (A&A in press; find it at
797 "http://arxiv.org/abs/astro-ph/0207407").
798
799 As a special case, you can pass in the boolean option "ignore_rgb"
800 (default 0), and if you pass in a 3-D FITS header in which the last
801 dimension has exactly 3 elements, it will be ignored in the output
802 transformation. That turns out to be handy for handling rgb images.
803
804 t_code
805 $f = t_code(<func>,[<inv>],[options]);
806
807 Transform implementing arbitrary perl code.
808
809 This is a way of getting quick-and-dirty new transforms. You pass in
810 anonymous (or otherwise) code refs pointing to subroutines that
811 implement the forward and, optionally, inverse transforms. The
812 subroutines should accept a data PDL followed by a parameter hash ref,
813 and return the transformed data PDL. The parameter hash ref can be set
814 via the options, if you want to.
815
816 Options that are accepted are:
817
818 p,params
819 The parameter hash that will be passed back to your code (defaults
820 to the empty hash).
821
822 n,name
823 The name of the transform (defaults to "code").
824
825 i, idim (default 2)
826 The number of input dimensions (additional ones should be passed
827 through unchanged)
828
829 o, odim (default 2)
830 The number of output dimensions
831
832 itype
833 The type of the input dimensions, in an array ref (optional and
834 advisiory)
835
836 otype
837 The type of the output dimension, in an array ref (optional and
838 advisory)
839
840 iunit
841 The units that are expected for the input dimensions (optional and
842 advisory)
843
844 ounit
845 The units that are returned in the output (optional and advisory).
846
847 The code variables are executable perl code, either as a code ref or as
848 a string that will be eval'ed to produce code refs. If you pass in a
849 string, it gets eval'ed at call time to get a code ref. If it compiles
850 OK but does not return a code ref, then it gets re-evaluated with "sub
851 { ... }" wrapped around it, to get a code ref.
852
853 Note that code callbacks like this can be used to do really weird
854 things and generate equally weird results -- caveat scriptor!
855
856 t_cylindrical
857 t_radial
858 $f = t_radial(<options>);
859
860 Convert Cartesian to radial/cylindrical coordinates. (2-D/3-D; with
861 inverse)
862
863 Converts 2-D Cartesian to radial (theta,r) coordinates. You can choose
864 direct or conformal conversion. Direct conversion preserves radial
865 distance from the origin; conformal conversion preserves local angles,
866 so that each small-enough part of the image only appears to be scaled
867 and rotated, not stretched. Conformal conversion puts the radius on a
868 logarithmic scale, so that scaling of the original image plane is
869 equivalent to a simple offset of the transformed image plane.
870
871 If you use three or more dimensions, the higher dimensions are ignored,
872 yielding a conversion from Cartesian to cylindrical coordinates, which
873 is why there are two aliases for the same transform. If you use higher
874 dimensionality than 2, you must manually specify the origin or you will
875 get dimension mismatch errors when you apply the transform.
876
877 Theta runs clockwise instead of the more usual counterclockwise; that
878 is to preserve the mirror sense of small structures.
879
880 OPTIONS:
881
882 d, direct, Direct
883 Generate (theta,r) coordinates out (this is the default);
884 incompatible with Conformal. Theta is in radians, and the radial
885 coordinate is in the units of distance in the input plane.
886
887 r0, c, conformal, Conformal
888 If defined, this floating-point value causes t_radial to generate
889 (theta, ln(r/r0)) coordinates out. Theta is in radians, and the
890 radial coordinate varies by 1 for each e-folding of the r0-scaled
891 distance from the input origin. The logarithmic scaling is useful
892 for viewing both large and small things at the same time, and for
893 keeping shapes of small things preserved in the image.
894
895 o, origin, Origin [default (0,0,0)]
896 This is the origin of the expansion. Pass in a PDL or an array ref.
897
898 u, unit, Unit [default 'radians']
899 This is the angular unit to be used for the azimuth.
900
901 EXAMPLES
902
903 These examples do transformations back into the same size image as they
904 started from; by suitable use of the "transform" option to unmap you
905 can send them to any size array you like.
906
907 Examine radial structure in M51: Here, we scale the output to stretch
908 2*pi radians out to the full image width in the horizontal direction,
909 and to stretch 1 radius out to a diameter in the vertical direction.
910
911 $a = rfits('m51.fits');
912 $ts = t_linear(s => [250/2.0/3.14159, 2]); # Scale to fill orig. image
913 $tu = t_radial(o => [130,130]); # Expand around galactic core
914 $b = $a->map($ts x $tu);
915
916 Examine radial structure in M51 (conformal): Here, we scale the output
917 to stretch 2*pi radians out to the full image width in the horizontal
918 direction, and scale the vertical direction by the exact same amount to
919 preserve conformality of the operation. Notice that each piece of the
920 image looks "natural" -- only scaled and not stretched.
921
922 $a = rfits('m51.fits')
923 $ts = t_linear(s=> 250/2.0/3.14159); # Note scalar (heh) scale.
924 $tu = t_radial(o=> [130,130], r0=>5); # 5 pix. radius -> bottom of image
925 $b = $ts->compose($tu)->unmap($a);
926
927 t_quadratic
928 $t = t_quadratic(<options>);
929
930 Quadratic scaling -- cylindrical pincushion (n-d; with inverse)
931
932 Quadratic scaling emulates pincushion in a cylindrical optical system:
933 separate quadratic scaling is applied to each axis. You can apply
934 separate distortion along any of the principal axes. If you want
935 different axes, use wrap and t_linear to rotate them to the correct
936 angle. The scaling options may be scalars or vectors; if they are
937 scalars then the expansion is isotropic.
938
939 The formula for the expansion is:
940
941 f(a) = ( <a> + <strength> * a^2/<L_0> ) / (abs(<strength>) + 1)
942
943 where <strength> is a scaling coefficient and <L_0> is a fundamental
944 length scale. Negative values of <strength> result in a pincushion
945 contraction.
946
947 OPTIONS
948
949 o,origin,Origin
950 The origin of the pincushion. (default is the, er, origin).
951
952 l,l0,length,Length,r0
953 The fundamental scale of the transformation -- the radius that
954 remains unchanged. (default=1)
955
956 s,str,strength,Strength
957 The relative strength of the pincushion. (default = 0.1)
958
959 d, dim, dims, Dims
960 The number of dimensions to quadratically scale (default is the
961 dimensionality of your input vectors)
962
963 t_spherical
964 $t = t_spherical(<options>);
965
966 Convert Cartesian to spherical coordinates. (3-D; with inverse)
967
968 Convert 3-D Cartesian to spherical (theta, phi, r) coordinates. Theta
969 is longitude, centered on 0, and phi is latitude, also centered on 0.
970 Unless you specify Euler angles, the pole points in the +Z direction
971 and the prime meridian is in the +X direction. The default is for
972 theta and phi to be in radians; you can select degrees if you want
973 them.
974
975 Just as the t_radial 2-D transform acts like a 3-D cylindrical
976 transform by ignoring third and higher dimensions, Spherical acts like
977 a hypercylindrical transform in four (or higher) dimensions. Also as
978 with t_radial, you must manually specify the origin if you want to use
979 more dimensions than 3.
980
981 To deal with latitude & longitude on the surface of a sphere (rather
982 than full 3-D coordinates), see t_unitsphere.
983
984 OPTIONS:
985
986 o, origin, Origin [default (0,0,0)]
987 This is the Cartesian origin of the spherical expansion. Pass in a
988 PDL or an array ref.
989
990 e, euler, Euler [default (0,0,0)]
991 This is a 3-vector containing Euler angles to change the angle of
992 the pole and ordinate. The first two numbers are the (theta, phi)
993 angles of the pole in a (+Z,+X) spherical expansion, and the last is
994 the angle that the new prime meridian makes with the meridian of a
995 simply tilted sphere. This is implemented by composing the output
996 transform with a PDL::Transform::Linear object.
997
998 u, unit, Unit (default radians)
999 This option sets the angular unit to be used. Acceptable values are
1000 "degrees","radians", or reasonable substrings thereof (e.g. "deg",
1001 and "rad", but "d" and "r" are deprecated). Once genuine unit
1002 processing comes online (a la Math::Units) any angular unit should
1003 be OK.
1004
1005 t_projective
1006 $t = t_projective(<options>);
1007
1008 Projective transformation
1009
1010 Projective transforms are simple quadratic, quasi-linear
1011 transformations. They are the simplest transformation that can
1012 continuously warp an image plane so that four arbitrarily chosen points
1013 exactly map to four other arbitrarily chosen points. They have the
1014 property that straight lines remain straight after transformation.
1015
1016 You can specify your projective transformation directly in homogeneous
1017 coordinates, or (in 2 dimensions only) as a set of four unique points
1018 that are mapped one to the other by the transformation.
1019
1020 Projective transforms are quasi-linear because they are most easily
1021 described as a linear transformation in homogeneous coordinates (e.g.
1022 (x',y',w) where w is a normalization factor: x = x'/w, etc.). In those
1023 coordinates, an N-D projective transformation is represented as simple
1024 multiplication of an N+1-vector by an N+1 x N+1 matrix, whose lower-
1025 right corner value is 1.
1026
1027 If the bottom row of the matrix consists of all zeroes, then the
1028 transformation reduces to a linear affine transformation (as in
1029 t_linear).
1030
1031 If the bottom row of the matrix contains nonzero elements, then the
1032 transformed x,y,z,etc. coordinates are related to the original
1033 coordinates by a quadratic polynomial, because the normalization factor
1034 'w' allows a second factor of x,y, and/or z to enter the equations.
1035
1036 OPTIONS:
1037
1038 m, mat, matrix, Matrix
1039 If specified, this is the homogeneous-coordinate matrix to use. It
1040 must be N+1 x N+1, for an N-dimensional transformation.
1041
1042 p, point, points, Points
1043 If specified, this is the set of four points that should be mapped
1044 one to the other. The homogeneous-coordinate matrix is calculated
1045 from them. You should feed in a 2x2x4 PDL, where the 0 dimension
1046 runs over coordinate, the 1 dimension runs between input and output,
1047 and the 2 dimension runs over point. For example, specifying
1048
1049 p=> pdl([ [[0,1],[0,1]], [[5,9],[5,8]], [[9,4],[9,3]], [[0,0],[0,0]] ])
1050
1051 maps the origin and the point (0,1) to themselves, the point (5,9)
1052 to (5,8), and the point (9,4) to (9,3).
1053
1054 This is similar to the behavior of fitwarp2d with a quadratic
1055 polynomial.
1056
1058 Copyright 2002, 2003 Craig DeForest. There is no warranty. You are
1059 allowed to redistribute this software under certain conditions. For
1060 details, see the file COPYING in the PDL distribution. If this file is
1061 separated from the PDL distribution, the copyright notice should be
1062 included in the file.
1063
1064
1065
1066perl v5.12.3 2011-03-31 Transform(3)