File Coverage

blib/lib/PDLA/Transform.pm

Criterion	Covered	Total	%
statement	395	883	44.7
branch	141	402	35.0
condition	51	166	30.7
subroutine	32	66	48.4
pod	17	24	70.8
total	636	1541	41.2

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