"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.080/GENERATED/PDL/Transform.pm" (28 May 2022, 104325 Bytes) of package /linux/misc/PDL-2.080.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Perl source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Transform.pm" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.079_vs_2.080.

A hint: This file contains one or more very long lines, so maybe it is better readable using the pure text view mode that shows the contents as wrapped lines within the browser window.


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