File Coverage

lib/PDL/Transform.pd

Criterion	Covered	Total	%
statement	382	391	97.7
branch	400	2034	19.6
condition			n/a
subroutine			n/a
pod			n/a
total	782	2425	32.2

pp_addhdr("#include \n") unless $^O =~ /MSWin32/i;
pp_addpm({At=>'Top'},<<'+======EOD======');

=head1 NAME

PDL::Transform - Coordinate transforms, image warping, and N-D functions

=head1 SYNOPSIS

use PDL::Transform;

 my $t = PDL::Transform::->new()

 $out = $t->apply($in)  # Apply transform to some N-vectors (Transform method)
 $out = $in->apply($t)  # Apply transform to some N-vectors (PDL method)

 $im1 = $t->map($im);   # Transform image coordinates (Transform method)
 $im1 = $im->map($t);   # Transform image coordinates (PDL method)

 $t2 = $t->compose($t1);  # compose two transforms
 $t2 = $t x $t1;          # compose two transforms (by analogy to matrix mult.)

 $t3 = $t2->inverse();    # invert a transform
 $t3 = !$t2;              # invert a transform (by analogy to logical "not")

=head1 DESCRIPTION

PDL::Transform is a convenient way to represent coordinate
transformations and resample images.  It embodies functions mapping
R^N -> R^M, both with and without inverses.  Provision exists for
parametrizing functions, and for composing them.  You can use this
part of the Transform object to keep track of arbitrary functions
mapping R^N -> R^M with or without inverses.

The simplest way to use a Transform object is to transform vector
data between coordinate systems.  The L method
accepts a PDL whose 0th dimension is coordinate index (all other
dimensions are broadcasted over) and transforms the vectors into the new
coordinate system.

Transform also includes image resampling, via the L method.
You define a coordinate transform using a Transform object, then use
it to remap an image PDL.  The output is a remapped, resampled image.

You can define and compose several transformations, then apply them
all at once to an image.  The image is interpolated only once, when
all the composed transformations are applied.

In keeping with standard practice, but somewhat counterintuitively,
the L engine uses the inverse transform to map coordinates
FROM the destination dataspace (or image plane) TO the source dataspace;
hence PDL::Transform keeps track of both the forward and inverse transform.

For terseness and convenience, most of the constructors are exported
into the current package with the name C<< t_ >>, so the following
(for example) are synonyms:

  $t = PDL::Transform::Radial->new;  # Long way

  $t = t_radial();                    # Short way

Several math operators are overloaded, so that you can compose and
invert functions with expression syntax instead of method syntax (see below).

=head1 EXAMPLE

Coordinate transformations and mappings are a little counterintuitive
at first.  Here are some examples of transforms in action:

   use PDL::Transform;
   use PDL::Graphics::Simple;
   $x = rfits('m51.fits');   # Substitute path if necessary!
   $ts = t_linear(Scale=>3); # Scaling transform

   $w = pgswin(xs);
   $w->plot(with=>'image', $x);

   ## Grow m51 by a factor of 3; origin is at lower left.
   $y = $ts->map($x,{pix=>1});    # pix option uses direct pixel coord system
   $w->plot(with=>'image', $y);

   ## Shrink m51 by a factor of 3; origin still at lower left.
   $c = $ts->unmap($x, {pix=>1});
   $w->plot(with=>'image', $c);

   ## Grow m51 by a factor of 3; origin is at scientific origin.
   $d = $ts->map($x,$x->hdr);    # FITS hdr template prevents autoscaling
   $w->plot(with=>'image', $d);

   ## Shrink m51 by a factor of 3; origin is still at sci. origin.
   $e = $ts->unmap($x,$x->hdr);
   $w->plot(with=>'image', $e);

   ## A no-op: shrink m51 by a factor of 3, then autoscale back to size
   $f = $ts->map($x);            # No template causes autoscaling of output

=head1 OPERATOR OVERLOADS

=over 3

=item '!'

The bang is a unary inversion operator.  It binds exactly as
tightly as the normal bang operator.

=item 'x'

By analogy to matrix multiplication, 'x' is the compose operator, so these
two expressions are equivalent:

  $f->inverse()->compose($g)->compose($f) # long way
  !$f x $g x $f                           # short way

Both of those expressions are equivalent to the mathematical expression
f^-1 o g o f, or f^-1(g(f(x))).

=item '**'

By analogy to numeric powers, you can apply an operator a positive
integer number of times with the ** operator:

  $f->compose($f)->compose($f)  # long way
  $f**3                         # short way

=back

=head1 INTERNALS

Transforms are perl hashes.  Here's a list of the meaning of each key:

=over 3

=item func

Ref to a subroutine that evaluates the transformed coordinates.  It's
called with the input coordinate, and the "params" hash.  This
springboarding is done via explicit ref rather than by subclassing,
for convenience both in coding new transforms (just add the
appropriate sub to the module) and in adding custom transforms at
run-time. Note that, if possible, new Cs should support
L operation to save memory when the data are flagged
inplace.  But C should always return its result even when
flagged to compute in-place.

