File Coverage

Transform.xs

Criterion	Covered	Total	%
statement	627	2679	23.4
branch	368	1766	20.8
condition			n/a
subroutine			n/a
pod			n/a
total	995	4445	22.3

/*
 * THIS FILE WAS GENERATED BY PDLA::PP! Do not modify!
 */

#define PDLA_COMMENT(comment)
PDLA_COMMENT("This preprocessor symbol is used to add commentary in the PDLA  ")
PDLA_COMMENT("autogenerated code. Normally, one would use typical C-style    ")
PDLA_COMMENT("multiline comments (i.e. /* comment */). However, because such ")
PDLA_COMMENT("comments do not nest, it's not possible for PDLA::PP users to   ")
PDLA_COMMENT("comment-out sections of code using multiline comments, as is   ")
PDLA_COMMENT("often the practice when debugging, for example. So, when you   ")
PDLA_COMMENT("see something like this:                                       ")
PDLA_COMMENT("                                                               ")
                PDLA_COMMENT("Memory access")
PDLA_COMMENT("                                                               ")
PDLA_COMMENT("just think of it as a C multiline comment like:                ")
PDLA_COMMENT("                                                               ")
PDLA_COMMENT("   /* Memory access */                                         ")

#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include "pdl.h"
#include "pdlcore.h"
static Core* PDLA; PDLA_COMMENT("Structure hold core C functions")
static int __pdl_debugging = 0;
static int __pdl_boundscheck = 0;
static SV* CoreSV;       PDLA_COMMENT("Gets pointer to perl var holding core structure")

#if ! 1
# define PP_INDTERM(max, at) at
#else
# define PP_INDTERM(max, at) (__pdl_boundscheck? PDLA->safe_indterm(max,at, __FILE__, __LINE__) : at)
#endif

#include


/*
 * 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.
 *
 */

void pdl_xform_svd(PDLA_Double *W, PDLA_Double *Z, int nRow, int nCol)
{
  int i, j, k, EstColRank, RotCount, SweepCount, slimit;
  PDLA_Double eps, e2, tol, vt, p, h2, x0, y0, q, r, c0, s0, c2, 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--;
    }
}

/*
 * PDLA_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 PDLA_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.
 *
 */

PDLA_Double PDLA_xform_aux ( pdl *map, PDLA_Indx *coords, PDLA_Double *tmp, PDLA_Double sv_min) {
  short ndims;
  PDLA_Long i, j, k;
  PDLA_Long offset;
  PDLA_Double det;
  PDLA_Double *jptr;
  PDLA_Double *svptr;
  PDLA_Double *aptr,*bptr;
  PDLA_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);
    PDLA_Double *ohi,*olo;
    PDLA_Long diff = map->dimincs[i+1];
    ohi = ((PDLA_Double *)map->data) + ( offset + ( top ? 0 : diff ));
    olo = ((PDLA_Double *)map->data) + ( offset - ( bot ? 0 : diff ));

    for( j=0; j
       PDLA_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;
}

typedef struct pdl_map_struct {
		PDLA_TRANS_START(1);
		pdl_thread  __pdlthread;
		SV  *in;SV  *out;SV  *map;SV  *boundary;SV  *method;SV  *big;SV  *blur;SV  *sv_min;SV  *flux;SV  *bv;
		char __ddone; PDLA_COMMENT("Dims done")
		} pdl_map_struct;

