File Coverage

blib/lib/PDL/Transform.pm

Criterion	Covered	Total	%
statement	9	9	100.0
branch			n/a
condition			n/a
subroutine	3	3	100.0
pod			n/a
total	12	12	100.0

line	stmt	sub	time	code
1				#
2				# GENERATED WITH PDL::PP from lib/PDL/Transform.pd! Don't modify!
3				#
4				package PDL::Transform;
5
6				our @EXPORT_OK = qw(apply invert map map unmap t_inverse t_compose t_wrap t_identity t_lookup t_linear t_scale t_offset t_rot t_fits t_code t_cylindrical t_radial t_quadratic t_cubic t_quadratic t_spherical t_projective );
7				our %EXPORT_TAGS = (Func=>\@EXPORT_OK);
8
9	1	1	884	use PDL::Core;
	1		1
	1		28
10	1	1	5	use PDL::Exporter;
	1		2
	1		4
11	1	1	4	use DynaLoader;
	1		1
	1		941
12
13
14
15				our @ISA = ( 'PDL::Exporter','DynaLoader' );
16				push @PDL::Core::PP, __PACKAGE__;
17				bootstrap PDL::Transform ;
18
19
20
21
22
23
24
25
26				#line 2 "lib/PDL/Transform.pd"
27
28				=head1 NAME
29
30				PDL::Transform - Coordinate transforms, image warping, and N-D functions
31
32				=head1 SYNOPSIS
33
34				use PDL::Transform;
35
36				my $t = PDL::Transform::->new()
37
38				$out = $t->apply($in) # Apply transform to some N-vectors (Transform method)
39				$out = $in->apply($t) # Apply transform to some N-vectors (PDL method)
40
41				$im1 = $t->map($im); # Transform image coordinates (Transform method)
42				$im1 = $im->map($t); # Transform image coordinates (PDL method)
43
44				$t2 = $t->compose($t1); # compose two transforms
45				$t2 = $t x $t1; # compose two transforms (by analogy to matrix mult.)
46
47				$t3 = $t2->inverse(); # invert a transform
48				$t3 = !$t2; # invert a transform (by analogy to logical "not")
49
50				=head1 DESCRIPTION
51
52				PDL::Transform is a convenient way to represent coordinate
53				transformations and resample images. It embodies functions mapping
54				R^N -> R^M, both with and without inverses. Provision exists for
55				parametrizing functions, and for composing them. You can use this
56				part of the Transform object to keep track of arbitrary functions
57				mapping R^N -> R^M with or without inverses.
58
59				The simplest way to use a Transform object is to transform vector
60				data between coordinate systems. The L method
61				accepts a PDL whose 0th dimension is coordinate index (all other
62				dimensions are broadcasted over) and transforms the vectors into the new
63				coordinate system.
64
65				Transform also includes image resampling, via the L method.
66				You define a coordinate transform using a Transform object, then use
67				it to remap an image PDL. The output is a remapped, resampled image.
68
69				You can define and compose several transformations, then apply them
70				all at once to an image. The image is interpolated only once, when
71				all the composed transformations are applied.
72
73				In keeping with standard practice, but somewhat counterintuitively,
74				the L engine uses the inverse transform to map coordinates
75				FROM the destination dataspace (or image plane) TO the source dataspace;
76				hence PDL::Transform keeps track of both the forward and inverse transform.
77
78				For terseness and convenience, most of the constructors are exported
79				into the current package with the name C<< t_ >>, so the following
80				(for example) are synonyms:
81
82				$t = PDL::Transform::Radial->new; # Long way
83
84				$t = t_radial(); # Short way
85
86				Several math operators are overloaded, so that you can compose and
87				invert functions with expression syntax instead of method syntax (see below).
88
89				=head1 EXAMPLE
90
91				Coordinate transformations and mappings are a little counterintuitive
92				at first. Here are some examples of transforms in action:
93
94				use PDL::Transform;
95				use PDL::Graphics::Simple;
96				$x = rfits('m51.fits'); # Substitute path if necessary!
97				$ts = t_linear(Scale=>3); # Scaling transform
98
99				$w = pgswin(xs);
100				$w->plot(with=>'image', $x);
101
102				## Grow m51 by a factor of 3; origin is at lower left.
103				$y = $ts->map($x,{pix=>1}); # pix option uses direct pixel coord system
104				$w->plot(with=>'image', $y);
105
106				## Shrink m51 by a factor of 3; origin still at lower left.
107				$c = $ts->unmap($x, {pix=>1});
108				$w->plot(with=>'image', $c);
109
110				## Grow m51 by a factor of 3; origin is at scientific origin.
111				$d = $ts->map($x,$x->hdr); # FITS hdr template prevents autoscaling
112				$w->plot(with=>'image', $d);
113
114				## Shrink m51 by a factor of 3; origin is still at sci. origin.
115				$e = $ts->unmap($x,$x->hdr);
116				$w->plot(with=>'image', $e);
117
118				## A no-op: shrink m51 by a factor of 3, then autoscale back to size
119				$f = $ts->map($x); # No template causes autoscaling of output
120
121				=head1 OPERATOR OVERLOADS
122
123				=over 3
124
125				=item '!'
126
127				The bang is a unary inversion operator. It binds exactly as
128				tightly as the normal bang operator.
129
130				=item 'x'
131
132				By analogy to matrix multiplication, 'x' is the compose operator, so these
133				two expressions are equivalent:
134
135				$f->inverse()->compose($g)->compose($f) # long way
136				!$f x $g x $f # short way
137
138				Both of those expressions are equivalent to the mathematical expression
139				f^-1 o g o f, or f^-1(g(f(x))).
140
141				=item '**'
142
143				By analogy to numeric powers, you can apply an operator a positive
144				integer number of times with the ** operator:
145
146				$f->compose($f)->compose($f) # long way
147				$f**3 # short way
148
149				=back
150
151				=head1 INTERNALS
152
153				Transforms are perl hashes. Here's a list of the meaning of each key:
154
155				=over 3
156
157				=item func
158
159				Ref to a subroutine that evaluates the transformed coordinates. It's
160				called with the input coordinate, and the "params" hash. This
161				springboarding is done via explicit ref rather than by subclassing,
162				for convenience both in coding new transforms (just add the
163				appropriate sub to the module) and in adding custom transforms at
164				run-time. Note that, if possible, new Cs should support
165				L operation to save memory when the data are flagged
166				inplace. But C should always return its result even when
167				flagged to compute in-place.
168
169				C should treat the 0th dimension of its input as a dimensional
170				index (running 0..N-1 for R^N operation) and broadcast over all other input
171				dimensions.
172
173				=item inv
174
175				Ref to an inverse method that reverses the transformation. It must
176				accept the same "params" hash that the forward method accepts. This
177				key can be left undefined in cases where there is no inverse.
178
179				=item idim, odim
180
181				Number of useful dimensions for indexing on the input and output sides
182				(ie the order of the 0th dimension of the coordinates to be fed in or
183				that come out). If this is set to 0, then as many are allocated as needed.
184
185				=item name
186
187				A shorthand name for the transformation (convenient for debugging).
188				You should plan on using UNIVERAL::isa to identify classes of
189				transformation, e.g. all linear transformations should be subclasses
190				of PDL::Transform::Linear. That makes it easier to add smarts to,
191				e.g., the compose() method.
192
193				=item itype
194
195				An array containing the name of the quantity that is expected from the
196				input ndarray for the transform, for each dimension. This field is advisory,
197				and can be left blank if there's no obvious quantity associated with
198				the transform. This is analogous to the CTYPEn field used in FITS headers.
199
200				=item oname
201
202				Same as itype, but reporting what quantity is delivered for each
203				dimension.
204
205				=item iunit
206
207				The units expected on input, if a specific unit (e.g. degrees) is expected.
208				This field is advisory, and can be left blank if there's no obvious
209				unit associated with the transform.
210
211				=item ounit
212
213				Same as iunit, but reporting what quantity is delivered for each dimension.
214
215				=item params
216
217				Hash ref containing relevant parameters or anything else the func needs to
218				work right.
219
220				=item is_inverse
221
222				Bit indicating whether the transform has been inverted. That is useful
223				for some stringifications (see the PDL::Transform::Linear
224				stringifier), and may be useful for other things.
225
226				=back
227
228				Transforms should be inplace-aware where possible, to prevent excessive
229				memory usage.
230
231				If you define a new type of transform, consider generating a new stringify
232				method for it. Just define the sub "stringify" in the subclass package.
233				It should call SUPER::stringify to generate the first line (though
234				the PDL::Transform::Composition bends this rule by tweaking the
235				top-level line), then output (indented) additional lines as necessary to
236				fully describe the transformation.
237
238				=head1 NOTES
239
240				Transforms have a mechanism for labeling the units and type of each
241				coordinate, but it is just advisory. A routine to identify and, if
242				necessary, modify units by scaling would be a good idea. Currently,
243				it just assumes that the coordinates are correct for (e.g.) FITS
244				scientific-to-pixel transformations.
245
246				Composition works OK but should probably be done in a more
247				sophisticated way so that, for example, linear transformations are
248				combined at the matrix level instead of just strung together
249				pixel-to-pixel.
250
251				=head1 MODULE INTERFACE
252
253				There are both operators and constructors. The constructors are all
254				exported, all begin with "t_", and all return objects that are subclasses
255				of PDL::Transform.
256
257				The L, L, L,
258				and L methods are also exported to the C package: they
259				are both Transform methods and PDL methods.
260
261				=cut
262
263				use strict;
264				use warnings;
265				#line 266 "lib/PDL/Transform.pm"
266
267
268				=head1 FUNCTIONS
269
270				=cut
271
272
273
274
275
276				#line 315 "lib/PDL/Transform.pd"
277
278				=head2 apply
279
280				=for sig
281
282				Signature: (data(); PDL::Transform t)
283
284				=for usage
285
286				$out = $data->apply($t);
287				$out = $t->apply($data);
288
289				=for ref
290
291				Apply a transformation to some input coordinates.
292
293				In the example, C<$t> is a PDL::Transform and C<$data> is a PDL to
294				be interpreted as a collection of N-vectors (with index in the 0th
295				dimension). The output is a similar but transformed PDL.
296
297				For convenience, this is both a PDL method and a Transform method.
298
299				=cut
300
301				use Carp;
302
303				*PDL::apply = \&apply;
304				sub apply {
305				my ($me, $from) = UNIVERSAL::isa($_[0],'PDL') ? @_[1,0] : @_;
306				croak "apply requires both a PDL and a PDL::Transform.\n"
307				unless UNIVERSAL::isa($me,'PDL::Transform') && UNIVERSAL::isa($from,'PDL');
308				croak "Applying a PDL::Transform with no func! Oops.\n"
309				unless(defined($me->{func}) and ref($me->{func}) eq 'CODE');
310				my $result = $me->{func}->($from,$me->{params});
311				$result->is_inplace(0); # clear inplace flag, just in case.
312				if (defined $from->gethdr && $from->hdr->{NAXIS}) {
313				my $nd = $me->{odim} \|\| $me->{idim} \|\| 2;
314				for my $d (1..$nd) {
315				$result->hdr->{"CTYPE$d"} = ( (!defined($me->{otype}) ? "" :
316				ref($me->{otype}) ne 'ARRAY' ? $me->{otype} :
317				$me->{otype}[$d-1])
318				\|\| $from->hdr->{"CTYPE$d"}
319				\|\| "");
320				$result->hdr->{"CUNIT$d"} = ( (!defined($me->{ounit}) ? "" :
321				ref($me->{ounit}) ne 'ARRAY' ? $me->{ounit} :
322				$me->{ounit}[$d-1])
323				\|\| $from->hdr->{"CUNIT$d"}
324				\|\| $from->hdr->{"CTYPE$d"}
325				\|\| "");
326				}
327				}
328				return $result;
329				}
330
331				#line 373 "lib/PDL/Transform.pd"
332
333				=head2 invert
334
335				=for sig
336
337				Signature: (data(); PDL::Transform t)
338
339				=for usage
340
341				$out = $t->invert($data);
342				$out = $data->invert($t);
343
344				=for ref
345
346				Apply an inverse transformation to some input coordinates.
347
348				In the example, C<$t> is a PDL::Transform and C<$data> is an ndarray to
349				be interpreted as a collection of N-vectors (with index in the 0th
350				dimension). The output is a similar ndarray.
351
352				For convenience this is both a PDL method and a PDL::Transform method.
353
354				=cut
355
356				*PDL::invert = \&invert;
357				sub invert {
358				my ($me, $from) = UNIVERSAL::isa($_[0],'PDL') ? @_[1,0] : @_;
359				croak "invert requires a PDL and a PDL::Transform (did you want 'inverse' instead?)\n"
360				unless UNIVERSAL::isa($me,'PDL::Transform') && UNIVERSAL::isa($from,'PDL');
361				croak "Inverting a PDL::Transform with no inverse! Oops.\n"
362				unless(defined($me->{inv}) and ref($me->{inv}) eq 'CODE');
363				my $result = $me->{inv}->($from, $me->{params});
364				$result->is_inplace(0); # make sure inplace flag is clear.
365				return $result;
366				}
367				#line 368 "lib/PDL/Transform.pm"
368
369
370				=head2 map
371
372				=for sig
373
374				Signature: (k0(); pdl in; pdl out; pdl map; SV boundary; SV *method;
375				long big; double blur; double sv_min; char flux; SV *bv)
376				Types: (sbyte byte short ushort long ulong indx ulonglong longlong
377				float double ldouble)
378
379				=head2 match
380
381				=for usage
382
383				$y = $x->match($c); # Match $c's header and size
384				$y = $x->match([100,200]); # Rescale to 100x200 pixels
385				$y = $x->match([100,200],{rect=>1}); # Rescale and remove rotation/skew.
386
387				=for ref
388
389				Resample a scientific image to the same coordinate system as another.
390
391				The example above is syntactic sugar for
392
393				$y = $x->map(t_identity, $c, ...);
394
395				it resamples the input PDL with the identity transformation in
396				scientific coordinates, and matches the pixel coordinate system to
397				$c's FITS header.
398
399				There is one difference between match and map: match makes the
400				C option to C default to 0, not 1. This only affects
401				matching where autoscaling is required (i.e. the array ref example
402				above). By default, that example simply scales $x to the new size and
403				maintains any rotation or skew in its scientific-to-pixel coordinate
404				transform.
405
406				=head2 map
407
408				=for usage
409
410				$y = $x->map($xform,[
411				$y = $x->map(t_identity,[$pdl],[\%opt]); # rescale $x to match $pdl's dims.
412
413				=for ref
414
415				Resample an image or N-D dataset using a coordinate transform.
416
417				The data are resampled so that the new pixel indices are proportional
418				to the transformed coordinates rather than the original ones.
419
420				The operation uses the inverse transform: each output pixel location
421				is inverse-transformed back to a location in the original dataset, and
422				the value is interpolated or sampled appropriately and copied into the
423				output domain. A variety of sampling options are available, trading
424				off speed and mathematical correctness.
425
426				For convenience, this is both a PDL method and a PDL::Transform method.
427
428				C is FITS-aware: if there is a FITS header in the input data,
429				then the coordinate transform acts on the scientific coordinate system
430				rather than the pixel coordinate system.
431
432				By default, the output coordinates are separated from pixel coordinates
433				by a single layer of indirection. You can specify the mapping between
434				output transform (scientific) coordinates to pixel coordinates using
435				the C and C options (see below), or by supplying a
436				FITS header in the template.
437
438				If you don't specify an output transform, then the output is
439				autoscaled: C transforms a few vectors in the forward direction
440				to generate a mapping that will put most of the data on the image
441				plane, for most transformations. The calculated mapping gets stuck in the
442				output's FITS header.
443
444				Autoscaling is especially useful for rescaling images -- if you specify
445				the identity transform and allow autoscaling, you duplicate the
446				functionality of L, but with more
447				options for interpolation.
448
449				You can operate in pixel space, and avoid autoscaling of the output,
450				by setting the C option (see below).
451
452				The output has the same data type as the input. This is a feature,
453				but it can lead to strange-looking banding behaviors if you use
454				interpolation on an integer input variable.
455
456				The C
457
458				=over 3
459
460				=item * a PDL
461
462				The PDL and its header are copied to the output array, which is then
463				populated with data. If the PDL has a FITS header, then the FITS
464				transform is automatically applied so that $t applies to the output
465				scientific coordinates and not to the output pixel coordinates. In
466				this case the NAXIS fields of the FITS header are ignored.
467
468				=item * a FITS header stored as a hash ref
469
470				The FITS NAXIS fields are used to define the output array, and the
471				FITS transformation is applied to the coordinates so that $t applies to the
472				output scientific coordinates.
473
474				=item * a list ref
475
476				This is a list of dimensions for the output array. The code estimates
477				appropriate pixel scaling factors to fill the available space. The
478				scaling factors are placed in the output FITS header.
479
480				=item * nothing
481
482				In this case, the input image size is used as a template, and scaling
483				is done as with the list ref case (above).
484
485				=back
486
487				OPTIONS:
488
489				The following options are interpreted:
490
491				=over 3
492
493				=item b, bound, boundary, Boundary (default = 'truncate')
494
495				This is the boundary condition to be applied to the input image; it is
496				passed verbatim to L or
497				L in the sampling or interpolating
498				stage. Other values are 'forbid','extend', and 'periodic'. You can
499				abbreviate this to a single letter. The default 'truncate' causes the
500				entire notional space outside the original image to be filled with 0.
501
502				=item pix, Pixel, nf, nofits, NoFITS (default = 0)
503
504				If you set this to a true value, then FITS headers and interpretation
505				are ignored; the transformation is treated as being in raw pixel coordinates.
506
507				=item j, J, just, justify, Justify (default = 0)
508
509				If you set this to 1, then output pixels are autoscaled to have unit
510				aspect ratio in the output coordinates. If you set it to a non-1
511				value, then it is the aspect ratio between the first dimension and all
512				subsequent dimensions -- or, for a 2-D transformation, the scientific
513				pixel aspect ratio. Values less than 1 shrink the scale in the first
514				dimension compared to the other dimensions; values greater than 1
515				enlarge it compared to the other dimensions.
516
517				=item ir, irange, input_range, Input_Range
518
519				This is a way to modify the autoscaling. It specifies the range of
520				input scientific (not necessarily pixel) coordinates that you want to be
521				mapped to the output image. It can be either a nested array ref or
522				an ndarray. The 0th dim (outside coordinate in the array ref) is
523				dimension index in the data; the 1st dim should have order 2.
524				For example, passing in either [[-1,2],[3,4]] or pdl([[-1,2],[3,4]])
525				limites the map to the quadrilateral in input space defined by the
526				four points (-1,3), (-1,4), (2,4), and (2,3).
527
528				As with plain autoscaling, the quadrilateral gets sparsely sampled by
529				the autoranger, so pathological transformations can give you strange
530				results.
531
532				This parameter is overridden by C, below.
533
534				=item or, orange, output_range, Output_Range
535
536				This sets the window of output space that is to be sampled onto the
537				output array. It works exactly like C, except that it specifies
538				a quadrilateral in output space. Since the output pixel array is itself
539				a quadrilateral, you get pretty much exactly what you asked for.
540
541				This parameter overrides C, if both are specified. It forces
542				rectification of the output (so that scientific axes align with the pixel
543				grid).
544
545				=item r, rect, rectify
546
547				This option defaults TRUE and controls how autoscaling is performed. If
548				it is true or undefined, then autoscaling adjusts so that pixel coordinates
549				in the output image are proportional to individual scientific coordinates.
550				If it is false, then autoscaling automatically applies the inverse of any
551				input FITS transformation before autoscaling the pixels. In the special
552				case of linear transformations, this preserves the rectangular shape of the
553				original pixel grid and makes output pixel coordinate proportional to input
554				coordinate.
555
556				=item m, method, Method
557
558				This option controls the interpolation method to be used.
559				Interpolation greatly affects both speed and quality of output. For
560				most cases the option is directly passed to
561				L for interpolation. Possible
562				options, in order from fastest to slowest, are:
563
564				=over 3
565
566				=item * s, sample (default for ints)
567
568				Pixel values in the output plane are sampled from the closest data value
569				in the input plane. This is very fast but not very accurate for either
570				magnification or decimation (shrinking). It is the default for templates
571				of integer type.
572
573				=item * l, linear (default for floats)
574
575				Pixel values are linearly interpolated from the closest data value in the
576				input plane. This is reasonably fast but only accurate for magnification.
577				Decimation (shrinking) of the image causes aliasing and loss of photometry
578				as features fall between the samples. It is the default for floating-point
579				templates.
580
581				=item * c, cubic
582
583				Pixel values are interpolated using an N-cubic scheme from a 4-pixel
584				N-cube around each coordinate value. As with linear interpolation,
585				this is only accurate for magnification.
586
587				=item * f, fft
588
589				Pixel values are interpolated using the term coefficients of the
590				Fourier transform of the original data. This is the most appropriate
591				technique for some kinds of data, but can yield undesired "ringing" for
592				expansion of normal images. Best suited to studying images with
593				repetitive or wavelike features.
594
595				=item * h, hanning
596
597				Pixel values are filtered through a spatially-variable filter tuned to
598				the computed Jacobian of the transformation, with hanning-window
599				(cosine) pixel rolloff in each dimension. This prevents aliasing in the
600				case where the image is distorted or shrunk, but allows small amounts
601				of aliasing at pixel edges wherever the image is enlarged.
602
603				=item * g, gaussian, j, jacobian
604
605				Pixel values are filtered through a spatially-variable filter tuned to
606				the computed Jacobian of the transformation, with radial Gaussian
607				rolloff. This is the most accurate resampling method, in the sense of
608				introducing the fewest artifacts into a properly sampled data set.
609				This method uses a lookup table to speed up calculation of the Gaussian
610				weighting.
611
612				=item * G
613
614				This method works exactly like 'g' (above), except that the Gaussian
615				values are computed explicitly for every sample -- which is considerably
616				slower.
617
618				=back
619
620				=item blur, Blur (default = 1.0)
621
622				This value scales the input-space footprint of each output pixel in
623				the gaussian and hanning methods. It's retained for historical
624				reasons. Larger values yield blurrier images; values significantly
625				smaller than unity cause aliasing. The parameter has slightly
626				different meanings for Gaussian and Hanning interpolation. For
627				Hanning interpolation, numbers smaller than unity control the
628				sharpness of the border at the edge of each pixel (so that blur=>0 is
629				equivalent to sampling for non-decimating transforms). For
630				Gaussian interpolation, the blur factor parameter controls the
631				width parameter of the Gaussian around the center of each pixel.
632
633				=item sv, SV (default = 1.0)
634
635				This value lets you set the lower limit of the transformation's
636				singular values in the hanning and gaussian methods, limiting the
637				minimum radius of influence associated with each output pixel. Large
638				numbers yield smoother interpolation in magnified parts of the image
639				but don't affect reduced parts of the image.
640
641				=item big, Big (default = 0.5)
642
643				This is the largest allowable input spot size which may be mapped to a
644				single output pixel by the hanning and gaussian methods, in units of
645				the largest non-broadcast input dimension. (i.e. the default won't let
646				you reduce the original image to less than 5 pixels across). This places
647				a limit on how long the processing can take for pathological transformations.
648				Smaller numbers keep the code from hanging for a long time; larger numbers
649				provide for photometric accuracy in more pathological cases. Numbers larger
650				than 1.0 are silly, because they allow the entire input array to be compressed
651				into a region smaller than a single pixel.
652
653				Wherever an output pixel would require averaging over an area that is too
654				big in input space, it instead gets NaN or the bad value.
655
656				=item phot, photometry, Photometry
657
658				This lets you set the style of photometric conversion to be used in the
659				hanning or gaussian methods. You may choose:
660
661				=over 3
662
663				=item * 0, s, surf, surface, Surface (default)
664
665				(this is the default): surface brightness is preserved over the transformation,
666				so features maintain their original intensity. This is what the sampling
667				and interpolation methods do.
668
669				=item * 1, f, flux, Flux
670
671				Total flux is preserved over the transformation, so that the brightness
672				integral over image regions is preserved. Parts of the image that are
673				shrunk wind up brighter; parts that are enlarged end up fainter.
674
675				=back
676
677				=back
678
679				VARIABLE FILTERING:
680
681				The 'hanning' and 'gaussian' methods of interpolation give
682				photometrically accurate resampling of the input data for arbitrary
683				transformations. At each pixel, the code generates a linear
684				approximation to the input transformation, and uses that linearization
685				to estimate the "footprint" of the output pixel in the input space.
686				The output value is a weighted average of the appropriate input spaces.
687
688				A caveat about these methods is that they assume the transformation is
689				continuous. Transformations that contain discontinuities will give
690				incorrect results near the discontinuity. In particular, the 180th
691				meridian isn't handled well in lat/lon mapping transformations (see
692				L) -- pixels along the 180th meridian get
693				the average value of everything along the parallel occupied by the
694				pixel. This flaw is inherent in the assumptions that underly creating
695				a Jacobian matrix. Maybe someone will write code to work around it.
696				Maybe that someone is you.
697
698				BAD VALUES:
699
700				C supports bad values in the data array. Bad values in the input
701				array are propagated to the output array. The 'g' and 'h' methods will
702				do some smoothing over bad values: if more than 1/3 of the weighted
703				input-array footprint of a given output pixel is bad, then the output
704				pixel gets marked bad.
705
706				=pod
707
708				Broadcasts over its inputs.
709
710				=for bad
711
712				C does not process bad values.
713				It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
714
715				=cut
716
717
718
719
720
721				#line 1537 "lib/PDL/Transform.pd"
722
723				sub PDL::match {
724				# Set default for rectification to 0 for simple matching...
725				push @_, {} if ref($_[-1]) ne 'HASH';
726				my @k = grep(m/^r(e(c(t)?)?)?/,sort keys %{$_[-1]});
727				#line 1542 "lib/PDL/Transform.pd"
728				unless(@k) {
729				$_[-1]->{rectify} = 0;
730				}
731				t_identity()->map(@_);
732				}
733
734				*PDL::map = \↦
735				sub map {
736				my ($me, $in, $tmp, $opt) = @_;
737				($in, $me) = ($me, $in)
738				if UNIVERSAL::isa($me,'PDL') && UNIVERSAL::isa($in,'PDL::Transform');
739				barf ("PDL::Transform::map: source is not defined or is not a PDL\n")
740				unless(defined $in and UNIVERSAL::isa($in,'PDL'));
741				barf ("PDL::Transform::map: source is empty\n")
742				unless($in->nelem);
743				# Check for options-but-no-template case
744				($opt, $tmp) = ($tmp, undef)
745				if ref $tmp eq 'HASH' && !defined($opt)
746				&& !defined($tmp->{NAXIS}); # FITS headers all have NAXIS.
747				croak("PDL::Transform::map: Option 'p' was ambiguous and has been removed. You probably want 'pix' or 'phot'.") if exists($opt->{'p'});
748				$tmp = [$in->dims] unless defined $tmp; # no //= because of Devel::Cover bug
749
750				# Generate an appropriate output ndarray for values to go in
751				my ($out, @odims, $ohdr);
752				if(UNIVERSAL::isa($tmp,'PDL')) {
753				@odims = $tmp->dims;
754				if (defined(my $x = $tmp->gethdr)) { $ohdr = {%$x} }
755				} elsif(ref $tmp eq 'HASH') {
756				# (must be a fits header -- or would be filtered above)
757				@odims = map $tmp->{"NAXIS$_"}, 1..$tmp->{NAXIS};
758				$ohdr = {%$tmp}; # copy FITS header into output
759				} elsif(ref $tmp eq 'ARRAY') {
760				@odims = @$tmp;
761				} else {
762				barf("map: confused about dimensions of the output array...\n");
763				}
764
765				if(scalar(@odims) < scalar($in->dims)) {
766				my @idims = $in->dims;
767				push(@odims, splice(@idims,scalar(@odims)));
768				}
769
770				$out = PDL->new_from_specification($in->type,@odims);
771				$out->sethdr($ohdr) if defined($ohdr);
772
773				# set badflag on output all the time if possible, to account for boundary violations
774				$out->badflag(1);
775
776				##############################
777				## Figure out the dimensionality of the
778				## transform itself (extra dimensions come along for the ride)
779				my $nd = $me->{odim} \|\| $me->{idim} \|\| 2;
780				my @dd = my @sizes = $out->dims;
781				splice @dd,$nd; # Cut out dimensions after the end
782
783				# Check that there are elements in the output fields...
784				barf "map: output has no dims!\n"
785				unless(@dd);
786				my $ddtotal = 1;
787				$ddtotal *= $_ for @dd;
788				barf "map: output has no elements (at least one dim is 0)!\n"
789				unless($ddtotal);
790
791				##############################
792				# If necessary, generate an appropriate FITS header for the output.
793
794				my $nofits = _opt($opt, ['nf','nofits','NoFITS','pix','pixel','Pixel']);
795
796				##############################
797				# Autoscale by transforming a subset of the input points' coordinates
798				# to the output range, and pick a FITS header that fits the output
799				# coordinates into the given template.
800				#
801				# Autoscaling always produces a simple, linear mapping in the FITS header.
802				# We support more complex mappings (via t_fits) but only to match a pre-existing
803				# FITS header (which doesn't use autoscaling).
804				#
805				# If the rectify option is set (the default) then the image is rectified
806				# in scientific coordinates; if it is clear, then the existing matrix
807				# is used, preserving any shear or rotation in the coordinate system.
808				# Since we eschew CROTA whenever possible, the CDi_j formalism is used instead.
809				my $f_in = (defined($in->hdr->{NAXIS}) ? t_fits($in,{ignore_rgb=>1}) : t_identity());
810
811				unless((defined $out->gethdr && $out->hdr->{NAXIS}) or $nofits) {
812				print "generating output FITS header..." if($PDL::Transform::debug);
813
814				$out->sethdr($in->hdr_copy) # Copy extraneous fields...
815				if(defined $in->hdr);
816
817				my $samp_ratio = 300;
818
819				my $orange = _opt($opt, ['or','orange','output_range','Output_Range'],
820				undef);
821
822				my $omin;
823				my $omax;
824				my $osize;
825
826				my $rectify = _opt($opt,['r','rect','rectify','Rectify'],1);
827
828				if (defined $orange) {
829				# orange always rectifies the coordinates -- the output scientific
830				# coordinates must align with the axes, or orange wouldn't make
831				# sense.
832				print "using user's orange..." if($PDL::Transform::debug);
833				$orange = pdl($orange) unless(UNIVERSAL::isa($orange,'PDL'));
834				barf "map: orange must be 2xN for an N-D transform"
835				unless ( (($orange->dim(1)) == $nd )
836				&& $orange->ndims == 2);
837
838				$omin = $orange->slice("(0)");
839				$omax = $orange->slice("(1)");
840				$osize = $omax - $omin;
841
842				$rectify = 1;
843
844				} else {
845
846				##############################
847				# Real autoscaling happens here.
848
849				if(!$rectify and ref( $f_in ) !~ /Linear/i) {
850				if( $f_in->{name} ne 'identity' ) {
851				print STDERR "Warning: map can't preserve nonlinear FITS distortions while autoscaling.\n";
852				}
853				$rectify=1;
854				}
855
856				my $f_tr = ( $rectify ?
857				$me x $f_in :
858				( ($me->{name} eq 'identity') ? # Simple optimization for match()
859				$me : # identity -- just matching
860				!$f_in x $me x $f_in # common case
861				)
862				);
863
864				my $samps = (pdl(($in->dims)[0..$nd-1]))->clip(0,$samp_ratio);
865
866				my $coords = ndcoords(($samps + 1)->list);
867
868				my $t;
869				my $irange = _opt($opt, ['ir','irange','input_range','Input_Range'],
870				undef);
871
872				# If input range is defined, sample that quadrilateral -- else
873				# sample the quad defined by the boundaries of the input image.
874				if(defined $irange) {
875				print "using user's irange..." if($PDL::Transform::debug);
876				$irange = pdl($irange) unless(UNIVERSAL::isa($irange,'PDL'));
877				barf "map: irange must be 2xN for an N-D transform"
878				unless ( (($irange->dim(1)) == $nd )
879				&& $irange->ndims == 2);
880
881				$coords *= ($irange->slice("(1)") - $irange->slice("(0)")) / $samps;
882				$coords += $irange->slice("(0)");
883				$coords -= 0.5; # offset to pixel corners...
884				$t = $me;
885				} else {
886				$coords *= pdl(($in->dims)[0..$nd-1]) / $samps;
887				$coords -= 0.5; # offset to pixel corners...
888				$t = $f_tr;
889				}
890				my $ocoords = $t->apply($coords)->mv(0,-1)->clump($nd);
891
892				# discard non-finite entries
893				my $oc2 = $ocoords->range(
894				which(
895				$ocoords->
896				transpose->
897				sumover->
898				isfinite
899				)
900				->dummy(0,1)
901				);
902
903				$omin = $oc2->minimum;
904				$omax = $oc2->maximum;
905
906				$osize = $omax - $omin;
907				my $tosize;
908				($tosize = $osize->where($osize == 0)) .= 1.0;
909				}
910
911				my ($scale) = $osize / pdl(($out->dims)[0..$nd-1]);
912
913				my $justify = _opt($opt,['j','J','just','justify','Justify'],0);
914				if($justify) {
915				$scale->slice("0") *= $justify;
916				$scale .= $scale->max;
917				$scale->slice("0") /= $justify;
918				}
919
920				print "done with autoscale. Making fits header....\n" if($PDL::Transform::debug);
921				if( $rectify ) {
922				# Rectified header generation -- make a simple coordinate header with no
923				# rotation or skew.
924				print "rectify\n" if($PDL::Transform::debug);
925				for my $d(1..$nd) {
926				$out->hdr->{"CRPIX$d"} = 1 + ($out->dim($d-1)-1)/2 ;
927				$out->hdr->{"CDELT$d"} = $scale->at($d-1);
928				$out->hdr->{"CRVAL$d"} = ( $omin->at($d-1) + $omax->at($d-1) ) /2 ;
929				$out->hdr->{"NAXIS$d"} = $out->dim($d-1);
930				$out->hdr->{"CTYPE$d"} = ( (!defined($me->{otype}) ? "" :
931				ref($me->{otype}) ne 'ARRAY' ? $me->{otype} :
932				$me->{otype}[$d-1])
933				\|\| $in->hdr->{"CTYPE$d"}
934				\|\| "");
935				$out->hdr->{"CUNIT$d"} = ( (!defined($me->{ounit}) ? "" :
936				ref($me->{ounit}) ne 'ARRAY' ? $me->{ounit} :
937				$me->{ounit}[$d-1])
938				\|\| $in->hdr->{"CUNIT$d"}
939				\|\| $in->hdr->{"CTYPE$d"}
940				\|\| "");
941				}
942				$out->hdr->{"NAXIS"} = $nd;
943
944				$out->hdr->{"SIMPLE"} = 'T';
945				$out->hdr->{"HISTORY"} .= "Header written by PDL::Transform::Cartography::map";
946
947				### Eliminate fancy newfangled output header pointing tags if they exist
948				### These are the CROTA, PCi_j, and CDi_j.
949				delete @{$out->hdr}{
950				grep /(^CROTA\d*$)\|(^(CD\|PC)\d+_\d+[A-Z]?$)/, keys %{$out->hdr}
951				#line 1768 "lib/PDL/Transform.pd"
952				};
953				} else {
954				# Non-rectified output -- generate a CDi_j matrix instead of the simple formalism.
955				# We have to deal with a linear transformation: we've got: (scaling) x !input x (t x input),
956				# where input is a linear transformation with offset and scaling is a simple scaling. We have
957				# the scaling parameters and the matrix for !input -- we have only to merge them and then we
958				# can write out the FITS header in CDi_j form.
959				print "non-rectify\n" if($PDL::Transform::debug);
960				my $midpoint_val = (pdl(($out->dims)[0..$nd-1])/2 * $scale)->apply( $f_in );
961				print "midpoint_val is $midpoint_val\n" if($PDL::Transform::debug);
962				# linear transformation
963				unless(ref($f_in) =~ m/Linear/) {
964				croak("Whups -- got a nonlinear t_fits transformation. Can't deal with it.");
965				}
966
967				my $inv_sc_mat = PDL::MatrixOps::stretcher($scale);
968				my $mat = $f_in->{params}->{matrix} x $inv_sc_mat;
969				print "scale is $scale; mat is $mat\n" if($PDL::Transform::debug);
970
971				print "looping dims 1..$nd: " if($PDL::Transform::debug);
972				for my $d(1..$nd) {
973				print "$d..." if($PDL::Transform::debug);
974				$out->hdr->{"CRPIX$d"} = 1 + ($out->dim($d-1)-1)/2;
975				$out->hdr->{"CRVAL$d"} = $midpoint_val->at($d-1);
976				$out->hdr->{"NAXIS$d"} = $out->dim($d-1);
977				$out->hdr->{"CTYPE$d"} = ( (!defined($me->{otype}) ? "" :
978				ref($me->{otype}) ne 'ARRAY' ? $me->{otype} :
979				$me->{otype}[$d-1])
980				\|\| $in->hdr->{"CTYPE$d"}
981				\|\| "");
982				$out->hdr->{"CUNIT$d"} = ( (!defined($me->{ounit}) ? "" :
983				ref($me->{ounit}) ne 'ARRAY' ? $me->{ounit} :
984				$me->{ounit}[$d-1])
985				\|\| $in->hdr->{"CUNIT$d"}
986				\|\| $in->hdr->{"CTYPE$d"}
987				\|\| "");
988				for my $e(1..$nd) {
989				$out->hdr->{"CD${d}_${e}"} = $mat->at($d-1,$e-1);
990				print "setting CD${d}_${e} to ".$mat->at($d-1,$e-1)."\n" if($PDL::Transform::debug);
991				}
992				}
993
994				## Eliminate competing header pointing tags if they exist
995				delete @{$out->hdr}{
996				grep /(^CROTA\d$)\|(^(PC)\d+_\d+[A-Z]?$)\|(CDELT\d$)/, keys %{$out->hdr}
997				#line 1813 "lib/PDL/Transform.pd"
998				};
999				}
1000				}
1001
1002				$out->hdrcpy(1);
1003
1004				##############################
1005				# Sandwich the transform between the input and output plane FITS headers.
1006				unless($nofits) {
1007				$me = !t_fits($out,{ignore_rgb=>1}) x $me x $f_in;
1008				}
1009
1010				##############################
1011				## Figure out the interpND options
1012				my $method = _opt($opt,['m','method','Method'],undef);
1013				my $bound = _opt($opt,['b','bound','boundary','Boundary'],'t');
1014
1015				##############################
1016				## Rubber meets the road: calculate the inverse transformed points.
1017				my $ndc = PDL::Basic::ndcoords(@dd);
1018				my $idx = $me->invert($ndc->double);
1019
1020				barf "map: Transformation had no inverse\n" unless defined($idx);
1021
1022				##############################
1023				## Integrate ? (Jacobian, Gaussian, Hanning)
1024				my $integrate = ($method =~ m/^[jghJGH]/) if defined($method);
1025
1026				##############################
1027				## Sampling code:
1028				## just transform and interpolate.
1029				## ( Kind of an anticlimax after all that, eh? )
1030				if(!$integrate) {
1031				my $x = $in->interpND($idx,{method=>$method, bound=>$bound});
1032				$out->slice(":") .= $x; # trivial slice prevents header overwrite...
1033				return $out;
1034				}
1035
1036				##############################
1037				## Anti-aliasing code:
1038				## Condition the input and call the pixelwise C interpolator.
1039				##
1040
1041				barf("PDL::Transform::map: Too many dims in transformation\n")
1042				if($in->ndims < $idx->ndims-1);
1043
1044				####################
1045				## Condition the broadcasting -- pixelwise interpolator only broadcasts
1046				## in 1 dimension, so squish all broadcast dimensions into 1, if necessary
1047				my @iddims = $idx->dims;
1048				my $in2 = $in->ndims == $#iddims
1049				? $in->dummy(-1,1)
1050				: $in->reorder($nd..$in->ndims-1, 0..$nd-1)
1051				->clump($in->ndims - $nd)
1052				->mv(0,-1);
1053
1054				####################
1055				# Allocate the output array
1056				my $o2 = PDL->new_from_specification($in2->type,
1057				@iddims[1..$#iddims],
1058				$in2->dim(-1)
1059				);
1060
1061				####################
1062				# Pack boundary string if necessary
1063				if(defined $bound) {
1064				if(ref $bound eq 'ARRAY') {
1065				my ($s,$el);
1066				foreach $el(@$bound) {
1067				barf "Illegal boundary value '$el' in range"
1068				unless( $el =~ m/^([0123fFtTeEpPmM])/ );
1069				$s .= $1;
1070				}
1071				$bound = $s;
1072				}
1073				elsif($bound !~ m/^[0123ftepx]+$/ && $bound =~ m/^([0123ftepx])/i ) {
1074				$bound = $1;
1075				}
1076				}
1077
1078				####################
1079				# Get the blur and minimum-sv values
1080				my $blur = _opt($opt,['blur','Blur'],1.0);
1081				my $svmin = _opt($opt,['sv','SV'],1.0);
1082				my $big = _opt($opt,['big','Big'],1.0);
1083				my $flux = _opt($opt,['phot','photometry'],0);
1084				my @idims = $in->dims;
1085
1086				$flux = ($flux =~ m/^[1fF]/);
1087				$big = $big * max(pdl(@idims[0..$nd]));
1088				$blur = $blur->at(0) if(ref $blur);
1089				$svmin = $svmin->at(0) if(ref $svmin);
1090
1091				my $bv = $in->badflag ? $in->badvalue->sclr : 0;
1092
1093				### The first argument is a dummy to set $GENERIC.
1094				$idx = double($idx) unless($idx->type == double);
1095				print "Calling _map_int...\n" if($PDL::Transform::debug);
1096				&PDL::_map_int( $in2->flat->index(0),
1097				$in2, $o2, $idx,
1098				$bound, $method, $big, $blur, $svmin, $flux, $bv);
1099
1100				my @rdims = (@iddims[1..$#iddims], @idims[$#iddims..$#idims]);
1101				{
1102				$out->slice(":") .= $o2->reshape(@rdims);
1103				}
1104				return $out;
1105				}
1106				#line 1107 "lib/PDL/Transform.pm"
1107
1108				*map = \&PDL::map;
1109
1110
1111
1112
1113
1114				#line 1928 "lib/PDL/Transform.pd"
1115
1116				######################################################################
1117
1118				=head2 unmap
1119
1120				=for sig
1121
1122				Signature: (data(); PDL::Transform a; template(); \%opt)
1123
1124				=for usage
1125
1126				$out_image = $in_image->unmap($t,[],[
1127				$out_image = $t->unmap($in_image,[],[
1128
1129				=for ref
1130
1131				Map an image or N-D dataset using the inverse as a coordinate transform.
1132
1133				This convenience function just inverts $t and calls L on
1134				the inverse; everything works the same otherwise. For convenience, it
1135				is both a PDL method and a PDL::Transform method.
1136
1137				=cut
1138
1139				*PDL::unmap = \&unmap;
1140				sub unmap {
1141				my($me) = shift;
1142				my($data) = shift;
1143				my(@params) = @_;
1144
1145				if(UNIVERSAL::isa($data,'PDL::Transform') && UNIVERSAL::isa($me,'PDL')) {
1146				my $x = $data;
1147				$data = $me;
1148				$me = $x;
1149				}
1150
1151				return $me->inverse->map($data,@params);
1152				}
1153
1154				#line 1971 "lib/PDL/Transform.pd"
1155
1156				=head2 t_inverse
1157
1158				=for usage
1159
1160				$t2 = t_inverse($t);
1161				$t2 = $t->inverse;
1162				$t2 = $t ** -1;
1163				$t2 = !$t;
1164
1165				=for ref
1166
1167				Return the inverse of a PDL::Transform. This just reverses the
1168				func/inv, idim/odim, itype/otype, and iunit/ounit pairs. Note that
1169				sometimes you end up with a transform that cannot be applied or
1170				mapped, because either the mathematical inverse doesn't exist or the
1171				inverse func isn't implemented.
1172
1173				You can invert a transform by raising it to a negative power, or by
1174				negating it with '!'.
1175
1176				The inverse transform remains connected to the main transform because
1177				they both point to the original parameters hash. That turns out to be
1178				useful.
1179
1180				=cut
1181
1182				*t_inverse = \&inverse;
1183
1184				sub inverse {
1185				my($me) = shift;
1186				barf("PDL::Transform::inverse: got a transform with no inverse")
1187				unless defined $me->{inv};
1188				bless {
1189				%$me, # force explicit copy of top-level
1190				inv => $me->{func},
1191				func => $me->{inv},
1192				(map +("o$_"=>$me->{"i$_"}, "i$_"=>$me->{"o$_"}), qw(dim type unit)),
1193				name => "(inverse ".$me->{name}.")",
1194				is_inverse => !($me->{is_inverse}),
1195				}, ref $me;
1196				}
1197
1198				#line 2018 "lib/PDL/Transform.pd"
1199
1200				=head2 t_compose
1201
1202				=for usage
1203
1204				$f2 = t_compose($f, $g,[...]);
1205				$f2 = $f->compose($g[,$h,$i,...]);
1206				$f2 = $f x $g x ...;
1207
1208				=for ref
1209
1210				Function composition: f(g(x)), f(g(h(x))), ...
1211
1212				You can also compose transforms using the overloaded matrix-multiplication
1213				(nee repeat) operator 'x'.
1214
1215				This is accomplished by inserting a splicing code ref into the C
1216				and C slots. It combines multiple compositions into a single
1217				list of transforms to be executed in order, from last to first (in
1218				keeping with standard mathematical notation). If one of the functions is
1219				itself a composition, it is interpolated into the list rather than left
1220				separate. Ultimately, linear transformations may also be combined within
1221				the list.
1222
1223				No checking is done that the itype/otype and iunit/ounit fields are
1224				compatible -- that may happen later, or you can implement it yourself
1225				if you like.
1226
1227				=cut
1228
1229				@PDL::Transform::Composition::ISA = ('PDL::Transform');
1230				sub PDL::Transform::Composition::stringify {
1231				package PDL::Transform::Composition;
1232				shift->SUPER::stringify;
1233				}
1234
1235				*t_compose = \&compose;
1236				sub compose {
1237				local($_);
1238				my(@funcs) = @_;
1239				my($me) = PDL::Transform->new;
1240				# No inputs case: return the identity function
1241				return $me if !@funcs;
1242				$me->{name} = "";
1243				my @clist;
1244				for my $f (@funcs) {
1245				$me->{idim} = $f->{idim} if $f->{idim} > $me->{idim};
1246				$me->{odim} = $f->{odim} if $f->{odim} > $me->{odim};
1247				if (UNIVERSAL::isa($f,"PDL::Transform::Composition")) {
1248				if($f->{is_inverse}) {
1249				for(reverse(@{$f->{params}->{clist}})) {
1250				push(@clist,$_->inverse);
1251				$me->{name} .= " o inverse ( ".$_->{name}." )";
1252				}
1253				} else {
1254				for(@{$f->{params}->{clist}}) {
1255				push(@clist,$_);
1256				$me->{name} .= " o ".$_->{name};
1257				}
1258				}
1259				} else { # Not a composition -- just push the transform onto the list.
1260				push(@clist,$f);
1261				$me->{name} .= " o ".$f->{name};
1262				}
1263				}
1264				$me->{name}=~ s/^ o //; # Get rid of leading composition mark
1265				$me->{otype} = $funcs[0]->{otype};
1266				$me->{ounit} = $funcs[0]->{ounit};
1267				$me->{itype} = $funcs[-1]->{itype};
1268				$me->{iunit} = $funcs[-1]->{iunit};
1269				$me->{params}->{clist} = \@clist;
1270				$me->{func} = sub {
1271				my ($data,$p) = @_;
1272				my ($ip) = $data->is_inplace;
1273				for my $t ( reverse @{$p->{clist}} ) {
1274				croak("Error: tried to apply a PDL::Transform with no function inside a composition!\n offending transform: $t\n")
1275				unless(defined($t->{func}) and ref($t->{func}) eq 'CODE');
1276				$data = $t->{func}($ip ? $data->inplace : $data, $t->{params});
1277				}
1278				$data->is_inplace(0); # clear inplace flag (avoid core bug with inplace)
1279				$data;
1280				};
1281				$me->{inv} = sub {
1282				my ($data,$p) = @_;
1283				my $ip = $data->is_inplace;
1284				for my $t ( @{$p->{clist}} ) {
1285				croak("Error: tried to invert a non-invertible PDL::Transform inside a composition!\n offending transform: $t\n")
1286				unless(defined($t->{inv}) and ref($t->{inv}) eq 'CODE');
1287				$data = &{$t->{inv}}($ip ? $data->inplace : $data, $t->{params});
1288				}
1289				$data;
1290				};
1291				return bless($me,'PDL::Transform::Composition');
1292				}
1293
1294				#line 2117 "lib/PDL/Transform.pd"
1295
1296				=head2 t_wrap
1297
1298				=for usage
1299
1300				$g1fg = $f->wrap($g);
1301				$g1fg = t_wrap($f,$g);
1302
1303				=for ref
1304
1305				Shift a transform into a different space by 'wrapping' it with a second.
1306
1307				This is just a convenience function for two
1308				L calls. C<< $x->wrap($y) >> is the same as
1309				C<(!$y) x $x x $y>: the resulting transform first hits the data with
1310				$y, then with $x, then with the inverse of $y.
1311
1312				For example, to shift the origin of rotation, do this:
1313
1314				$im = rfits('m51.fits');
1315				$tf = t_fits($im);
1316				$tr = t_linear({rot=>30});
1317				$im1 = $tr->map($tr); # Rotate around pixel origin
1318				$im2 = $tr->map($tr->wrap($tf)); # Rotate round FITS scientific origin
1319
1320				=cut
1321
1322				*t_wrap = \&wrap;
1323
1324				sub wrap {
1325				my($f) = shift;
1326				my($g) = shift;
1327
1328				return $g->inverse->compose($f,$g);
1329				}
1330
1331				######################################################################
1332
1333				# Composition operator -- handles 'x'.
1334				sub _compose_op {
1335				my($x,$y,$c) = @_;
1336				$c ? compose($y,$x) : compose($x,$y);
1337				}
1338
1339				# Raise-to-power operator -- handles '**'.
1340
1341				sub _pow_op {
1342				my($x,$y,$c) = @_;
1343
1344				barf("Can't raise anything to the power of a transform")
1345				if($c \|\| UNIVERSAL::isa($y,'PDL::Transform')) ;
1346
1347				$x = $x->inverse
1348				if($y < 0);
1349
1350				return $x if(abs($y) == 1);
1351				return new PDL::Transform if(abs($y) == 0);
1352
1353				my(@l);
1354				for my $i(1..abs($y)) {
1355				push(@l,$x);
1356				}
1357
1358				t_compose(@l);
1359				}
1360
1361				#line 2190 "lib/PDL/Transform.pd"
1362
1363				=head2 t_identity
1364
1365				=for usage
1366
1367				my $xform = t_identity
1368				my $xform = PDL::Transform->new;
1369
1370				=for ref
1371
1372				Generic constructor generates the identity transform.
1373
1374				This constructor really is trivial -- it is mainly used by the other transform
1375				constructors. It takes no parameters and returns the identity transform.
1376
1377				=cut
1378
1379				sub _identity { return shift; }
1380				sub t_identity { new PDL::Transform(@_) };
1381
1382				sub new {
1383				my($class) = shift;
1384				my $me = {name=>'identity',
1385				idim => 0,
1386				odim => 0,
1387				func=>\&PDL::Transform::_identity,
1388				inv=>\&PDL::Transform::_identity,
1389				params=>{}
1390				};
1391
1392				return bless $me,$class;
1393				}
1394
1395				#line 2228 "lib/PDL/Transform.pd"
1396
1397				=head2 t_lookup
1398
1399				=for usage
1400
1401				$f = t_lookup($lookup, {});
1402
1403				=for ref
1404
1405				Transform by lookup into an explicit table.
1406
1407				You specify an N+1-D PDL that is interpreted as an N-D lookup table of
1408				column vectors (vector index comes last). The last dimension has
1409				order equal to the output dimensionality of the transform.
1410
1411				For added flexibility in data space, You can specify pre-lookup linear
1412				scaling and offset of the data. Of course you can specify the
1413				interpolation method to be used. The linear scaling stuff is a little
1414				primitive; if you want more, try composing the linear transform with
1415				this one.
1416
1417				The prescribed values in the lookup table are treated as
1418				pixel-centered: that is, if your input array has N elements per row
1419				then valid data exist between the locations (-0.5) and (N-0.5) in
1420				lookup pixel space, because the pixels (which are numbered from 0 to
1421				N-1) are centered on their locations.
1422
1423				Lookup is done using L, so the boundary conditions
1424				and broadcasting behaviour follow from that.
1425
1426				The indexed-over dimensions come first in the table, followed by a
1427				single dimension containing the column vector to be output for each
1428				set of other dimensions -- ie to output 2-vectors from 2 input
1429				parameters, each of which can range from 0 to 49, you want an index
1430				that has dimension list (50,50,2). For the identity lookup table
1431				you could use C.
1432
1433				If you want to output a single value per input vector, you still need
1434				that last index broadcasting dimension -- if necessary, use C.
1435
1436				The lookup index scaling is: out = lookup[ (scale * data) + offset ].
1437
1438				A simplistic table inversion routine is included. This means that
1439				you can (for example) use the C method with C transformations.
1440				But the table inversion is exceedingly slow, and not practical for tables
1441				larger than about 100x100. The inversion table is calculated in its
1442				entirety the first time it is needed, and then cached until the object is
1443				destroyed.
1444
1445				Options are listed below; there are several synonyms for each.
1446
1447				=over 3
1448
1449				=item s, scale, Scale
1450
1451				(default 1.0) Specifies the linear amount of scaling to be done before
1452				lookup. You can feed in a scalar or an N-vector; other values may cause
1453				trouble. If you want to save space in your table, then specify smaller
1454				scale numbers.
1455
1456				=item o, offset, Offset
1457
1458				(default 0.0) Specifies the linear amount of offset before lookup.
1459				This is only a scalar, because it is intended to let you switch to
1460				corner-centered coordinates if you want to (just feed in o=-0.25).
1461
1462				=item b, bound, boundary, Boundary
1463
1464				Boundary condition to be fed to L
1465
1466				=item m, method, Method
1467
1468				Interpolation method to be fed to L
1469
1470				=back
1471
1472				EXAMPLE
1473
1474				To scale logarithmically the Y axis of m51, try:
1475
1476				$x = float rfits('m51.fits');
1477				$lookup = xvals(128,128) -> cat( 10*(yvals(128,128)/50) 256/10**2.55 );
1478				$t = t_lookup($lookup);
1479				$y = $t->map($x);
1480
1481				To do the same thing but with a smaller lookup table, try:
1482
1483				$lookup = 16 * xvals(17,17)->cat(10*(yvals(17,17)/(100/16)) 16/10**2.55);
1484				$t = t_lookup($lookup,{scale=>1/16.0});
1485				$y = $t->map($x);
1486
1487				(Notice that, although the lookup table coordinates are is divided by 16,
1488				it is a 17x17 -- so linear interpolation works right to the edge of the original
1489				domain.)
1490
1491				NOTES
1492
1493				Inverses are "only weakly supported" -- the best way to do it might be by
1494				judicious use of map() on the forward transformation.
1495
1496				the type/unit fields are ignored.
1497
1498				=cut
1499
1500				sub t_lookup {
1501				my($class) = 'PDL::Transform';
1502				my($source)= shift;
1503				my($o) = shift;
1504
1505				if(!defined($o) && ((ref $source) eq 'HASH')) {
1506				Carp::cluck("lookup transform called as sub not method; using 'PDL::Transform' as class...\n");
1507				$o = $source;
1508				$source = $class;
1509				$class = "PDL::Transform";
1510				}
1511
1512				$o = {} unless(ref $o eq 'HASH');
1513
1514				my($me) = $class->new;
1515
1516				my($bound) = _opt($o,['b','bound','boundary','Boundary']);
1517				my($method)= _opt($o,['m','meth','method','Method']);
1518
1519				$me->{idim} = $source->ndims - 1;
1520				$me->{odim} = $source->dim($source->ndims-1);
1521
1522				$me->{params} = {
1523				table => $source,
1524				scale => _opt($o,['s','scale','Scale'],1.0),
1525				offset => _opt($o,['o','off','offset','Offset'],0.0),
1526				interpND_opt => {
1527				method => $method,
1528				bound => $bound,
1529				bad => _opt($o,['bad'],0)
1530				}
1531				};
1532
1533				my $lookup_func = sub {
1534				my($data,$p,$table,$scale,$offset) = @_;
1535
1536				$data = pdl($data)
1537				unless ((ref $data) && (UNIVERSAL::isa($data,'PDL')));
1538				$main::foo = $data;
1539
1540				if($data->dim(0) > $me->{idim}) {
1541				barf("Too many dims (".$data->dim(0).") for your table (".$me->{idim}.")\n");
1542				};
1543
1544				my($x)= ($table
1545				->interpND(float($data) * $scale + $offset,
1546				$p->{interpND_opt}
1547				)
1548				);
1549
1550				# Put the index dimension (and broadcasted indices) back at the front of
1551				# the dimension list.
1552				my($dnd) = $data->ndims - 1;
1553				return ($x -> ndims > $data->ndims - 1) ?
1554				($x->reorder( $dnd..($dnd + $table->ndims - $data->dim(0)-1)
1555				, 0..$data->ndims-2
1556				)
1557				) : $x;
1558				};
1559
1560				$me->{func} = sub {my($data,$p) = @_; &$lookup_func($data,$p,$p->{table},$p->{scale},$p->{offset})};
1561
1562				#######
1563				## Lazy inverse -- find it if and only if we need it...
1564				$me->{inv} = sub {
1565				my $data = shift;
1566				my $p = shift;
1567				if(!defined($p->{'itable'})) {
1568				if($me->{idim} != $me->{odim}) {
1569				barf("t_lookup: can't calculate an inverse of a projection operation! (idim != odim)");
1570				}
1571				print "t_lookup: Warning, table inversion is only weakly supported. \n I'll try to invert it using a pretty boneheaded algorithm...\n (If it takes too long, consider writing a faster algorithm!)\n Calculating inverse table (will be cached)...\n" if($PDL::verbose \|\| $PDL::debug \|\| $PDL::Transform::debug);
1572				my $itable = zeroes($p->{table});
1573				my $minvals = $p->{table}->clump($me->{idim})->minimum;
1574				my $maxvals = $p->{table}->clump($me->{idim})->maximum;
1575
1576				# Scale so that the range runs from 0 through the top pixel in the table
1577				my $scale = ( $itable->shape->slice("0:-2")-1 ) /
1578				(($maxvals - $minvals)+ (($maxvals-$minvals) == 0));
1579				my $offset = - ($minvals * $scale);
1580
1581				$p->{iscale} = $scale;
1582				$p->{ioffset} = $offset;
1583
1584				my $enl_scale = $p->{'enl_scale'} \|\| 10;
1585				my $smallcoords = ndcoords((pdl($enl_scale * 2 + 1)->at(0)) x $me->{idim})/ $enl_scale - 1;
1586
1587				# $oloop runs over (point, index) for all points in the output table, in
1588				# $p->{table} output space
1589				my $oloop = ndcoords($itable->mv(-1,0)->slice("(0)"))->
1590				double->
1591				mv(0,-1)->
1592				clump($itable->ndims-1); # oloop: (pixel, index)
1593				{
1594				$oloop->mv(-1,0) -= $offset;
1595				$oloop->mv(-1,0) /= $scale;
1596				}
1597
1598				# The Right Thing to do here is to take the outer product of $itable and $otable, then collapse
1599				# to find minimum distance. But memory demands for that would be HUGE. Instead we do an
1600				# elementwise calculation.
1601
1602				print "t_lookup: inverting ".$oloop->dim(0)." points...\n" if($PDL::verbose \|\| $PDL::debug \|\| $PDL::Transform::debug);
1603				my $pt = $p->{table}->mv(-1,0); # pt runs (index, x,y,...)
1604
1605				my $itable_flattened = zeroes($oloop);
1606
1607				for my $i (0..$oloop->dim(0)-1) {
1608
1609				my $olp = $oloop->slice("($i)"); # olp runs (index)
1610				my $diff = ($pt - $olp); # diff runs (index, x, y, ...)
1611				my $r2 = ($diff * $diff)->sumover; # r2 runs (x,y,...)
1612				my $c = whichND($r2==$r2->min)->slice(":,(0)"); # c runs (index)
1613
1614				# Now zero in on the neighborhood around the point of closest approach.
1615				my $neighborhood = $p->{table}->interpND($c + $smallcoords,{method=>'linear',bound=>'t'})->
1616				mv(-1,0); # neighborhood runs (index, dx, dy,...)
1617				$diff = $neighborhood - $olp; # diff runs (index, dx, dy, ...)
1618				$r2 = ($diff * $diff)->sumover; # r2 runs (dx,dy,...)
1619				my $dc = $smallcoords->mv(0,-1)->indexND(0+whichND($r2==$r2->min)->slice(":,(0)"));
1620
1621				my $coord = $c + $dc;
1622				# At last, we've found the best-fit to an enl_scale'th of an input-table pixel.
1623				# Back it out to input-science coordinates, and stuff it in the inverse table.
1624				$itable_flattened->slice("($i)") .= $coord;
1625
1626				print " $i..." if( ($i%1000==0) && ( $PDL::verbose \|\| $PDL::debug \|\| $PDL::Transform::debug));
1627				}
1628
1629				{
1630				$itable->clump($itable->ndims-1) .= $itable_flattened;
1631				$itable->mv(-1,0) -= $p->{offset};
1632				$itable->mv(-1,0) /= $p->{scale};
1633				}
1634
1635				$p->{itable} = $itable;
1636				}
1637				&$lookup_func($data,$p, $p->{itable},$p->{iscale},$p->{ioffset}) ;
1638				};
1639
1640				$me->{name} = 'Lookup';
1641
1642				return $me;
1643				}
1644
1645				#line 2485 "lib/PDL/Transform.pd"
1646
1647				=head2 t_linear
1648
1649				=for usage
1650
1651				$f = t_linear({options});
1652
1653				=for ref
1654
1655				Linear (affine) transformations with optional offset
1656
1657				t_linear implements simple matrix multiplication with offset,
1658				also known as the affine transformations.
1659
1660				You specify the linear transformation with pre-offset, a mixing
1661				matrix, and a post-offset. That overspecifies the transformation, so
1662				you can choose your favorite method to specify the transform you want.
1663				The inverse transform is automagically generated, provided that it
1664				actually exists (the transform matrix is invertible). Otherwise, the
1665				inverse transform just croaks.
1666
1667				Extra dimensions in the input vector are ignored, so if you pass a
1668				3xN vector into a 3-D linear transformation, the final dimension is passed
1669				through unchanged.
1670
1671				The options you can usefully pass in are:
1672
1673				=over 3
1674
1675				=item s, scale, Scale
1676
1677				A scaling scalar (heh), vector, or matrix. If you specify a vector
1678				it is treated as a diagonal matrix (for convenience). It gets
1679				left-multiplied with the transformation matrix you specify (or the
1680				identity), so that if you specify both a scale and a matrix the
1681				scaling is done after the rotation or skewing or whatever.
1682
1683				=item r, rot, rota, rotation, Rotation
1684
1685				A rotation angle in degrees -- useful for 2-D and 3-D data only. If
1686				you pass in a scalar, it specifies a rotation from the 0th axis toward
1687				the 1st axis. If you pass in a 3-vector as either a PDL or an array
1688				ref (as in "rot=>[3,4,5]"), then it is treated as a set of Euler
1689				angles in three dimensions, and a rotation matrix is generated that
1690				does the following, in order:
1691
1692				=over 3
1693
1694				=item * Rotate by rot->(2) degrees from 0th to 1st axis
1695
1696				=item * Rotate by rot->(1) degrees from the 2nd to the 0th axis
1697
1698				=item * Rotate by rot->(0) degrees from the 1st to the 2nd axis
1699
1700				=back
1701
1702				The rotation matrix is left-multiplied with the transformation matrix
1703				you specify, so that if you specify both rotation and a general matrix
1704				the rotation happens after the more general operation -- though that is
1705				deprecated.
1706
1707				Of course, you can duplicate this functionality -- and get more
1708				general -- by generating your own rotation matrix and feeding it in
1709				with the C option.
1710
1711				=item m, matrix, Matrix
1712
1713				The transformation matrix. It does not even have to be square, if you want
1714				to change the dimensionality of your input. If it is invertible (note:
1715				must be square for that), then you automagically get an inverse transform too.
1716
1717				=item pre, preoffset, offset, Offset
1718
1719				The vector to be added to the data before they get multiplied by the matrix
1720				(equivalent of CRVAL in FITS, if you are converting from scientific to
1721				pixel units).
1722
1723				=item post, postoffset, shift, Shift
1724
1725				The vector to be added to the data after it gets multiplied by the matrix
1726				(equivalent of CRPIX-1 in FITS, if youre converting from scientific to pixel
1727				units).
1728
1729				=item d, dim, dims, Dims
1730
1731				Most of the time it is obvious how many dimensions you want to deal
1732				with: if you supply a matrix, it defines the transformation; if you
1733				input offset vectors in the C and C options, those define
1734				the number of dimensions. But if you only supply scalars, there is no way
1735				to tell and the default number of dimensions is 2. This provides a way
1736				to do, e.g., 3-D scaling: just set C<{s=>, dims=>3}> and
1737				you are on your way.
1738
1739				=item idim
1740
1741				=item odim
1742
1743				=item iunit
1744
1745				=item ounit
1746
1747				=item itype
1748
1749				=item otype
1750
1751				These are incorporated in the object, though not much is done with them
1752				yet as of 2.094.
1753
1754				=back
1755
1756				=cut
1757
1758				{ package PDL::Transform::Linear;
1759				our @ISA = ('PDL::Transform');
1760				*_opt = \&PDL::Transform::_opt;
1761				*barf = \&PDL::Transform::barf;
1762				*identity = \&PDL::MatrixOps::identity;
1763
1764				sub PDL::Transform::t_linear { PDL::Transform::Linear->new(@_); }
1765
1766				sub new {
1767				my($class) = shift;
1768				my($o) = $_[0];
1769				pop @_ if (($#_ % 2 ==0) && !defined($_[-1]));
1770				#suppresses a warning if @_ has an odd number of elements and the
1771				#last is undef
1772
1773				if(!(ref $o)) {
1774				$o = {@_};
1775				}
1776
1777				my($me) = __PACKAGE__->SUPER::new;
1778
1779				$me->{name} = "linear";
1780
1781				$me->{params}{pre} = _opt($o,['pre','Pre','preoffset','offset',
1782				'Offset','PreOffset','Preoffset'],0);
1783				$me->{params}{pre} = PDL->pdl($me->{params}{pre})
1784				if(defined $me->{params}{pre});
1785
1786				$me->{params}{post} = _opt($o,['post','Post','postoffset','PostOffset',
1787				'shift','Shift'],0);
1788				$me->{params}{post} = PDL->pdl($me->{params}{post})
1789				if(defined $me->{params}{post});
1790
1791				$me->{params}{matrix} = _opt($o,['m','matrix','Matrix','mat','Mat']);
1792				$me->{params}{matrix} = PDL->pdl($me->{params}{matrix})
1793				if(defined $me->{params}{matrix});
1794
1795				$me->{params}{rot} = _opt($o,['r','rot','rota','rotation','Rotation'], 0);
1796				$me->{params}{rot} = PDL->pdl($me->{params}{rot});
1797
1798				my $o_dims = _opt($o,['d','dim','dims','Dims']);
1799				$o_dims = PDL->pdl($o_dims) if defined($o_dims);
1800
1801				my $scale = _opt($o,['s','scale','Scale']);
1802				$scale = PDL->pdl($scale) if defined($scale);
1803
1804				# Figure out the number of dimensions to transform, and,
1805				# if necessary, generate a new matrix.
1806
1807				my $rot = $me->{params}{rot};
1808				if(defined($me->{params}{matrix})) {
1809				my $mat = $me->{params}{matrix} = $me->{params}{matrix}->slice(":,:");
1810				$me->{idim} = $mat->dim(0);
1811				$me->{odim} = $mat->dim(1);
1812
1813				} else {
1814				if(defined($me->{params}->{rot}) &&
1815				UNIVERSAL::isa($me->{params}->{rot},'PDL')) {
1816				$me->{idim} = $me->{odim} = 2 if($me->{params}{rot}->nelem == 1);
1817				$me->{idim} = $me->{odim} = 3 if($me->{params}{rot}->nelem == 3);
1818				}
1819
1820				if(defined($scale) &&
1821				UNIVERSAL::isa($scale,'PDL') &&
1822				$scale->getndims > 0) {
1823				$me->{idim} = $me->{odim} = $scale->dim(0);
1824				$me->{odim} = $scale->dim(0);
1825
1826				} elsif(defined($me->{params}{pre}) &&
1827				UNIVERSAL::isa($me->{params}{pre},'PDL') &&
1828				$me->{params}{pre}->getndims > 0) {
1829				$me->{idim} = $me->{odim} = $me->{params}{pre}->dim(0);
1830
1831				} elsif(defined($me->{params}{post}) &&
1832				UNIVERSAL::isa($me->{params}{post},'PDL') &&
1833				$me->{params}{post}->getndims > 0) {
1834				$me->{idim} = $me->{odim} = $me->{params}{post}->dim(0);
1835				} elsif(defined($o_dims)) {
1836				$me->{idim} = $me->{odim} = $o_dims;
1837				} elsif (UNIVERSAL::isa($rot,'PDL') && (my $nrots = $rot->nelem) > 1) {
1838				$me->{idim} = $me->{odim} = $nrots;
1839				} else {
1840				print "PDL::Transform::Linear: assuming 2-D transform (set dims option to change)\n" if($PDL::Transform::debug);
1841				$me->{idim} = $me->{odim} = 2;
1842				}
1843
1844				$me->{params}->{matrix} = PDL->zeroes($me->{idim},$me->{odim});
1845				$me->{params}->{matrix}->diagonal(0,1) .= 1;
1846				}
1847
1848				### Handle rotation option
1849				if(defined($rot)) {
1850				# Subrotation closure -- rotates from axis $d->(0) --> $d->(1).
1851				my $subrot = sub {
1852				my($d,$angle,$m)=@_;
1853				my($i) = identity($m->dim(0));
1854				my($subm) = $i->dice($d,$d);
1855
1856				$angle = $angle->at(0)
1857				if(UNIVERSAL::isa($angle,'PDL'));
1858
1859				my($x) = $angle * $PDL::Transform::DEG2RAD;
1860				$subm .= $subm x PDL->pdl([cos($x),-sin($x)],[sin($x),cos($x)]);
1861				$m .= $m x $i;
1862				};
1863
1864				if(UNIVERSAL::isa($rot,'PDL') && $rot->nelem > 1) {
1865				if($rot->ndims == 2) {
1866				$me->{params}{matrix} x= $rot;
1867				} elsif($rot->nelem == 3) {
1868				my $rm = identity(3);
1869
1870				# Do these in reverse order to make it more like
1871				# function composition!
1872				&$subrot(PDL->pdl(0,1),$rot->at(2),$rm);
1873				&$subrot(PDL->pdl(2,0),$rot->at(1),$rm);
1874				&$subrot(PDL->pdl(1,2),$rot->at(0),$rm);
1875
1876				$me->{params}{matrix} .= $me->{params}{matrix} x $rm;
1877				} else {
1878				barf("PDL::Transform::Linear: Got a strange rot option -- giving up.\n");
1879				}
1880				} else {
1881				if($rot != 0 and $me->{params}{matrix}->dim(0)>1) {
1882				&$subrot(PDL->pdl(0,1),$rot,$me->{params}{matrix});
1883				}
1884				}
1885				}
1886
1887				# Apply scaling
1888				$me->{params}{matrix} = $me->{params}{matrix}->slice(":,:");
1889				$me->{params}{matrix}->diagonal(0,1) *= $scale
1890				if defined($scale);
1891
1892				# Check for an inverse and apply it if possible
1893				my($o2);
1894				if($me->{params}{matrix}->det($o2 = {lu=>undef})) {
1895				$me->{params}{inverse} = $me->{params}{matrix}->inv($o2);
1896				} else {
1897				delete $me->{params}{inverse};
1898				}
1899
1900				$me->{params}{idim} = $me->{idim};
1901				$me->{params}{odim} = $me->{odim};
1902
1903				for (qw(iunit ounit itype otype)) {
1904				next unless defined (my $val = _opt($o,[$_]));
1905				$me->{$_} = $val;
1906				}
1907
1908				##############################
1909				# The meat -- just shift, matrix-multiply, and shift again.
1910				$me->{func} = sub {
1911				my ($in,$opt) = @_;
1912				my $d = $opt->{matrix}->dim(0)-1;
1913				barf("Linear transform: transform is @{[$d+1]}-D; data is ".$in->dim(0))
1914				if $in->dim(0) <= $d;
1915				my $x = $in->slice("0:$d")->copy + $opt->{pre};
1916				my $out = $in->is_inplace ? $in : $in->copy;
1917				$out->slice("0:$d") .= $x x $opt->{matrix} + $opt->{post};
1918				return $out;
1919				};
1920
1921				$me->{inv} = (defined $me->{params}{inverse}) ? sub {
1922				my ($in,$opt) = @_;
1923				my $d = $opt->{inverse}->dim(0)-1;
1924				barf("Linear transform: transform is @{[$d+1]}-D; data is ".$in->dim(0))
1925				if $in->dim(0) <= $d;
1926				my $x = $in->slice("0:$d")->copy - $opt->{post};
1927				my $out = $in->is_inplace ? $in : $in->copy;
1928				$out->slice("0:$d") .= $x x $opt->{inverse} - $opt->{pre};
1929				$out;
1930				} : undef;
1931
1932				return $me;
1933				}
1934
1935				sub stringify {
1936				my($me) = shift; my($out) = SUPER::stringify $me;
1937				my $mp = $me->{params};
1938
1939				if(!($me->{is_inverse})){
1940				$out .= "Pre-add: ".($mp->{pre})."\n"
1941				if(defined $mp->{pre});
1942
1943				$out .= "Post-add: ".($mp->{post})."\n"
1944				if(defined $mp->{post});
1945
1946				$out .= "Forward matrix:".($mp->{matrix})
1947				if(defined $mp->{matrix});
1948
1949				$out .= "Inverse matrix:".($mp->{inverse})
1950				if(defined $mp->{inverse});
1951				} else {
1952				$out .= "Pre-add: ".(-$mp->{post})."\n"
1953				if(defined $mp->{post});
1954
1955				$out .= "Post-add: ".(-$mp->{pre})."\n"
1956				if(defined $mp->{pre});
1957
1958				$out .= "Forward matrix:".($mp->{inverse})
1959				if(defined $mp->{inverse});
1960
1961				$out .= "Inverse matrix:".($mp->{matrix})
1962				if(defined $mp->{matrix});
1963				}
1964
1965				$out =~ s/\n/\n /go;
1966				$out;
1967				}
1968				}
1969
1970				#line 2813 "lib/PDL/Transform.pd"
1971
1972				=head2 t_scale
1973
1974				=for usage
1975
1976				$f = t_scale()
1977
1978				=for ref
1979
1980				Convenience interface to L.
1981
1982				t_scale produces a transform that scales around the origin by a fixed
1983				amount. It acts exactly the same as C<<< t_linear(Scale=>$scale) >>>.
1984
1985				=cut
1986
1987				sub t_scale {
1988				my($scale) = shift;
1989				my($y) = shift;
1990				return t_linear(scale=>$scale,%{$y})
1991				#line 2834 "lib/PDL/Transform.pd"
1992				if(ref $y eq 'HASH');
1993				t_linear(Scale=>$scale,$y,@_);
1994				}
1995
1996				#line 2841 "lib/PDL/Transform.pd"
1997
1998				=head2 t_offset
1999
2000				=for usage
2001
2002				$f = t_offset()
2003
2004				=for ref
2005
2006				Convenience interface to L.
2007
2008				t_offset produces a transform that shifts the origin to a new location.
2009				It acts exactly the same as C<<< t_linear(Pre=>$shift) >>>.
2010
2011				=cut
2012
2013				sub t_offset {
2014				my($pre) = shift;
2015				my($y) = shift;
2016				return t_linear(pre=>$pre,%{$y})
2017				#line 2862 "lib/PDL/Transform.pd"
2018				if(ref $y eq 'HASH');
2019
2020				t_linear(pre=>$pre,$y,@_);
2021				}
2022
2023				#line 2870 "lib/PDL/Transform.pd"
2024
2025				=head2 t_rot
2026
2027				=for usage
2028
2029				$f = t_rot(\@rotation_in_degrees)
2030
2031				=for ref
2032
2033				Convenience interface to L.
2034
2035				t_rot produces a rotation transform in 2-D (scalar), 3-D (3-vector), or
2036				N-D (matrix). It acts exactly the same as Cshift)>.
2037
2038				=cut
2039
2040				*t_rot = \&t_rotate;
2041				sub t_rotate {
2042				my $rot = shift;
2043				my($y) = shift;
2044				return t_linear(rot=>$rot,%{$y})
2045				#line 2892 "lib/PDL/Transform.pd"
2046				if(ref $y eq 'HASH');
2047
2048				t_linear(rot=>$rot,$y,@_);
2049				}
2050
2051				#line 2902 "lib/PDL/Transform.pd"
2052
2053				=head2 t_fits
2054
2055				=for usage
2056
2057				$f = t_fits($fits,[option]);
2058
2059				=for ref
2060
2061				FITS pixel-to-scientific transformation with inverse
2062
2063				You feed in a hash ref or a PDL with one of those as a header, and you
2064				get back a transform that converts 0-originated, pixel-centered
2065				coordinates into scientific coordinates via the transformation in the
2066				FITS header. For most FITS headers, the transform is reversible, so
2067				applying the inverse goes the other way. This is just an instance
2068				of PDL::Transform::Linear, but with unit/type support
2069				using the FITS header you supply.
2070
2071				For now, this transform is rather limited -- it really ought to
2072				accept units differences and stuff like that, but they are just
2073				ignored for now. Probably that would require putting units into
2074				the whole transform framework.
2075
2076				This transform implements the linear transform part of the WCS FITS
2077				standard outlined in Greisen & Calabata 2002 (A&A in press; find it at
2078				L).
2079
2080				As a special case, you can pass in the boolean option "ignore_rgb"
2081				(default 0), and if you pass in a 3-D FITS header in which the last
2082				dimension has exactly 3 elements, it will be ignored in the output
2083				transformation. That turns out to be handy for handling rgb images.
2084
2085				=cut
2086
2087				sub t_fits {
2088				my($hdr) = shift;
2089				my($opt) = shift;
2090
2091				if(ref $opt ne 'HASH') {
2092				$opt = defined $opt ? {$opt,@_} : {} ;
2093				}
2094
2095				$hdr = $hdr->gethdr
2096				if(defined $hdr && UNIVERSAL::isa($hdr,'PDL'));
2097
2098				croak('t_fits requires a FITS header hash\n')
2099				if(!defined $hdr \|\| ref $hdr ne 'HASH' \|\| !defined($hdr->{NAXIS}));
2100
2101				my($n) = $hdr->{NAXIS}; $n = $n->at(0) if(UNIVERSAL::isa($n,'PDL'));
2102
2103				$n = 2
2104				if($opt->{ignore_rgb} && $n==3 && $hdr->{NAXIS3} == 3);
2105
2106				my($matrix) = PDL->zeroes($n,$n);
2107				my($pre) = PDL->zeroes($n);
2108				my($post) = PDL->zeroes($n);
2109
2110				##############################
2111				# Scaling: Use CDi_j formalism if present; otherwise use the
2112				# older PCi_j + CDELTi formalism.
2113
2114				my(@hgrab);
2115
2116				if(@hgrab = grep(m/^CD\d{1,3}_\d{1,3}$/,sort keys %$hdr)) { # assignment
2117				#
2118				# CDi_j formalism
2119				#
2120				for my $h(@hgrab) {
2121				$h =~ m/CD(\d{1,3})_(\d{1,3})/; # Should always match
2122				$matrix->slice("(".($1-1)."),(".($2-1).")") .= $hdr->{$h};
2123				}
2124				print "t_fits: Detected CDi_j matrix: \n",$matrix,"\n"
2125				if($PDL::Transform::debug);
2126
2127				} else {
2128
2129				#
2130				# PCi_j + CDELTi formalism
2131				# If PCi_j aren't present, and N=2, then try using CROTA or
2132				# CROTA2 to generate a rotation matrix instead.
2133				#
2134
2135				my($cdm) = PDL->zeroes($n,$n);
2136				my($cd) = $cdm->diagonal(0,1);
2137
2138				my($cpm) = PDL->zeroes($n,$n);
2139				$cpm->diagonal(0,1) .= 1; # PC: diagonal defaults to unity
2140				$cd .= 1;
2141
2142				if( @hgrab = grep(m/^PC\d{1,3}_\d{1,3}$/,sort keys %$hdr) ) { # assignment
2143
2144				for my $h(@hgrab) {
2145				$h =~ m/PC(\d{1,3})_(\d{1,3})$/ \|\| die "t_fits - match failed! This should never happen!";
2146				$cpm->slice("(".($1-1)."),(".($2-1).")") .= $hdr->{$h};
2147				}
2148				print "t_fits: Detected PCi_j matrix: \n",$cpm,"\n"
2149				if($PDL::Transform::debug && @hgrab);
2150
2151				} elsif ($n==2 && grep defined $hdr->{$_}, qw(CROTA CROTA1 CROTA2)) {
2152
2153				## CROTA is deprecated; CROTA1 was used for a while but is unofficial;
2154				## CROTA2 is encouraged instead.
2155				my $cr = $hdr->{CROTA2} // $hdr->{CROTA} // $hdr->{CROTA1};
2156
2157				$cr *= $PDL::Transform::DEG2RAD;
2158				# Rotation matrix rotates counterclockwise to get from sci to pixel coords
2159				# (detector has been rotated ccw, according to FITS standard)
2160				$cpm .= pdl( [cos($cr), sin($cr)],[-sin($cr),cos($cr)] );
2161
2162				}
2163
2164				for my $i(1..$n) {
2165				$cd->slice("(".($i-1).")") .= $hdr->{"CDELT$i"};
2166				}
2167				#If there are no CDELTs, then we assume they are all 1.0
2168				$cd+=1 if (all($cd==0));
2169
2170				$matrix = $cdm x $cpm;
2171				}
2172
2173				my($i1) = 0;
2174				for my $i(1..$n) {
2175				$pre->slice("($i1)") .= 1 - $hdr->{"CRPIX$i"};
2176				$post->slice("($i1)") .= $hdr->{"CRVAL$i"};
2177				$i1++;
2178				}
2179
2180				my $me = PDL::Transform::Linear->new({
2181				pre=>$pre, post=>$post, matrix=>$matrix
2182				});
2183				$me->{name} = 'FITS';
2184
2185				my (@otype,@ounit,@itype,@iunit);
2186				my @names = qw(X Y Z);
2187				for my $i(1..$hdr->{NAXIS}) {
2188				push(@otype,$hdr->{"CTYPE$i"});
2189				push(@ounit,$hdr->{"CUNIT$i"});
2190				push(@itype,"Image ". ( ($i-1<=$#names) ? $names[$i-1] : "${i}th dim" ));
2191				push(@iunit,"Pixels");
2192				}
2193				$me->{otype} = \@otype;
2194				$me->{itype} = \@itype;
2195				$me->{ounit} = \@ounit;
2196				$me->{iunit} = \@iunit;
2197
2198				# Check for nonlinear projection...
2199				# if($hdr->{CTYPE1} =~ m/(\w\w\w\w)\-(\w\w\w)/) {
2200				# print "Nonlinear transformation found... ignoring nonlinear part...\n";
2201				# }
2202
2203				return $me;
2204				}
2205
2206				#line 3065 "lib/PDL/Transform.pd"
2207
2208				=head2 t_code
2209
2210				=for usage
2211
2212				$f = t_code(,[],[options]);
2213
2214				=for ref
2215
2216				Transform implementing arbitrary perl code.
2217
2218				This is a way of getting quick-and-dirty new transforms. You pass in
2219				anonymous (or otherwise) code refs pointing to subroutines that
2220				implement the forward and, optionally, inverse transforms. The
2221				subroutines should accept a data PDL followed by a parameter hash ref,
2222				and return the transformed data PDL. The parameter hash ref can be
2223				set via the options, if you want to.
2224
2225				Options that are accepted are:
2226
2227				=over 3
2228
2229				=item p,params
2230
2231				The parameter hash that will be passed back to your code (defaults to the
2232				empty hash).
2233
2234				=item n,name
2235
2236				The name of the transform (defaults to "code").
2237
2238				=item i, idim (default 2)
2239
2240				The number of input dimensions (additional ones should be passed through
2241				unchanged)
2242
2243				=item o, odim (default 2)
2244
2245				The number of output dimensions
2246
2247				=item itype
2248
2249				The type of the input dimensions, in an array ref (optional and advisiory)
2250
2251				=item otype
2252
2253				The type of the output dimension, in an array ref (optional and advisory)
2254
2255				=item iunit
2256
2257				The units that are expected for the input dimensions (optional and advisory)
2258
2259				=item ounit
2260
2261				The units that are returned in the output (optional and advisory).
2262
2263				=back
2264
2265				The code variables are executable perl code, either as a code ref or
2266				as a string that will be eval'ed to produce code refs. If you pass in
2267				a string, it gets eval'ed at call time to get a code ref. If it compiles
2268				OK but does not return a code ref, then it gets re-evaluated with "sub {
2269				... }" wrapped around it, to get a code ref.
2270
2271				Note that code callbacks like this can be used to do really weird
2272				things and generate equally weird results -- caveat scriptor!
2273
2274				=cut
2275
2276				sub t_code {
2277				my($class) = 'PDL::Transform';
2278				my($func, $inv, $o) = @_;
2279				if(ref $inv eq 'HASH') {
2280				$o = $inv;
2281				$inv = undef;
2282				}
2283
2284				my($me) = $class->new;
2285				$me->{name} = _opt($o,['n','name','Name']) \|\| "code";
2286				$me->{func} = $func;
2287				$me->{inv} = $inv;
2288				$me->{params} = _opt($o,['p','params','Params']) \|\| {};
2289				$me->{idim} = _opt($o,['i','idim']) \|\| 2;
2290				$me->{odim} = _opt($o,['o','odim']) \|\| 2;
2291				$me->{itype} = _opt($o,['itype']) \|\| [];
2292				$me->{otype} = _opt($o,['otype']) \|\| [];
2293				$me->{iunit} = _opt($o,['iunit']) \|\| [];
2294				$me->{ounit} = _opt($o,['ounit']) \|\| [];
2295
2296				$me;
2297				}
2298
2299				#line 3164 "lib/PDL/Transform.pd"
2300
2301				=head2 t_cylindrical
2302
2303				C is an alias for C
2304
2305				=head2 t_radial
2306
2307				=for usage
2308
2309				$f = t_radial();
2310
2311				=for ref
2312
2313				Convert Cartesian to radial/cylindrical coordinates. (2-D/3-D; with inverse)
2314
2315				Converts 2-D Cartesian to radial (theta,r) coordinates. You can choose
2316				direct or conformal conversion. Direct conversion preserves radial
2317				distance from the origin; conformal conversion preserves local angles,
2318				so that each small-enough part of the image only appears to be scaled
2319				and rotated, not stretched. Conformal conversion puts the radius on a
2320				logarithmic scale, so that scaling of the original image plane is
2321				equivalent to a simple offset of the transformed image plane.
2322
2323				If you use three or more dimensions, the higher dimensions are ignored,
2324				yielding a conversion from Cartesian to cylindrical coordinates, which
2325				is why there are two aliases for the same transform. If you use higher
2326				dimensionality than 2, you must manually specify the origin or you will
2327				get dimension mismatch errors when you apply the transform.
2328
2329				Theta runs B instead of the more usual counterclockwise; that is
2330				to preserve the mirror sense of small structures.
2331
2332				OPTIONS:
2333
2334				=over 3
2335
2336				=item d, direct, Direct
2337
2338				Generate (theta,r) coordinates out (this is the default); incompatible
2339				with Conformal. Theta is in radians, and the radial coordinate is
2340				in the units of distance in the input plane.
2341
2342				=item r0, c, conformal, Conformal
2343
2344				If defined, this floating-point value causes t_radial to generate
2345				(theta, ln(r/r0)) coordinates out. Theta is in radians, and the
2346				radial coordinate varies by 1 for each e-folding of the r0-scaled
2347				distance from the input origin. The logarithmic scaling is useful for
2348				viewing both large and small things at the same time, and for keeping
2349				shapes of small things preserved in the image.
2350
2351				=item o, origin, Origin [default (0,0,0)]
2352
2353				This is the origin of the expansion. Pass in a PDL or an array ref.
2354
2355				=item u, unit, Unit [default 'radians']
2356
2357				This is the angular unit to be used for the azimuth.
2358
2359				=back
2360
2361				EXAMPLES
2362
2363				These examples do transformations back into the same size image as they
2364				started from; by suitable use of the "transform" option to
2365				L you can send them to any size array you like.
2366
2367				Examine radial structure in M51:
2368				Here, we scale the output to stretch 2*pi radians out to the
2369				full image width in the horizontal direction, and to stretch 1 radius out
2370				to a diameter in the vertical direction.
2371
2372				$x = rfits('m51.fits');
2373				$ts = t_linear(s => [250/2.0/3.14159, 2]); # Scale to fill orig. image
2374				$tu = t_radial(o => [130,130]); # Expand around galactic core
2375				$y = $x->map($ts x $tu);
2376
2377				Examine radial structure in M51 (conformal):
2378				Here, we scale the output to stretch 2*pi radians out to the full image width
2379				in the horizontal direction, and scale the vertical direction by the exact
2380				same amount to preserve conformality of the operation. Notice that
2381				each piece of the image looks "natural" -- only scaled and not stretched.
2382
2383				$x = rfits('m51.fits')
2384				$ts = t_linear(s=> 250/2.0/3.14159); # Note scalar (heh) scale.
2385				$tu = t_radial(o=> [130,130], r0=>5); # 5 pix. radius -> bottom of image
2386				$y = $ts->compose($tu)->unmap($x);
2387
2388				=cut
2389
2390				*t_cylindrical = \&t_radial;
2391				sub t_radial {
2392				my($class) = 'PDL::Transform';
2393				my($o) = $_[0];
2394				if(ref $o ne 'HASH') {
2395				$o = { @_ };
2396				}
2397
2398				my($me) = $class->new;
2399
2400				$me->{params}{origin} = _opt($o,['o','origin','Origin']);
2401				$me->{params}{origin} = pdl(0,0)
2402				unless defined($me->{params}{origin});
2403				$me->{params}{origin} = PDL->pdl($me->{params}{origin});
2404
2405				$me->{params}{r0} = _opt($o,['r0','R0','c','conformal','Conformal']);
2406				$me->{params}{origin} = PDL->pdl($me->{params}{origin});
2407
2408				$me->{params}{u} = _opt($o,['u','unit','Unit'],'radians');
2409				### Replace this kludge with a units call
2410				$me->{params}{angunit} = ($me->{params}{u} =~ m/^d/i) ? $PDL::Transform::RAD2DEG : 1.0;
2411				print "radial: conversion is $me->{params}{angunit}\n" if($PDL::Transform::debug);
2412
2413				$me->{name} = "radial (direct)";
2414
2415				$me->{idim} = 2;
2416				$me->{odim} = 2;
2417
2418				if($me->{params}{r0}) {
2419				$me->{otype} = ["Azimuth", "Ln radius" . ($me->{params}->{r0} != 1.0 ? "/$me->{params}->{r0}" : "")];
2420				$me->{ounit} = [$me->{params}->{u},'']; # true-but-null prevents copying
2421				} else {
2422				$me->{otype} = ["Azimuth","Radius"];
2423				$me->{ounit} = [$me->{params}->{u},'']; # false value copies prev. unit
2424				}
2425
2426				$me->{func} = sub {
2427
2428				my($data,$o) = @_;
2429
2430				my($out) = ($data->new_or_inplace);
2431
2432				my($d) = $data->copy;
2433				$d->slice("0:1") -= $o->{origin};
2434
2435				my($d0) = $d->slice("(0)");
2436				my($d1) = $d->slice("(1)");
2437
2438				# (mod operator on atan2 puts everything in the interval [0,2*PI).)
2439				$out->slice("(0)") .= (atan2(-$d1,$d0) % (2$PDL::Transform::PI)) $me->{params}->{angunit};
2440
2441				$out->slice("(1)") .= (defined $o->{r0}) ?
2442				0.5 * log( ($d1$d1 + $d0 $d0) / ($o->{r0} * $o->{r0}) ) :
2443				sqrt($d1$d1 + $d0$d0);
2444
2445				$out;
2446				};
2447
2448				$me->{inv} = sub {
2449
2450				my($d,$o) = @_;
2451				my($d0,$d1,$out)=
2452				( ($d->is_inplace) ?
2453				($d->slice("(0)")->copy, $d->slice("(1)")->copy->dummy(0,2), $d) :
2454				($d->slice("(0)"), $d->slice("(1)")->dummy(0,2), $d->copy)
2455				);
2456
2457				$d0 /= $me->{params}->{angunit};
2458
2459				my($os) = $out->slice("0:1");
2460				$os .= append(cos($d0)->dummy(0,1),-sin($d0)->dummy(0,1));
2461				$os = defined $o->{r0} ? ($o->{r0} exp($d1)) : $d1;
2462				$os += $o->{origin};
2463
2464				$out;
2465				};
2466
2467				$me;
2468				}
2469
2470				#line 3341 "lib/PDL/Transform.pd"
2471
2472				=head2 t_quadratic
2473
2474				=for usage
2475
2476				$t = t_quadratic();
2477
2478				=for ref
2479
2480				Quadratic scaling -- cylindrical pincushion (n-d; with inverse)
2481
2482				Quadratic scaling emulates pincushion in a cylindrical optical system:
2483				separate quadratic scaling is applied to each axis. You can apply
2484				separate distortion along any of the principal axes. If you want
2485				different axes, use L and L to rotate
2486				them to the correct angle. The scaling options may be scalars or
2487				vectors; if they are scalars then the expansion is isotropic.
2488
2489				The formula for the expansion is:
2490
2491				f(a) = ( + * a^2/ ) / (abs() + 1)
2492
2493				where is a scaling coefficient and is a fundamental
2494				length scale. Negative values of result in a pincushion
2495				contraction.
2496
2497				Note that, because quadratic scaling does not have a strict inverse for
2498				coordinate systems that cross the origin, we cheat slightly and use
2499				$x * abs($x) rather than $x**2. This does the Right thing for pincushion
2500				and barrel distortion, but means that t_quadratic does not behave exactly
2501				like t_cubic with a null cubic strength coefficient.
2502
2503				OPTIONS
2504
2505				=over 3
2506
2507				=item o,origin,Origin
2508
2509				The origin of the pincushion. (default is the, er, origin).
2510
2511				=item l,l0,length,Length,r0
2512
2513				The fundamental scale of the transformation -- the radius that remains
2514				unchanged. (default=1)
2515
2516				=item s,str,strength,Strength
2517
2518				The relative strength of the pincushion. (default = 0.1)
2519
2520				=item honest (default=0)
2521
2522				Sets whether this is a true quadratic coordinate transform. The more
2523				common form is pincushion or cylindrical distortion, which switches
2524				branches of the square root at the origin (for symmetric expansion).
2525				Setting honest to a non-false value forces true quadratic behavior,
2526				which is not mirror-symmetric about the origin.
2527
2528				=item d, dim, dims, Dims
2529
2530				The number of dimensions to quadratically scale (default is the
2531				dimensionality of your input vectors)
2532
2533				=back
2534
2535				=cut
2536
2537				sub t_quadratic {
2538				my($class) = 'PDL::Transform';
2539				my($o) = $_[0];
2540				if(ref $o ne 'HASH') {
2541				$o = {@_};
2542				}
2543				my($me) = $class->new;
2544
2545				$me->{params}->{origin} = _opt($o,['o','origin','Origin'],pdl(0));
2546				$me->{params}->{l0} = _opt($o,['r0','l','l0','length','Length'],pdl(1));
2547				$me->{params}->{str} = _opt($o,['s','str','strength','Strength'],pdl(0.1));
2548				$me->{params}->{dim} = _opt($o,['d','dim','dims','Dims']);
2549				$me->{name} = "quadratic";
2550
2551				$me->{func} = sub {
2552				my($data,$o) = @_;
2553				my($dd) = $data->copy - $o->{origin};
2554				my($d) = (defined $o->{dim}) ? $dd->slice("0:".($o->{dim}-1)) : $dd;
2555				$d += $o->{str} * ($d * abs($d)) / $o->{l0};
2556				$d /= (abs($o->{str}) + 1);
2557				$d += $o->{origin};
2558				if($data->is_inplace) {
2559				$data .= $dd;
2560				return $data;
2561				}
2562				$dd;
2563				};
2564
2565				$me->{inv} = sub {
2566				my($data,$opt) = @_;
2567				my($dd) = $data->copy ;
2568				my($d) = (defined $opt->{dim}) ? $dd->slice("0:".($opt->{dim}-1)) : $dd;
2569				my($o) = $opt->{origin};
2570				my($s) = $opt->{str};
2571				my($l) = $opt->{l0};
2572
2573				$d .= ((-1 + sqrt(1 + 4 * $s/$l * abs($d-$o) * (1+abs($s))))
2574				/ 2 / $s * $l) * (1 - 2*($d < $o));
2575				$d += $o;
2576				if($data->is_inplace) {
2577				$data .= $dd;
2578				return $data;
2579				}
2580				$dd;
2581				};
2582				$me;
2583				}
2584
2585				#line 3460 "lib/PDL/Transform.pd"
2586
2587				=head2 t_cubic
2588
2589				=for usage
2590
2591				$t = t_cubic();
2592
2593				=for ref
2594
2595				Cubic scaling - cubic pincushion (n-d; with inverse)
2596
2597				Cubic scaling is a generalization of t_quadratic to a purely
2598				cubic expansion.
2599
2600				The formula for the expansion is:
2601
2602				f(a) = ( a' + st * a'^3/L_0^2 ) / (1 + abs(st)) + origin
2603
2604				where a'=(a-origin). That is a simple pincushion
2605				expansion/contraction that is fixed at a distance of L_0 from the
2606				origin.
2607
2608				Because there is no quadratic term the result is always invertible with
2609				one real root, and there is no mucking about with complex numbers or
2610				multivalued solutions.
2611
2612				OPTIONS
2613
2614				=over 3
2615
2616				=item o,origin,Origin
2617
2618				The origin of the pincushion. (default is the, er, origin).
2619
2620				=item l,l0,length,Length,r0
2621
2622				The fundamental scale of the transformation -- the radius that remains
2623				unchanged. (default=1)
2624
2625				=item d, dim, dims, Dims
2626
2627				The number of dimensions to treat (default is the
2628				dimensionality of your input vectors)
2629
2630				=back
2631
2632				=cut
2633
2634				sub t_cubic {
2635				my ($class) = 'PDL::Transform';
2636				my($o) = $_[0];
2637				if(ref $o ne 'HASH') {
2638				$o = {@_};
2639				}
2640				my($me) = $class->new;
2641
2642				$me->{params}->{dim} = _opt($o,['d','dim','dims','Dims'],undef);
2643				$me->{params}->{origin} = _opt($o,['o','origin','Origin'],pdl(0));
2644				$me->{params}->{l0} = _opt($o,['r0','l','l0','length','Length'],pdl(1));
2645				$me->{params}->{st} = _opt($o,['s','st','str'],pdl(0));
2646				$me->{name} = "cubic";
2647
2648				$me->{params}->{cuberoot} = sub {
2649				my $x = shift;
2650				my $as = 1 - 2*($x<0);
2651				return $as * ( abs($x) ** (1/3) );
2652				};
2653
2654				$me->{func} = sub {
2655				my($data,$o) = @_;
2656				my($dd) = $data->copy;
2657				my($d) = (defined $o->{dim}) ? $dd->slice("0:".($o->{dim}-1)) : $dd;
2658
2659				$d -= $o->{origin};
2660
2661				my $dl0 = $d / $o->{l0};
2662
2663				# f(x) = x + x*3 ($o->{st} / $o->{l0}**2):
2664				# A = ($o->{st}/$dl0**2)
2665				# B = 0
2666				# C = 1
2667				# D = -f(x)
2668				$d += $o->{st} * $d * $dl0 * $dl0;
2669				$d /= ($o->{st}**2 + 1);
2670
2671				$d += $o->{origin};
2672
2673				if($data->is_inplace) {
2674				$data .= $dd;
2675				return $data;
2676				}
2677				return $dd;
2678				};
2679
2680				$me->{inv} = sub {
2681				my($data,$o) = @_;
2682				my($l) = $o->{l0};
2683
2684				my($dd) = $data->copy;
2685				my($d) = (defined $o->{dim}) ? $dd->slice("0:".($o->{dim}-1)) : $dd;
2686
2687				$d -= $o->{origin};
2688				$d *= ($o->{st}+1);
2689
2690				# Now we have to solve:
2691				# A x^3 + B X^2 + C x + D dd = 0
2692				# with the assignments above for A,B,C,D.
2693				# Since B is zero, this is greatly simplified - the discriminant is always negative,
2694				# so there is always exactly one real root.
2695				#
2696				# The only real hassle is creating a symmetric cube root; for convenience
2697				# is stashed in the params hash.
2698
2699				# First: create coefficients for mnemonics.
2700				my ($A, $C, $D) = ( $o->{st} / $l / $l, 1, - $d );
2701				my $alpha = 27 * $A * $A * $D;
2702				my $beta = 3 * $A * $C;
2703
2704				my $inner_root = sqrt( $alpha * $alpha + 4 * $beta * $beta * $beta );
2705
2706				$d .= (-1 / (3 * $A)) *
2707				(
2708				+ &{$o->{cuberoot}}( 0.5 * ( $alpha + $inner_root ) )
2709				+ &{$o->{cuberoot}}( 0.5 * ( $alpha - $inner_root ) )
2710				);
2711
2712				$d += $me->{params}{origin};
2713
2714				if($data->is_inplace) {
2715				$data .= $dd;
2716				return $data;
2717				} else {
2718				return $dd;
2719				}
2720				};
2721
2722				$me;
2723				}
2724
2725				#line 3606 "lib/PDL/Transform.pd"
2726
2727				=head2 t_quartic
2728
2729				=for usage
2730
2731				$t = t_quartic();
2732
2733				=for ref
2734
2735				Quartic scaling -- cylindrical pincushion (n-d; with inverse)
2736
2737				Quartic scaling is a generalization of t_quadratic to a quartic
2738				expansion. Only even powers of the input coordinates are retained,
2739				and (as with t_quadratic) sign is preserved, making it an odd function
2740				although a true quartic transformation would be an even function.
2741
2742				You can apply separate distortion along any of the principal axes. If
2743				you want different axes, use L and
2744				L to rotate them to the correct angle. The
2745				scaling options may be scalars or vectors; if they are scalars then
2746				the expansion is isotropic.
2747
2748				The formula for the expansion is:
2749
2750				f(a) = ( + * a^2/ ) / (abs() + 1)
2751
2752				where is a scaling coefficient and is a fundamental
2753				length scale. Negative values of result in a pincushion
2754				contraction.
2755
2756				Note that, because quadratic scaling does not have a strict inverse for
2757				coordinate systems that cross the origin, we cheat slightly and use
2758				$x * abs($x) rather than $x**2. This does the Right thing for pincushion
2759				and barrel distortion, but means that t_quadratic does not behave exactly
2760				like t_cubic with a null cubic strength coefficient.
2761
2762				OPTIONS
2763
2764				=over 3
2765
2766				=item o,origin,Origin
2767
2768				The origin of the pincushion. (default is the, er, origin).
2769
2770				=item l,l0,length,Length,r0
2771
2772				The fundamental scale of the transformation -- the radius that remains
2773				unchanged. (default=1)
2774
2775				=item s,str,strength,Strength
2776
2777				The relative strength of the pincushion. (default = 0.1)
2778
2779				=item honest (default=0)
2780
2781				Sets whether this is a true quadratic coordinate transform. The more
2782				common form is pincushion or cylindrical distortion, which switches
2783				branches of the square root at the origin (for symmetric expansion).
2784				Setting honest to a non-false value forces true quadratic behavior,
2785				which is not mirror-symmetric about the origin.
2786
2787				=item d, dim, dims, Dims
2788
2789				The number of dimensions to quadratically scale (default is the
2790				dimensionality of your input vectors)
2791
2792				=back
2793
2794				=cut
2795
2796				sub t_quartic {
2797				my($class) = 'PDL::Transform';
2798				my($o) = $_[0];
2799				if(ref $o ne 'HASH') {
2800				$o = {@_};
2801				}
2802				my($me) = $class->new;
2803
2804				$me->{params}->{origin} = _opt($o,['o','origin','Origin'],pdl(0));
2805				$me->{params}->{l0} = _opt($o,['r0','l','l0','length','Length'],pdl(1));
2806				$me->{params}->{str} = _opt($o,['s','str','strength','Strength'],pdl(0.1));
2807				$me->{params}->{dim} = _opt($o,['d','dim','dims','Dims']);
2808				$me->{name} = "quadratic";
2809
2810				$me->{func} = sub {
2811				my($data,$o) = @_;
2812				my($dd) = $data->copy - $o->{origin};
2813				my($d) = (defined $o->{dim}) ? $dd->slice("0:".($o->{dim}-1)) : $dd;
2814				$d += $o->{str} * ($d * abs($d)) / $o->{l0};
2815				$d /= (abs($o->{str}) + 1);
2816				$d += $o->{origin};
2817				if($data->is_inplace) {
2818				$data .= $dd;
2819				return $data;
2820				}
2821				$dd;
2822				};
2823
2824				$me->{inv} = sub {
2825				my($data,$opt) = @_;
2826				my($dd) = $data->copy ;
2827				my($d) = (defined $opt->{dim}) ? $dd->slice("0:".($opt->{dim}-1)) : $dd;
2828				my($o) = $opt->{origin};
2829				my($s) = $opt->{str};
2830				my($l) = $opt->{l0};
2831
2832				$d .= ((-1 + sqrt(1 + 4 * $s/$l * abs($d-$o) * (1+abs($s))))
2833				/ 2 / $s * $l) * (1 - 2*($d < $o));
2834				$d += $o;
2835				if($data->is_inplace) {
2836				$data .= $dd;
2837				return $data;
2838				}
2839				$dd;
2840				};
2841				$me;
2842				}
2843
2844				#line 3729 "lib/PDL/Transform.pd"
2845
2846				=head2 t_spherical
2847
2848				=for usage
2849
2850				$t = t_spherical();
2851
2852				=for ref
2853
2854				Convert Cartesian to spherical coordinates. (3-D; with inverse)
2855
2856				Convert 3-D Cartesian to spherical (theta, phi, r) coordinates. Theta
2857				is longitude, centered on 0, and phi is latitude, also centered on 0.
2858				Unless you specify Euler angles, the pole points in the +Z direction
2859				and the prime meridian is in the +X direction. The default is for
2860				theta and phi to be in radians; you can select degrees if you want
2861				them.
2862
2863				Just as the L 2-D transform acts like a 3-D
2864				cylindrical transform by ignoring third and higher dimensions,
2865				Spherical acts like a hypercylindrical transform in four (or higher)
2866				dimensions. Also as with L, you must manually specify
2867				the origin if you want to use more dimensions than 3.
2868
2869				To deal with latitude & longitude on the surface of a sphere (rather
2870				than full 3-D coordinates), see
2871				L.
2872
2873				OPTIONS:
2874
2875				=over 3
2876
2877				=item o, origin, Origin [default (0,0,0)]
2878
2879				This is the Cartesian origin of the spherical expansion. Pass in a PDL
2880				or an array ref.
2881
2882				=item e, euler, Euler [default (0,0,0)]
2883
2884				This is a 3-vector containing Euler angles to change the angle of the
2885				pole and ordinate. The first two numbers are the (theta, phi) angles
2886				of the pole in a (+Z,+X) spherical expansion, and the last is the
2887				angle that the new prime meridian makes with the meridian of a simply
2888				tilted sphere. This is implemented by composing the output transform
2889				with a PDL::Transform::Linear object.
2890
2891				=item u, unit, Unit (default radians)
2892
2893				This option sets the angular unit to be used. Acceptable values are
2894				"degrees","radians", or reasonable substrings thereof (e.g. "deg", and
2895				"rad", but "d" and "r" are deprecated). Once genuine unit processing
2896				comes online (a la Math::Units) any angular unit should be OK.
2897
2898				=back
2899
2900				=cut
2901
2902				sub t_spherical {
2903				my($class) = 'PDL::Transform';
2904				my($o) = $_[0];
2905				if(ref $o ne 'HASH') {
2906				$o = { @_ } ;
2907				}
2908
2909				my($me) = $class->new;
2910
2911				$me->{idim}=3;
2912				$me->{odim}=3;
2913
2914				$me->{params}->{origin} = _opt($o,['o','origin','Origin']);
2915				$me->{params}->{origin} //= PDL->zeroes(3);
2916				$me->{params}->{origin} = PDL->pdl($me->{params}->{origin});
2917
2918				$me->{params}->{deg} = _opt($o,['d','degrees','Degrees']);
2919
2920				my $unit = _opt($o,['u','unit','Unit']);
2921				$me->{params}{angunit} = ($unit && $unit =~ m/^d/i) ?
2922				$PDL::Transform::DEG2RAD : undef;
2923
2924				$me->{name} = "spherical";
2925
2926				$me->{func} = sub {
2927				my($data,$o) = @_;
2928				my($d) = $data->copy - $o->{origin};
2929
2930				my($d0,$d1,$d2) = ($d->slice("(0)"),$d->slice("(1)"),$d->slice("(2)"));
2931				my($out) = ($d->is_inplace) ? $data : $data->copy;
2932
2933				$out->slice("(0)") .= PDL::atan2($d1, $d0);
2934				$out->slice("(2)") .= PDL::sqrt($d0$d0 + $d1$d1 + $d2*$d2);
2935				$out->slice("(1)") .= PDL::asin($d2 / $out->slice("(2)"));
2936
2937				$out->slice("0:1") /= $o->{angunit}
2938				if(defined $o->{angunit});
2939
2940				$out;
2941				};
2942
2943				$me->{inv} = sub {
2944				my($d,$o) = @_;
2945
2946				my($theta,$phi,$r,$out) =
2947				( ($d->is_inplace) ?
2948				($d->slice("(0)")->copy, $d->slice("(1)")->copy, $d->slice("(2)")->copy, $d) :
2949				($d->slice("(0)"), $d->slice("(1)"), $d->slice("(2)"), $d->copy)
2950				);
2951
2952				my($x,$y,$z) =
2953				($out->slice("(0)"),$out->slice("(1)"),$out->slice("(2)"));
2954
2955				my($ph,$th);
2956				if(defined $o->{angunit}){
2957				$ph = $o->{angunit} * $phi;
2958				$th = $o->{angunit} * $theta;
2959				} else {
2960				$ph = $phi;
2961				$th = $theta;
2962				}
2963
2964				$z .= $r * sin($ph);
2965				$x .= $r * cos($ph);
2966				$y .= $x * sin($th);
2967				$x *= cos($th);
2968				$out += $o->{origin};
2969
2970				$out;
2971				};
2972
2973				$me;
2974				}
2975
2976				#line 3865 "lib/PDL/Transform.pd"
2977
2978				=head2 t_projective
2979
2980				=for usage
2981
2982				$t = t_projective();
2983
2984				=for ref
2985
2986				Projective transformation
2987
2988				Projective transforms are simple quadratic, quasi-linear
2989				transformations. They are the simplest transformation that
2990				can continuously warp an image plane so that four arbitrarily chosen
2991				points exactly map to four other arbitrarily chosen points. They
2992				have the property that straight lines remain straight after transformation.
2993
2994				You can specify your projective transformation directly in homogeneous
2995				coordinates, or (in 2 dimensions only) as a set of four unique points that
2996				are mapped one to the other by the transformation.
2997
2998				Projective transforms are quasi-linear because they are most easily
2999				described as a linear transformation in homogeneous coordinates
3000				(e.g. (x',y',w) where w is a normalization factor: x = x'/w, etc.).
3001				In those coordinates, an N-D projective transformation is represented
3002				as simple multiplication of an N+1-vector by an N+1 x N+1 matrix,
3003				whose lower-right corner value is 1.
3004
3005				If the bottom row of the matrix consists of all zeroes, then the
3006				transformation reduces to a linear affine transformation (as in
3007				L).
3008
3009				If the bottom row of the matrix contains nonzero elements, then the
3010				transformed x,y,z,etc. coordinates are related to the original coordinates
3011				by a quadratic polynomial, because the normalization factor 'w' allows
3012				a second factor of x,y, and/or z to enter the equations.
3013
3014				OPTIONS:
3015
3016				=over 3
3017
3018				=item m, mat, matrix, Matrix
3019
3020				If specified, this is the homogeneous-coordinate matrix to use. It must
3021				be N+1 x N+1, for an N-dimensional transformation.
3022
3023				=item p, point, points, Points
3024
3025				If specified, this is the set of four points that should be mapped one to the other.
3026				The homogeneous-coordinate matrix is calculated from them. You should feed in a
3027				2x2x4 PDL, where the 0 dimension runs over coordinate, the 1 dimension runs between input
3028				and output, and the 2 dimension runs over point. For example, specifying
3029
3030				p=> pdl([ [[0,1],[0,1]], [[5,9],[5,8]], [[9,4],[9,3]], [[0,0],[0,0]] ])
3031
3032				maps the origin and the point (0,1) to themselves, the point (5,9) to (5,8), and
3033				the point (9,4) to (9,3).
3034
3035				This is similar to the behavior of fitwarp2d with a quadratic polynomial.
3036
3037				=back
3038
3039				=cut
3040
3041				sub t_projective {
3042				my($class) = 'PDL::Transform';
3043				my($o) = $_[0];
3044				if(ref $o ne 'HASH') {
3045				$o = { @_ };
3046				}
3047
3048				my($me) = $class->new;
3049
3050				$me->{name} = "projective";
3051
3052				### Set options...
3053
3054				$me->{params}->{idim} = $me->{idim} = _opt($o,['d','dim','Dim']);
3055
3056				my $matrix;
3057				if(defined ($matrix=_opt($o,['m','matrix','Matrix']))) {
3058				$matrix = pdl($matrix);
3059				die "t_projective: needs a square matrix"
3060				if($matrix->dims != 2 \|\| $matrix->dim(0) != $matrix->dim(1));
3061
3062				$me->{params}->{idim} = $matrix->dim(0)-1
3063				unless(defined($me->{params}->{idim}));
3064
3065				$me->{idim} = $me->{params}->{idim};
3066
3067				die "t_projective: matrix not compatible with given dimension (should be N+1xN+1)\n"
3068				unless($me->{params}->{idim}==$matrix->dim(0)-1);
3069
3070				my $inv = $matrix->inv;
3071				print STDERR "t_projective: warning - received non-invertible matrix\n"
3072				unless(all($inv*0 == 0));
3073
3074				$me->{params}->{matrix} = $matrix;
3075				$me->{params}->{matinv} = $inv;
3076
3077				} elsif(defined (my $p=pdl(_opt($o,['p','point','points','Point','Points'])))) {
3078				die "t_projective: points array should be 2(x,y) x 2(in,out) x 4(point)\n\t(only 2-D points spec is available just now, sorry)\n"
3079				unless($p->dims==3 && all(pdl($p->dims)==pdl(2,2,4)));
3080
3081				# Solve the system of N equations to find the projective transform
3082				my ($p0,$p1,$p2,$p3) = ( $p->slice(":,(0),(0)"), $p->slice(":,(0),(1)"), $p->slice(":,(0),(2)"), $p->slice(":,(0),(3)") );
3083				my ($P0,$P1,$P2,$P3) = ( $p->slice(":,(1),(0)"), $p->slice(":,(1),(1)"), $p->slice(":,(1),(2)"), $p->slice(":,(1),(3)") );
3084				#print "declaring PDL...\n" ;
3085				my $M = pdl( [ [$p0->at(0), $p0->at(1), 1, 0, 0, 0, -$P0->at(0)$p0->at(0), -$P0->at(0)$p0->at(1)],
3086				[ 0, 0, 0, $p0->at(0), $p0->at(1), 1, -$P0->at(1)$p0->at(0), -$P0->at(1)$p0->at(1)],
3087				[$p1->at(0), $p1->at(1), 1, 0, 0, 0, -$P1->at(0)$p1->at(0), -$P1->at(0)$p1->at(1)],
3088				[ 0, 0, 0, $p1->at(0), $p1->at(1), 1, -$P1->at(1)$p1->at(0), -$P1->at(1)$p1->at(1)],
3089				[$p2->at(0), $p2->at(1), 1, 0, 0, 0, -$P2->at(0)$p2->at(0), -$P2->at(0)$p2->at(1)],
3090				[ 0, 0, 0, $p2->at(0), $p2->at(1), 1, -$P2->at(1)$p2->at(0), -$P2->at(1)$p2->at(1)],
3091				[$p3->at(0), $p3->at(1), 1, 0, 0, 0, -$P3->at(0)$p3->at(0), -$P3->at(0)$p3->at(1)],
3092				[ 0, 0, 0, $p3->at(0), $p3->at(1), 1, -$P3->at(1)$p3->at(0), -$P3->at(1)$p3->at(1)]
3093				] );
3094				#print "ok. Finding inverse...\n";
3095				my $h = ($M->inv x $p->slice(":,(1),:")->flat->slice("*1"))->slice("(0)");
3096				# print "ok\n";
3097				my $matrix = ones(3,3);
3098				$matrix->flat->slice("0:7") .= $h;
3099				$me->{params}->{matrix} = $matrix;
3100
3101				$me->{params}->{matinv} = $matrix->inv;
3102				}
3103
3104				$me->{params}->{idim} = 2 unless defined $me->{params}->{idim};
3105				$me->{params}->{odim} = $me->{params}->{idim};
3106				$me->{idim} = $me->{params}->{idim};
3107				$me->{odim} = $me->{params}->{odim};
3108
3109				$me->{func} = sub {
3110				my($data,$o) = @_;
3111				my($id) = $data->dim(0);
3112				my($d) = $data->glue(0,ones($data->slice("0")));
3113				my($out) = ($o->{matrix} x $d->slice("*1"))->slice("(0)");
3114				return ($out->slice("0:".($id-1))/$out->slice("$id"));
3115				};
3116
3117				$me->{inv} = sub {
3118				my($data,$o) = @_;
3119				my($id) = $data->dim(0);
3120				my($d) = $data->glue(0,ones($data->slice("0")));
3121				my($out) = ($o->{matinv} x $d->slice("*1"))->slice("(0)");
3122				return ($out->slice("0:".($id-1))/$out->slice("$id"));
3123				};
3124
3125				$me;
3126				}
3127
3128				#line 245 "lib/PDL/Transform.pd"
3129
3130				=head1 AUTHOR
3131
3132				Copyright 2002, 2003 Craig DeForest. There is no warranty. You are allowed
3133				to redistribute this software under certain conditions. For details,
3134				see the file COPYING in the PDL distribution. If this file is
3135				separated from the PDL distribution, the copyright notice should be
3136				included in the file.
3137
3138				=cut
3139
3140				package PDL::Transform;
3141				use Carp;
3142				use overload '""' => \&_strval;
3143				use overload 'x' => \&_compose_op;
3144				use overload '**' => \&_pow_op;
3145				use overload '!' => \&t_inverse;
3146
3147				use PDL::LiteF;
3148				use PDL::MatrixOps;
3149
3150				our $PI = 3.1415926535897932384626;
3151				our $DEG2RAD = $PI / 180;
3152				our $RAD2DEG = 180/$PI;
3153				our $E = exp(1);
3154
3155				#### little helper kludge parses a list of synonyms
3156				sub _opt {
3157				my($hash) = shift;
3158				my($synonyms) = shift;
3159				my($alt) = shift; # default is undef -- ok.
3160				local($_);
3161				foreach $_(@$synonyms){
3162				return (UNIVERSAL::isa($alt,'PDL')) ? PDL->pdl($hash->{$_}) : $hash->{$_}
3163				if defined($hash->{$_}) ;
3164				}
3165				return $alt;
3166				}
3167
3168				######################################################################
3169				#
3170				# Stringification hack. _strval just does a method search on stringify
3171				# for the object itself. This gets around the fact that stringification
3172				# overload is a subroutine call, not a method search.
3173				#
3174
3175				sub _strval {
3176				my($me) = shift;
3177				$me->stringify();
3178				}
3179
3180				######################################################################
3181				#
3182				# PDL::Transform overall stringifier. Subclassed stringifiers should
3183				# call this routine first then append auxiliary information.
3184				#
3185				sub stringify {
3186				my($me) = shift;
3187				my($mestr) = (ref $me);
3188				$mestr =~ s/PDL::Transform:://;
3189				my $out = $mestr . " (" . $me->{name} . "): ";
3190				$out .= "fwd ". ((defined ($me->{func})) ? ( (ref($me->{func}) eq 'CODE') ? "ok" : "non-CODE(!!)" ): "missing")."; ";
3191				$out .= "inv ". ((defined ($me->{inv})) ? ( (ref($me->{inv}) eq 'CODE') ? "ok" : "non-CODE(!!)" ):"missing").".\n";
3192				}
3193				#line 3194 "lib/PDL/Transform.pm"
3194
3195				# Exit with OK status
3196
3197				1;