C should treat the 0th dimension of its input as a dimensional
index (running 0..N-1 for R^N operation) and broadcast over all other input
dimensions.

=item inv

Ref to an inverse method that reverses the transformation.  It must
accept the same "params" hash that the forward method accepts.  This
key can be left undefined in cases where there is no inverse.

=item idim, odim

Number of useful dimensions for indexing on the input and output sides
(ie the order of the 0th dimension of the coordinates to be fed in or
that come out).  If this is set to 0, then as many are allocated as needed.

=item name

A shorthand name for the transformation (convenient for debugging).
You should plan on using UNIVERAL::isa to identify classes of
transformation, e.g. all linear transformations should be subclasses
of PDL::Transform::Linear.  That makes it easier to add smarts to,
e.g., the compose() method.

=item itype

An array containing the name of the quantity that is expected from the
input ndarray for the transform, for each dimension.  This field is advisory,
and can be left blank if there's no obvious quantity associated with
the transform.  This is analogous to the CTYPEn field used in FITS headers.

=item oname

Same as itype, but reporting what quantity is delivered for each
dimension.

=item iunit

The units expected on input, if a specific unit (e.g. degrees) is expected.
This field is advisory, and can be left blank if there's no obvious
unit associated with the transform.

=item ounit

Same as iunit, but reporting what quantity is delivered for each dimension.

=item params

Hash ref containing relevant parameters or anything else the func needs to
work right.

=item is_inverse

Bit indicating whether the transform has been inverted.  That is useful
for some stringifications (see the PDL::Transform::Linear
stringifier), and may be useful for other things.

=back

Transforms should be inplace-aware where possible, to prevent excessive
memory usage.

If you define a new type of transform, consider generating a new stringify
method for it.  Just define the sub "stringify" in the subclass package.
It should call SUPER::stringify to generate the first line (though
the PDL::Transform::Composition bends this rule by tweaking the
top-level line), then output (indented) additional lines as necessary to
fully describe the transformation.

=head1 NOTES

Transforms have a mechanism for labeling the units and type of each
coordinate, but it is just advisory.  A routine to identify and, if
necessary, modify units by scaling would be a good idea.  Currently,
it just assumes that the coordinates are correct for (e.g.)  FITS
scientific-to-pixel transformations.

Composition works OK but should probably be done in a more
sophisticated way so that, for example, linear transformations are
combined at the matrix level instead of just strung together
pixel-to-pixel.

=head1 MODULE INTERFACE

There are both operators and constructors.  The constructors are all
exported, all begin with "t_", and all return objects that are subclasses
of PDL::Transform.

The L, L, L,
and L methods are also exported to the C package: they
are both Transform methods and PDL methods.

=cut

use strict;
use warnings;

+======EOD======


pp_addpm({At=>'Bot'},<<'+======EOD======');

=head1 AUTHOR

Copyright 2002, 2003 Craig DeForest.  There is no warranty.  You are allowed
to redistribute this software under certain conditions.  For details,
see the file COPYING in the PDL distribution.  If this file is
separated from the PDL distribution, the copyright notice should be
included in the file.

=cut

package PDL::Transform;
use Carp;
use overload '""' => \&_strval;
use overload 'x' => \&_compose_op;
use overload '**' => \&_pow_op;
use overload '!'  => \&t_inverse;

use PDL::LiteF;
use PDL::MatrixOps;

our $PI = 3.1415926535897932384626;
our $DEG2RAD = $PI / 180;
our $RAD2DEG = 180/$PI;
our $E  = exp(1);


#### little helper kludge parses a list of synonyms
sub _opt {
  my($hash) = shift;
  my($synonyms) = shift;
  my($alt) = shift;  # default is undef -- ok.
  local($_);
  foreach $_(@$synonyms){
    return (UNIVERSAL::isa($alt,'PDL')) ? PDL->pdl($hash->{$_}) : $hash->{$_}
    if defined($hash->{$_}) ;
  }
  return $alt;
}

######################################################################
#
# Stringification hack.  _strval just does a method search on stringify
# for the object itself.  This gets around the fact that stringification
# overload is a subroutine call, not a method search.
#

sub _strval {
  my($me) = shift;
  $me->stringify();
}

######################################################################
#
# PDL::Transform overall stringifier.  Subclassed stringifiers should
# call this routine first then append auxiliary information.
#
sub stringify {
  my($me) = shift;
  my($mestr) = (ref $me);
  $mestr =~ s/PDL::Transform:://;
  my $out = $mestr . " (" . $me->{name} . "): ";
  $out .= "fwd ". ((defined ($me->{func})) ? ( (ref($me->{func}) eq 'CODE') ? "ok" : "non-CODE(!!)" ): "missing")."; ";
  $out .= "inv ". ((defined ($me->{inv})) ?  ( (ref($me->{inv}) eq 'CODE') ? "ok" : "non-CODE(!!)" ):"missing").".\n";
}

+======EOD======

pp_add_exported('apply');
pp_addpm(<<'+======EOD_apply======');

=head2 apply

=for sig

  Signature: (data(); PDL::Transform t)