void pdl_map_redodims(pdl_trans *__tr ) {
	int __dim;
	pdl_map_struct *__privtrans = (pdl_map_struct *) __tr;

	{
	    PDLA_Indx __creating[1];
__creating[0] = 0;
 {
{PDLA_COMMENT("Start generic loop")

	switch(__privtrans->__datatype) { case -42: PDLA_COMMENT("Warning eater") {(void)1;
	} break; case PDLA_B: {
	PDLA_Byte * k0_datap = ((PDLA_Byte *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Byte * k0_physdatap = ((PDLA_Byte *)((__privtrans->pdls[0])->data));

{
PDLA_COMMENT("none")
}	} break; case PDLA_S: {
	PDLA_Short * k0_datap = ((PDLA_Short *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Short * k0_physdatap = ((PDLA_Short *)((__privtrans->pdls[0])->data));

{
PDLA_COMMENT("none")
}	} break; case PDLA_US: {
	PDLA_Ushort * k0_datap = ((PDLA_Ushort *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Ushort * k0_physdatap = ((PDLA_Ushort *)((__privtrans->pdls[0])->data));

{
PDLA_COMMENT("none")
}	} break; case PDLA_L: {
	PDLA_Long * k0_datap = ((PDLA_Long *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Long * k0_physdatap = ((PDLA_Long *)((__privtrans->pdls[0])->data));

{
PDLA_COMMENT("none")
}	} break; case PDLA_IND: {
	PDLA_Indx * k0_datap = ((PDLA_Indx *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Indx * k0_physdatap = ((PDLA_Indx *)((__privtrans->pdls[0])->data));

{
PDLA_COMMENT("none")
}	} break; case PDLA_LL: {
	PDLA_LongLong * k0_datap = ((PDLA_LongLong *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_LongLong * k0_physdatap = ((PDLA_LongLong *)((__privtrans->pdls[0])->data));

{
PDLA_COMMENT("none")
}	} break; case PDLA_F: {
	PDLA_Float * k0_datap = ((PDLA_Float *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Float * k0_physdatap = ((PDLA_Float *)((__privtrans->pdls[0])->data));

{
PDLA_COMMENT("none")
}	} break; case PDLA_D: {
	PDLA_Double * k0_datap = ((PDLA_Double *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Double * k0_physdatap = ((PDLA_Double *)((__privtrans->pdls[0])->data));

{
PDLA_COMMENT("none")
}	break;}
	default:barf("PP INTERNAL ERROR! PLEASE MAKE A BUG REPORT\n");}
}
}
 {
 static char *__parnames[] = {"k0"};
		static PDLA_Indx __realdims[] = {0};
		static char __funcname[] = "PDLA::Transform::map";
		static pdl_errorinfo __einfo = {
			__funcname, __parnames, 1
		};

		 PDLA->initthreadstruct(2,__privtrans->pdls,
			__realdims,__creating,1,
                      &__einfo,&(__privtrans->__pdlthread),
                        __privtrans->vtable->per_pdl_flags,
			0 );
		}

{ PDLA_COMMENT("convenience block")
  void *hdrp = NULL;
  char propagate_hdrcpy = 0;
  SV *hdr_copy = NULL;
      if(!hdrp &&
	      __privtrans->pdls[0]->hdrsv &&
	 (__privtrans->pdls[0]->state & PDLA_HDRCPY)
	 ) {
	hdrp = __privtrans->pdls[0]->hdrsv;
	propagate_hdrcpy = ((__privtrans->pdls[0]->state & PDLA_HDRCPY) != 0);
      }
if (hdrp) {
  if(hdrp == &PL_sv_undef)
    hdr_copy = &PL_sv_undef;
  else  {  PDLA_COMMENT("Call the perl routine _hdr_copy...")
    int count;
    PDLA_COMMENT("Call the perl routine PDLA::_hdr_copy(hdrp)")
    dSP;
    ENTER ;
    SAVETMPS ;
    PUSHMARK(SP) ;
    XPUSHs( hdrp );
    PUTBACK ;
    count = call_pv("PDLA::_hdr_copy",G_SCALAR);
    SPAGAIN ;
    if(count != 1)
	croak("PDLA::_hdr_copy didn't return a single value - please report this bug (A).");

    hdr_copy = (SV *)POPs;

    if(hdr_copy && hdr_copy != &PL_sv_undef) {
       (void)SvREFCNT_inc(hdr_copy); PDLA_COMMENT("Keep hdr_copy from vanishing during FREETMPS")
    }

    FREETMPS ;
    LEAVE ;


  } PDLA_COMMENT("end of callback  block")


         if(hdr_copy != &PL_sv_undef)
            SvREFCNT_dec(hdr_copy); PDLA_COMMENT("make hdr_copy mortal again")
      } PDLA_COMMENT("end of if(hdrp) block")
   } PDLA_COMMENT("end of conv. block")
 __privtrans->__ddone = 1;
	}
    }


pdl_trans * pdl_map_copy(pdl_trans *__tr ) {
	int __dim;
	pdl_map_struct *__privtrans = (pdl_map_struct *) __tr;

	{
	    pdl_map_struct *__copy = malloc(sizeof(pdl_map_struct));
PDLA_THR_CLRMAGIC(&__copy->__pdlthread);			PDLA_TR_CLRMAGIC(__copy);
                        __copy->has_badvalue = __privtrans->has_badvalue;
                        __copy->badvalue = __privtrans->badvalue;
			__copy->flags = __privtrans->flags;
			__copy->vtable = __privtrans->vtable;
			__copy->__datatype = __privtrans->__datatype;
			__copy->freeproc = NULL;
			__copy->__ddone = __privtrans->__ddone;
			{int i;
			 for(i=0; i<__copy->vtable->npdls; i++)
				__copy->pdls[i] = __privtrans->pdls[i];
			}
			(__copy->in) = newSVsv(__privtrans->in);;(__copy->out) = newSVsv(__privtrans->out);;(__copy->map) = newSVsv(__privtrans->map);;(__copy->boundary) = newSVsv(__privtrans->boundary);;(__copy->method) = newSVsv(__privtrans->method);;(__copy->big) = newSVsv(__privtrans->big);;(__copy->blur) = newSVsv(__privtrans->blur);;(__copy->sv_min) = newSVsv(__privtrans->sv_min);;(__copy->flux) = newSVsv(__privtrans->flux);;(__copy->bv) = newSVsv(__privtrans->bv);;
			if(__copy->__ddone) {
				PDLA->thread_copy(&(__privtrans->__pdlthread),&(__copy->__pdlthread));
			}
			return (pdl_trans*)__copy;
	}
    }


void pdl_map_readdata(pdl_trans *__tr ) {
	int __dim;
	pdl_map_struct *__privtrans = (pdl_map_struct *) __tr;

	{
	    {PDLA_COMMENT("Start generic loop")

	switch(__privtrans->__datatype) { case -42: PDLA_COMMENT("Warning eater") {(void)1;
	} break; case PDLA_B: {
	PDLA_Byte * k0_datap = ((PDLA_Byte *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Byte * k0_physdatap = ((PDLA_Byte *)((__privtrans->pdls[0])->data));


	PDLA_COMMENT("THREADLOOPBEGIN")
	if ( PDLA->startthreadloop(&(__privtrans->__pdlthread),__privtrans->vtable->readdata, __tr) ) return;
           do { register PDLA_Indx __tind1=0,__tind2=0;
                register PDLA_Indx __tnpdls = __privtrans->__pdlthread.npdls;
                register PDLA_Indx __tdims1 = __privtrans->__pdlthread.dims[1];
                register PDLA_Indx __tdims0 = __privtrans->__pdlthread.dims[0];
                register PDLA_Indx *__offsp = PDLA->get_threadoffsp(&__privtrans->__pdlthread);
	register PDLA_Indx __tinc0_0 = __privtrans->__pdlthread.incs[0];
	register PDLA_Indx __tinc1_0 = __privtrans->__pdlthread.incs[__tnpdls+0];
	k0_datap += __offsp[0];
	for( __tind2 = 0 ;
                __tind2 < __tdims1 ;
                __tind2++
			,k0_datap += __tinc1_0 - __tinc0_0 * __tdims0
	)
	{
	     for( __tind1 = 0 ;
                  __tind1 < __tdims0 ;
                  __tind1++
			,k0_datap += __tinc0_0
	)
	{  PDLA_COMMENT("This is the tightest threadloop. Make sure inside is optimal."){

/*
 * 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 threadloops; threading
 * is handled internally.  To simplify the threading business, any
 * thread dimensions should all be collapsed to a single one by the
 * perl front-end.
 *
 */

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

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

 pdl *in  = PDLA->SvPDLAV(__privtrans->in);
 pdl *out = PDLA->SvPDLAV(__privtrans->out);
 pdl *map = PDLA->SvPDLAV(__privtrans->map);

 PDLA->make_physical(in);
 PDLA->make_physical(out);
 PDLA->make_physical(map);

 ndims = map->ndims -1;

 /*
  * Allocate all our dynamic workspaces at once
  * */
 ovec = (PDLA_Indx *)(PDLA->smalloc( (STRLEN)
                                   ( + sizeof(PDLA_Indx)    * 3 * ndims                                                    + sizeof(PDLA_Double) * (3*ndims)                                                     + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * 3 * ndims*ndims + ndims                                       + sizeof(char) * ndims                                                               + sizeof(PDLA_Indx)   * ndims                                                         )
                                   )
                );
 ivec   =                   &(ovec[ndims]);
 ibvec  =                   &(ivec[ndims]);
 dvec   =    (PDLA_Double *)(&(ibvec[ndims]));
 tvec   =                   &(dvec[ndims]);
 acc    =                   &(tvec[ndims]);
 wgt    =                   &(acc[in->dims[ndims]]);   wgt2   =                   &(wgt[in->dims[ndims]]);   tmp    =                   &(wgt2[in->dims[ndims]]);
 bounds =          (char *)(&(tmp [3*ndims*ndims+ndims]));
 index_stash = (PDLA_Indx *) &(bounds[ndims]);


/***
 * Fill in the boundary condition array
 */
 {
  char *bstr;
  STRLEN blen;
  bstr = SvPV(__privtrans->boundary,blen);

  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    int i;
    for (i=0;i
      bounds[i] = 1;
  } else {
    int 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:
        {
          char buf[BUFSIZ];
          sprintf(buf,"Error in map: Unknown boundary condition '%c'",bstr[i]);
          barf("%s", buf);
        }
        break;
      }
    }
  }
 }

/***
 * Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
 */
 big    = labs((PDLA_Long)   (SvNV(__privtrans->big)));
 if(big <= 0)
   barf("%s","map: 'big' parameter must be >0");

 blur   = fabs((PDLA_Double) (SvNV(__privtrans->blur)));
 if(blur < 0)
   barf("%s","map: 'blur' parameter must be >= 0");

 sv_min = fabs((PDLA_Double) (SvNV(__privtrans->sv_min)));
 if(sv_min < 0)
   barf("%s","map: 'sv_min' parameter must be >= 0");

 flux   = (SvNV(__privtrans->flux) != 0);

 {
  char *mstr;
  STRLEN mlen;
  mstr = SvPV(__privtrans->method,mlen);

  if(mlen==0)
        method = 'h';
  else switch(*mstr) {
    case 'H':                       method='H'; break;
    case 'h':                       method = 'h';
      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;
      hanning_offset = (blur >= 1) ?
                0 :
                0.5 * (1.0 - blur);
      break;

    case 'g': case 'j':             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':             method = 'G'; break;
    default:
      {
        char err[80];
        sprintf(err,"Bug in map: unknown method '%c'",*mstr);
        barf("%s", err);
      }
      break;
    }
 }



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

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

 map_ptr = (PDLA_Double *)(map->data);


/* Main pixel loop (iterates over pixels in the output plane) */
 do {
   PDLA_Indx psize;    PDLA_Indx i_off;    PDLA_Indx j;        char t_vio;        char carry;
   /* Prefrobnicate the transformation matrix */
   psize = (PDLA_Long)(blur * PDLA_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */

#ifdef DEBUG_MAP
 {
     int k; PDLA_Indx foo = 0;
     printf("ovec: [");
     for(k=0;k
       foo += ovec[k] * map->dimincs[k+1];
       printf(" %2d ",(int)(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)(((PDLA_Byte *)(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.
      */
     {
       PDLA_Double *mp = map_ptr;
       for (i=0;i
        dvec[i]  = *mp - ( ibvec[i] = (PDLA_Long)(*mp + 0.5));   /* assignment */
        mp += map->dimincs[0];
       }
     }

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

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

     }
     {
       PDLA_Double *wg = wgt;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
     }
     {
       PDLA_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 */
           barf("%s","index out-of-bounds in map");
           break;
         case 1: /* truncation */
           t_vio++;
           /* fall through */
         case 2: /* extension -- crop */
           if(j<0)
             j=0;
           else              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:
           barf("%s", "Unknown boundary condition in map -- bug alert!");
           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 */
         PDLA_Double *cp;
         PDLA_Double alpha;
       /* 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 ) {

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

         switch(method) {
           PDLA_Double dd;
         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
             int lodex;
             int hidex;
             PDLA_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;                   if(lodex > HANNING_LOOKUP_SIZE)
                 lodex = HANNING_LOOKUP_SIZE;
                 hidex = lodex+1;
                                  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. */
	  {
	    PDLA_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 {
	      int lodex,hidex;
	      PDLA_Double beta = fabs(zeta * sum);

	      lodex = beta;
	      beta -= lodex; 	      hidex = lodex+1;
	                    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. */
              {
                 PDLA_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:
           {
             char buf[80];
             sprintf(buf,"This can't happen: method='%c'",method);
             barf("%s", buf);
           }
         }

         {  /* 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 threadloop here, so we accumulate each threaded element separately. */
            PDLA_Byte *dat = ((PDLA_Byte *)(in->data)) + i_off;
            PDLA_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;              }
          }
       } /* 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 */
           int 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) */

     {
       PDLA_Double *ac = acc;
       PDLA_Double *wg = wgt;
       PDLA_Double *wg2 = wgt2;
       PDLA_Byte *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++) {
	   if(*wg && (*wg2 / *wg) < 1.5 ) {
              *dat = *(ac++) / *(wg++);
	      wg2++;
           } else {
	      *dat = badval;
	      ac++; wg++; wg2++;
	   }
              dat += out->dimincs[ndims];
         }
       } else {
         /* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
         PDLA_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 */
     PDLA_Byte *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 */
   {
     for(i=0;
             (i
             (map_ptr += map->dimincs[i+1]) &&  /* Funky pre-test increment */
             (++(ovec[i]) >= out->dims[i]);     /* Actual carry test */
         i++) {
        ovec[i] = 0;
        map_ptr -= out->dims[i] * map->dimincs[i+1];
     }
    }
  } while(i



}
	PDLA_COMMENT("THREADLOOPEND")
	}
	}
	k0_datap -= __tinc1_0 * __tdims1 + __offsp[0];
	} while(PDLA->iterthreadloop(&__privtrans->__pdlthread,2));	} break; case PDLA_S: {
	PDLA_Short * k0_datap = ((PDLA_Short *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Short * k0_physdatap = ((PDLA_Short *)((__privtrans->pdls[0])->data));


	PDLA_COMMENT("THREADLOOPBEGIN")
	if ( PDLA->startthreadloop(&(__privtrans->__pdlthread),__privtrans->vtable->readdata, __tr) ) return;
           do { register PDLA_Indx __tind1=0,__tind2=0;
                register PDLA_Indx __tnpdls = __privtrans->__pdlthread.npdls;
                register PDLA_Indx __tdims1 = __privtrans->__pdlthread.dims[1];
                register PDLA_Indx __tdims0 = __privtrans->__pdlthread.dims[0];
                register PDLA_Indx *__offsp = PDLA->get_threadoffsp(&__privtrans->__pdlthread);
	register PDLA_Indx __tinc0_0 = __privtrans->__pdlthread.incs[0];
	register PDLA_Indx __tinc1_0 = __privtrans->__pdlthread.incs[__tnpdls+0];
	k0_datap += __offsp[0];
	for( __tind2 = 0 ;
                __tind2 < __tdims1 ;
                __tind2++
			,k0_datap += __tinc1_0 - __tinc0_0 * __tdims0
	)
	{
	     for( __tind1 = 0 ;
                  __tind1 < __tdims0 ;
                  __tind1++
			,k0_datap += __tinc0_0
	)
	{  PDLA_COMMENT("This is the tightest threadloop. Make sure inside is optimal."){

/*
 * 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 threadloops; threading
 * is handled internally.  To simplify the threading business, any
 * thread dimensions should all be collapsed to a single one by the
 * perl front-end.
 *
 */

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

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

 pdl *in  = PDLA->SvPDLAV(__privtrans->in);
 pdl *out = PDLA->SvPDLAV(__privtrans->out);
 pdl *map = PDLA->SvPDLAV(__privtrans->map);

 PDLA->make_physical(in);
 PDLA->make_physical(out);
 PDLA->make_physical(map);

 ndims = map->ndims -1;

 /*
  * Allocate all our dynamic workspaces at once
  * */
 ovec = (PDLA_Indx *)(PDLA->smalloc( (STRLEN)
                                   ( + sizeof(PDLA_Indx)    * 3 * ndims                                                    + sizeof(PDLA_Double) * (3*ndims)                                                     + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * 3 * ndims*ndims + ndims                                       + sizeof(char) * ndims                                                               + sizeof(PDLA_Indx)   * ndims                                                         )
                                   )
                );
 ivec   =                   &(ovec[ndims]);
 ibvec  =                   &(ivec[ndims]);
 dvec   =    (PDLA_Double *)(&(ibvec[ndims]));
 tvec   =                   &(dvec[ndims]);
 acc    =                   &(tvec[ndims]);
 wgt    =                   &(acc[in->dims[ndims]]);   wgt2   =                   &(wgt[in->dims[ndims]]);   tmp    =                   &(wgt2[in->dims[ndims]]);
 bounds =          (char *)(&(tmp [3*ndims*ndims+ndims]));
 index_stash = (PDLA_Indx *) &(bounds[ndims]);


/***
 * Fill in the boundary condition array
 */
 {
  char *bstr;
  STRLEN blen;
  bstr = SvPV(__privtrans->boundary,blen);

  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    int i;
    for (i=0;i
      bounds[i] = 1;
  } else {
    int 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:
        {
          char buf[BUFSIZ];
          sprintf(buf,"Error in map: Unknown boundary condition '%c'",bstr[i]);
          barf("%s", buf);
        }
        break;
      }
    }
  }
 }

/***
 * Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
 */
 big    = labs((PDLA_Long)   (SvNV(__privtrans->big)));
 if(big <= 0)
   barf("%s","map: 'big' parameter must be >0");

 blur   = fabs((PDLA_Double) (SvNV(__privtrans->blur)));
 if(blur < 0)
   barf("%s","map: 'blur' parameter must be >= 0");

 sv_min = fabs((PDLA_Double) (SvNV(__privtrans->sv_min)));
 if(sv_min < 0)
   barf("%s","map: 'sv_min' parameter must be >= 0");

 flux   = (SvNV(__privtrans->flux) != 0);

 {
  char *mstr;
  STRLEN mlen;
  mstr = SvPV(__privtrans->method,mlen);

  if(mlen==0)
        method = 'h';
  else switch(*mstr) {
    case 'H':                       method='H'; break;
    case 'h':                       method = 'h';
      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;
      hanning_offset = (blur >= 1) ?
                0 :
                0.5 * (1.0 - blur);
      break;

    case 'g': case 'j':             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':             method = 'G'; break;
    default:
      {
        char err[80];
        sprintf(err,"Bug in map: unknown method '%c'",*mstr);
        barf("%s", err);
      }
      break;
    }
 }



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

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

 map_ptr = (PDLA_Double *)(map->data);


/* Main pixel loop (iterates over pixels in the output plane) */
 do {
   PDLA_Indx psize;    PDLA_Indx i_off;    PDLA_Indx j;        char t_vio;        char carry;
   /* Prefrobnicate the transformation matrix */
   psize = (PDLA_Long)(blur * PDLA_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */

#ifdef DEBUG_MAP
 {
     int k; PDLA_Indx foo = 0;
     printf("ovec: [");
     for(k=0;k
       foo += ovec[k] * map->dimincs[k+1];
       printf(" %2d ",(int)(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)(((PDLA_Short *)(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.
      */
     {
       PDLA_Double *mp = map_ptr;
       for (i=0;i
        dvec[i]  = *mp - ( ibvec[i] = (PDLA_Long)(*mp + 0.5));   /* assignment */
        mp += map->dimincs[0];
       }
     }

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

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

     }
     {
       PDLA_Double *wg = wgt;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
     }
     {
       PDLA_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 */
           barf("%s","index out-of-bounds in map");
           break;
         case 1: /* truncation */
           t_vio++;
           /* fall through */
         case 2: /* extension -- crop */
           if(j<0)
             j=0;
           else              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:
           barf("%s", "Unknown boundary condition in map -- bug alert!");
           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 */
         PDLA_Double *cp;
         PDLA_Double alpha;
       /* 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 ) {

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

         switch(method) {
           PDLA_Double dd;
         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
             int lodex;
             int hidex;
             PDLA_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;                   if(lodex > HANNING_LOOKUP_SIZE)
                 lodex = HANNING_LOOKUP_SIZE;
                 hidex = lodex+1;
                                  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. */
	  {
	    PDLA_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 {
	      int lodex,hidex;
	      PDLA_Double beta = fabs(zeta * sum);

	      lodex = beta;
	      beta -= lodex; 	      hidex = lodex+1;
	                    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. */
              {
                 PDLA_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:
           {
             char buf[80];
             sprintf(buf,"This can't happen: method='%c'",method);
             barf("%s", buf);
           }
         }

         {  /* 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 threadloop here, so we accumulate each threaded element separately. */
            PDLA_Short *dat = ((PDLA_Short *)(in->data)) + i_off;
            PDLA_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;              }
          }
       } /* 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 */
           int 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) */

     {
       PDLA_Double *ac = acc;
       PDLA_Double *wg = wgt;
       PDLA_Double *wg2 = wgt2;
       PDLA_Short *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++) {
	   if(*wg && (*wg2 / *wg) < 1.5 ) {
              *dat = *(ac++) / *(wg++);
	      wg2++;
           } else {
	      *dat = badval;
	      ac++; wg++; wg2++;
	   }
              dat += out->dimincs[ndims];
         }
       } else {
         /* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
         PDLA_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 */
     PDLA_Short *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 */
   {
     for(i=0;
             (i
             (map_ptr += map->dimincs[i+1]) &&  /* Funky pre-test increment */
             (++(ovec[i]) >= out->dims[i]);     /* Actual carry test */
         i++) {
        ovec[i] = 0;
        map_ptr -= out->dims[i] * map->dimincs[i+1];
     }
    }
  } while(i



}
	PDLA_COMMENT("THREADLOOPEND")
	}
	}
	k0_datap -= __tinc1_0 * __tdims1 + __offsp[0];
	} while(PDLA->iterthreadloop(&__privtrans->__pdlthread,2));	} break; case PDLA_US: {
	PDLA_Ushort * k0_datap = ((PDLA_Ushort *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Ushort * k0_physdatap = ((PDLA_Ushort *)((__privtrans->pdls[0])->data));


	PDLA_COMMENT("THREADLOOPBEGIN")
	if ( PDLA->startthreadloop(&(__privtrans->__pdlthread),__privtrans->vtable->readdata, __tr) ) return;
           do { register PDLA_Indx __tind1=0,__tind2=0;
                register PDLA_Indx __tnpdls = __privtrans->__pdlthread.npdls;
                register PDLA_Indx __tdims1 = __privtrans->__pdlthread.dims[1];
                register PDLA_Indx __tdims0 = __privtrans->__pdlthread.dims[0];
                register PDLA_Indx *__offsp = PDLA->get_threadoffsp(&__privtrans->__pdlthread);
	register PDLA_Indx __tinc0_0 = __privtrans->__pdlthread.incs[0];
	register PDLA_Indx __tinc1_0 = __privtrans->__pdlthread.incs[__tnpdls+0];
	k0_datap += __offsp[0];
	for( __tind2 = 0 ;
                __tind2 < __tdims1 ;
                __tind2++
			,k0_datap += __tinc1_0 - __tinc0_0 * __tdims0
	)
	{
	     for( __tind1 = 0 ;
                  __tind1 < __tdims0 ;
                  __tind1++
			,k0_datap += __tinc0_0
	)
	{  PDLA_COMMENT("This is the tightest threadloop. Make sure inside is optimal."){

/*
 * 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 threadloops; threading
 * is handled internally.  To simplify the threading business, any
 * thread dimensions should all be collapsed to a single one by the
 * perl front-end.
 *
 */

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

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

 pdl *in  = PDLA->SvPDLAV(__privtrans->in);
 pdl *out = PDLA->SvPDLAV(__privtrans->out);
 pdl *map = PDLA->SvPDLAV(__privtrans->map);

 PDLA->make_physical(in);
 PDLA->make_physical(out);
 PDLA->make_physical(map);

 ndims = map->ndims -1;

 /*
  * Allocate all our dynamic workspaces at once
  * */
 ovec = (PDLA_Indx *)(PDLA->smalloc( (STRLEN)
                                   ( + sizeof(PDLA_Indx)    * 3 * ndims                                                    + sizeof(PDLA_Double) * (3*ndims)                                                     + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * 3 * ndims*ndims + ndims                                       + sizeof(char) * ndims                                                               + sizeof(PDLA_Indx)   * ndims                                                         )
                                   )
                );
 ivec   =                   &(ovec[ndims]);
 ibvec  =                   &(ivec[ndims]);
 dvec   =    (PDLA_Double *)(&(ibvec[ndims]));
 tvec   =                   &(dvec[ndims]);
 acc    =                   &(tvec[ndims]);
 wgt    =                   &(acc[in->dims[ndims]]);   wgt2   =                   &(wgt[in->dims[ndims]]);   tmp    =                   &(wgt2[in->dims[ndims]]);
 bounds =          (char *)(&(tmp [3*ndims*ndims+ndims]));
 index_stash = (PDLA_Indx *) &(bounds[ndims]);


/***
 * Fill in the boundary condition array
 */
 {
  char *bstr;
  STRLEN blen;
  bstr = SvPV(__privtrans->boundary,blen);

  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    int i;
    for (i=0;i
      bounds[i] = 1;
  } else {
    int 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:
        {
          char buf[BUFSIZ];
          sprintf(buf,"Error in map: Unknown boundary condition '%c'",bstr[i]);
          barf("%s", buf);
        }
        break;
      }
    }
  }
 }

/***
 * Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
 */
 big    = labs((PDLA_Long)   (SvNV(__privtrans->big)));
 if(big <= 0)
   barf("%s","map: 'big' parameter must be >0");

 blur   = fabs((PDLA_Double) (SvNV(__privtrans->blur)));
 if(blur < 0)
   barf("%s","map: 'blur' parameter must be >= 0");

 sv_min = fabs((PDLA_Double) (SvNV(__privtrans->sv_min)));
 if(sv_min < 0)
   barf("%s","map: 'sv_min' parameter must be >= 0");

 flux   = (SvNV(__privtrans->flux) != 0);

 {
  char *mstr;
  STRLEN mlen;
  mstr = SvPV(__privtrans->method,mlen);

  if(mlen==0)
        method = 'h';
  else switch(*mstr) {
    case 'H':                       method='H'; break;
    case 'h':                       method = 'h';
      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;
      hanning_offset = (blur >= 1) ?
                0 :
                0.5 * (1.0 - blur);
      break;

    case 'g': case 'j':             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':             method = 'G'; break;
    default:
      {
        char err[80];
        sprintf(err,"Bug in map: unknown method '%c'",*mstr);
        barf("%s", err);
      }
      break;
    }
 }



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

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

 map_ptr = (PDLA_Double *)(map->data);


/* Main pixel loop (iterates over pixels in the output plane) */
 do {
   PDLA_Indx psize;    PDLA_Indx i_off;    PDLA_Indx j;        char t_vio;        char carry;
   /* Prefrobnicate the transformation matrix */
   psize = (PDLA_Long)(blur * PDLA_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */

#ifdef DEBUG_MAP
 {
     int k; PDLA_Indx foo = 0;
     printf("ovec: [");
     for(k=0;k
       foo += ovec[k] * map->dimincs[k+1];
       printf(" %2d ",(int)(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)(((PDLA_Ushort *)(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.
      */
     {
       PDLA_Double *mp = map_ptr;
       for (i=0;i
        dvec[i]  = *mp - ( ibvec[i] = (PDLA_Long)(*mp + 0.5));   /* assignment */
        mp += map->dimincs[0];
       }
     }

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

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

     }
     {
       PDLA_Double *wg = wgt;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
     }
     {
       PDLA_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 */
           barf("%s","index out-of-bounds in map");
           break;
         case 1: /* truncation */
           t_vio++;
           /* fall through */
         case 2: /* extension -- crop */
           if(j<0)
             j=0;
           else              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:
           barf("%s", "Unknown boundary condition in map -- bug alert!");
           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 */
         PDLA_Double *cp;
         PDLA_Double alpha;
       /* 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 ) {

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

         switch(method) {
           PDLA_Double dd;
         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
             int lodex;
             int hidex;
             PDLA_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;                   if(lodex > HANNING_LOOKUP_SIZE)
                 lodex = HANNING_LOOKUP_SIZE;
                 hidex = lodex+1;
                                  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. */
	  {
	    PDLA_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 {
	      int lodex,hidex;
	      PDLA_Double beta = fabs(zeta * sum);

	      lodex = beta;
	      beta -= lodex; 	      hidex = lodex+1;
	                    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. */
              {
                 PDLA_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:
           {
             char buf[80];
             sprintf(buf,"This can't happen: method='%c'",method);
             barf("%s", buf);
           }
         }

         {  /* 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 threadloop here, so we accumulate each threaded element separately. */
            PDLA_Ushort *dat = ((PDLA_Ushort *)(in->data)) + i_off;
            PDLA_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;              }
          }
       } /* 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 */
           int 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) */

     {
       PDLA_Double *ac = acc;
       PDLA_Double *wg = wgt;
       PDLA_Double *wg2 = wgt2;
       PDLA_Ushort *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++) {
	   if(*wg && (*wg2 / *wg) < 1.5 ) {
              *dat = *(ac++) / *(wg++);
	      wg2++;
           } else {
	      *dat = badval;
	      ac++; wg++; wg2++;
	   }
              dat += out->dimincs[ndims];
         }
       } else {
         /* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
         PDLA_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 */
     PDLA_Ushort *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 */
   {
     for(i=0;
             (i
             (map_ptr += map->dimincs[i+1]) &&  /* Funky pre-test increment */
             (++(ovec[i]) >= out->dims[i]);     /* Actual carry test */
         i++) {
        ovec[i] = 0;
        map_ptr -= out->dims[i] * map->dimincs[i+1];
     }
    }
  } while(i



}
	PDLA_COMMENT("THREADLOOPEND")
	}
	}
	k0_datap -= __tinc1_0 * __tdims1 + __offsp[0];
	} while(PDLA->iterthreadloop(&__privtrans->__pdlthread,2));	} break; case PDLA_L: {
	PDLA_Long * k0_datap = ((PDLA_Long *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Long * k0_physdatap = ((PDLA_Long *)((__privtrans->pdls[0])->data));


	PDLA_COMMENT("THREADLOOPBEGIN")
	if ( PDLA->startthreadloop(&(__privtrans->__pdlthread),__privtrans->vtable->readdata, __tr) ) return;
           do { register PDLA_Indx __tind1=0,__tind2=0;
                register PDLA_Indx __tnpdls = __privtrans->__pdlthread.npdls;
                register PDLA_Indx __tdims1 = __privtrans->__pdlthread.dims[1];
                register PDLA_Indx __tdims0 = __privtrans->__pdlthread.dims[0];
                register PDLA_Indx *__offsp = PDLA->get_threadoffsp(&__privtrans->__pdlthread);
	register PDLA_Indx __tinc0_0 = __privtrans->__pdlthread.incs[0];
	register PDLA_Indx __tinc1_0 = __privtrans->__pdlthread.incs[__tnpdls+0];
	k0_datap += __offsp[0];
	for( __tind2 = 0 ;
                __tind2 < __tdims1 ;
                __tind2++
			,k0_datap += __tinc1_0 - __tinc0_0 * __tdims0
	)
	{
	     for( __tind1 = 0 ;
                  __tind1 < __tdims0 ;
                  __tind1++
			,k0_datap += __tinc0_0
	)
	{  PDLA_COMMENT("This is the tightest threadloop. Make sure inside is optimal."){

/*
 * 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 threadloops; threading
 * is handled internally.  To simplify the threading business, any
 * thread dimensions should all be collapsed to a single one by the
 * perl front-end.
 *
 */

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

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

 pdl *in  = PDLA->SvPDLAV(__privtrans->in);
 pdl *out = PDLA->SvPDLAV(__privtrans->out);
 pdl *map = PDLA->SvPDLAV(__privtrans->map);

 PDLA->make_physical(in);
 PDLA->make_physical(out);
 PDLA->make_physical(map);

 ndims = map->ndims -1;

 /*
  * Allocate all our dynamic workspaces at once
  * */
 ovec = (PDLA_Indx *)(PDLA->smalloc( (STRLEN)
                                   ( + sizeof(PDLA_Indx)    * 3 * ndims                                                    + sizeof(PDLA_Double) * (3*ndims)                                                     + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * 3 * ndims*ndims + ndims                                       + sizeof(char) * ndims                                                               + sizeof(PDLA_Indx)   * ndims                                                         )
                                   )
                );
 ivec   =                   &(ovec[ndims]);
 ibvec  =                   &(ivec[ndims]);
 dvec   =    (PDLA_Double *)(&(ibvec[ndims]));
 tvec   =                   &(dvec[ndims]);
 acc    =                   &(tvec[ndims]);
 wgt    =                   &(acc[in->dims[ndims]]);   wgt2   =                   &(wgt[in->dims[ndims]]);   tmp    =                   &(wgt2[in->dims[ndims]]);
 bounds =          (char *)(&(tmp [3*ndims*ndims+ndims]));
 index_stash = (PDLA_Indx *) &(bounds[ndims]);


/***
 * Fill in the boundary condition array
 */
 {
  char *bstr;
  STRLEN blen;
  bstr = SvPV(__privtrans->boundary,blen);

  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    int i;
    for (i=0;i
      bounds[i] = 1;
  } else {
    int 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:
        {
          char buf[BUFSIZ];
          sprintf(buf,"Error in map: Unknown boundary condition '%c'",bstr[i]);
          barf("%s", buf);
        }
        break;
      }
    }
  }
 }

/***
 * Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
 */
 big    = labs((PDLA_Long)   (SvNV(__privtrans->big)));
 if(big <= 0)
   barf("%s","map: 'big' parameter must be >0");

 blur   = fabs((PDLA_Double) (SvNV(__privtrans->blur)));
 if(blur < 0)
   barf("%s","map: 'blur' parameter must be >= 0");

 sv_min = fabs((PDLA_Double) (SvNV(__privtrans->sv_min)));
 if(sv_min < 0)
   barf("%s","map: 'sv_min' parameter must be >= 0");

 flux   = (SvNV(__privtrans->flux) != 0);

 {
  char *mstr;
  STRLEN mlen;
  mstr = SvPV(__privtrans->method,mlen);

  if(mlen==0)
        method = 'h';
  else switch(*mstr) {
    case 'H':                       method='H'; break;
    case 'h':                       method = 'h';
      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;
      hanning_offset = (blur >= 1) ?
                0 :
                0.5 * (1.0 - blur);
      break;

    case 'g': case 'j':             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':             method = 'G'; break;
    default:
      {
        char err[80];
        sprintf(err,"Bug in map: unknown method '%c'",*mstr);
        barf("%s", err);
      }
      break;
    }
 }



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

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

 map_ptr = (PDLA_Double *)(map->data);


/* Main pixel loop (iterates over pixels in the output plane) */
 do {
   PDLA_Indx psize;    PDLA_Indx i_off;    PDLA_Indx j;        char t_vio;        char carry;
   /* Prefrobnicate the transformation matrix */
   psize = (PDLA_Long)(blur * PDLA_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */

#ifdef DEBUG_MAP
 {
     int k; PDLA_Indx foo = 0;
     printf("ovec: [");
     for(k=0;k
       foo += ovec[k] * map->dimincs[k+1];
       printf(" %2d ",(int)(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)(((PDLA_Long *)(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.
      */
     {
       PDLA_Double *mp = map_ptr;
       for (i=0;i
        dvec[i]  = *mp - ( ibvec[i] = (PDLA_Long)(*mp + 0.5));   /* assignment */
        mp += map->dimincs[0];
       }
     }

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

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

     }
     {
       PDLA_Double *wg = wgt;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
     }
     {
       PDLA_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 */
           barf("%s","index out-of-bounds in map");
           break;
         case 1: /* truncation */
           t_vio++;
           /* fall through */
         case 2: /* extension -- crop */
           if(j<0)
             j=0;
           else              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:
           barf("%s", "Unknown boundary condition in map -- bug alert!");
           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 */
         PDLA_Double *cp;
         PDLA_Double alpha;
       /* 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 ) {

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

         switch(method) {
           PDLA_Double dd;
         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
             int lodex;
             int hidex;
             PDLA_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;                   if(lodex > HANNING_LOOKUP_SIZE)
                 lodex = HANNING_LOOKUP_SIZE;
                 hidex = lodex+1;
                                  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. */
	  {
	    PDLA_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 {
	      int lodex,hidex;
	      PDLA_Double beta = fabs(zeta * sum);

	      lodex = beta;
	      beta -= lodex; 	      hidex = lodex+1;
	                    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. */
              {
                 PDLA_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:
           {
             char buf[80];
             sprintf(buf,"This can't happen: method='%c'",method);
             barf("%s", buf);
           }
         }

         {  /* 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 threadloop here, so we accumulate each threaded element separately. */
            PDLA_Long *dat = ((PDLA_Long *)(in->data)) + i_off;
            PDLA_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;              }
          }
       } /* 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 */
           int 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) */

     {
       PDLA_Double *ac = acc;
       PDLA_Double *wg = wgt;
       PDLA_Double *wg2 = wgt2;
       PDLA_Long *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++) {
	   if(*wg && (*wg2 / *wg) < 1.5 ) {
              *dat = *(ac++) / *(wg++);
	      wg2++;
           } else {
	      *dat = badval;
	      ac++; wg++; wg2++;
	   }
              dat += out->dimincs[ndims];
         }
       } else {
         /* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
         PDLA_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 */
     PDLA_Long *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 */
   {
     for(i=0;
             (i
             (map_ptr += map->dimincs[i+1]) &&  /* Funky pre-test increment */
             (++(ovec[i]) >= out->dims[i]);     /* Actual carry test */
         i++) {
        ovec[i] = 0;
        map_ptr -= out->dims[i] * map->dimincs[i+1];
     }
    }
  } while(i



}
	PDLA_COMMENT("THREADLOOPEND")
	}
	}
	k0_datap -= __tinc1_0 * __tdims1 + __offsp[0];
	} while(PDLA->iterthreadloop(&__privtrans->__pdlthread,2));	} break; case PDLA_IND: {
	PDLA_Indx * k0_datap = ((PDLA_Indx *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Indx * k0_physdatap = ((PDLA_Indx *)((__privtrans->pdls[0])->data));


	PDLA_COMMENT("THREADLOOPBEGIN")
	if ( PDLA->startthreadloop(&(__privtrans->__pdlthread),__privtrans->vtable->readdata, __tr) ) return;
           do { register PDLA_Indx __tind1=0,__tind2=0;
                register PDLA_Indx __tnpdls = __privtrans->__pdlthread.npdls;
                register PDLA_Indx __tdims1 = __privtrans->__pdlthread.dims[1];
                register PDLA_Indx __tdims0 = __privtrans->__pdlthread.dims[0];
                register PDLA_Indx *__offsp = PDLA->get_threadoffsp(&__privtrans->__pdlthread);
	register PDLA_Indx __tinc0_0 = __privtrans->__pdlthread.incs[0];
	register PDLA_Indx __tinc1_0 = __privtrans->__pdlthread.incs[__tnpdls+0];
	k0_datap += __offsp[0];
	for( __tind2 = 0 ;
                __tind2 < __tdims1 ;
                __tind2++
			,k0_datap += __tinc1_0 - __tinc0_0 * __tdims0
	)
	{
	     for( __tind1 = 0 ;
                  __tind1 < __tdims0 ;
                  __tind1++
			,k0_datap += __tinc0_0
	)
	{  PDLA_COMMENT("This is the tightest threadloop. Make sure inside is optimal."){

/*
 * 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 threadloops; threading
 * is handled internally.  To simplify the threading business, any
 * thread dimensions should all be collapsed to a single one by the
 * perl front-end.
 *
 */

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

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

 pdl *in  = PDLA->SvPDLAV(__privtrans->in);
 pdl *out = PDLA->SvPDLAV(__privtrans->out);
 pdl *map = PDLA->SvPDLAV(__privtrans->map);

 PDLA->make_physical(in);
 PDLA->make_physical(out);
 PDLA->make_physical(map);

 ndims = map->ndims -1;

 /*
  * Allocate all our dynamic workspaces at once
  * */
 ovec = (PDLA_Indx *)(PDLA->smalloc( (STRLEN)
                                   ( + sizeof(PDLA_Indx)    * 3 * ndims                                                    + sizeof(PDLA_Double) * (3*ndims)                                                     + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * 3 * ndims*ndims + ndims                                       + sizeof(char) * ndims                                                               + sizeof(PDLA_Indx)   * ndims                                                         )
                                   )
                );
 ivec   =                   &(ovec[ndims]);
 ibvec  =                   &(ivec[ndims]);
 dvec   =    (PDLA_Double *)(&(ibvec[ndims]));
 tvec   =                   &(dvec[ndims]);
 acc    =                   &(tvec[ndims]);
 wgt    =                   &(acc[in->dims[ndims]]);   wgt2   =                   &(wgt[in->dims[ndims]]);   tmp    =                   &(wgt2[in->dims[ndims]]);
 bounds =          (char *)(&(tmp [3*ndims*ndims+ndims]));
 index_stash = (PDLA_Indx *) &(bounds[ndims]);


/***
 * Fill in the boundary condition array
 */
 {
  char *bstr;
  STRLEN blen;
  bstr = SvPV(__privtrans->boundary,blen);

  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    int i;
    for (i=0;i
      bounds[i] = 1;
  } else {
    int 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:
        {
          char buf[BUFSIZ];
          sprintf(buf,"Error in map: Unknown boundary condition '%c'",bstr[i]);
          barf("%s", buf);
        }
        break;
      }
    }
  }
 }

/***
 * Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
 */
 big    = labs((PDLA_Long)   (SvNV(__privtrans->big)));
 if(big <= 0)
   barf("%s","map: 'big' parameter must be >0");

 blur   = fabs((PDLA_Double) (SvNV(__privtrans->blur)));
 if(blur < 0)
   barf("%s","map: 'blur' parameter must be >= 0");

 sv_min = fabs((PDLA_Double) (SvNV(__privtrans->sv_min)));
 if(sv_min < 0)
   barf("%s","map: 'sv_min' parameter must be >= 0");

 flux   = (SvNV(__privtrans->flux) != 0);

 {
  char *mstr;
  STRLEN mlen;
  mstr = SvPV(__privtrans->method,mlen);

  if(mlen==0)
        method = 'h';
  else switch(*mstr) {
    case 'H':                       method='H'; break;
    case 'h':                       method = 'h';
      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;
      hanning_offset = (blur >= 1) ?
                0 :
                0.5 * (1.0 - blur);
      break;

    case 'g': case 'j':             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':             method = 'G'; break;
    default:
      {
        char err[80];
        sprintf(err,"Bug in map: unknown method '%c'",*mstr);
        barf("%s", err);
      }
      break;
    }
 }



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

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

 map_ptr = (PDLA_Double *)(map->data);


/* Main pixel loop (iterates over pixels in the output plane) */
 do {
   PDLA_Indx psize;    PDLA_Indx i_off;    PDLA_Indx j;        char t_vio;        char carry;
   /* Prefrobnicate the transformation matrix */
   psize = (PDLA_Long)(blur * PDLA_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */

#ifdef DEBUG_MAP
 {
     int k; PDLA_Indx foo = 0;
     printf("ovec: [");
     for(k=0;k
       foo += ovec[k] * map->dimincs[k+1];
       printf(" %2d ",(int)(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)(((PDLA_Indx *)(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.
      */
     {
       PDLA_Double *mp = map_ptr;
       for (i=0;i
        dvec[i]  = *mp - ( ibvec[i] = (PDLA_Long)(*mp + 0.5));   /* assignment */
        mp += map->dimincs[0];
       }
     }

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

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

     }
     {
       PDLA_Double *wg = wgt;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
     }
     {
       PDLA_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 */
           barf("%s","index out-of-bounds in map");
           break;
         case 1: /* truncation */
           t_vio++;
           /* fall through */
         case 2: /* extension -- crop */
           if(j<0)
             j=0;
           else              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:
           barf("%s", "Unknown boundary condition in map -- bug alert!");
           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 */
         PDLA_Double *cp;
         PDLA_Double alpha;
       /* 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 ) {

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

         switch(method) {
           PDLA_Double dd;
         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
             int lodex;
             int hidex;
             PDLA_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;                   if(lodex > HANNING_LOOKUP_SIZE)
                 lodex = HANNING_LOOKUP_SIZE;
                 hidex = lodex+1;
                                  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. */
	  {
	    PDLA_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 {
	      int lodex,hidex;
	      PDLA_Double beta = fabs(zeta * sum);

	      lodex = beta;
	      beta -= lodex; 	      hidex = lodex+1;
	                    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. */
              {
                 PDLA_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:
           {
             char buf[80];
             sprintf(buf,"This can't happen: method='%c'",method);
             barf("%s", buf);
           }
         }

         {  /* 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 threadloop here, so we accumulate each threaded element separately. */
            PDLA_Indx *dat = ((PDLA_Indx *)(in->data)) + i_off;
            PDLA_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;              }
          }
       } /* 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 */
           int 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) */

     {
       PDLA_Double *ac = acc;
       PDLA_Double *wg = wgt;
       PDLA_Double *wg2 = wgt2;
       PDLA_Indx *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++) {
	   if(*wg && (*wg2 / *wg) < 1.5 ) {
              *dat = *(ac++) / *(wg++);
	      wg2++;
           } else {
	      *dat = badval;
	      ac++; wg++; wg2++;
	   }
              dat += out->dimincs[ndims];
         }
       } else {
         /* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
         PDLA_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 */
     PDLA_Indx *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 */
   {
     for(i=0;
             (i
             (map_ptr += map->dimincs[i+1]) &&  /* Funky pre-test increment */
             (++(ovec[i]) >= out->dims[i]);     /* Actual carry test */
         i++) {
        ovec[i] = 0;
        map_ptr -= out->dims[i] * map->dimincs[i+1];
     }
    }
  } while(i



}
	PDLA_COMMENT("THREADLOOPEND")
	}
	}
	k0_datap -= __tinc1_0 * __tdims1 + __offsp[0];
	} while(PDLA->iterthreadloop(&__privtrans->__pdlthread,2));	} break; case PDLA_LL: {
	PDLA_LongLong * k0_datap = ((PDLA_LongLong *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_LongLong * k0_physdatap = ((PDLA_LongLong *)((__privtrans->pdls[0])->data));


	PDLA_COMMENT("THREADLOOPBEGIN")
	if ( PDLA->startthreadloop(&(__privtrans->__pdlthread),__privtrans->vtable->readdata, __tr) ) return;
           do { register PDLA_Indx __tind1=0,__tind2=0;
                register PDLA_Indx __tnpdls = __privtrans->__pdlthread.npdls;
                register PDLA_Indx __tdims1 = __privtrans->__pdlthread.dims[1];
                register PDLA_Indx __tdims0 = __privtrans->__pdlthread.dims[0];
                register PDLA_Indx *__offsp = PDLA->get_threadoffsp(&__privtrans->__pdlthread);
	register PDLA_Indx __tinc0_0 = __privtrans->__pdlthread.incs[0];
	register PDLA_Indx __tinc1_0 = __privtrans->__pdlthread.incs[__tnpdls+0];
	k0_datap += __offsp[0];
	for( __tind2 = 0 ;
                __tind2 < __tdims1 ;
                __tind2++
			,k0_datap += __tinc1_0 - __tinc0_0 * __tdims0
	)
	{
	     for( __tind1 = 0 ;
                  __tind1 < __tdims0 ;
                  __tind1++
			,k0_datap += __tinc0_0
	)
	{  PDLA_COMMENT("This is the tightest threadloop. Make sure inside is optimal."){

/*
 * 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 threadloops; threading
 * is handled internally.  To simplify the threading business, any
 * thread dimensions should all be collapsed to a single one by the
 * perl front-end.
 *
 */

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

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

 pdl *in  = PDLA->SvPDLAV(__privtrans->in);
 pdl *out = PDLA->SvPDLAV(__privtrans->out);
 pdl *map = PDLA->SvPDLAV(__privtrans->map);

 PDLA->make_physical(in);
 PDLA->make_physical(out);
 PDLA->make_physical(map);

 ndims = map->ndims -1;

 /*
  * Allocate all our dynamic workspaces at once
  * */
 ovec = (PDLA_Indx *)(PDLA->smalloc( (STRLEN)
                                   ( + sizeof(PDLA_Indx)    * 3 * ndims                                                    + sizeof(PDLA_Double) * (3*ndims)                                                     + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * 3 * ndims*ndims + ndims                                       + sizeof(char) * ndims                                                               + sizeof(PDLA_Indx)   * ndims                                                         )
                                   )
                );
 ivec   =                   &(ovec[ndims]);
 ibvec  =                   &(ivec[ndims]);
 dvec   =    (PDLA_Double *)(&(ibvec[ndims]));
 tvec   =                   &(dvec[ndims]);
 acc    =                   &(tvec[ndims]);
 wgt    =                   &(acc[in->dims[ndims]]);   wgt2   =                   &(wgt[in->dims[ndims]]);   tmp    =                   &(wgt2[in->dims[ndims]]);
 bounds =          (char *)(&(tmp [3*ndims*ndims+ndims]));
 index_stash = (PDLA_Indx *) &(bounds[ndims]);


/***
 * Fill in the boundary condition array
 */
 {
  char *bstr;
  STRLEN blen;
  bstr = SvPV(__privtrans->boundary,blen);

  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    int i;
    for (i=0;i
      bounds[i] = 1;
  } else {
    int 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:
        {
          char buf[BUFSIZ];
          sprintf(buf,"Error in map: Unknown boundary condition '%c'",bstr[i]);
          barf("%s", buf);
        }
        break;
      }
    }
  }
 }

/***
 * Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
 */
 big    = labs((PDLA_Long)   (SvNV(__privtrans->big)));
 if(big <= 0)
   barf("%s","map: 'big' parameter must be >0");

 blur   = fabs((PDLA_Double) (SvNV(__privtrans->blur)));
 if(blur < 0)
   barf("%s","map: 'blur' parameter must be >= 0");

 sv_min = fabs((PDLA_Double) (SvNV(__privtrans->sv_min)));
 if(sv_min < 0)
   barf("%s","map: 'sv_min' parameter must be >= 0");

 flux   = (SvNV(__privtrans->flux) != 0);

 {
  char *mstr;
  STRLEN mlen;
  mstr = SvPV(__privtrans->method,mlen);

  if(mlen==0)
        method = 'h';
  else switch(*mstr) {
    case 'H':                       method='H'; break;
    case 'h':                       method = 'h';
      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;
      hanning_offset = (blur >= 1) ?
                0 :
                0.5 * (1.0 - blur);
      break;

    case 'g': case 'j':             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':             method = 'G'; break;
    default:
      {
        char err[80];
        sprintf(err,"Bug in map: unknown method '%c'",*mstr);
        barf("%s", err);
      }
      break;
    }
 }



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

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

 map_ptr = (PDLA_Double *)(map->data);


/* Main pixel loop (iterates over pixels in the output plane) */
 do {
   PDLA_Indx psize;    PDLA_Indx i_off;    PDLA_Indx j;        char t_vio;        char carry;
   /* Prefrobnicate the transformation matrix */
   psize = (PDLA_Long)(blur * PDLA_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */

#ifdef DEBUG_MAP
 {
     int k; PDLA_Indx foo = 0;
     printf("ovec: [");
     for(k=0;k
       foo += ovec[k] * map->dimincs[k+1];
       printf(" %2d ",(int)(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)(((PDLA_LongLong *)(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.
      */
     {
       PDLA_Double *mp = map_ptr;
       for (i=0;i
        dvec[i]  = *mp - ( ibvec[i] = (PDLA_Long)(*mp + 0.5));   /* assignment */
        mp += map->dimincs[0];
       }
     }

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

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

     }
     {
       PDLA_Double *wg = wgt;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
     }
     {
       PDLA_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 */
           barf("%s","index out-of-bounds in map");
           break;
         case 1: /* truncation */
           t_vio++;
           /* fall through */
         case 2: /* extension -- crop */
           if(j<0)
             j=0;
           else              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:
           barf("%s", "Unknown boundary condition in map -- bug alert!");
           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 */
         PDLA_Double *cp;
         PDLA_Double alpha;
       /* 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 ) {

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

         switch(method) {
           PDLA_Double dd;
         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
             int lodex;
             int hidex;
             PDLA_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;                   if(lodex > HANNING_LOOKUP_SIZE)
                 lodex = HANNING_LOOKUP_SIZE;
                 hidex = lodex+1;
                                  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. */
	  {
	    PDLA_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 {
	      int lodex,hidex;
	      PDLA_Double beta = fabs(zeta * sum);

	      lodex = beta;
	      beta -= lodex; 	      hidex = lodex+1;
	                    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. */
              {
                 PDLA_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:
           {
             char buf[80];
             sprintf(buf,"This can't happen: method='%c'",method);
             barf("%s", buf);
           }
         }

         {  /* 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 threadloop here, so we accumulate each threaded element separately. */
            PDLA_LongLong *dat = ((PDLA_LongLong *)(in->data)) + i_off;
            PDLA_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;              }
          }
       } /* 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 */
           int 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) */

     {
       PDLA_Double *ac = acc;
       PDLA_Double *wg = wgt;
       PDLA_Double *wg2 = wgt2;
       PDLA_LongLong *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++) {
	   if(*wg && (*wg2 / *wg) < 1.5 ) {
              *dat = *(ac++) / *(wg++);
	      wg2++;
           } else {
	      *dat = badval;
	      ac++; wg++; wg2++;
	   }
              dat += out->dimincs[ndims];
         }
       } else {
         /* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
         PDLA_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 */
     PDLA_LongLong *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 */
   {
     for(i=0;
             (i
             (map_ptr += map->dimincs[i+1]) &&  /* Funky pre-test increment */
             (++(ovec[i]) >= out->dims[i]);     /* Actual carry test */
         i++) {
        ovec[i] = 0;
        map_ptr -= out->dims[i] * map->dimincs[i+1];
     }
    }
  } while(i



}
	PDLA_COMMENT("THREADLOOPEND")
	}
	}
	k0_datap -= __tinc1_0 * __tdims1 + __offsp[0];
	} while(PDLA->iterthreadloop(&__privtrans->__pdlthread,2));	} break; case PDLA_F: {
	PDLA_Float * k0_datap = ((PDLA_Float *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Float * k0_physdatap = ((PDLA_Float *)((__privtrans->pdls[0])->data));


	PDLA_COMMENT("THREADLOOPBEGIN")
	if ( PDLA->startthreadloop(&(__privtrans->__pdlthread),__privtrans->vtable->readdata, __tr) ) return;
           do { register PDLA_Indx __tind1=0,__tind2=0;
                register PDLA_Indx __tnpdls = __privtrans->__pdlthread.npdls;
                register PDLA_Indx __tdims1 = __privtrans->__pdlthread.dims[1];
                register PDLA_Indx __tdims0 = __privtrans->__pdlthread.dims[0];
                register PDLA_Indx *__offsp = PDLA->get_threadoffsp(&__privtrans->__pdlthread);
	register PDLA_Indx __tinc0_0 = __privtrans->__pdlthread.incs[0];
	register PDLA_Indx __tinc1_0 = __privtrans->__pdlthread.incs[__tnpdls+0];
	k0_datap += __offsp[0];
	for( __tind2 = 0 ;
                __tind2 < __tdims1 ;
                __tind2++
			,k0_datap += __tinc1_0 - __tinc0_0 * __tdims0
	)
	{
	     for( __tind1 = 0 ;
                  __tind1 < __tdims0 ;
                  __tind1++
			,k0_datap += __tinc0_0
	)
	{  PDLA_COMMENT("This is the tightest threadloop. Make sure inside is optimal."){

/*
 * 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 threadloops; threading
 * is handled internally.  To simplify the threading business, any
 * thread dimensions should all be collapsed to a single one by the
 * perl front-end.
 *
 */

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

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

 pdl *in  = PDLA->SvPDLAV(__privtrans->in);
 pdl *out = PDLA->SvPDLAV(__privtrans->out);
 pdl *map = PDLA->SvPDLAV(__privtrans->map);

 PDLA->make_physical(in);
 PDLA->make_physical(out);
 PDLA->make_physical(map);

 ndims = map->ndims -1;

 /*
  * Allocate all our dynamic workspaces at once
  * */
 ovec = (PDLA_Indx *)(PDLA->smalloc( (STRLEN)
                                   ( + sizeof(PDLA_Indx)    * 3 * ndims                                                    + sizeof(PDLA_Double) * (3*ndims)                                                     + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * 3 * ndims*ndims + ndims                                       + sizeof(char) * ndims                                                               + sizeof(PDLA_Indx)   * ndims                                                         )
                                   )
                );
 ivec   =                   &(ovec[ndims]);
 ibvec  =                   &(ivec[ndims]);
 dvec   =    (PDLA_Double *)(&(ibvec[ndims]));
 tvec   =                   &(dvec[ndims]);
 acc    =                   &(tvec[ndims]);
 wgt    =                   &(acc[in->dims[ndims]]);   wgt2   =                   &(wgt[in->dims[ndims]]);   tmp    =                   &(wgt2[in->dims[ndims]]);
 bounds =          (char *)(&(tmp [3*ndims*ndims+ndims]));
 index_stash = (PDLA_Indx *) &(bounds[ndims]);


/***
 * Fill in the boundary condition array
 */
 {
  char *bstr;
  STRLEN blen;
  bstr = SvPV(__privtrans->boundary,blen);

  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    int i;
    for (i=0;i
      bounds[i] = 1;
  } else {
    int 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:
        {
          char buf[BUFSIZ];
          sprintf(buf,"Error in map: Unknown boundary condition '%c'",bstr[i]);
          barf("%s", buf);
        }
        break;
      }
    }
  }
 }

/***
 * Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
 */
 big    = labs((PDLA_Long)   (SvNV(__privtrans->big)));
 if(big <= 0)
   barf("%s","map: 'big' parameter must be >0");

 blur   = fabs((PDLA_Double) (SvNV(__privtrans->blur)));
 if(blur < 0)
   barf("%s","map: 'blur' parameter must be >= 0");

 sv_min = fabs((PDLA_Double) (SvNV(__privtrans->sv_min)));
 if(sv_min < 0)
   barf("%s","map: 'sv_min' parameter must be >= 0");

 flux   = (SvNV(__privtrans->flux) != 0);

 {
  char *mstr;
  STRLEN mlen;
  mstr = SvPV(__privtrans->method,mlen);

  if(mlen==0)
        method = 'h';
  else switch(*mstr) {
    case 'H':                       method='H'; break;
    case 'h':                       method = 'h';
      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;
      hanning_offset = (blur >= 1) ?
                0 :
                0.5 * (1.0 - blur);
      break;

    case 'g': case 'j':             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':             method = 'G'; break;
    default:
      {
        char err[80];
        sprintf(err,"Bug in map: unknown method '%c'",*mstr);
        barf("%s", err);
      }
      break;
    }
 }



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

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

 map_ptr = (PDLA_Double *)(map->data);


/* Main pixel loop (iterates over pixels in the output plane) */
 do {
   PDLA_Indx psize;    PDLA_Indx i_off;    PDLA_Indx j;        char t_vio;        char carry;
   /* Prefrobnicate the transformation matrix */
   psize = (PDLA_Long)(blur * PDLA_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */

#ifdef DEBUG_MAP
 {
     int k; PDLA_Indx foo = 0;
     printf("ovec: [");
     for(k=0;k
       foo += ovec[k] * map->dimincs[k+1];
       printf(" %2d ",(int)(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)(((PDLA_Float *)(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.
      */
     {
       PDLA_Double *mp = map_ptr;
       for (i=0;i
        dvec[i]  = *mp - ( ibvec[i] = (PDLA_Long)(*mp + 0.5));   /* assignment */
        mp += map->dimincs[0];
       }
     }

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

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

     }
     {
       PDLA_Double *wg = wgt;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
     }
     {
       PDLA_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 */
           barf("%s","index out-of-bounds in map");
           break;
         case 1: /* truncation */
           t_vio++;
           /* fall through */
         case 2: /* extension -- crop */
           if(j<0)
             j=0;
           else              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:
           barf("%s", "Unknown boundary condition in map -- bug alert!");
           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 */
         PDLA_Double *cp;
         PDLA_Double alpha;
       /* 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 ) {

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

         switch(method) {
           PDLA_Double dd;
         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
             int lodex;
             int hidex;
             PDLA_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;                   if(lodex > HANNING_LOOKUP_SIZE)
                 lodex = HANNING_LOOKUP_SIZE;
                 hidex = lodex+1;
                                  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. */
	  {
	    PDLA_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 {
	      int lodex,hidex;
	      PDLA_Double beta = fabs(zeta * sum);

	      lodex = beta;
	      beta -= lodex; 	      hidex = lodex+1;
	                    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. */
              {
                 PDLA_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:
           {
             char buf[80];
             sprintf(buf,"This can't happen: method='%c'",method);
             barf("%s", buf);
           }
         }

         {  /* 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 threadloop here, so we accumulate each threaded element separately. */
            PDLA_Float *dat = ((PDLA_Float *)(in->data)) + i_off;
            PDLA_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;              }
          }
       } /* 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 */
           int 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) */

     {
       PDLA_Double *ac = acc;
       PDLA_Double *wg = wgt;
       PDLA_Double *wg2 = wgt2;
       PDLA_Float *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++) {
	   if(*wg && (*wg2 / *wg) < 1.5 ) {
              *dat = *(ac++) / *(wg++);
	      wg2++;
           } else {
	      *dat = badval;
	      ac++; wg++; wg2++;
	   }
              dat += out->dimincs[ndims];
         }
       } else {
         /* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
         PDLA_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 */
     PDLA_Float *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 */
   {
     for(i=0;
             (i
             (map_ptr += map->dimincs[i+1]) &&  /* Funky pre-test increment */
             (++(ovec[i]) >= out->dims[i]);     /* Actual carry test */
         i++) {
        ovec[i] = 0;
        map_ptr -= out->dims[i] * map->dimincs[i+1];
     }
    }
  } while(i



}
	PDLA_COMMENT("THREADLOOPEND")
	}
	}
	k0_datap -= __tinc1_0 * __tdims1 + __offsp[0];
	} while(PDLA->iterthreadloop(&__privtrans->__pdlthread,2));	} break; case PDLA_D: {
	PDLA_Double * k0_datap = ((PDLA_Double *)(PDLA_REPRP_TRANS((__privtrans->pdls[0]),(__privtrans->vtable->per_pdl_flags[0]))));
	PDLA_Double * k0_physdatap = ((PDLA_Double *)((__privtrans->pdls[0])->data));


	PDLA_COMMENT("THREADLOOPBEGIN")
	if ( PDLA->startthreadloop(&(__privtrans->__pdlthread),__privtrans->vtable->readdata, __tr) ) return;
           do { register PDLA_Indx __tind1=0,__tind2=0;
                register PDLA_Indx __tnpdls = __privtrans->__pdlthread.npdls;
                register PDLA_Indx __tdims1 = __privtrans->__pdlthread.dims[1];
                register PDLA_Indx __tdims0 = __privtrans->__pdlthread.dims[0];
                register PDLA_Indx *__offsp = PDLA->get_threadoffsp(&__privtrans->__pdlthread);
	register PDLA_Indx __tinc0_0 = __privtrans->__pdlthread.incs[0];
	register PDLA_Indx __tinc1_0 = __privtrans->__pdlthread.incs[__tnpdls+0];
	k0_datap += __offsp[0];
	for( __tind2 = 0 ;
                __tind2 < __tdims1 ;
                __tind2++
			,k0_datap += __tinc1_0 - __tinc0_0 * __tdims0
	)
	{
	     for( __tind1 = 0 ;
                  __tind1 < __tdims0 ;
                  __tind1++
			,k0_datap += __tinc0_0
	)
	{  PDLA_COMMENT("This is the tightest threadloop. Make sure inside is optimal."){

/*
 * 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 threadloops; threading
 * is handled internally.  To simplify the threading business, any
 * thread dimensions should all be collapsed to a single one by the
 * perl front-end.
 *
 */

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

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

 pdl *in  = PDLA->SvPDLAV(__privtrans->in);
 pdl *out = PDLA->SvPDLAV(__privtrans->out);
 pdl *map = PDLA->SvPDLAV(__privtrans->map);

 PDLA->make_physical(in);
 PDLA->make_physical(out);
 PDLA->make_physical(map);

 ndims = map->ndims -1;

 /*
  * Allocate all our dynamic workspaces at once
  * */
 ovec = (PDLA_Indx *)(PDLA->smalloc( (STRLEN)
                                   ( + sizeof(PDLA_Indx)    * 3 * ndims                                                    + sizeof(PDLA_Double) * (3*ndims)                                                     + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * in->dims[ndims]                                               + sizeof(PDLA_Double) * 3 * ndims*ndims + ndims                                       + sizeof(char) * ndims                                                               + sizeof(PDLA_Indx)   * ndims                                                         )
                                   )
                );
 ivec   =                   &(ovec[ndims]);
 ibvec  =                   &(ivec[ndims]);
 dvec   =    (PDLA_Double *)(&(ibvec[ndims]));
 tvec   =                   &(dvec[ndims]);
 acc    =                   &(tvec[ndims]);
 wgt    =                   &(acc[in->dims[ndims]]);   wgt2   =                   &(wgt[in->dims[ndims]]);   tmp    =                   &(wgt2[in->dims[ndims]]);
 bounds =          (char *)(&(tmp [3*ndims*ndims+ndims]));
 index_stash = (PDLA_Indx *) &(bounds[ndims]);


/***
 * Fill in the boundary condition array
 */
 {
  char *bstr;
  STRLEN blen;
  bstr = SvPV(__privtrans->boundary,blen);

  if(blen == 0) {
    /* If no boundary is specified then every dim gets truncated */
    int i;
    for (i=0;i
      bounds[i] = 1;
  } else {
    int 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:
        {
          char buf[BUFSIZ];
          sprintf(buf,"Error in map: Unknown boundary condition '%c'",bstr[i]);
          barf("%s", buf);
        }
        break;
      }
    }
  }
 }

/***
 * Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
 */
 big    = labs((PDLA_Long)   (SvNV(__privtrans->big)));
 if(big <= 0)
   barf("%s","map: 'big' parameter must be >0");

 blur   = fabs((PDLA_Double) (SvNV(__privtrans->blur)));
 if(blur < 0)
   barf("%s","map: 'blur' parameter must be >= 0");

 sv_min = fabs((PDLA_Double) (SvNV(__privtrans->sv_min)));
 if(sv_min < 0)
   barf("%s","map: 'sv_min' parameter must be >= 0");

 flux   = (SvNV(__privtrans->flux) != 0);

 {
  char *mstr;
  STRLEN mlen;
  mstr = SvPV(__privtrans->method,mlen);

  if(mlen==0)
        method = 'h';
  else switch(*mstr) {
    case 'H':                       method='H'; break;
    case 'h':                       method = 'h';
      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;
      hanning_offset = (blur >= 1) ?
                0 :
                0.5 * (1.0 - blur);
      break;

    case 'g': case 'j':             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':             method = 'G'; break;
    default:
      {
        char err[80];
        sprintf(err,"Bug in map: unknown method '%c'",*mstr);
        barf("%s", err);
      }
      break;
    }
 }



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

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

 map_ptr = (PDLA_Double *)(map->data);


/* Main pixel loop (iterates over pixels in the output plane) */
 do {
   PDLA_Indx psize;    PDLA_Indx i_off;    PDLA_Indx j;        char t_vio;        char carry;
   /* Prefrobnicate the transformation matrix */
   psize = (PDLA_Long)(blur * PDLA_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */

#ifdef DEBUG_MAP
 {
     int k; PDLA_Indx foo = 0;
     printf("ovec: [");
     for(k=0;k
       foo += ovec[k] * map->dimincs[k+1];
       printf(" %2d ",(int)(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)(((PDLA_Double *)(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.
      */
     {
       PDLA_Double *mp = map_ptr;
       for (i=0;i
        dvec[i]  = *mp - ( ibvec[i] = (PDLA_Long)(*mp + 0.5));   /* assignment */
        mp += map->dimincs[0];
       }
     }

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

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

     }
     {
       PDLA_Double *wg = wgt;
       for(i=0;i < in->dims[ndims]; i++)
         *(wg++) = 0.0;
     }
     {
       PDLA_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 */
           barf("%s","index out-of-bounds in map");
           break;
         case 1: /* truncation */
           t_vio++;
           /* fall through */
         case 2: /* extension -- crop */
           if(j<0)
             j=0;
           else              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:
           barf("%s", "Unknown boundary condition in map -- bug alert!");
           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 */
         PDLA_Double *cp;
         PDLA_Double alpha;
       /* 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 ) {

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

         switch(method) {
           PDLA_Double dd;
         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
             int lodex;
             int hidex;
             PDLA_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;                   if(lodex > HANNING_LOOKUP_SIZE)
                 lodex = HANNING_LOOKUP_SIZE;
                 hidex = lodex+1;
                                  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. */
	  {
	    PDLA_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 {
	      int lodex,hidex;
	      PDLA_Double beta = fabs(zeta * sum);

	      lodex = beta;
	      beta -= lodex; 	      hidex = lodex+1;
	                    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. */
              {
                 PDLA_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:
           {
             char buf[80];
             sprintf(buf,"This can't happen: method='%c'",method);
             barf("%s", buf);
           }
         }

         {  /* 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 threadloop here, so we accumulate each threaded element separately. */
            PDLA_Double *dat = ((PDLA_Double *)(in->data)) + i_off;
            PDLA_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;              }
          }
       } /* 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 */
           int 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) */

     {
       PDLA_Double *ac = acc;
       PDLA_Double *wg = wgt;
       PDLA_Double *wg2 = wgt2;
       PDLA_Double *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++) {
	   if(*wg && (*wg2 / *wg) < 1.5 ) {
              *dat = *(ac++) / *(wg++);
	      wg2++;
           } else {
	      *dat = badval;
	      ac++; wg++; wg2++;
	   }
              dat += out->dimincs[ndims];
         }
       } else {
         /* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
         PDLA_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 */
     PDLA_Double *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 */
   {
     for(i=0;
             (i
             (map_ptr += map->dimincs[i+1]) &&  /* Funky pre-test increment */
             (++(ovec[i]) >= out->dims[i]);     /* Actual carry test */
         i++) {
        ovec[i] = 0;
        map_ptr -= out->dims[i] * map->dimincs[i+1];
     }
    }
  } while(i



}
	PDLA_COMMENT("THREADLOOPEND")
	}
	}
	k0_datap -= __tinc1_0 * __tdims1 + __offsp[0];
	} while(PDLA->iterthreadloop(&__privtrans->__pdlthread,2));	break;}
	default:barf("PP INTERNAL ERROR! PLEASE MAKE A BUG REPORT\n");}
}
	}
    }




void pdl_map_free(pdl_trans *__tr ) {
	int __dim;
	pdl_map_struct *__privtrans = (pdl_map_struct *) __tr;

	{

			PDLA_TR_CLRMAGIC(__privtrans);
			SvREFCNT_dec(__privtrans->in);;SvREFCNT_dec(__privtrans->out);;SvREFCNT_dec(__privtrans->map);;SvREFCNT_dec(__privtrans->boundary);;SvREFCNT_dec(__privtrans->method);;SvREFCNT_dec(__privtrans->big);;SvREFCNT_dec(__privtrans->blur);;SvREFCNT_dec(__privtrans->sv_min);;SvREFCNT_dec(__privtrans->flux);;SvREFCNT_dec(__privtrans->bv);;
			if(__privtrans->__ddone) {
				PDLA->freethreadloop(&(__privtrans->__pdlthread));
				;
			}

	}
    }




static char pdl_map_vtable_flags[] =
	 	{ PDLA_TPDLA_VAFFINE_OK};
	 pdl_transvtable pdl_map_vtable = {
		0,0, 1, 1, pdl_map_vtable_flags,
		pdl_map_redodims, pdl_map_readdata, NULL,
		pdl_map_free,NULL,NULL,pdl_map_copy,
		sizeof(pdl_map_struct),"pdl_map_vtable"
	 };



MODULE = PDLA::Transform PACKAGE = PDLA::Transform

PROTOTYPES: ENABLE

int
set_debugging(i)
	int i;
	CODE:
	RETVAL = __pdl_debugging;
	__pdl_debugging = i;
	OUTPUT:
	RETVAL

int
set_boundscheck(i)
       int i;
       CODE:
       if (! 1)
         warn("Bounds checking is disabled for PDLA::Transform");
       RETVAL = __pdl_boundscheck;
       __pdl_boundscheck = i;
       OUTPUT:
       RETVAL


MODULE = PDLA::Transform PACKAGE = PDLA


void
_map_int(k0,in,out,map,boundary,method,big,blur,sv_min,flux,bv)
	pdl  *k0
	SV  *in
	SV  *out
	SV  *map
	SV  *boundary
	SV  *method
	SV  *big
	SV  *blur
	SV  *sv_min
	SV  *flux
	SV  *bv
 CODE:
 { pdl_map_struct  *__privtrans;
  int badflag_cache = 0;
   __privtrans = malloc(sizeof(*__privtrans));
   PDLA_THR_CLRMAGIC(&__privtrans->__pdlthread);
   PDLA_TR_SETMAGIC(__privtrans);
   __privtrans->flags = 0;
   __privtrans->__ddone = 0;
   __privtrans->vtable = &pdl_map_vtable;
   __privtrans->freeproc = PDLA->trans_mallocfreeproc;
__privtrans->bvalflag = 0;
  badflag_cache = ((k0->state & PDLA_BADVAL) > 0);
  if (badflag_cache) __privtrans->bvalflag = 1;
__privtrans->__datatype = 0;if(__privtrans->__datatype < k0->datatype) {
		 	__privtrans->__datatype = k0->datatype;
		    }
if(__privtrans->__datatype == PDLA_B) {}
else if(__privtrans->__datatype == PDLA_S) {}
else if(__privtrans->__datatype == PDLA_U) {}
else if(__privtrans->__datatype == PDLA_L) {}
else if(__privtrans->__datatype == PDLA_N) {}
else if(__privtrans->__datatype == PDLA_Q) {}
else if(__privtrans->__datatype == PDLA_F) {}
else if(__privtrans->__datatype == PDLA_D) {}
else __privtrans->__datatype = PDLA_D;
if(__privtrans->__datatype != k0->datatype) {
	     k0 = PDLA->get_convertedpdl(k0,__privtrans->__datatype);
	  }{(__privtrans->in) = newSVsv(in);(__privtrans->out) = newSVsv(out);(__privtrans->map) = newSVsv(map);(__privtrans->boundary) = newSVsv(boundary);(__privtrans->method) = newSVsv(method);(__privtrans->big) = newSVsv(big);(__privtrans->blur) = newSVsv(blur);(__privtrans->sv_min) = newSVsv(sv_min);(__privtrans->flux) = newSVsv(flux);(__privtrans->bv) = newSVsv(bv);}PDLA_COMMENT("No flow")__privtrans->__pdlthread.inds = 0;__privtrans->pdls[0] = k0;
PDLA->make_trans_mutual((pdl_trans *)__privtrans);
if (badflag_cache) {
}
 XSRETURN(0);
}


BOOT:

   PDLA_COMMENT("Get pointer to structure of core shared C routines")
   PDLA_COMMENT("make sure PDLA::Core is loaded")

   perl_require_pv ("PDLA/Core.pm"); /* make sure PDLA::Core is loaded */
#ifndef aTHX_
#define aTHX_
#endif
   if (SvTRUE (ERRSV)) Perl_croak(aTHX_ "%s",SvPV_nolen (ERRSV));
   CoreSV = perl_get_sv("PDLA::SHARE",FALSE);  /* SV* value */
   if (CoreSV==NULL)
     Perl_croak(aTHX_ "We require the PDLA::Core module, which was not found");
   PDLA = INT2PTR(Core*,SvIV( CoreSV ));  /* Core* value */
   if (PDLA->Version != PDLA_CORE_VERSION)
     Perl_croak(aTHX_ "[PDLA->Version: %d PDLA_CORE_VERSION: %d XS_VERSION: %s] PDLA::Transform needs to be recompiled against the newly installed PDLA", PDLA->Version, PDLA_CORE_VERSION, XS_VERSION);