"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.080/Libtmp/Transform/transform.pd" (7 May 2022, 130880 Bytes) of package /linux/misc/PDL-2.080.tar.gz:


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