=for usage

  $out = $data->apply($t);
  $out = $t->apply($data);

=for ref

Apply a transformation to some input coordinates.

In the example, C<$t> is a PDL::Transform and C<$data> is a PDL to
be interpreted as a collection of N-vectors (with index in the 0th
dimension).  The output is a similar but transformed PDL.

For convenience, this is both a PDL method and a Transform method.

=cut

use Carp;

*PDL::apply = \&apply;
sub apply {
  my ($me, $from) = UNIVERSAL::isa($_[0],'PDL') ? @_[1,0] : @_;
  croak "apply requires both a PDL and a PDL::Transform.\n"
    unless UNIVERSAL::isa($me,'PDL::Transform') && UNIVERSAL::isa($from,'PDL');
  croak "Applying a PDL::Transform with no func! Oops.\n"
    unless(defined($me->{func}) and ref($me->{func}) eq 'CODE');
  my $result = $me->{func}->($from,$me->{params});
  $result->is_inplace(0); # clear inplace flag, just in case.
  if (defined $from->gethdr && $from->hdr->{NAXIS}) {
    my $nd = $me->{odim} || $me->{idim} || 2;
    for my $d (1..$nd) {
      $result->hdr->{"CTYPE$d"} = ( (!defined($me->{otype}) ? "" :
          ref($me->{otype}) ne 'ARRAY' ? $me->{otype} :
          $me->{otype}[$d-1])
         || $from->hdr->{"CTYPE$d"}
         || "");
      $result->hdr->{"CUNIT$d"} = ( (!defined($me->{ounit}) ? "" :
          ref($me->{ounit}) ne 'ARRAY' ? $me->{ounit} :
          $me->{ounit}[$d-1])
         || $from->hdr->{"CUNIT$d"}
         || $from->hdr->{"CTYPE$d"}
         || "");
    }
  }
  return $result;
}

+======EOD_apply======

pp_add_exported('invert');
pp_addpm(<<'+======EOD_invert======');

=head2 invert

=for sig

  Signature: (data(); PDL::Transform t)

=for usage

  $out = $t->invert($data);
  $out = $data->invert($t);

=for ref

Apply an inverse transformation to some input coordinates.

In the example, C<$t> is a PDL::Transform and C<$data> is an ndarray to
be interpreted as a collection of N-vectors (with index in the 0th
dimension).  The output is a similar ndarray.

For convenience this is both a PDL method and a PDL::Transform method.

=cut

*PDL::invert = \&invert;
sub invert {
  my ($me, $from) = UNIVERSAL::isa($_[0],'PDL') ? @_[1,0] : @_;
  croak "invert requires a PDL and a PDL::Transform (did you want 'inverse' instead?)\n"
    unless UNIVERSAL::isa($me,'PDL::Transform') && UNIVERSAL::isa($from,'PDL');
  croak "Inverting a PDL::Transform with no inverse! Oops.\n"
    unless(defined($me->{inv}) and ref($me->{inv}) eq 'CODE');
  my $result = $me->{inv}->($from, $me->{params});
  $result->is_inplace(0);  # make sure inplace flag is clear.
  return $result;
}

+======EOD_invert======

pp_add_exported('map');
pp_def('map',
        Pars=>'k0()',  # Dummy to set type (should match the type of "in").
        OtherPars=>'pdl *in; pdl *out; pdl *map; SV *boundary; SV *method;
                    long big; double blur; double sv_min; char flux; SV *bv',
        CHeader => pp_line_numbers(__LINE__, <<'+==EOD_map_auxiliary=='),
/*
 * Singular-value decomposition code is borrowed from
 * MatrixOps -- cut-and-pasted here because of linker trouble.
 * It's used by the auxiliary matrix manipulation code, below.
 */
static inline void pdl_xform_svd(PDL_Double *W, PDL_Double *Z, PDL_Indx nRow, PDL_Indx nCol)
{
  int i, j, k, EstColRank, RotCount, SweepCount, slimit;
  PDL_Double eps, e2, tol, vt, p, x0, y0, q, r, c0, s0, d1, d2;
  eps = 1e-6;
  slimit = nCol/4;
  if (slimit < 6.0)
    slimit = 6;
  SweepCount = 0;
  e2 = 10.0*nRow*eps*eps;
  tol = eps*.1;
  EstColRank = nCol;
  for (i=0; i
    for (j=0; j
      W[nCol*(nRow+i)+j] = 0.0;
    }
    W[nCol*(nRow+i)+i] = 1.0;  /* rjrw 7/7/99: moved this line out of j loop */
  }
  RotCount = EstColRank*(EstColRank-1)/2;
  while (RotCount != 0 && SweepCount <= slimit)
    {
      RotCount = EstColRank*(EstColRank-1)/2;
      SweepCount++;
      for (j=0; j
        {
          for (k=j+1; k
            {
              p = q = r = 0.0;
              for (i=0; i
                {
                  x0 = W[nCol*i+j]; y0 = W[nCol*i+k];
                  p += x0*y0; q += x0*x0; r += y0*y0;
                }
              Z[j] = q; Z[k] = r;
              if (q >= r)
                {
                  if (q<=e2*Z[0] || fabs(p)<=tol*q) RotCount--;
                  else
                    {
                      p /= q; r = 1 - r/q; vt = sqrt(4*p*p+r*r);
                      c0 = sqrt(fabs(.5*(1+r/vt))); s0 = p/(vt*c0);
                      for (i=0; i
                        {
                          d1 = W[nCol*i+j]; d2 = W[nCol*i+k];
                          W[nCol*i+j] = d1*c0+d2*s0; W[nCol*i+k] = -d1*s0+d2*c0;
                        }
                    }
                }
              else
                {
                  p /= r; q = q/r-1; vt = sqrt(4*p*p+q*q);
                  s0 = sqrt(fabs(.5*(1-q/vt)));
                  if (p<0) s0 = -s0;
                  c0 = p/(vt*s0);
                  for (i=0; i
                    {
                      d1 = W[nCol*i+j]; d2 = W[nCol*i+k];
                      W[nCol*i+j] = d1*c0+d2*s0; W[nCol*i+k] = -d1*s0+d2*c0;
                    }
                }
            }
        }
      while (EstColRank>=3 && Z[(EstColRank-1)]<=Z[0]*tol+tol*tol)
        EstColRank--;
    }
}

/*
 * PDL_xform_aux:
 *  This handles the matrix manipulation part of the Jacobian filtered
 *  mapping code.  It's separate from the main code because it's
 *  independent of the data type of the original arrays.
 *
 *Given a pre-allocated workspace  and
 * an integer set of coordinates, generate the discrete Jacobian
 * from the map, pad the singular values, and return the inverse
 * Jacobian, the largest singular value of the Jacobian itself, and
 * the determinant of the original Jacobian.  Boundary values use the
 * asymmetric discrete Jacobian; others use the symmetric discrete Jacobian.
 *
 * The map and workspace must be of type PDL_D.  If the dimensionality is
 * d, then the workspace must have at least 3*n^2+n elements.  The
 * inverse of the padded Jacobian is returned in the first n^2 elements.
 * The determinant of the original Jacobian gets stuffed into the n^2
 * element of the workspace. The largest padded singular value is returned.
 *
 */

static inline PDL_Double PDL_xform_aux ( pdl *map, PDL_Indx *coords, PDL_Double *tmp, PDL_Double sv_min) {
  PDL_Indx ndims;
  PDL_Indx i, j, k;
  PDL_Indx offset;
  PDL_Double det;
  PDL_Double *jptr;
  PDL_Double *svptr;
  PDL_Double *aptr,*bptr;
  PDL_Double max_sv = 0.0;

  ndims = map->ndims-1;

  /****** Accumulate the Jacobian */
  /* Accumulate the offset into the map array */
  for( i=offset=0; i
    offset += coords[i] * map->dimincs[i+1];

  jptr = tmp + ndims*ndims;

  for( i=0; i
    char bot = (coords[i] <=0);
    char top = (coords[i] >= map->dims[i+1]-1);
    char symmetric = !(bot || top);
    PDL_Double *ohi,*olo;
    PDL_Indx diff = map->dimincs[i+1];
    ohi = ((PDL_Double *)map->data) + ( offset + ( top ? 0 : diff ));
    olo = ((PDL_Double *)map->data) + ( offset - ( bot ? 0 : diff ));

    for( j=0; j
       PDL_Double jel =   *ohi - *olo;

       ohi += map->dimincs[0];
       olo += map->dimincs[0];

       if(symmetric)
         jel /= 2;
       *(jptr++) = jel;
    }
  }


  /****** Singular-value decompose the Jacobian
   * The svd routine produces the squares of the singular values,
   * and requires normalization for one of the rotation matrices.
   */

  jptr =   tmp + ndims*ndims;
  svptr = tmp + 3*ndims*ndims;
  pdl_xform_svd(jptr,svptr,ndims,ndims);

  aptr = svptr;
  for (i=0;i
    *aptr = sqrt(*aptr);


  /* fix-up matrices here */
  bptr = jptr;
  for(i=0; i
    aptr = svptr;
    for(j=0;j
     *(bptr++) /= *(aptr++);
  }

  /****** Store the determinant, and pad the singular values as necessary.*/
  aptr = svptr;
  det = 1.0;
  for(i=0;i
    det *= *aptr;
    if(*aptr < sv_min)
      *aptr = sv_min;
    if(*aptr > max_sv )
        max_sv = *aptr;
    aptr++;
  }

  /****** Generate the inverse matrix */
  /* Multiply B-transpose times 1/S times A-transpose. */
  /* since S is diagonal we just divide by the appropriate element. */
  /* */
  aptr = tmp  + ndims*ndims;
  bptr = aptr + ndims*ndims;
  jptr=  tmp;

  for(i=0;i
    for(j=0;j
        *jptr = 0;

        for(k=0;k
          *jptr += aptr[j*ndims + k] * bptr[k*ndims + i] / *svptr;

        jptr++;
    }
    svptr++;
  }

 *jptr = det;
  return max_sv;
}
+==EOD_map_auxiliary==
        Code => pp_line_numbers(__LINE__, <<'+==EOD_map_c_code=='),

/*
 * Pixel interpolation & averaging code
 *
 * Calls a common coordinate-transformation block (see following hdr)
 * that isn't dependent on the type of the input variable.
 *
 * The inputs are SVs to avoid hassling with broadcastloops; broadcasting
 * is handled internally.  To simplify the broadcasting business, any
 * broadcast dimensions should all be collapsed to a single one by the
 * perl front-end.
 *
 */

 pdl *in  = $COMP(in);
 pdl *out = $COMP(out);
 pdl *map = $COMP(map);

 PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
 PDL_Double tmp[3*ndims*ndims+ndims]; /* Workspace for prefrobnication */
 PDL_Indx ovec[ndims];      /* output pixel loop vector                 */
 PDL_Indx ivec[ndims];      /* input pixel loop vector                  */
 PDL_Indx ibvec[ndims];     /* input pixel base offset vector           */
 PDL_Double dvec[ndims];    /* Residual vector for linearization        */
 PDL_Double tvec[ndims];    /* Temporary floating-point vector          */
 PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
 PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
 PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
 char bounds[ndims];        /* Boundary condition packed string */
 PDL_Indx  index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
 char   method;             /* Method identifier (gets one of 'h','g')  */
 PDL_Long   big;            /* Max size of input footprint for each pix */
 PDL_Double blur;           /* Scaling of filter                        */
 PDL_Double sv_min;         /* minimum singular value                   */
 char flux;                 /* Flag to indicate flux conservation       */
 PDL_Indx i;
 $GENERIC() badval = SvNV($COMP(bv));
 #define HANNING_LOOKUP_SIZE 2500
 static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
 static int needs_hanning_calc = 1;
 PDL_Double zeta = 0;
 PDL_Double hanning_offset;

 #define GAUSSIAN_LOOKUP_SIZE 4000
 #define GAUSSIAN_MAXVAL 6.25  /* 2.5 HWHMs (square it) */
 static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
 static int needs_gaussian_calc = 1;

 PDL_RETERROR(PDL_err, PDL->make_physical(in));
 PDL_RETERROR(PDL_err, PDL->make_physical(out));
 PDL_RETERROR(PDL_err, PDL->make_physical(map));

/***
 * Fill in the boundary condition array
 */
 {
  char *bstr;
  STRLEN blen;
  bstr = SvPV($COMP(boundary),blen);

  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    PDL_Indx i;
    for (i=0;i
      bounds[i] = 1;
  } else {
    PDL_Indx i;
    for(i=0;i
      switch(bstr[i < blen ? i : blen-1 ]) {
      case '0': case 'f': case 'F':               /* forbid */
        bounds[i] = 0;
        break;
      case '1': case 't': case 'T':               /* truncate */
        bounds[i] = 1;
        break;
      case '2': case 'e': case 'E':               /* extend */
        bounds[i] = 2;
        break;
      case '3': case 'p': case 'P':               /* periodic */
        bounds[i] = 3;
        break;
      case '4': case 'm': case 'M':               /* mirror */
        bounds[i] = 4;
        break;
      default:
        $CROAK("Error in map: Unknown boundary condition '%c'",bstr[i]);
        break;
      }
    }
  }
 }

/***
 * Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
 */
 big    = labs($COMP(big));
 if(big <= 0)
   $CROAK("map: 'big' parameter must be >0");

 blur   = fabs($COMP(blur));
 if(blur < 0)
   $CROAK("map: 'blur' parameter must be >= 0");

 sv_min = fabs($COMP(sv_min));
 if(sv_min < 0)
   $CROAK("map: 'sv_min' parameter must be >= 0");

 flux   = $COMP(flux) != 0;

 {
  char *mstr;
  STRLEN mlen;
  mstr = SvPV($COMP(method),mlen);

  if(mlen==0)
        method = 'h';
  switch(*mstr) {
    case 'h':                 // h - lookup-table hanning window
      if( needs_hanning_calc ) {
         int i;
         for(i=0;i
           hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
         }
         hanning_lookup[HANNING_LOOKUP_SIZE]   = 0;
         hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
         needs_hanning_calc = 0;
      }
      zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
    case 'H':                 // H - rigorous hanning window
      method = *mstr;
      hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
      break;

    case 'g': case 'j':       // Gaussian and/or Jacobian, using lookup table
      method = 'g';
      zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;

      if( needs_gaussian_calc ) {
        int i;
	for(i=0;i
	  gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
	}
	gaussian_lookup[GAUSSIAN_LOOKUP_SIZE]   = 0;
	gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
	needs_gaussian_calc = 0;
      }
      break;

    case 'G': case 'J':       // Jacobian -- same thing, with direct calculation
      method = 'G'; break;
    default:
      $CROAK("Bug in map: unknown method '%c'",*mstr);
      break;
    }
 }



 /* End of initialization */
 /*************************************************************/
 /* Start of Real Work */

 /* Initialize coordinate vector and map offset
  */
 for(i=0;i
   ovec[i] = 0;

 PDL_Double *map_ptr = (PDL_Double *)(map->data);
 PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
 for (i=0; i < npdls; i++) offs[i] = 0;
 for (i=0; i < ndims; i++) {
   incs[i*npdls + 0] = map->dimincs[i+1];
 }

/* Main pixel loop (iterates over pixels in the output plane) */
 do {
   PDL_Indx psize; // Size of the region to accumulate over for this pixel
   PDL_Indx i_off; // Offset into the data array of the source PDL
   PDL_Indx j;     // Generic loop index
   char t_vio;     // counter used for truncation boundary violations
   char carry;     // flag used for multidimensional loop iteration

   /* Prefrobnicate the transformation matrix */
   psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */

#ifdef DEBUG_MAP
 {
     PDL_Indx k; PDL_Indx foo = 0;
     printf("ovec: [");
     for(k=0;k
       foo += ovec[k] * map->dimincs[k+1];
       printf(" %2td ",(PDL_Indx)(ovec[k]));
     }
     printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
     for(k=0;k
       printf("%8.2f",(double)((($GENERIC() *)(map->data))[foo + k*map->dimincs[0]]));
     }
     printf("]\n");
 }
#endif

   /* Don't bother accumulating output if psize is too large */
   if(psize <= big) {
     /* Use the prefrobnicated matrix to generate a local linearization.
      * dvec gets the delta; ibvec gets the base.
      */
     {
       PDL_Double *mp = map_ptr + offs[0];
       for (i=0;i
        dvec[i]  = *mp - ( ibvec[i] = (PDL_Long)(*mp + 0.5));   /* assignment */
        mp += map->dimincs[0];
       }
     }

     /* Initialize input delta vector  */
     for(i=0;i
       ivec[i] = -psize;

     /* Initialize accumulators */
     {
       PDL_Double *ac = acc;
       for(i=0; i < in->dims[ndims]; i++)
         *(ac++) = 0.0;

     }
     {
       PDL_Double *wg = wgt;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
     }
     {
       PDL_Double *wg = wgt2;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
      }


     /*
      * Calculate the original offset into the data array, to enable
      * delta calculations in the pixel loop
      *
      * i runs over dims; j holds the working integer index in the
      * current dim.
      *
      * This code matches the incrementation code at the bottom of the accumulation loop
      */

     t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
     i_off = 0;
     for(i=0;i
       j = ivec[i] + ibvec[i];
       if(j<0 || j >= in->dims[i]) {
         switch(bounds[i]) {
         case 0: /* no breakage allowed */
           $CROAK("index out-of-bounds in map");
           break;
         case 1: /* truncation */
           t_vio++;
           /* fall through */
         case 2: /* extension -- crop */
           if(j<0)
             j=0;
           else //if(j>=in->dims[i])
             j = in->dims[i] - 1;
           break;
         case 3: /* periodic -- mod it */
           j %= in->dims[i];
           if(j<0)
             j += in->dims[i];
           break;
         case 4: /* mirror -- reflect off the edges */
           j += in->dims[i];
           j %= (in->dims[i]*2);
           if(j<0)
             j += in->dims[i]*2;
           j -= in->dims[i];
           if(j<0) {
             j *= -1;
             j -= 1;
           }
           break;
         default:
           $CROAK("Unknown boundary condition %td in map -- bug alert!", (PDL_Indx)bounds[i]);
           break;
         }
       }
       i_off += in->dimincs[i] * j;
     }

     /* Initialize index stashes for later reference as we scan the footprint */
     /* It's a pain in the ass to deal with boundaries, and doubly so at the  */
     /* end of a dimensional scan.  So we stash the index location at the     */
     /* start of each dimensional scan here.  When we finish incrementing     */
     /* through a particular dim, we pull its value back out of the stash.    */
     for(i=0;i
       index_stash[i] = i_off;
    }

     /* The input accumulation loop is the hotspot for the entire operation.        */
     /* We loop over pixels in the region of interest (+/- psize in each dimension) */
     /* in the input array, use the linearized transform to bring each pixel center */
     /* forward to the output plane, and calculate a weighting based on the chosen  */
     /* filter function.   'h' is a fast Hanning window rolloff using a lookup      */
     /* table that is initialized the first time through the code.  'H' is the      */
     /* same process, but explicitly calculated for each interation (~2x slower).   */
     /* 'g' uses a radial Gaussian filter.  Rather than calculate the array offset  */
     /* into the input array fresh from the current input array vector each time,   */
     /* we walk through the array using dimincs and the old offset.  This saves     */
     /* about half of the time spent on index calculation.                          */

     do {  /* Input accumulation loop */
         PDL_Double *cp;
         PDL_Double alpha; // weighting coefficient for the current point

       /* Calculate the weight of the current input point.  Don't bother if we're
        * violating any truncation boundaries (in that case our value is zero, but
        * for the interpolation we also set the weight to zero).
        */
       if( !t_vio ) {

         PDL_Double *ap = tvec;
         PDL_Double *bp = dvec;
         PDL_Indx   *ip = ivec;
         for(i=0; i
           *(ap++) = *(ip++) - *(bp++);

         switch(method) {
           PDL_Double dd; // tmp space for calculations in filters

         case 'h':
           /* This is the Hanning window rolloff.  It is a product of a simple */
           /* cos^2(theta) rolloff in each dimension.  Using a lookup table    */
           /* is about 2x faster than using cos(theta) directly in each        */
           /* weighting calculation, so we do.  Using 2500 entries and linear  */
           /* interpolation is accurate to about 10^-7, and should preserve    */
           /* the contents of cache pretty well.                               */
           alpha = 1;
           cp = tmp;
           for (i=0; i
             PDL_Indx lodex;
             PDL_Indx hidex;
             PDL_Double beta;
             dd = 0;
             ap = tvec;
             /* Get the matrix-multiply element for this dimension */
             for (j=0;j
               dd += *(ap++) * *(cp++);

             /* Do linear interpolation from the table */
             /* The table captures a hanning window centered 0.5 pixel from center. */
             /* We scale the filter by the blur parameter -- but if blur is less    */
             /* than unity, we shrink the hanning blur window while keeping the 0.5 */
             /* value on the pixel edge at 0.5.  For blur greater than unity, we    */
             /* scale simply. */
             beta = fabs(dd) - hanning_offset;
             if (beta > 0) {
               if (beta >= blur) {
                 alpha = 0;
                 i = ndims;
               } else {
                 beta *= zeta;
                 lodex = beta;
                 beta -= lodex;  //lodex now has the integer part, beta has the fractional part
                 if (lodex > HANNING_LOOKUP_SIZE)
                 lodex = HANNING_LOOKUP_SIZE;
                 hidex = lodex+1;
                 //do a linear interpolation between the two nearest values
                 alpha *= hanning_lookup[hidex]*beta + hanning_lookup[lodex]*(1-beta);
               } /* end of interpolation branch */
             } /* end of beta > 0 branch */
           } /* end of dimension loop */
         break;

        case 'H':
          /* This is the Hanning window rolloff with explicit calculation, preserved */
          /* in case someone actually wants the slower longer method. */
               alpha = 1;
               cp = tmp;
               for (i=0; i
                 dd = 0;
                 ap = tvec;
                 for (j=0;j
                     dd += *(ap++) * *(cp++);
                 dd = (fabs(dd) - hanning_offset) / blur;
                 if ( dd > 1 ) {
                     alpha = 0;
                     i = ndims;
                 } else
                     alpha *= (0.5 + 0.5 * cos( dd * 3.1415926536 ));
               }
              break;

	case 'g':
	  /* This is the Gaussian rolloff.  It does lookup into a precalculated exponential. */
	  {
	    PDL_Double sum = 0;
	    cp = tmp;
	    for(i=0; i
	      dd = 0;
              ap = tvec;
	      for(j=0;j
	        dd += *(ap++) * *(cp++);
	      dd /= blur;
	      sum += dd * dd;
	      if(sum > GAUSSIAN_MAXVAL) {
		  i = ndims; /* exit early if we're too far out */
		  alpha = 0;
              }
	    }
	    if( sum > GAUSSIAN_MAXVAL || !isfinite(sum) || isnan(sum) ) {
	      alpha = 0;
	    } else {
	      PDL_Indx lodex,hidex;
	      PDL_Double beta = fabs(zeta * sum);

	      lodex = beta;
	      beta -= lodex; // lodex now has the integer part, beta has the fractional part
	      hidex = lodex+1;
	      //do a linear interpolation between the two nearest values
              alpha = gaussian_lookup[hidex]*beta + gaussian_lookup[lodex]*(1 - beta);

            }
          }
          break;

        case 'G':
          /* This is the Gaussian rolloff with explicit calculation, preserved */
          /* in case someone actually wants the slower longer method. */
              {
                 PDL_Double sum = 0;
                 cp = tmp;
                 for(i=0; i
                   dd = 0;
                   ap = tvec;
                   for(j=0;j
                     dd += *(ap++) * *(cp++);
                   dd /= blur;
                   sum += dd * dd;
                   if(sum > 4) /* 2 pixels -- four half-widths */
                    i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
                 }

                 if(sum > GAUSSIAN_MAXVAL)
                   alpha = 0;
                 else
                   alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
              }
              break;
         default:
           $CROAK("This can't happen: method='%c'",method);
         }

         {  /* convenience block -- accumulate the current point into the weighted sum.  */
            /* This is more than simple assignment because we have our own explicit poor */
            /* man's broadcastloop here, so we accumulate each broadcasted element separately. */
            $GENERIC() *dat = (($GENERIC() *)(in->data)) + i_off;
            PDL_Indx max = out->dims[ndims];
            for( i=0; i < max; i++ ) {
	      if( (badval==0) || (*dat != badval) ) {
                 acc[i] += *dat * alpha;
                 dat += in->dimincs[ndims];
	         wgt[i] += alpha;
	      }
	      wgt2[i] += alpha;  // Accumulate what weight we would have with no bad values
            }
          }
       } /* end of t_vio check (i.e. of input accumulation) */


       /* Advance input accumulation loop.  */
       /* We both increment the total vector and also advance the index. */
       carry = 1;
       for (i=0; i
         /* Advance the current element of the offset vector */
         ivec[i]++;
         j = ivec[i] + ibvec[i];

         /* Advance the offset into the data array */
         if( j > 0 && j <= in->dims[i]-1 ) {
            /* Normal case -- just advance the input vector */
            i_off += in->dimincs[i];
         } else {
            /* Busted a boundary - either before or after.  */
            switch(bounds[i]){
            case 0: /* no breakage allowed -- treat as truncation for interpolation */
            case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
              if( j == 0 )
                t_vio--;
              else if( j == in->dims[i] )
                t_vio++;
              break;
             case 2: /* extension -- do nothing (so the same input point is re-used) */
               break;
             case 3: /* periodic  -- advance and mod into the allowed range */
               if((j % in->dims[i]) == 0) {
                 i_off -= in->dimincs[i] * (in->dims[i]-1);
               } else {
                 i_off += in->dimincs[i];
               }
               break;
            case 4: /* mirror -- advance or retreat depending on phase */
               j += in->dims[i];
               j %= (in->dims[i]*2);
               j -= in->dims[i];
               if( j!=0 && j!= -in->dims[i] ) {
                 if(j<0)
                   i_off -= in->dimincs[i];
                 else
                   i_off += in->dimincs[i];
               }
               break;
            }
         }

         /* Now check for carry */
         if (ivec[i] <= psize) {
           /* Normal case -- copy the current offset to the faster-running dim stashes */
           PDL_Indx k;
           for(k=0;k
             index_stash[k] = i_off;
           }
           carry = 0;

         } else { /* End of this scan -- recover the last position, and mark carry */
           i_off = index_stash[i];
           if(bounds[i]==1) {
             j = ivec[i] + ibvec[i];
             if( j < 0 || j >= in->dims[i] )
               t_vio--;
             ivec[i] = -psize;
             j = ivec[i] + ibvec[i];
             if( j < 0 || j >= in->dims[i] )
               t_vio++;
             carry = 1;
           } else {
             ivec[i] = -psize;
           }
         }
       } /* End of counter-advance loop */
     } while (carry==0);  /* end of total data accumulation loop (termination condition has carry on last dim) */

     {
       PDL_Double *ac = acc;
       PDL_Double *wg = wgt;
       PDL_Double *wg2 = wgt2;
       $GENERIC() *dat = out->data;

       /* Calculate output vector offset */
       for(i=0;i
         dat += out->dimincs[i] * ovec[i];

       if (!flux) {
         /* Flux flag is NOT set -- normal case.  Copy the weighted accumulated data. */
         for (i=0; i < out->dims[ndims]; i++) {
           *dat = (*wg && (*wg2 / *wg) < 1.5) ? *ac / *wg : badval;
           ac++; wg++; wg2++;
           dat += out->dimincs[ndims];
         }
       } else {
         /* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
         PDL_Double det = tmp[ndims*ndims];
         for(i=0; i < out->dims[ndims]; i++) {
	   if(*wg && (*wg2 / *wg) < 1.5 ) {
              *dat = *(ac++) / *(wg++) * det;
	      wg2++;
           } else {
	      *dat = badval;
	      ac++; wg++; wg2++;
	   }
           dat += out->dimincs[ndims];
         } /* end of for loop */
       } /* end of flux flag set conditional */
     } /* end of convenience block */

     /* End of code for normal pixels */
   } else {
     /* The pixel was ludicrously huge -- just set this pixel to nan */
     $GENERIC() *dat = out->data;
      for(i=0;i
        dat += out->dimincs[i] * ovec[i];
      for(i=0;idims[ndims];i++) {
        *dat = badval;          /* Should handle bad values too -- not yet */
        dat += out->dimincs[ndims];
      }
   }

   /* Increment the pixel counter */
   if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
  } while (1);
  PDL->changed(out, PDL_PARENTDATACHANGED, 0);
+==EOD_map_c_code==

        Doc=><<'+==EOD_map_doc==',

=head2 match

=for usage

  $y = $x->match($c);                  # Match $c's header and size
  $y = $x->match([100,200]);           # Rescale to 100x200 pixels
  $y = $x->match([100,200],{rect=>1}); # Rescale and remove rotation/skew.

=for ref

Resample a scientific image to the same coordinate system as another.

The example above is syntactic sugar for

 $y = $x->map(t_identity, $c, ...);

it resamples the input PDL with the identity transformation in
scientific coordinates, and matches the pixel coordinate system to
$c's FITS header.

There is one difference between match and map: match makes the
C option to C default to 0, not 1.  This only affects
matching where autoscaling is required (i.e. the array ref example
above).  By default, that example simply scales $x to the new size and
maintains any rotation or skew in its scientific-to-pixel coordinate
transform.

=head2 map

=for usage

  $y = $x->map($xform,[