File Coverage

lib/PDL/Transform-pp-map.c

Criterion	Covered	Total	%
statement	40	91	43.9
branch	34	300	11.3
condition			n/a
subroutine			n/a
pod			n/a
total	74	391	18.9

line	stmt	bran	code
1
2			#line 453 "lib/PDL/PP.pm"
3			/*
4			* THIS FILE WAS GENERATED BY PDL::PP from lib/PDL/Transform.pd! Do not modify!
5			*/
6
7			#define PDL_FREE_CODE(trans, destroy, comp_free_code, ntpriv_free_code) \
8			if (destroy) { \
9			comp_free_code \
10			} \
11			if ((trans)->dims_redone) { \
12			ntpriv_free_code \
13			}
14
15			#include "EXTERN.h"
16			#include "perl.h"
17			#include "XSUB.h"
18			#include "pdl.h"
19			#include "pdlcore.h"
20			#define PDL PDL_Transform
21			extern Core* PDL; /* Structure hold core C functions */
22			#line 23 "lib/PDL/Transform-pp-map.c"
23			#include
24
25			#line 418 "lib/PDL/Transform.pd"
26			/*
27			* Singular-value decomposition code is borrowed from
28			* MatrixOps -- cut-and-pasted here because of linker trouble.
29			* It's used by the auxiliary matrix manipulation code, below.
30			*/
31			static inline void pdl_xform_svd(PDL_Double W, PDL_Double Z, PDL_Indx nRow, PDL_Indx nCol)
32			{
33			int i, j, k, EstColRank, RotCount, SweepCount, slimit;
34			PDL_Double eps, e2, tol, vt, p, x0, y0, q, r, c0, s0, d1, d2;
35			eps = 1e-6;
36			slimit = nCol/4;
37			if (slimit < 6.0)
38			slimit = 6;
39			SweepCount = 0;
40			e2 = 10.0nRoweps*eps;
41			tol = eps*.1;
42			EstColRank = nCol;
43			for (i=0; i
44			for (j=0; j
45			W[nCol*(nRow+i)+j] = 0.0;
46			}
47			W[nCol(nRow+i)+i] = 1.0; / rjrw 7/7/99: moved this line out of j loop */
48			}
49			RotCount = EstColRank*(EstColRank-1)/2;
50			while (RotCount != 0 && SweepCount <= slimit)
51			{
52			RotCount = EstColRank*(EstColRank-1)/2;
53			SweepCount++;
54			for (j=0; j
55			{
56			for (k=j+1; k
57			{
58			p = q = r = 0.0;
59			for (i=0; i
60			{
61			x0 = W[nColi+j]; y0 = W[nColi+k];
62			p += x0y0; q += x0x0; r += y0*y0;
63			}
64			Z[j] = q; Z[k] = r;
65			if (q >= r)
66			{
67			if (q<=e2Z[0] \|\| fabs(p)<=tolq) RotCount--;
68			else
69			{
70			p /= q; r = 1 - r/q; vt = sqrt(4pp+r*r);
71			c0 = sqrt(fabs(.5(1+r/vt))); s0 = p/(vtc0);
72			for (i=0; i
73			{
74			d1 = W[nColi+j]; d2 = W[nColi+k];
75			W[nColi+j] = d1c0+d2s0; W[nColi+k] = -d1s0+d2c0;
76			}
77			}
78			}
79			else
80			{
81			p /= r; q = q/r-1; vt = sqrt(4pp+q*q);
82			s0 = sqrt(fabs(.5*(1-q/vt)));
83			if (p<0) s0 = -s0;
84			c0 = p/(vt*s0);
85			for (i=0; i
86			{
87			d1 = W[nColi+j]; d2 = W[nColi+k];
88			W[nColi+j] = d1c0+d2s0; W[nColi+k] = -d1s0+d2c0;
89			}
90			}
91			}
92			}
93			while (EstColRank>=3 && Z[(EstColRank-1)]<=Z[0]tol+toltol)
94			EstColRank--;
95			}
96			}
97
98			/*
99			* PDL_xform_aux:
100			* This handles the matrix manipulation part of the Jacobian filtered
101			* mapping code. It's separate from the main code because it's
102			* independent of the data type of the original arrays.
103			*
104			*Given a pre-allocated workspace and
105			* an integer set of coordinates, generate the discrete Jacobian
106			* from the map, pad the singular values, and return the inverse
107			* Jacobian, the largest singular value of the Jacobian itself, and
108			* the determinant of the original Jacobian. Boundary values use the
109			* asymmetric discrete Jacobian; others use the symmetric discrete Jacobian.
110			*
111			* The map and workspace must be of type PDL_D. If the dimensionality is
112			* d, then the workspace must have at least 3*n^2+n elements. The
113			* inverse of the padded Jacobian is returned in the first n^2 elements.
114			* The determinant of the original Jacobian gets stuffed into the n^2
115			* element of the workspace. The largest padded singular value is returned.
116			*
117			*/
118
119			static inline PDL_Double PDL_xform_aux ( pdl map, PDL_Indx coords, PDL_Double *tmp, PDL_Double sv_min) {
120			PDL_Indx ndims;
121			PDL_Indx i, j, k;
122			PDL_Indx offset;
123			PDL_Double det;
124			PDL_Double *jptr;
125			PDL_Double *svptr;
126			PDL_Double aptr,bptr;
127			PDL_Double max_sv = 0.0;
128
129			ndims = map->ndims-1;
130
131			/****** Accumulate the Jacobian */
132			/* Accumulate the offset into the map array */
133			for( i=offset=0; i
134			offset += coords[i] * map->dimincs[i+1];
135
136			jptr = tmp + ndims*ndims;
137
138			for( i=0; i
139			char bot = (coords[i] <=0);
140			char top = (coords[i] >= map->dims[i+1]-1);
141			char symmetric = !(bot \|\| top);
142			PDL_Double ohi,olo;
143			PDL_Indx diff = map->dimincs[i+1];
144			ohi = ((PDL_Double *)map->data) + ( offset + ( top ? 0 : diff ));
145			olo = ((PDL_Double *)map->data) + ( offset - ( bot ? 0 : diff ));
146
147			for( j=0; j
148			PDL_Double jel = ohi - olo;
149
150			ohi += map->dimincs[0];
151			olo += map->dimincs[0];
152
153			if(symmetric)
154			jel /= 2;
155			*(jptr++) = jel;
156			}
157			}
158
159
160			/****** Singular-value decompose the Jacobian
161			* The svd routine produces the squares of the singular values,
162			* and requires normalization for one of the rotation matrices.
163			*/
164
165			jptr = tmp + ndims*ndims;
166			svptr = tmp + 3ndimsndims;
167			pdl_xform_svd(jptr,svptr,ndims,ndims);
168
169			aptr = svptr;
170			for (i=0;i
171			aptr = sqrt(aptr);
172
173
174			/* fix-up matrices here */
175			bptr = jptr;
176			for(i=0; i
177			aptr = svptr;
178			for(j=0;j
179			(bptr++) /= (aptr++);
180			}
181
182			/****** Store the determinant, and pad the singular values as necessary.*/
183			aptr = svptr;
184			det = 1.0;
185			for(i=0;i
186			det = aptr;
187			if(*aptr < sv_min)
188			*aptr = sv_min;
189			if(*aptr > max_sv )
190			max_sv = *aptr;
191			aptr++;
192			}
193
194			/****** Generate the inverse matrix */
195			/* Multiply B-transpose times 1/S times A-transpose. */
196			/* since S is diagonal we just divide by the appropriate element. */
197			/* */
198			aptr = tmp + ndims*ndims;
199			bptr = aptr + ndims*ndims;
200			jptr= tmp;
201
202			for(i=0;i
203			for(j=0;j
204			*jptr = 0;
205
206			for(k=0;k
207			jptr += aptr[jndims + k] * bptr[kndims + i] / svptr;
208
209			jptr++;
210			}
211			svptr++;
212			}
213
214			*jptr = det;
215			return max_sv;
216			}
217
218			#line 1846 "lib/PDL/PP.pm"
219			typedef struct pdl_params_map {
220			#line 221 "lib/PDL/Transform-pp-map.c"
221			pdl *in;
222			pdl *out;
223			pdl *map;
224			SV *boundary;
225			SV *method;
226			long big;
227			double blur;
228			double sv_min;
229			char flux;
230			SV *bv;
231			} pdl_params_map;
232
233
234			#line 1857 "lib/PDL/PP.pm"
235			pdl_error pdl_map_readdata(pdl_trans *__privtrans) {
236			pdl_error PDL_err = {0, NULL, 0};
237			#line 238 "lib/PDL/Transform-pp-map.c"
238	12		pdl_params_map *__params = __privtrans->params; (void)__params;
239	12	50	if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in map:" "broadcast.incs NULL");
240			/* broadcastloop declarations */
241			int __brcloopval;
242			register PDL_Indx __tind0,__tind1; /* counters along dim */
243	12		register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
244			/* dims here are how many steps along those dims */
245	12		register PDL_Indx __tinc0_k0 = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
246	12		register PDL_Indx __tinc1_k0 = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
247			#define PDL_BROADCASTLOOP_START_map_readdata PDL_BROADCASTLOOP_START( \
248			readdata, \
249			__privtrans->broadcast, \
250			__privtrans->vtable, \
251			k0_datap += __offsp[0]; \
252			, \
253			( ,k0_datap += __tinc1_k0 - __tinc0_k0 * __tdims0 \
254			), \
255			( ,k0_datap += __tinc0_k0 \
256			) \
257			)
258			#define PDL_BROADCASTLOOP_END_map_readdata PDL_BROADCASTLOOP_END( \
259			__privtrans->broadcast, \
260			k0_datap -= __tinc1_k0 * __tdims1 + __offsp[0]; \
261			)
262			#ifndef PDL_DECLARE_PARAMS_map_1
263			#define PDL_DECLARE_PARAMS_map_1(PDL_TYPE_OP,PDL_PPSYM_OP) \
264			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, k0, (__privtrans->pdls[0]), 1, PDL_PPSYM_OP)
265			#endif
266			#define PDL_IF_BAD(t,f) f
267	12		switch (__privtrans->__datatype) { /* Start generic switch */
268	0		case PDL_SB: {
269	0	0	PDL_DECLARE_PARAMS_map_1(PDL_SByte,A)
		0
		0
270	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
271			#line 611 "lib/PDL/Transform.pd"
272
273			/*
274			* Pixel interpolation & averaging code
275			*
276			* Calls a common coordinate-transformation block (see following hdr)
277			* that isn't dependent on the type of the input variable.
278			*
279			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
280			* is handled internally. To simplify the broadcasting business, any
281			* broadcast dimensions should all be collapsed to a single one by the
282			* perl front-end.
283			*
284			*/
285
286			pdl *in = __params->in;
287			pdl *out = __params->out;
288			pdl *map = __params->map;
289
290			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
291			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
292			PDL_Indx ovec[ndims]; /* output pixel loop vector */
293			PDL_Indx ivec[ndims]; /* input pixel loop vector */
294			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
295			PDL_Double dvec[ndims]; /* Residual vector for linearization */
296			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
297			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
298			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
299			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
300			char bounds[ndims]; /* Boundary condition packed string */
301			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
302			char method; /* Method identifier (gets one of 'h','g') */
303			PDL_Long big; /* Max size of input footprint for each pix */
304			PDL_Double blur; /* Scaling of filter */
305			PDL_Double sv_min; /* minimum singular value */
306			char flux; /* Flag to indicate flux conservation */
307			PDL_Indx i;
308			PDL_SByte badval = SvNV(__params->bv);
309			#define HANNING_LOOKUP_SIZE 2500
310			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
311			static int needs_hanning_calc = 1;
312			PDL_Double zeta = 0;
313			PDL_Double hanning_offset;
314
315			#define GAUSSIAN_LOOKUP_SIZE 4000
316			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
317			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
318			static int needs_gaussian_calc = 1;
319
320			PDL_RETERROR(PDL_err, PDL->make_physical(in));
321			PDL_RETERROR(PDL_err, PDL->make_physical(out));
322			PDL_RETERROR(PDL_err, PDL->make_physical(map));
323
324			/***
325			* Fill in the boundary condition array
326			*/
327			{
328			char *bstr;
329			STRLEN blen;
330			bstr = SvPV(__params->boundary,blen);
331
332			if(blen == 0) {
333			/* If no boundary is specified then every dim gets truncated */
334			PDL_Indx i;
335			for (i=0;i
336			bounds[i] = 1;
337			} else {
338			PDL_Indx i;
339			for(i=0;i
340			switch(bstr[i < blen ? i : blen-1 ]) {
341			case '0': case 'f': case 'F': /* forbid */
342			bounds[i] = 0;
343			break;
344			case '1': case 't': case 'T': /* truncate */
345			bounds[i] = 1;
346			break;
347			case '2': case 'e': case 'E': /* extend */
348			bounds[i] = 2;
349			break;
350			case '3': case 'p': case 'P': /* periodic */
351			bounds[i] = 3;
352			break;
353			case '4': case 'm': case 'M': /* mirror */
354			bounds[i] = 4;
355			break;
356			default:
357			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
358			break;
359			}
360			}
361			}
362			}
363
364			/***
365			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
366			*/
367			big = labs(__params->big);
368			if(big <= 0)
369			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
370
371			blur = fabs(__params->blur);
372			if(blur < 0)
373			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
374
375			sv_min = fabs(__params->sv_min);
376			if(sv_min < 0)
377			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
378
379			flux = __params->flux != 0;
380
381			{
382			char *mstr;
383			STRLEN mlen;
384			mstr = SvPV(__params->method,mlen);
385
386			if(mlen==0)
387			method = 'h';
388			switch(*mstr) {
389			case 'h': if( needs_hanning_calc ) {
390			int i;
391			for(i=0;i
392			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
393			}
394			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
395			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
396			needs_hanning_calc = 0;
397			}
398			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
399			case 'H': method = *mstr;
400			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
401			break;
402
403			case 'g': case 'j': method = 'g';
404			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
405
406			if( needs_gaussian_calc ) {
407			int i;
408			for(i=0;i
409			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
410			}
411			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
412			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
413			needs_gaussian_calc = 0;
414			}
415			break;
416
417			case 'G': case 'J': method = 'G'; break;
418			default:
419			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
420			break;
421			}
422			}
423
424
425
426			/* End of initialization */
427			/*************************************************************/
428			/* Start of Real Work */
429
430			/* Initialize coordinate vector and map offset
431			*/
432			for(i=0;i
433			ovec[i] = 0;
434
435			PDL_Double map_ptr = (PDL_Double )(map->data);
436			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
437			for (i=0; i < npdls; i++) offs[i] = 0;
438			for (i=0; i < ndims; i++) {
439			incs[i*npdls + 0] = map->dimincs[i+1];
440			}
441
442			/* Main pixel loop (iterates over pixels in the output plane) */
443			do {
444			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
445			/* Prefrobnicate the transformation matrix */
446			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
447
448			#ifdef DEBUG_MAP
449			{
450			PDL_Indx k; PDL_Indx foo = 0;
451			printf("ovec: [");
452			for(k=0;k
453			foo += ovec[k] * map->dimincs[k+1];
454			printf(" %2td ",(PDL_Indx)(ovec[k]));
455			}
456			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
457			for(k=0;k
458			printf("%8.2f",(double)(((PDL_SByte )(map->data))[foo + kmap->dimincs[0]]));
459			}
460			printf("]\n");
461			}
462			#endif
463
464			/* Don't bother accumulating output if psize is too large */
465			if(psize <= big) {
466			/* Use the prefrobnicated matrix to generate a local linearization.
467			* dvec gets the delta; ibvec gets the base.
468			*/
469			{
470			PDL_Double *mp = map_ptr + offs[0];
471			for (i=0;i
472			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
473			mp += map->dimincs[0];
474			}
475			}
476
477			/* Initialize input delta vector */
478			for(i=0;i
479			ivec[i] = -psize;
480
481			/* Initialize accumulators */
482			{
483			PDL_Double *ac = acc;
484			for(i=0; i < in->dims[ndims]; i++)
485			*(ac++) = 0.0;
486
487			}
488			{
489			PDL_Double *wg = wgt;
490			for(i=0;i < in->dims[ndims]; i++)
491			*(wg++) = 0.0;
492			}
493			{
494			PDL_Double *wg = wgt2;
495			for(i=0;i < in->dims[ndims]; i++)
496			*(wg++) = 0.0;
497			}
498
499
500			/*
501			* Calculate the original offset into the data array, to enable
502			* delta calculations in the pixel loop
503			*
504			* i runs over dims; j holds the working integer index in the
505			* current dim.
506			*
507			* This code matches the incrementation code at the bottom of the accumulation loop
508			*/
509
510			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
511			i_off = 0;
512			for(i=0;i
513			j = ivec[i] + ibvec[i];
514			if(j<0 \|\| j >= in->dims[i]) {
515			switch(bounds[i]) {
516			case 0: /* no breakage allowed */
517			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
518			break;
519			case 1: /* truncation */
520			t_vio++;
521			/* fall through */
522			case 2: /* extension -- crop */
523			if(j<0)
524			j=0;
525			else j = in->dims[i] - 1;
526			break;
527			case 3: /* periodic -- mod it */
528			j %= in->dims[i];
529			if(j<0)
530			j += in->dims[i];
531			break;
532			case 4: /* mirror -- reflect off the edges */
533			j += in->dims[i];
534			j %= (in->dims[i]*2);
535			if(j<0)
536			j += in->dims[i]*2;
537			j -= in->dims[i];
538			if(j<0) {
539			j *= -1;
540			j -= 1;
541			}
542			break;
543			default:
544			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
545			break;
546			}
547			}
548			i_off += in->dimincs[i] * j;
549			}
550
551			/* Initialize index stashes for later reference as we scan the footprint */
552			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
553			/* end of a dimensional scan. So we stash the index location at the */
554			/* start of each dimensional scan here. When we finish incrementing */
555			/* through a particular dim, we pull its value back out of the stash. */
556			for(i=0;i
557			index_stash[i] = i_off;
558			}
559
560			/* The input accumulation loop is the hotspot for the entire operation. */
561			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
562			/* in the input array, use the linearized transform to bring each pixel center */
563			/* forward to the output plane, and calculate a weighting based on the chosen */
564			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
565			/* table that is initialized the first time through the code. 'H' is the */
566			/* same process, but explicitly calculated for each interation (~2x slower). */
567			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
568			/* into the input array fresh from the current input array vector each time, */
569			/* we walk through the array using dimincs and the old offset. This saves */
570			/* about half of the time spent on index calculation. */
571
572			do { /* Input accumulation loop */
573			PDL_Double *cp;
574			PDL_Double alpha;
575			/* Calculate the weight of the current input point. Don't bother if we're
576			* violating any truncation boundaries (in that case our value is zero, but
577			* for the interpolation we also set the weight to zero).
578			*/
579			if( !t_vio ) {
580
581			PDL_Double *ap = tvec;
582			PDL_Double *bp = dvec;
583			PDL_Indx *ip = ivec;
584			for(i=0; i
585			(ap++) = (ip++) - *(bp++);
586
587			switch(method) {
588			PDL_Double dd;
589			case 'h':
590			/* This is the Hanning window rolloff. It is a product of a simple */
591			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
592			/* is about 2x faster than using cos(theta) directly in each */
593			/* weighting calculation, so we do. Using 2500 entries and linear */
594			/* interpolation is accurate to about 10^-7, and should preserve */
595			/* the contents of cache pretty well. */
596			alpha = 1;
597			cp = tmp;
598			for (i=0; i
599			PDL_Indx lodex;
600			PDL_Indx hidex;
601			PDL_Double beta;
602			dd = 0;
603			ap = tvec;
604			/* Get the matrix-multiply element for this dimension */
605			for (j=0;j
606			dd += (ap++) *(cp++);
607
608			/* Do linear interpolation from the table */
609			/* The table captures a hanning window centered 0.5 pixel from center. */
610			/* We scale the filter by the blur parameter -- but if blur is less */
611			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
612			/* value on the pixel edge at 0.5. For blur greater than unity, we */
613			/* scale simply. */
614			beta = fabs(dd) - hanning_offset;
615			if (beta > 0) {
616			if (beta >= blur) {
617			alpha = 0;
618			i = ndims;
619			} else {
620			beta *= zeta;
621			lodex = beta;
622			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
623			lodex = HANNING_LOOKUP_SIZE;
624			hidex = lodex+1;
625			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
626			} /* end of interpolation branch */
627			} /* end of beta > 0 branch */
628			} /* end of dimension loop */
629			break;
630
631			case 'H':
632			/* This is the Hanning window rolloff with explicit calculation, preserved */
633			/* in case someone actually wants the slower longer method. */
634			alpha = 1;
635			cp = tmp;
636			for (i=0; i
637			dd = 0;
638			ap = tvec;
639			for (j=0;j
640			dd += (ap++) *(cp++);
641			dd = (fabs(dd) - hanning_offset) / blur;
642			if ( dd > 1 ) {
643			alpha = 0;
644			i = ndims;
645			} else
646			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
647			}
648			break;
649
650			case 'g':
651			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
652			{
653			PDL_Double sum = 0;
654			cp = tmp;
655			for(i=0; i
656			dd = 0;
657			ap = tvec;
658			for(j=0;j
659			dd += (ap++) *(cp++);
660			dd /= blur;
661			sum += dd * dd;
662			if(sum > GAUSSIAN_MAXVAL) {
663			i = ndims; /* exit early if we're too far out */
664			alpha = 0;
665			}
666			}
667			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
668			alpha = 0;
669			} else {
670			PDL_Indx lodex,hidex;
671			PDL_Double beta = fabs(zeta * sum);
672
673			lodex = beta;
674			beta -= lodex; hidex = lodex+1;
675			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
676
677			}
678			}
679			break;
680
681			case 'G':
682			/* This is the Gaussian rolloff with explicit calculation, preserved */
683			/* in case someone actually wants the slower longer method. */
684			{
685			PDL_Double sum = 0;
686			cp = tmp;
687			for(i=0; i
688			dd = 0;
689			ap = tvec;
690			for(j=0;j
691			dd += (ap++) *(cp++);
692			dd /= blur;
693			sum += dd * dd;
694			if(sum > 4) /* 2 pixels -- four half-widths */
695			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
696			}
697
698			if(sum > GAUSSIAN_MAXVAL)
699			alpha = 0;
700			else
701			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
702			}
703			break;
704			default:
705			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
706			}
707
708			{ /* convenience block -- accumulate the current point into the weighted sum. */
709			/* This is more than simple assignment because we have our own explicit poor */
710			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
711			PDL_SByte dat = ((PDL_SByte )(in->data)) + i_off;
712			PDL_Indx max = out->dims[ndims];
713			for( i=0; i < max; i++ ) {
714			if( (badval==0) \|\| (*dat != badval) ) {
715			acc[i] += dat alpha;
716			dat += in->dimincs[ndims];
717			wgt[i] += alpha;
718			}
719			wgt2[i] += alpha; }
720			}
721			} /* end of t_vio check (i.e. of input accumulation) */
722
723
724			/* Advance input accumulation loop. */
725			/* We both increment the total vector and also advance the index. */
726			carry = 1;
727			for (i=0; i
728			/* Advance the current element of the offset vector */
729			ivec[i]++;
730			j = ivec[i] + ibvec[i];
731
732			/* Advance the offset into the data array */
733			if( j > 0 && j <= in->dims[i]-1 ) {
734			/* Normal case -- just advance the input vector */
735			i_off += in->dimincs[i];
736			} else {
737			/* Busted a boundary - either before or after. */
738			switch(bounds[i]){
739			case 0: /* no breakage allowed -- treat as truncation for interpolation */
740			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
741			if( j == 0 )
742			t_vio--;
743			else if( j == in->dims[i] )
744			t_vio++;
745			break;
746			case 2: /* extension -- do nothing (so the same input point is re-used) */
747			break;
748			case 3: /* periodic -- advance and mod into the allowed range */
749			if((j % in->dims[i]) == 0) {
750			i_off -= in->dimincs[i] * (in->dims[i]-1);
751			} else {
752			i_off += in->dimincs[i];
753			}
754			break;
755			case 4: /* mirror -- advance or retreat depending on phase */
756			j += in->dims[i];
757			j %= (in->dims[i]*2);
758			j -= in->dims[i];
759			if( j!=0 && j!= -in->dims[i] ) {
760			if(j<0)
761			i_off -= in->dimincs[i];
762			else
763			i_off += in->dimincs[i];
764			}
765			break;
766			}
767			}
768
769			/* Now check for carry */
770			if (ivec[i] <= psize) {
771			/* Normal case -- copy the current offset to the faster-running dim stashes */
772			PDL_Indx k;
773			for(k=0;k
774			index_stash[k] = i_off;
775			}
776			carry = 0;
777
778			} else { /* End of this scan -- recover the last position, and mark carry */
779			i_off = index_stash[i];
780			if(bounds[i]==1) {
781			j = ivec[i] + ibvec[i];
782			if( j < 0 \|\| j >= in->dims[i] )
783			t_vio--;
784			ivec[i] = -psize;
785			j = ivec[i] + ibvec[i];
786			if( j < 0 \|\| j >= in->dims[i] )
787			t_vio++;
788			carry = 1;
789			} else {
790			ivec[i] = -psize;
791			}
792			}
793			} /* End of counter-advance loop */
794			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
795
796			{
797			PDL_Double *ac = acc;
798			PDL_Double *wg = wgt;
799			PDL_Double *wg2 = wgt2;
800			PDL_SByte *dat = out->data;
801
802			/* Calculate output vector offset */
803			for(i=0;i
804			dat += out->dimincs[i] * ovec[i];
805
806			if (!flux) {
807			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
808			for (i=0; i < out->dims[ndims]; i++) {
809			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
810			ac++; wg++; wg2++;
811			dat += out->dimincs[ndims];
812			}
813			} else {
814			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
815			PDL_Double det = tmp[ndims*ndims];
816			for(i=0; i < out->dims[ndims]; i++) {
817			if(wg && (wg2 / *wg) < 1.5 ) {
818			dat = (ac++) / (wg++) det;
819			wg2++;
820			} else {
821			*dat = badval;
822			ac++; wg++; wg2++;
823			}
824			dat += out->dimincs[ndims];
825			} /* end of for loop */
826			} /* end of flux flag set conditional */
827			} /* end of convenience block */
828
829			/* End of code for normal pixels */
830			} else {
831			/* The pixel was ludicrously huge -- just set this pixel to nan */
832			PDL_SByte *dat = out->data;
833			for(i=0;i
834			dat += out->dimincs[i] * ovec[i];
835			for(i=0;idims[ndims];i++) {
836			dat = badval; / Should handle bad values too -- not yet */
837			dat += out->dimincs[ndims];
838			}
839			}
840
841			/* Increment the pixel counter */
842			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
843			} while (1);
844			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
845			#line 846 "lib/PDL/Transform-pp-map.c"
846	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
847	0		} break;
848	0		case PDL_B: {
849	0	0	PDL_DECLARE_PARAMS_map_1(PDL_Byte,B)
		0
		0
850	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
851			#line 611 "lib/PDL/Transform.pd"
852
853			/*
854			* Pixel interpolation & averaging code
855			*
856			* Calls a common coordinate-transformation block (see following hdr)
857			* that isn't dependent on the type of the input variable.
858			*
859			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
860			* is handled internally. To simplify the broadcasting business, any
861			* broadcast dimensions should all be collapsed to a single one by the
862			* perl front-end.
863			*
864			*/
865
866			pdl *in = __params->in;
867			pdl *out = __params->out;
868			pdl *map = __params->map;
869
870			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
871			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
872			PDL_Indx ovec[ndims]; /* output pixel loop vector */
873			PDL_Indx ivec[ndims]; /* input pixel loop vector */
874			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
875			PDL_Double dvec[ndims]; /* Residual vector for linearization */
876			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
877			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
878			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
879			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
880			char bounds[ndims]; /* Boundary condition packed string */
881			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
882			char method; /* Method identifier (gets one of 'h','g') */
883			PDL_Long big; /* Max size of input footprint for each pix */
884			PDL_Double blur; /* Scaling of filter */
885			PDL_Double sv_min; /* minimum singular value */
886			char flux; /* Flag to indicate flux conservation */
887			PDL_Indx i;
888			PDL_Byte badval = SvNV(__params->bv);
889			#define HANNING_LOOKUP_SIZE 2500
890			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
891			static int needs_hanning_calc = 1;
892			PDL_Double zeta = 0;
893			PDL_Double hanning_offset;
894
895			#define GAUSSIAN_LOOKUP_SIZE 4000
896			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
897			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
898			static int needs_gaussian_calc = 1;
899
900			PDL_RETERROR(PDL_err, PDL->make_physical(in));
901			PDL_RETERROR(PDL_err, PDL->make_physical(out));
902			PDL_RETERROR(PDL_err, PDL->make_physical(map));
903
904			/***
905			* Fill in the boundary condition array
906			*/
907			{
908			char *bstr;
909			STRLEN blen;
910			bstr = SvPV(__params->boundary,blen);
911
912			if(blen == 0) {
913			/* If no boundary is specified then every dim gets truncated */
914			PDL_Indx i;
915			for (i=0;i
916			bounds[i] = 1;
917			} else {
918			PDL_Indx i;
919			for(i=0;i
920			switch(bstr[i < blen ? i : blen-1 ]) {
921			case '0': case 'f': case 'F': /* forbid */
922			bounds[i] = 0;
923			break;
924			case '1': case 't': case 'T': /* truncate */
925			bounds[i] = 1;
926			break;
927			case '2': case 'e': case 'E': /* extend */
928			bounds[i] = 2;
929			break;
930			case '3': case 'p': case 'P': /* periodic */
931			bounds[i] = 3;
932			break;
933			case '4': case 'm': case 'M': /* mirror */
934			bounds[i] = 4;
935			break;
936			default:
937			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
938			break;
939			}
940			}
941			}
942			}
943
944			/***
945			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
946			*/
947			big = labs(__params->big);
948			if(big <= 0)
949			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
950
951			blur = fabs(__params->blur);
952			if(blur < 0)
953			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
954
955			sv_min = fabs(__params->sv_min);
956			if(sv_min < 0)
957			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
958
959			flux = __params->flux != 0;
960
961			{
962			char *mstr;
963			STRLEN mlen;
964			mstr = SvPV(__params->method,mlen);
965
966			if(mlen==0)
967			method = 'h';
968			switch(*mstr) {
969			case 'h': if( needs_hanning_calc ) {
970			int i;
971			for(i=0;i
972			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
973			}
974			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
975			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
976			needs_hanning_calc = 0;
977			}
978			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
979			case 'H': method = *mstr;
980			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
981			break;
982
983			case 'g': case 'j': method = 'g';
984			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
985
986			if( needs_gaussian_calc ) {
987			int i;
988			for(i=0;i
989			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
990			}
991			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
992			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
993			needs_gaussian_calc = 0;
994			}
995			break;
996
997			case 'G': case 'J': method = 'G'; break;
998			default:
999			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
1000			break;
1001			}
1002			}
1003
1004
1005
1006			/* End of initialization */
1007			/*************************************************************/
1008			/* Start of Real Work */
1009
1010			/* Initialize coordinate vector and map offset
1011			*/
1012			for(i=0;i
1013			ovec[i] = 0;
1014
1015			PDL_Double map_ptr = (PDL_Double )(map->data);
1016			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
1017			for (i=0; i < npdls; i++) offs[i] = 0;
1018			for (i=0; i < ndims; i++) {
1019			incs[i*npdls + 0] = map->dimincs[i+1];
1020			}
1021
1022			/* Main pixel loop (iterates over pixels in the output plane) */
1023			do {
1024			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
1025			/* Prefrobnicate the transformation matrix */
1026			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
1027
1028			#ifdef DEBUG_MAP
1029			{
1030			PDL_Indx k; PDL_Indx foo = 0;
1031			printf("ovec: [");
1032			for(k=0;k
1033			foo += ovec[k] * map->dimincs[k+1];
1034			printf(" %2td ",(PDL_Indx)(ovec[k]));
1035			}
1036			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
1037			for(k=0;k
1038			printf("%8.2f",(double)(((PDL_Byte )(map->data))[foo + kmap->dimincs[0]]));
1039			}
1040			printf("]\n");
1041			}
1042			#endif
1043
1044			/* Don't bother accumulating output if psize is too large */
1045			if(psize <= big) {
1046			/* Use the prefrobnicated matrix to generate a local linearization.
1047			* dvec gets the delta; ibvec gets the base.
1048			*/
1049			{
1050			PDL_Double *mp = map_ptr + offs[0];
1051			for (i=0;i
1052			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
1053			mp += map->dimincs[0];
1054			}
1055			}
1056
1057			/* Initialize input delta vector */
1058			for(i=0;i
1059			ivec[i] = -psize;
1060
1061			/* Initialize accumulators */
1062			{
1063			PDL_Double *ac = acc;
1064			for(i=0; i < in->dims[ndims]; i++)
1065			*(ac++) = 0.0;
1066
1067			}
1068			{
1069			PDL_Double *wg = wgt;
1070			for(i=0;i < in->dims[ndims]; i++)
1071			*(wg++) = 0.0;
1072			}
1073			{
1074			PDL_Double *wg = wgt2;
1075			for(i=0;i < in->dims[ndims]; i++)
1076			*(wg++) = 0.0;
1077			}
1078
1079
1080			/*
1081			* Calculate the original offset into the data array, to enable
1082			* delta calculations in the pixel loop
1083			*
1084			* i runs over dims; j holds the working integer index in the
1085			* current dim.
1086			*
1087			* This code matches the incrementation code at the bottom of the accumulation loop
1088			*/
1089
1090			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
1091			i_off = 0;
1092			for(i=0;i
1093			j = ivec[i] + ibvec[i];
1094			if(j<0 \|\| j >= in->dims[i]) {
1095			switch(bounds[i]) {
1096			case 0: /* no breakage allowed */
1097			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
1098			break;
1099			case 1: /* truncation */
1100			t_vio++;
1101			/* fall through */
1102			case 2: /* extension -- crop */
1103			if(j<0)
1104			j=0;
1105			else j = in->dims[i] - 1;
1106			break;
1107			case 3: /* periodic -- mod it */
1108			j %= in->dims[i];
1109			if(j<0)
1110			j += in->dims[i];
1111			break;
1112			case 4: /* mirror -- reflect off the edges */
1113			j += in->dims[i];
1114			j %= (in->dims[i]*2);
1115			if(j<0)
1116			j += in->dims[i]*2;
1117			j -= in->dims[i];
1118			if(j<0) {
1119			j *= -1;
1120			j -= 1;
1121			}
1122			break;
1123			default:
1124			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
1125			break;
1126			}
1127			}
1128			i_off += in->dimincs[i] * j;
1129			}
1130
1131			/* Initialize index stashes for later reference as we scan the footprint */
1132			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
1133			/* end of a dimensional scan. So we stash the index location at the */
1134			/* start of each dimensional scan here. When we finish incrementing */
1135			/* through a particular dim, we pull its value back out of the stash. */
1136			for(i=0;i
1137			index_stash[i] = i_off;
1138			}
1139
1140			/* The input accumulation loop is the hotspot for the entire operation. */
1141			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
1142			/* in the input array, use the linearized transform to bring each pixel center */
1143			/* forward to the output plane, and calculate a weighting based on the chosen */
1144			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
1145			/* table that is initialized the first time through the code. 'H' is the */
1146			/* same process, but explicitly calculated for each interation (~2x slower). */
1147			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
1148			/* into the input array fresh from the current input array vector each time, */
1149			/* we walk through the array using dimincs and the old offset. This saves */
1150			/* about half of the time spent on index calculation. */
1151
1152			do { /* Input accumulation loop */
1153			PDL_Double *cp;
1154			PDL_Double alpha;
1155			/* Calculate the weight of the current input point. Don't bother if we're
1156			* violating any truncation boundaries (in that case our value is zero, but
1157			* for the interpolation we also set the weight to zero).
1158			*/
1159			if( !t_vio ) {
1160
1161			PDL_Double *ap = tvec;
1162			PDL_Double *bp = dvec;
1163			PDL_Indx *ip = ivec;
1164			for(i=0; i
1165			(ap++) = (ip++) - *(bp++);
1166
1167			switch(method) {
1168			PDL_Double dd;
1169			case 'h':
1170			/* This is the Hanning window rolloff. It is a product of a simple */
1171			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
1172			/* is about 2x faster than using cos(theta) directly in each */
1173			/* weighting calculation, so we do. Using 2500 entries and linear */
1174			/* interpolation is accurate to about 10^-7, and should preserve */
1175			/* the contents of cache pretty well. */
1176			alpha = 1;
1177			cp = tmp;
1178			for (i=0; i
1179			PDL_Indx lodex;
1180			PDL_Indx hidex;
1181			PDL_Double beta;
1182			dd = 0;
1183			ap = tvec;
1184			/* Get the matrix-multiply element for this dimension */
1185			for (j=0;j
1186			dd += (ap++) *(cp++);
1187
1188			/* Do linear interpolation from the table */
1189			/* The table captures a hanning window centered 0.5 pixel from center. */
1190			/* We scale the filter by the blur parameter -- but if blur is less */
1191			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
1192			/* value on the pixel edge at 0.5. For blur greater than unity, we */
1193			/* scale simply. */
1194			beta = fabs(dd) - hanning_offset;
1195			if (beta > 0) {
1196			if (beta >= blur) {
1197			alpha = 0;
1198			i = ndims;
1199			} else {
1200			beta *= zeta;
1201			lodex = beta;
1202			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
1203			lodex = HANNING_LOOKUP_SIZE;
1204			hidex = lodex+1;
1205			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
1206			} /* end of interpolation branch */
1207			} /* end of beta > 0 branch */
1208			} /* end of dimension loop */
1209			break;
1210
1211			case 'H':
1212			/* This is the Hanning window rolloff with explicit calculation, preserved */
1213			/* in case someone actually wants the slower longer method. */
1214			alpha = 1;
1215			cp = tmp;
1216			for (i=0; i
1217			dd = 0;
1218			ap = tvec;
1219			for (j=0;j
1220			dd += (ap++) *(cp++);
1221			dd = (fabs(dd) - hanning_offset) / blur;
1222			if ( dd > 1 ) {
1223			alpha = 0;
1224			i = ndims;
1225			} else
1226			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
1227			}
1228			break;
1229
1230			case 'g':
1231			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
1232			{
1233			PDL_Double sum = 0;
1234			cp = tmp;
1235			for(i=0; i
1236			dd = 0;
1237			ap = tvec;
1238			for(j=0;j
1239			dd += (ap++) *(cp++);
1240			dd /= blur;
1241			sum += dd * dd;
1242			if(sum > GAUSSIAN_MAXVAL) {
1243			i = ndims; /* exit early if we're too far out */
1244			alpha = 0;
1245			}
1246			}
1247			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
1248			alpha = 0;
1249			} else {
1250			PDL_Indx lodex,hidex;
1251			PDL_Double beta = fabs(zeta * sum);
1252
1253			lodex = beta;
1254			beta -= lodex; hidex = lodex+1;
1255			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
1256
1257			}
1258			}
1259			break;
1260
1261			case 'G':
1262			/* This is the Gaussian rolloff with explicit calculation, preserved */
1263			/* in case someone actually wants the slower longer method. */
1264			{
1265			PDL_Double sum = 0;
1266			cp = tmp;
1267			for(i=0; i
1268			dd = 0;
1269			ap = tvec;
1270			for(j=0;j
1271			dd += (ap++) *(cp++);
1272			dd /= blur;
1273			sum += dd * dd;
1274			if(sum > 4) /* 2 pixels -- four half-widths */
1275			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
1276			}
1277
1278			if(sum > GAUSSIAN_MAXVAL)
1279			alpha = 0;
1280			else
1281			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
1282			}
1283			break;
1284			default:
1285			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
1286			}
1287
1288			{ /* convenience block -- accumulate the current point into the weighted sum. */
1289			/* This is more than simple assignment because we have our own explicit poor */
1290			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
1291			PDL_Byte dat = ((PDL_Byte )(in->data)) + i_off;
1292			PDL_Indx max = out->dims[ndims];
1293			for( i=0; i < max; i++ ) {
1294			if( (badval==0) \|\| (*dat != badval) ) {
1295			acc[i] += dat alpha;
1296			dat += in->dimincs[ndims];
1297			wgt[i] += alpha;
1298			}
1299			wgt2[i] += alpha; }
1300			}
1301			} /* end of t_vio check (i.e. of input accumulation) */
1302
1303
1304			/* Advance input accumulation loop. */
1305			/* We both increment the total vector and also advance the index. */
1306			carry = 1;
1307			for (i=0; i
1308			/* Advance the current element of the offset vector */
1309			ivec[i]++;
1310			j = ivec[i] + ibvec[i];
1311
1312			/* Advance the offset into the data array */
1313			if( j > 0 && j <= in->dims[i]-1 ) {
1314			/* Normal case -- just advance the input vector */
1315			i_off += in->dimincs[i];
1316			} else {
1317			/* Busted a boundary - either before or after. */
1318			switch(bounds[i]){
1319			case 0: /* no breakage allowed -- treat as truncation for interpolation */
1320			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
1321			if( j == 0 )
1322			t_vio--;
1323			else if( j == in->dims[i] )
1324			t_vio++;
1325			break;
1326			case 2: /* extension -- do nothing (so the same input point is re-used) */
1327			break;
1328			case 3: /* periodic -- advance and mod into the allowed range */
1329			if((j % in->dims[i]) == 0) {
1330			i_off -= in->dimincs[i] * (in->dims[i]-1);
1331			} else {
1332			i_off += in->dimincs[i];
1333			}
1334			break;
1335			case 4: /* mirror -- advance or retreat depending on phase */
1336			j += in->dims[i];
1337			j %= (in->dims[i]*2);
1338			j -= in->dims[i];
1339			if( j!=0 && j!= -in->dims[i] ) {
1340			if(j<0)
1341			i_off -= in->dimincs[i];
1342			else
1343			i_off += in->dimincs[i];
1344			}
1345			break;
1346			}
1347			}
1348
1349			/* Now check for carry */
1350			if (ivec[i] <= psize) {
1351			/* Normal case -- copy the current offset to the faster-running dim stashes */
1352			PDL_Indx k;
1353			for(k=0;k
1354			index_stash[k] = i_off;
1355			}
1356			carry = 0;
1357
1358			} else { /* End of this scan -- recover the last position, and mark carry */
1359			i_off = index_stash[i];
1360			if(bounds[i]==1) {
1361			j = ivec[i] + ibvec[i];
1362			if( j < 0 \|\| j >= in->dims[i] )
1363			t_vio--;
1364			ivec[i] = -psize;
1365			j = ivec[i] + ibvec[i];
1366			if( j < 0 \|\| j >= in->dims[i] )
1367			t_vio++;
1368			carry = 1;
1369			} else {
1370			ivec[i] = -psize;
1371			}
1372			}
1373			} /* End of counter-advance loop */
1374			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
1375
1376			{
1377			PDL_Double *ac = acc;
1378			PDL_Double *wg = wgt;
1379			PDL_Double *wg2 = wgt2;
1380			PDL_Byte *dat = out->data;
1381
1382			/* Calculate output vector offset */
1383			for(i=0;i
1384			dat += out->dimincs[i] * ovec[i];
1385
1386			if (!flux) {
1387			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
1388			for (i=0; i < out->dims[ndims]; i++) {
1389			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
1390			ac++; wg++; wg2++;
1391			dat += out->dimincs[ndims];
1392			}
1393			} else {
1394			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
1395			PDL_Double det = tmp[ndims*ndims];
1396			for(i=0; i < out->dims[ndims]; i++) {
1397			if(wg && (wg2 / *wg) < 1.5 ) {
1398			dat = (ac++) / (wg++) det;
1399			wg2++;
1400			} else {
1401			*dat = badval;
1402			ac++; wg++; wg2++;
1403			}
1404			dat += out->dimincs[ndims];
1405			} /* end of for loop */
1406			} /* end of flux flag set conditional */
1407			} /* end of convenience block */
1408
1409			/* End of code for normal pixels */
1410			} else {
1411			/* The pixel was ludicrously huge -- just set this pixel to nan */
1412			PDL_Byte *dat = out->data;
1413			for(i=0;i
1414			dat += out->dimincs[i] * ovec[i];
1415			for(i=0;idims[ndims];i++) {
1416			dat = badval; / Should handle bad values too -- not yet */
1417			dat += out->dimincs[ndims];
1418			}
1419			}
1420
1421			/* Increment the pixel counter */
1422			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
1423			} while (1);
1424			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
1425			#line 1426 "lib/PDL/Transform-pp-map.c"
1426	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
1427	0		} break;
1428	0		case PDL_S: {
1429	0	0	PDL_DECLARE_PARAMS_map_1(PDL_Short,S)
		0
		0
1430	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
1431			#line 611 "lib/PDL/Transform.pd"
1432
1433			/*
1434			* Pixel interpolation & averaging code
1435			*
1436			* Calls a common coordinate-transformation block (see following hdr)
1437			* that isn't dependent on the type of the input variable.
1438			*
1439			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
1440			* is handled internally. To simplify the broadcasting business, any
1441			* broadcast dimensions should all be collapsed to a single one by the
1442			* perl front-end.
1443			*
1444			*/
1445
1446			pdl *in = __params->in;
1447			pdl *out = __params->out;
1448			pdl *map = __params->map;
1449
1450			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
1451			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
1452			PDL_Indx ovec[ndims]; /* output pixel loop vector */
1453			PDL_Indx ivec[ndims]; /* input pixel loop vector */
1454			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
1455			PDL_Double dvec[ndims]; /* Residual vector for linearization */
1456			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
1457			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
1458			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
1459			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
1460			char bounds[ndims]; /* Boundary condition packed string */
1461			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
1462			char method; /* Method identifier (gets one of 'h','g') */
1463			PDL_Long big; /* Max size of input footprint for each pix */
1464			PDL_Double blur; /* Scaling of filter */
1465			PDL_Double sv_min; /* minimum singular value */
1466			char flux; /* Flag to indicate flux conservation */
1467			PDL_Indx i;
1468			PDL_Short badval = SvNV(__params->bv);
1469			#define HANNING_LOOKUP_SIZE 2500
1470			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
1471			static int needs_hanning_calc = 1;
1472			PDL_Double zeta = 0;
1473			PDL_Double hanning_offset;
1474
1475			#define GAUSSIAN_LOOKUP_SIZE 4000
1476			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
1477			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
1478			static int needs_gaussian_calc = 1;
1479
1480			PDL_RETERROR(PDL_err, PDL->make_physical(in));
1481			PDL_RETERROR(PDL_err, PDL->make_physical(out));
1482			PDL_RETERROR(PDL_err, PDL->make_physical(map));
1483
1484			/***
1485			* Fill in the boundary condition array
1486			*/
1487			{
1488			char *bstr;
1489			STRLEN blen;
1490			bstr = SvPV(__params->boundary,blen);
1491
1492			if(blen == 0) {
1493			/* If no boundary is specified then every dim gets truncated */
1494			PDL_Indx i;
1495			for (i=0;i
1496			bounds[i] = 1;
1497			} else {
1498			PDL_Indx i;
1499			for(i=0;i
1500			switch(bstr[i < blen ? i : blen-1 ]) {
1501			case '0': case 'f': case 'F': /* forbid */
1502			bounds[i] = 0;
1503			break;
1504			case '1': case 't': case 'T': /* truncate */
1505			bounds[i] = 1;
1506			break;
1507			case '2': case 'e': case 'E': /* extend */
1508			bounds[i] = 2;
1509			break;
1510			case '3': case 'p': case 'P': /* periodic */
1511			bounds[i] = 3;
1512			break;
1513			case '4': case 'm': case 'M': /* mirror */
1514			bounds[i] = 4;
1515			break;
1516			default:
1517			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
1518			break;
1519			}
1520			}
1521			}
1522			}
1523
1524			/***
1525			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
1526			*/
1527			big = labs(__params->big);
1528			if(big <= 0)
1529			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
1530
1531			blur = fabs(__params->blur);
1532			if(blur < 0)
1533			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
1534
1535			sv_min = fabs(__params->sv_min);
1536			if(sv_min < 0)
1537			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
1538
1539			flux = __params->flux != 0;
1540
1541			{
1542			char *mstr;
1543			STRLEN mlen;
1544			mstr = SvPV(__params->method,mlen);
1545
1546			if(mlen==0)
1547			method = 'h';
1548			switch(*mstr) {
1549			case 'h': if( needs_hanning_calc ) {
1550			int i;
1551			for(i=0;i
1552			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
1553			}
1554			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
1555			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
1556			needs_hanning_calc = 0;
1557			}
1558			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
1559			case 'H': method = *mstr;
1560			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
1561			break;
1562
1563			case 'g': case 'j': method = 'g';
1564			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
1565
1566			if( needs_gaussian_calc ) {
1567			int i;
1568			for(i=0;i
1569			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
1570			}
1571			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
1572			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
1573			needs_gaussian_calc = 0;
1574			}
1575			break;
1576
1577			case 'G': case 'J': method = 'G'; break;
1578			default:
1579			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
1580			break;
1581			}
1582			}
1583
1584
1585
1586			/* End of initialization */
1587			/*************************************************************/
1588			/* Start of Real Work */
1589
1590			/* Initialize coordinate vector and map offset
1591			*/
1592			for(i=0;i
1593			ovec[i] = 0;
1594
1595			PDL_Double map_ptr = (PDL_Double )(map->data);
1596			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
1597			for (i=0; i < npdls; i++) offs[i] = 0;
1598			for (i=0; i < ndims; i++) {
1599			incs[i*npdls + 0] = map->dimincs[i+1];
1600			}
1601
1602			/* Main pixel loop (iterates over pixels in the output plane) */
1603			do {
1604			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
1605			/* Prefrobnicate the transformation matrix */
1606			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
1607
1608			#ifdef DEBUG_MAP
1609			{
1610			PDL_Indx k; PDL_Indx foo = 0;
1611			printf("ovec: [");
1612			for(k=0;k
1613			foo += ovec[k] * map->dimincs[k+1];
1614			printf(" %2td ",(PDL_Indx)(ovec[k]));
1615			}
1616			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
1617			for(k=0;k
1618			printf("%8.2f",(double)(((PDL_Short )(map->data))[foo + kmap->dimincs[0]]));
1619			}
1620			printf("]\n");
1621			}
1622			#endif
1623
1624			/* Don't bother accumulating output if psize is too large */
1625			if(psize <= big) {
1626			/* Use the prefrobnicated matrix to generate a local linearization.
1627			* dvec gets the delta; ibvec gets the base.
1628			*/
1629			{
1630			PDL_Double *mp = map_ptr + offs[0];
1631			for (i=0;i
1632			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
1633			mp += map->dimincs[0];
1634			}
1635			}
1636
1637			/* Initialize input delta vector */
1638			for(i=0;i
1639			ivec[i] = -psize;
1640
1641			/* Initialize accumulators */
1642			{
1643			PDL_Double *ac = acc;
1644			for(i=0; i < in->dims[ndims]; i++)
1645			*(ac++) = 0.0;
1646
1647			}
1648			{
1649			PDL_Double *wg = wgt;
1650			for(i=0;i < in->dims[ndims]; i++)
1651			*(wg++) = 0.0;
1652			}
1653			{
1654			PDL_Double *wg = wgt2;
1655			for(i=0;i < in->dims[ndims]; i++)
1656			*(wg++) = 0.0;
1657			}
1658
1659
1660			/*
1661			* Calculate the original offset into the data array, to enable
1662			* delta calculations in the pixel loop
1663			*
1664			* i runs over dims; j holds the working integer index in the
1665			* current dim.
1666			*
1667			* This code matches the incrementation code at the bottom of the accumulation loop
1668			*/
1669
1670			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
1671			i_off = 0;
1672			for(i=0;i
1673			j = ivec[i] + ibvec[i];
1674			if(j<0 \|\| j >= in->dims[i]) {
1675			switch(bounds[i]) {
1676			case 0: /* no breakage allowed */
1677			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
1678			break;
1679			case 1: /* truncation */
1680			t_vio++;
1681			/* fall through */
1682			case 2: /* extension -- crop */
1683			if(j<0)
1684			j=0;
1685			else j = in->dims[i] - 1;
1686			break;
1687			case 3: /* periodic -- mod it */
1688			j %= in->dims[i];
1689			if(j<0)
1690			j += in->dims[i];
1691			break;
1692			case 4: /* mirror -- reflect off the edges */
1693			j += in->dims[i];
1694			j %= (in->dims[i]*2);
1695			if(j<0)
1696			j += in->dims[i]*2;
1697			j -= in->dims[i];
1698			if(j<0) {
1699			j *= -1;
1700			j -= 1;
1701			}
1702			break;
1703			default:
1704			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
1705			break;
1706			}
1707			}
1708			i_off += in->dimincs[i] * j;
1709			}
1710
1711			/* Initialize index stashes for later reference as we scan the footprint */
1712			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
1713			/* end of a dimensional scan. So we stash the index location at the */
1714			/* start of each dimensional scan here. When we finish incrementing */
1715			/* through a particular dim, we pull its value back out of the stash. */
1716			for(i=0;i
1717			index_stash[i] = i_off;
1718			}
1719
1720			/* The input accumulation loop is the hotspot for the entire operation. */
1721			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
1722			/* in the input array, use the linearized transform to bring each pixel center */
1723			/* forward to the output plane, and calculate a weighting based on the chosen */
1724			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
1725			/* table that is initialized the first time through the code. 'H' is the */
1726			/* same process, but explicitly calculated for each interation (~2x slower). */
1727			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
1728			/* into the input array fresh from the current input array vector each time, */
1729			/* we walk through the array using dimincs and the old offset. This saves */
1730			/* about half of the time spent on index calculation. */
1731
1732			do { /* Input accumulation loop */
1733			PDL_Double *cp;
1734			PDL_Double alpha;
1735			/* Calculate the weight of the current input point. Don't bother if we're
1736			* violating any truncation boundaries (in that case our value is zero, but
1737			* for the interpolation we also set the weight to zero).
1738			*/
1739			if( !t_vio ) {
1740
1741			PDL_Double *ap = tvec;
1742			PDL_Double *bp = dvec;
1743			PDL_Indx *ip = ivec;
1744			for(i=0; i
1745			(ap++) = (ip++) - *(bp++);
1746
1747			switch(method) {
1748			PDL_Double dd;
1749			case 'h':
1750			/* This is the Hanning window rolloff. It is a product of a simple */
1751			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
1752			/* is about 2x faster than using cos(theta) directly in each */
1753			/* weighting calculation, so we do. Using 2500 entries and linear */
1754			/* interpolation is accurate to about 10^-7, and should preserve */
1755			/* the contents of cache pretty well. */
1756			alpha = 1;
1757			cp = tmp;
1758			for (i=0; i
1759			PDL_Indx lodex;
1760			PDL_Indx hidex;
1761			PDL_Double beta;
1762			dd = 0;
1763			ap = tvec;
1764			/* Get the matrix-multiply element for this dimension */
1765			for (j=0;j
1766			dd += (ap++) *(cp++);
1767
1768			/* Do linear interpolation from the table */
1769			/* The table captures a hanning window centered 0.5 pixel from center. */
1770			/* We scale the filter by the blur parameter -- but if blur is less */
1771			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
1772			/* value on the pixel edge at 0.5. For blur greater than unity, we */
1773			/* scale simply. */
1774			beta = fabs(dd) - hanning_offset;
1775			if (beta > 0) {
1776			if (beta >= blur) {
1777			alpha = 0;
1778			i = ndims;
1779			} else {
1780			beta *= zeta;
1781			lodex = beta;
1782			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
1783			lodex = HANNING_LOOKUP_SIZE;
1784			hidex = lodex+1;
1785			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
1786			} /* end of interpolation branch */
1787			} /* end of beta > 0 branch */
1788			} /* end of dimension loop */
1789			break;
1790
1791			case 'H':
1792			/* This is the Hanning window rolloff with explicit calculation, preserved */
1793			/* in case someone actually wants the slower longer method. */
1794			alpha = 1;
1795			cp = tmp;
1796			for (i=0; i
1797			dd = 0;
1798			ap = tvec;
1799			for (j=0;j
1800			dd += (ap++) *(cp++);
1801			dd = (fabs(dd) - hanning_offset) / blur;
1802			if ( dd > 1 ) {
1803			alpha = 0;
1804			i = ndims;
1805			} else
1806			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
1807			}
1808			break;
1809
1810			case 'g':
1811			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
1812			{
1813			PDL_Double sum = 0;
1814			cp = tmp;
1815			for(i=0; i
1816			dd = 0;
1817			ap = tvec;
1818			for(j=0;j
1819			dd += (ap++) *(cp++);
1820			dd /= blur;
1821			sum += dd * dd;
1822			if(sum > GAUSSIAN_MAXVAL) {
1823			i = ndims; /* exit early if we're too far out */
1824			alpha = 0;
1825			}
1826			}
1827			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
1828			alpha = 0;
1829			} else {
1830			PDL_Indx lodex,hidex;
1831			PDL_Double beta = fabs(zeta * sum);
1832
1833			lodex = beta;
1834			beta -= lodex; hidex = lodex+1;
1835			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
1836
1837			}
1838			}
1839			break;
1840
1841			case 'G':
1842			/* This is the Gaussian rolloff with explicit calculation, preserved */
1843			/* in case someone actually wants the slower longer method. */
1844			{
1845			PDL_Double sum = 0;
1846			cp = tmp;
1847			for(i=0; i
1848			dd = 0;
1849			ap = tvec;
1850			for(j=0;j
1851			dd += (ap++) *(cp++);
1852			dd /= blur;
1853			sum += dd * dd;
1854			if(sum > 4) /* 2 pixels -- four half-widths */
1855			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
1856			}
1857
1858			if(sum > GAUSSIAN_MAXVAL)
1859			alpha = 0;
1860			else
1861			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
1862			}
1863			break;
1864			default:
1865			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
1866			}
1867
1868			{ /* convenience block -- accumulate the current point into the weighted sum. */
1869			/* This is more than simple assignment because we have our own explicit poor */
1870			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
1871			PDL_Short dat = ((PDL_Short )(in->data)) + i_off;
1872			PDL_Indx max = out->dims[ndims];
1873			for( i=0; i < max; i++ ) {
1874			if( (badval==0) \|\| (*dat != badval) ) {
1875			acc[i] += dat alpha;
1876			dat += in->dimincs[ndims];
1877			wgt[i] += alpha;
1878			}
1879			wgt2[i] += alpha; }
1880			}
1881			} /* end of t_vio check (i.e. of input accumulation) */
1882
1883
1884			/* Advance input accumulation loop. */
1885			/* We both increment the total vector and also advance the index. */
1886			carry = 1;
1887			for (i=0; i
1888			/* Advance the current element of the offset vector */
1889			ivec[i]++;
1890			j = ivec[i] + ibvec[i];
1891
1892			/* Advance the offset into the data array */
1893			if( j > 0 && j <= in->dims[i]-1 ) {
1894			/* Normal case -- just advance the input vector */
1895			i_off += in->dimincs[i];
1896			} else {
1897			/* Busted a boundary - either before or after. */
1898			switch(bounds[i]){
1899			case 0: /* no breakage allowed -- treat as truncation for interpolation */
1900			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
1901			if( j == 0 )
1902			t_vio--;
1903			else if( j == in->dims[i] )
1904			t_vio++;
1905			break;
1906			case 2: /* extension -- do nothing (so the same input point is re-used) */
1907			break;
1908			case 3: /* periodic -- advance and mod into the allowed range */
1909			if((j % in->dims[i]) == 0) {
1910			i_off -= in->dimincs[i] * (in->dims[i]-1);
1911			} else {
1912			i_off += in->dimincs[i];
1913			}
1914			break;
1915			case 4: /* mirror -- advance or retreat depending on phase */
1916			j += in->dims[i];
1917			j %= (in->dims[i]*2);
1918			j -= in->dims[i];
1919			if( j!=0 && j!= -in->dims[i] ) {
1920			if(j<0)
1921			i_off -= in->dimincs[i];
1922			else
1923			i_off += in->dimincs[i];
1924			}
1925			break;
1926			}
1927			}
1928
1929			/* Now check for carry */
1930			if (ivec[i] <= psize) {
1931			/* Normal case -- copy the current offset to the faster-running dim stashes */
1932			PDL_Indx k;
1933			for(k=0;k
1934			index_stash[k] = i_off;
1935			}
1936			carry = 0;
1937
1938			} else { /* End of this scan -- recover the last position, and mark carry */
1939			i_off = index_stash[i];
1940			if(bounds[i]==1) {
1941			j = ivec[i] + ibvec[i];
1942			if( j < 0 \|\| j >= in->dims[i] )
1943			t_vio--;
1944			ivec[i] = -psize;
1945			j = ivec[i] + ibvec[i];
1946			if( j < 0 \|\| j >= in->dims[i] )
1947			t_vio++;
1948			carry = 1;
1949			} else {
1950			ivec[i] = -psize;
1951			}
1952			}
1953			} /* End of counter-advance loop */
1954			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
1955
1956			{
1957			PDL_Double *ac = acc;
1958			PDL_Double *wg = wgt;
1959			PDL_Double *wg2 = wgt2;
1960			PDL_Short *dat = out->data;
1961
1962			/* Calculate output vector offset */
1963			for(i=0;i
1964			dat += out->dimincs[i] * ovec[i];
1965
1966			if (!flux) {
1967			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
1968			for (i=0; i < out->dims[ndims]; i++) {
1969			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
1970			ac++; wg++; wg2++;
1971			dat += out->dimincs[ndims];
1972			}
1973			} else {
1974			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
1975			PDL_Double det = tmp[ndims*ndims];
1976			for(i=0; i < out->dims[ndims]; i++) {
1977			if(wg && (wg2 / *wg) < 1.5 ) {
1978			dat = (ac++) / (wg++) det;
1979			wg2++;
1980			} else {
1981			*dat = badval;
1982			ac++; wg++; wg2++;
1983			}
1984			dat += out->dimincs[ndims];
1985			} /* end of for loop */
1986			} /* end of flux flag set conditional */
1987			} /* end of convenience block */
1988
1989			/* End of code for normal pixels */
1990			} else {
1991			/* The pixel was ludicrously huge -- just set this pixel to nan */
1992			PDL_Short *dat = out->data;
1993			for(i=0;i
1994			dat += out->dimincs[i] * ovec[i];
1995			for(i=0;idims[ndims];i++) {
1996			dat = badval; / Should handle bad values too -- not yet */
1997			dat += out->dimincs[ndims];
1998			}
1999			}
2000
2001			/* Increment the pixel counter */
2002			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
2003			} while (1);
2004			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
2005			#line 2006 "lib/PDL/Transform-pp-map.c"
2006	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
2007	0		} break;
2008	0		case PDL_US: {
2009	0	0	PDL_DECLARE_PARAMS_map_1(PDL_Ushort,U)
		0
		0
2010	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
2011			#line 611 "lib/PDL/Transform.pd"
2012
2013			/*
2014			* Pixel interpolation & averaging code
2015			*
2016			* Calls a common coordinate-transformation block (see following hdr)
2017			* that isn't dependent on the type of the input variable.
2018			*
2019			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
2020			* is handled internally. To simplify the broadcasting business, any
2021			* broadcast dimensions should all be collapsed to a single one by the
2022			* perl front-end.
2023			*
2024			*/
2025
2026			pdl *in = __params->in;
2027			pdl *out = __params->out;
2028			pdl *map = __params->map;
2029
2030			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
2031			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
2032			PDL_Indx ovec[ndims]; /* output pixel loop vector */
2033			PDL_Indx ivec[ndims]; /* input pixel loop vector */
2034			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
2035			PDL_Double dvec[ndims]; /* Residual vector for linearization */
2036			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
2037			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
2038			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
2039			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
2040			char bounds[ndims]; /* Boundary condition packed string */
2041			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
2042			char method; /* Method identifier (gets one of 'h','g') */
2043			PDL_Long big; /* Max size of input footprint for each pix */
2044			PDL_Double blur; /* Scaling of filter */
2045			PDL_Double sv_min; /* minimum singular value */
2046			char flux; /* Flag to indicate flux conservation */
2047			PDL_Indx i;
2048			PDL_Ushort badval = SvNV(__params->bv);
2049			#define HANNING_LOOKUP_SIZE 2500
2050			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
2051			static int needs_hanning_calc = 1;
2052			PDL_Double zeta = 0;
2053			PDL_Double hanning_offset;
2054
2055			#define GAUSSIAN_LOOKUP_SIZE 4000
2056			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
2057			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
2058			static int needs_gaussian_calc = 1;
2059
2060			PDL_RETERROR(PDL_err, PDL->make_physical(in));
2061			PDL_RETERROR(PDL_err, PDL->make_physical(out));
2062			PDL_RETERROR(PDL_err, PDL->make_physical(map));
2063
2064			/***
2065			* Fill in the boundary condition array
2066			*/
2067			{
2068			char *bstr;
2069			STRLEN blen;
2070			bstr = SvPV(__params->boundary,blen);
2071
2072			if(blen == 0) {
2073			/* If no boundary is specified then every dim gets truncated */
2074			PDL_Indx i;
2075			for (i=0;i
2076			bounds[i] = 1;
2077			} else {
2078			PDL_Indx i;
2079			for(i=0;i
2080			switch(bstr[i < blen ? i : blen-1 ]) {
2081			case '0': case 'f': case 'F': /* forbid */
2082			bounds[i] = 0;
2083			break;
2084			case '1': case 't': case 'T': /* truncate */
2085			bounds[i] = 1;
2086			break;
2087			case '2': case 'e': case 'E': /* extend */
2088			bounds[i] = 2;
2089			break;
2090			case '3': case 'p': case 'P': /* periodic */
2091			bounds[i] = 3;
2092			break;
2093			case '4': case 'm': case 'M': /* mirror */
2094			bounds[i] = 4;
2095			break;
2096			default:
2097			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
2098			break;
2099			}
2100			}
2101			}
2102			}
2103
2104			/***
2105			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
2106			*/
2107			big = labs(__params->big);
2108			if(big <= 0)
2109			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
2110
2111			blur = fabs(__params->blur);
2112			if(blur < 0)
2113			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
2114
2115			sv_min = fabs(__params->sv_min);
2116			if(sv_min < 0)
2117			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
2118
2119			flux = __params->flux != 0;
2120
2121			{
2122			char *mstr;
2123			STRLEN mlen;
2124			mstr = SvPV(__params->method,mlen);
2125
2126			if(mlen==0)
2127			method = 'h';
2128			switch(*mstr) {
2129			case 'h': if( needs_hanning_calc ) {
2130			int i;
2131			for(i=0;i
2132			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
2133			}
2134			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
2135			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
2136			needs_hanning_calc = 0;
2137			}
2138			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
2139			case 'H': method = *mstr;
2140			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
2141			break;
2142
2143			case 'g': case 'j': method = 'g';
2144			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
2145
2146			if( needs_gaussian_calc ) {
2147			int i;
2148			for(i=0;i
2149			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
2150			}
2151			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
2152			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
2153			needs_gaussian_calc = 0;
2154			}
2155			break;
2156
2157			case 'G': case 'J': method = 'G'; break;
2158			default:
2159			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
2160			break;
2161			}
2162			}
2163
2164
2165
2166			/* End of initialization */
2167			/*************************************************************/
2168			/* Start of Real Work */
2169
2170			/* Initialize coordinate vector and map offset
2171			*/
2172			for(i=0;i
2173			ovec[i] = 0;
2174
2175			PDL_Double map_ptr = (PDL_Double )(map->data);
2176			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
2177			for (i=0; i < npdls; i++) offs[i] = 0;
2178			for (i=0; i < ndims; i++) {
2179			incs[i*npdls + 0] = map->dimincs[i+1];
2180			}
2181
2182			/* Main pixel loop (iterates over pixels in the output plane) */
2183			do {
2184			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
2185			/* Prefrobnicate the transformation matrix */
2186			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
2187
2188			#ifdef DEBUG_MAP
2189			{
2190			PDL_Indx k; PDL_Indx foo = 0;
2191			printf("ovec: [");
2192			for(k=0;k
2193			foo += ovec[k] * map->dimincs[k+1];
2194			printf(" %2td ",(PDL_Indx)(ovec[k]));
2195			}
2196			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
2197			for(k=0;k
2198			printf("%8.2f",(double)(((PDL_Ushort )(map->data))[foo + kmap->dimincs[0]]));
2199			}
2200			printf("]\n");
2201			}
2202			#endif
2203
2204			/* Don't bother accumulating output if psize is too large */
2205			if(psize <= big) {
2206			/* Use the prefrobnicated matrix to generate a local linearization.
2207			* dvec gets the delta; ibvec gets the base.
2208			*/
2209			{
2210			PDL_Double *mp = map_ptr + offs[0];
2211			for (i=0;i
2212			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
2213			mp += map->dimincs[0];
2214			}
2215			}
2216
2217			/* Initialize input delta vector */
2218			for(i=0;i
2219			ivec[i] = -psize;
2220
2221			/* Initialize accumulators */
2222			{
2223			PDL_Double *ac = acc;
2224			for(i=0; i < in->dims[ndims]; i++)
2225			*(ac++) = 0.0;
2226
2227			}
2228			{
2229			PDL_Double *wg = wgt;
2230			for(i=0;i < in->dims[ndims]; i++)
2231			*(wg++) = 0.0;
2232			}
2233			{
2234			PDL_Double *wg = wgt2;
2235			for(i=0;i < in->dims[ndims]; i++)
2236			*(wg++) = 0.0;
2237			}
2238
2239
2240			/*
2241			* Calculate the original offset into the data array, to enable
2242			* delta calculations in the pixel loop
2243			*
2244			* i runs over dims; j holds the working integer index in the
2245			* current dim.
2246			*
2247			* This code matches the incrementation code at the bottom of the accumulation loop
2248			*/
2249
2250			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
2251			i_off = 0;
2252			for(i=0;i
2253			j = ivec[i] + ibvec[i];
2254			if(j<0 \|\| j >= in->dims[i]) {
2255			switch(bounds[i]) {
2256			case 0: /* no breakage allowed */
2257			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
2258			break;
2259			case 1: /* truncation */
2260			t_vio++;
2261			/* fall through */
2262			case 2: /* extension -- crop */
2263			if(j<0)
2264			j=0;
2265			else j = in->dims[i] - 1;
2266			break;
2267			case 3: /* periodic -- mod it */
2268			j %= in->dims[i];
2269			if(j<0)
2270			j += in->dims[i];
2271			break;
2272			case 4: /* mirror -- reflect off the edges */
2273			j += in->dims[i];
2274			j %= (in->dims[i]*2);
2275			if(j<0)
2276			j += in->dims[i]*2;
2277			j -= in->dims[i];
2278			if(j<0) {
2279			j *= -1;
2280			j -= 1;
2281			}
2282			break;
2283			default:
2284			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
2285			break;
2286			}
2287			}
2288			i_off += in->dimincs[i] * j;
2289			}
2290
2291			/* Initialize index stashes for later reference as we scan the footprint */
2292			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
2293			/* end of a dimensional scan. So we stash the index location at the */
2294			/* start of each dimensional scan here. When we finish incrementing */
2295			/* through a particular dim, we pull its value back out of the stash. */
2296			for(i=0;i
2297			index_stash[i] = i_off;
2298			}
2299
2300			/* The input accumulation loop is the hotspot for the entire operation. */
2301			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
2302			/* in the input array, use the linearized transform to bring each pixel center */
2303			/* forward to the output plane, and calculate a weighting based on the chosen */
2304			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
2305			/* table that is initialized the first time through the code. 'H' is the */
2306			/* same process, but explicitly calculated for each interation (~2x slower). */
2307			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
2308			/* into the input array fresh from the current input array vector each time, */
2309			/* we walk through the array using dimincs and the old offset. This saves */
2310			/* about half of the time spent on index calculation. */
2311
2312			do { /* Input accumulation loop */
2313			PDL_Double *cp;
2314			PDL_Double alpha;
2315			/* Calculate the weight of the current input point. Don't bother if we're
2316			* violating any truncation boundaries (in that case our value is zero, but
2317			* for the interpolation we also set the weight to zero).
2318			*/
2319			if( !t_vio ) {
2320
2321			PDL_Double *ap = tvec;
2322			PDL_Double *bp = dvec;
2323			PDL_Indx *ip = ivec;
2324			for(i=0; i
2325			(ap++) = (ip++) - *(bp++);
2326
2327			switch(method) {
2328			PDL_Double dd;
2329			case 'h':
2330			/* This is the Hanning window rolloff. It is a product of a simple */
2331			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
2332			/* is about 2x faster than using cos(theta) directly in each */
2333			/* weighting calculation, so we do. Using 2500 entries and linear */
2334			/* interpolation is accurate to about 10^-7, and should preserve */
2335			/* the contents of cache pretty well. */
2336			alpha = 1;
2337			cp = tmp;
2338			for (i=0; i
2339			PDL_Indx lodex;
2340			PDL_Indx hidex;
2341			PDL_Double beta;
2342			dd = 0;
2343			ap = tvec;
2344			/* Get the matrix-multiply element for this dimension */
2345			for (j=0;j
2346			dd += (ap++) *(cp++);
2347
2348			/* Do linear interpolation from the table */
2349			/* The table captures a hanning window centered 0.5 pixel from center. */
2350			/* We scale the filter by the blur parameter -- but if blur is less */
2351			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
2352			/* value on the pixel edge at 0.5. For blur greater than unity, we */
2353			/* scale simply. */
2354			beta = fabs(dd) - hanning_offset;
2355			if (beta > 0) {
2356			if (beta >= blur) {
2357			alpha = 0;
2358			i = ndims;
2359			} else {
2360			beta *= zeta;
2361			lodex = beta;
2362			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
2363			lodex = HANNING_LOOKUP_SIZE;
2364			hidex = lodex+1;
2365			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
2366			} /* end of interpolation branch */
2367			} /* end of beta > 0 branch */
2368			} /* end of dimension loop */
2369			break;
2370
2371			case 'H':
2372			/* This is the Hanning window rolloff with explicit calculation, preserved */
2373			/* in case someone actually wants the slower longer method. */
2374			alpha = 1;
2375			cp = tmp;
2376			for (i=0; i
2377			dd = 0;
2378			ap = tvec;
2379			for (j=0;j
2380			dd += (ap++) *(cp++);
2381			dd = (fabs(dd) - hanning_offset) / blur;
2382			if ( dd > 1 ) {
2383			alpha = 0;
2384			i = ndims;
2385			} else
2386			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
2387			}
2388			break;
2389
2390			case 'g':
2391			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
2392			{
2393			PDL_Double sum = 0;
2394			cp = tmp;
2395			for(i=0; i
2396			dd = 0;
2397			ap = tvec;
2398			for(j=0;j
2399			dd += (ap++) *(cp++);
2400			dd /= blur;
2401			sum += dd * dd;
2402			if(sum > GAUSSIAN_MAXVAL) {
2403			i = ndims; /* exit early if we're too far out */
2404			alpha = 0;
2405			}
2406			}
2407			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
2408			alpha = 0;
2409			} else {
2410			PDL_Indx lodex,hidex;
2411			PDL_Double beta = fabs(zeta * sum);
2412
2413			lodex = beta;
2414			beta -= lodex; hidex = lodex+1;
2415			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
2416
2417			}
2418			}
2419			break;
2420
2421			case 'G':
2422			/* This is the Gaussian rolloff with explicit calculation, preserved */
2423			/* in case someone actually wants the slower longer method. */
2424			{
2425			PDL_Double sum = 0;
2426			cp = tmp;
2427			for(i=0; i
2428			dd = 0;
2429			ap = tvec;
2430			for(j=0;j
2431			dd += (ap++) *(cp++);
2432			dd /= blur;
2433			sum += dd * dd;
2434			if(sum > 4) /* 2 pixels -- four half-widths */
2435			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
2436			}
2437
2438			if(sum > GAUSSIAN_MAXVAL)
2439			alpha = 0;
2440			else
2441			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
2442			}
2443			break;
2444			default:
2445			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
2446			}
2447
2448			{ /* convenience block -- accumulate the current point into the weighted sum. */
2449			/* This is more than simple assignment because we have our own explicit poor */
2450			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
2451			PDL_Ushort dat = ((PDL_Ushort )(in->data)) + i_off;
2452			PDL_Indx max = out->dims[ndims];
2453			for( i=0; i < max; i++ ) {
2454			if( (badval==0) \|\| (*dat != badval) ) {
2455			acc[i] += dat alpha;
2456			dat += in->dimincs[ndims];
2457			wgt[i] += alpha;
2458			}
2459			wgt2[i] += alpha; }
2460			}
2461			} /* end of t_vio check (i.e. of input accumulation) */
2462
2463
2464			/* Advance input accumulation loop. */
2465			/* We both increment the total vector and also advance the index. */
2466			carry = 1;
2467			for (i=0; i
2468			/* Advance the current element of the offset vector */
2469			ivec[i]++;
2470			j = ivec[i] + ibvec[i];
2471
2472			/* Advance the offset into the data array */
2473			if( j > 0 && j <= in->dims[i]-1 ) {
2474			/* Normal case -- just advance the input vector */
2475			i_off += in->dimincs[i];
2476			} else {
2477			/* Busted a boundary - either before or after. */
2478			switch(bounds[i]){
2479			case 0: /* no breakage allowed -- treat as truncation for interpolation */
2480			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
2481			if( j == 0 )
2482			t_vio--;
2483			else if( j == in->dims[i] )
2484			t_vio++;
2485			break;
2486			case 2: /* extension -- do nothing (so the same input point is re-used) */
2487			break;
2488			case 3: /* periodic -- advance and mod into the allowed range */
2489			if((j % in->dims[i]) == 0) {
2490			i_off -= in->dimincs[i] * (in->dims[i]-1);
2491			} else {
2492			i_off += in->dimincs[i];
2493			}
2494			break;
2495			case 4: /* mirror -- advance or retreat depending on phase */
2496			j += in->dims[i];
2497			j %= (in->dims[i]*2);
2498			j -= in->dims[i];
2499			if( j!=0 && j!= -in->dims[i] ) {
2500			if(j<0)
2501			i_off -= in->dimincs[i];
2502			else
2503			i_off += in->dimincs[i];
2504			}
2505			break;
2506			}
2507			}
2508
2509			/* Now check for carry */
2510			if (ivec[i] <= psize) {
2511			/* Normal case -- copy the current offset to the faster-running dim stashes */
2512			PDL_Indx k;
2513			for(k=0;k
2514			index_stash[k] = i_off;
2515			}
2516			carry = 0;
2517
2518			} else { /* End of this scan -- recover the last position, and mark carry */
2519			i_off = index_stash[i];
2520			if(bounds[i]==1) {
2521			j = ivec[i] + ibvec[i];
2522			if( j < 0 \|\| j >= in->dims[i] )
2523			t_vio--;
2524			ivec[i] = -psize;
2525			j = ivec[i] + ibvec[i];
2526			if( j < 0 \|\| j >= in->dims[i] )
2527			t_vio++;
2528			carry = 1;
2529			} else {
2530			ivec[i] = -psize;
2531			}
2532			}
2533			} /* End of counter-advance loop */
2534			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
2535
2536			{
2537			PDL_Double *ac = acc;
2538			PDL_Double *wg = wgt;
2539			PDL_Double *wg2 = wgt2;
2540			PDL_Ushort *dat = out->data;
2541
2542			/* Calculate output vector offset */
2543			for(i=0;i
2544			dat += out->dimincs[i] * ovec[i];
2545
2546			if (!flux) {
2547			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
2548			for (i=0; i < out->dims[ndims]; i++) {
2549			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
2550			ac++; wg++; wg2++;
2551			dat += out->dimincs[ndims];
2552			}
2553			} else {
2554			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
2555			PDL_Double det = tmp[ndims*ndims];
2556			for(i=0; i < out->dims[ndims]; i++) {
2557			if(wg && (wg2 / *wg) < 1.5 ) {
2558			dat = (ac++) / (wg++) det;
2559			wg2++;
2560			} else {
2561			*dat = badval;
2562			ac++; wg++; wg2++;
2563			}
2564			dat += out->dimincs[ndims];
2565			} /* end of for loop */
2566			} /* end of flux flag set conditional */
2567			} /* end of convenience block */
2568
2569			/* End of code for normal pixels */
2570			} else {
2571			/* The pixel was ludicrously huge -- just set this pixel to nan */
2572			PDL_Ushort *dat = out->data;
2573			for(i=0;i
2574			dat += out->dimincs[i] * ovec[i];
2575			for(i=0;idims[ndims];i++) {
2576			dat = badval; / Should handle bad values too -- not yet */
2577			dat += out->dimincs[ndims];
2578			}
2579			}
2580
2581			/* Increment the pixel counter */
2582			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
2583			} while (1);
2584			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
2585			#line 2586 "lib/PDL/Transform-pp-map.c"
2586	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
2587	0		} break;
2588	5		case PDL_L: {
2589	5	50	PDL_DECLARE_PARAMS_map_1(PDL_Long,L)
		50
		50
2590	15	50	PDL_BROADCASTLOOP_START_map_readdata {
		50
		50
		50
		50
		100
		100
2591			#line 611 "lib/PDL/Transform.pd"
2592
2593			/*
2594			* Pixel interpolation & averaging code
2595			*
2596			* Calls a common coordinate-transformation block (see following hdr)
2597			* that isn't dependent on the type of the input variable.
2598			*
2599			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
2600			* is handled internally. To simplify the broadcasting business, any
2601			* broadcast dimensions should all be collapsed to a single one by the
2602			* perl front-end.
2603			*
2604			*/
2605
2606			pdl *in = __params->in;
2607			pdl *out = __params->out;
2608			pdl *map = __params->map;
2609
2610			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
2611			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
2612			PDL_Indx ovec[ndims]; /* output pixel loop vector */
2613			PDL_Indx ivec[ndims]; /* input pixel loop vector */
2614			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
2615			PDL_Double dvec[ndims]; /* Residual vector for linearization */
2616			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
2617			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
2618			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
2619			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
2620			char bounds[ndims]; /* Boundary condition packed string */
2621			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
2622			char method; /* Method identifier (gets one of 'h','g') */
2623			PDL_Long big; /* Max size of input footprint for each pix */
2624			PDL_Double blur; /* Scaling of filter */
2625			PDL_Double sv_min; /* minimum singular value */
2626			char flux; /* Flag to indicate flux conservation */
2627			PDL_Indx i;
2628			PDL_Long badval = SvNV(__params->bv);
2629			#define HANNING_LOOKUP_SIZE 2500
2630			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
2631			static int needs_hanning_calc = 1;
2632			PDL_Double zeta = 0;
2633			PDL_Double hanning_offset;
2634
2635			#define GAUSSIAN_LOOKUP_SIZE 4000
2636			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
2637			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
2638			static int needs_gaussian_calc = 1;
2639
2640			PDL_RETERROR(PDL_err, PDL->make_physical(in));
2641			PDL_RETERROR(PDL_err, PDL->make_physical(out));
2642			PDL_RETERROR(PDL_err, PDL->make_physical(map));
2643
2644			/***
2645			* Fill in the boundary condition array
2646			*/
2647			{
2648			char *bstr;
2649			STRLEN blen;
2650			bstr = SvPV(__params->boundary,blen);
2651
2652			if(blen == 0) {
2653			/* If no boundary is specified then every dim gets truncated */
2654			PDL_Indx i;
2655			for (i=0;i
2656			bounds[i] = 1;
2657			} else {
2658			PDL_Indx i;
2659			for(i=0;i
2660			switch(bstr[i < blen ? i : blen-1 ]) {
2661			case '0': case 'f': case 'F': /* forbid */
2662			bounds[i] = 0;
2663			break;
2664			case '1': case 't': case 'T': /* truncate */
2665			bounds[i] = 1;
2666			break;
2667			case '2': case 'e': case 'E': /* extend */
2668			bounds[i] = 2;
2669			break;
2670			case '3': case 'p': case 'P': /* periodic */
2671			bounds[i] = 3;
2672			break;
2673			case '4': case 'm': case 'M': /* mirror */
2674			bounds[i] = 4;
2675			break;
2676			default:
2677			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
2678			break;
2679			}
2680			}
2681			}
2682			}
2683
2684			/***
2685			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
2686			*/
2687			big = labs(__params->big);
2688			if(big <= 0)
2689			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
2690
2691			blur = fabs(__params->blur);
2692			if(blur < 0)
2693			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
2694
2695			sv_min = fabs(__params->sv_min);
2696			if(sv_min < 0)
2697			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
2698
2699			flux = __params->flux != 0;
2700
2701			{
2702			char *mstr;
2703			STRLEN mlen;
2704			mstr = SvPV(__params->method,mlen);
2705
2706			if(mlen==0)
2707			method = 'h';
2708			switch(*mstr) {
2709			case 'h': if( needs_hanning_calc ) {
2710			int i;
2711			for(i=0;i
2712			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
2713			}
2714			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
2715			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
2716			needs_hanning_calc = 0;
2717			}
2718			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
2719			case 'H': method = *mstr;
2720			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
2721			break;
2722
2723			case 'g': case 'j': method = 'g';
2724			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
2725
2726			if( needs_gaussian_calc ) {
2727			int i;
2728			for(i=0;i
2729			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
2730			}
2731			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
2732			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
2733			needs_gaussian_calc = 0;
2734			}
2735			break;
2736
2737			case 'G': case 'J': method = 'G'; break;
2738			default:
2739			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
2740			break;
2741			}
2742			}
2743
2744
2745
2746			/* End of initialization */
2747			/*************************************************************/
2748			/* Start of Real Work */
2749
2750			/* Initialize coordinate vector and map offset
2751			*/
2752			for(i=0;i
2753			ovec[i] = 0;
2754
2755			PDL_Double map_ptr = (PDL_Double )(map->data);
2756			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
2757			for (i=0; i < npdls; i++) offs[i] = 0;
2758			for (i=0; i < ndims; i++) {
2759			incs[i*npdls + 0] = map->dimincs[i+1];
2760			}
2761
2762			/* Main pixel loop (iterates over pixels in the output plane) */
2763			do {
2764			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
2765			/* Prefrobnicate the transformation matrix */
2766			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
2767
2768			#ifdef DEBUG_MAP
2769			{
2770			PDL_Indx k; PDL_Indx foo = 0;
2771			printf("ovec: [");
2772			for(k=0;k
2773			foo += ovec[k] * map->dimincs[k+1];
2774			printf(" %2td ",(PDL_Indx)(ovec[k]));
2775			}
2776			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
2777			for(k=0;k
2778			printf("%8.2f",(double)(((PDL_Long )(map->data))[foo + kmap->dimincs[0]]));
2779			}
2780			printf("]\n");
2781			}
2782			#endif
2783
2784			/* Don't bother accumulating output if psize is too large */
2785			if(psize <= big) {
2786			/* Use the prefrobnicated matrix to generate a local linearization.
2787			* dvec gets the delta; ibvec gets the base.
2788			*/
2789			{
2790			PDL_Double *mp = map_ptr + offs[0];
2791			for (i=0;i
2792			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
2793			mp += map->dimincs[0];
2794			}
2795			}
2796
2797			/* Initialize input delta vector */
2798			for(i=0;i
2799			ivec[i] = -psize;
2800
2801			/* Initialize accumulators */
2802			{
2803			PDL_Double *ac = acc;
2804			for(i=0; i < in->dims[ndims]; i++)
2805			*(ac++) = 0.0;
2806
2807			}
2808			{
2809			PDL_Double *wg = wgt;
2810			for(i=0;i < in->dims[ndims]; i++)
2811			*(wg++) = 0.0;
2812			}
2813			{
2814			PDL_Double *wg = wgt2;
2815			for(i=0;i < in->dims[ndims]; i++)
2816			*(wg++) = 0.0;
2817			}
2818
2819
2820			/*
2821			* Calculate the original offset into the data array, to enable
2822			* delta calculations in the pixel loop
2823			*
2824			* i runs over dims; j holds the working integer index in the
2825			* current dim.
2826			*
2827			* This code matches the incrementation code at the bottom of the accumulation loop
2828			*/
2829
2830			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
2831			i_off = 0;
2832			for(i=0;i
2833			j = ivec[i] + ibvec[i];
2834			if(j<0 \|\| j >= in->dims[i]) {
2835			switch(bounds[i]) {
2836			case 0: /* no breakage allowed */
2837			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
2838			break;
2839			case 1: /* truncation */
2840			t_vio++;
2841			/* fall through */
2842			case 2: /* extension -- crop */
2843			if(j<0)
2844			j=0;
2845			else j = in->dims[i] - 1;
2846			break;
2847			case 3: /* periodic -- mod it */
2848			j %= in->dims[i];
2849			if(j<0)
2850			j += in->dims[i];
2851			break;
2852			case 4: /* mirror -- reflect off the edges */
2853			j += in->dims[i];
2854			j %= (in->dims[i]*2);
2855			if(j<0)
2856			j += in->dims[i]*2;
2857			j -= in->dims[i];
2858			if(j<0) {
2859			j *= -1;
2860			j -= 1;
2861			}
2862			break;
2863			default:
2864			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
2865			break;
2866			}
2867			}
2868			i_off += in->dimincs[i] * j;
2869			}
2870
2871			/* Initialize index stashes for later reference as we scan the footprint */
2872			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
2873			/* end of a dimensional scan. So we stash the index location at the */
2874			/* start of each dimensional scan here. When we finish incrementing */
2875			/* through a particular dim, we pull its value back out of the stash. */
2876			for(i=0;i
2877			index_stash[i] = i_off;
2878			}
2879
2880			/* The input accumulation loop is the hotspot for the entire operation. */
2881			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
2882			/* in the input array, use the linearized transform to bring each pixel center */
2883			/* forward to the output plane, and calculate a weighting based on the chosen */
2884			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
2885			/* table that is initialized the first time through the code. 'H' is the */
2886			/* same process, but explicitly calculated for each interation (~2x slower). */
2887			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
2888			/* into the input array fresh from the current input array vector each time, */
2889			/* we walk through the array using dimincs and the old offset. This saves */
2890			/* about half of the time spent on index calculation. */
2891
2892			do { /* Input accumulation loop */
2893			PDL_Double *cp;
2894			PDL_Double alpha;
2895			/* Calculate the weight of the current input point. Don't bother if we're
2896			* violating any truncation boundaries (in that case our value is zero, but
2897			* for the interpolation we also set the weight to zero).
2898			*/
2899			if( !t_vio ) {
2900
2901			PDL_Double *ap = tvec;
2902			PDL_Double *bp = dvec;
2903			PDL_Indx *ip = ivec;
2904			for(i=0; i
2905			(ap++) = (ip++) - *(bp++);
2906
2907			switch(method) {
2908			PDL_Double dd;
2909			case 'h':
2910			/* This is the Hanning window rolloff. It is a product of a simple */
2911			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
2912			/* is about 2x faster than using cos(theta) directly in each */
2913			/* weighting calculation, so we do. Using 2500 entries and linear */
2914			/* interpolation is accurate to about 10^-7, and should preserve */
2915			/* the contents of cache pretty well. */
2916			alpha = 1;
2917			cp = tmp;
2918			for (i=0; i
2919			PDL_Indx lodex;
2920			PDL_Indx hidex;
2921			PDL_Double beta;
2922			dd = 0;
2923			ap = tvec;
2924			/* Get the matrix-multiply element for this dimension */
2925			for (j=0;j
2926			dd += (ap++) *(cp++);
2927
2928			/* Do linear interpolation from the table */
2929			/* The table captures a hanning window centered 0.5 pixel from center. */
2930			/* We scale the filter by the blur parameter -- but if blur is less */
2931			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
2932			/* value on the pixel edge at 0.5. For blur greater than unity, we */
2933			/* scale simply. */
2934			beta = fabs(dd) - hanning_offset;
2935			if (beta > 0) {
2936			if (beta >= blur) {
2937			alpha = 0;
2938			i = ndims;
2939			} else {
2940			beta *= zeta;
2941			lodex = beta;
2942			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
2943			lodex = HANNING_LOOKUP_SIZE;
2944			hidex = lodex+1;
2945			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
2946			} /* end of interpolation branch */
2947			} /* end of beta > 0 branch */
2948			} /* end of dimension loop */
2949			break;
2950
2951			case 'H':
2952			/* This is the Hanning window rolloff with explicit calculation, preserved */
2953			/* in case someone actually wants the slower longer method. */
2954			alpha = 1;
2955			cp = tmp;
2956			for (i=0; i
2957			dd = 0;
2958			ap = tvec;
2959			for (j=0;j
2960			dd += (ap++) *(cp++);
2961			dd = (fabs(dd) - hanning_offset) / blur;
2962			if ( dd > 1 ) {
2963			alpha = 0;
2964			i = ndims;
2965			} else
2966			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
2967			}
2968			break;
2969
2970			case 'g':
2971			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
2972			{
2973			PDL_Double sum = 0;
2974			cp = tmp;
2975			for(i=0; i
2976			dd = 0;
2977			ap = tvec;
2978			for(j=0;j
2979			dd += (ap++) *(cp++);
2980			dd /= blur;
2981			sum += dd * dd;
2982			if(sum > GAUSSIAN_MAXVAL) {
2983			i = ndims; /* exit early if we're too far out */
2984			alpha = 0;
2985			}
2986			}
2987			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
2988			alpha = 0;
2989			} else {
2990			PDL_Indx lodex,hidex;
2991			PDL_Double beta = fabs(zeta * sum);
2992
2993			lodex = beta;
2994			beta -= lodex; hidex = lodex+1;
2995			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
2996
2997			}
2998			}
2999			break;
3000
3001			case 'G':
3002			/* This is the Gaussian rolloff with explicit calculation, preserved */
3003			/* in case someone actually wants the slower longer method. */
3004			{
3005			PDL_Double sum = 0;
3006			cp = tmp;
3007			for(i=0; i
3008			dd = 0;
3009			ap = tvec;
3010			for(j=0;j
3011			dd += (ap++) *(cp++);
3012			dd /= blur;
3013			sum += dd * dd;
3014			if(sum > 4) /* 2 pixels -- four half-widths */
3015			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
3016			}
3017
3018			if(sum > GAUSSIAN_MAXVAL)
3019			alpha = 0;
3020			else
3021			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
3022			}
3023			break;
3024			default:
3025			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
3026			}
3027
3028			{ /* convenience block -- accumulate the current point into the weighted sum. */
3029			/* This is more than simple assignment because we have our own explicit poor */
3030			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
3031			PDL_Long dat = ((PDL_Long )(in->data)) + i_off;
3032			PDL_Indx max = out->dims[ndims];
3033			for( i=0; i < max; i++ ) {
3034			if( (badval==0) \|\| (*dat != badval) ) {
3035			acc[i] += dat alpha;
3036			dat += in->dimincs[ndims];
3037			wgt[i] += alpha;
3038			}
3039			wgt2[i] += alpha; }
3040			}
3041			} /* end of t_vio check (i.e. of input accumulation) */
3042
3043
3044			/* Advance input accumulation loop. */
3045			/* We both increment the total vector and also advance the index. */
3046			carry = 1;
3047			for (i=0; i
3048			/* Advance the current element of the offset vector */
3049			ivec[i]++;
3050			j = ivec[i] + ibvec[i];
3051
3052			/* Advance the offset into the data array */
3053			if( j > 0 && j <= in->dims[i]-1 ) {
3054			/* Normal case -- just advance the input vector */
3055			i_off += in->dimincs[i];
3056			} else {
3057			/* Busted a boundary - either before or after. */
3058			switch(bounds[i]){
3059			case 0: /* no breakage allowed -- treat as truncation for interpolation */
3060			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
3061			if( j == 0 )
3062			t_vio--;
3063			else if( j == in->dims[i] )
3064			t_vio++;
3065			break;
3066			case 2: /* extension -- do nothing (so the same input point is re-used) */
3067			break;
3068			case 3: /* periodic -- advance and mod into the allowed range */
3069			if((j % in->dims[i]) == 0) {
3070			i_off -= in->dimincs[i] * (in->dims[i]-1);
3071			} else {
3072			i_off += in->dimincs[i];
3073			}
3074			break;
3075			case 4: /* mirror -- advance or retreat depending on phase */
3076			j += in->dims[i];
3077			j %= (in->dims[i]*2);
3078			j -= in->dims[i];
3079			if( j!=0 && j!= -in->dims[i] ) {
3080			if(j<0)
3081			i_off -= in->dimincs[i];
3082			else
3083			i_off += in->dimincs[i];
3084			}
3085			break;
3086			}
3087			}
3088
3089			/* Now check for carry */
3090			if (ivec[i] <= psize) {
3091			/* Normal case -- copy the current offset to the faster-running dim stashes */
3092			PDL_Indx k;
3093			for(k=0;k
3094			index_stash[k] = i_off;
3095			}
3096			carry = 0;
3097
3098			} else { /* End of this scan -- recover the last position, and mark carry */
3099			i_off = index_stash[i];
3100			if(bounds[i]==1) {
3101			j = ivec[i] + ibvec[i];
3102			if( j < 0 \|\| j >= in->dims[i] )
3103			t_vio--;
3104			ivec[i] = -psize;
3105			j = ivec[i] + ibvec[i];
3106			if( j < 0 \|\| j >= in->dims[i] )
3107			t_vio++;
3108			carry = 1;
3109			} else {
3110			ivec[i] = -psize;
3111			}
3112			}
3113			} /* End of counter-advance loop */
3114			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
3115
3116			{
3117			PDL_Double *ac = acc;
3118			PDL_Double *wg = wgt;
3119			PDL_Double *wg2 = wgt2;
3120			PDL_Long *dat = out->data;
3121
3122			/* Calculate output vector offset */
3123			for(i=0;i
3124			dat += out->dimincs[i] * ovec[i];
3125
3126			if (!flux) {
3127			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
3128			for (i=0; i < out->dims[ndims]; i++) {
3129			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
3130			ac++; wg++; wg2++;
3131			dat += out->dimincs[ndims];
3132			}
3133			} else {
3134			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
3135			PDL_Double det = tmp[ndims*ndims];
3136			for(i=0; i < out->dims[ndims]; i++) {
3137			if(wg && (wg2 / *wg) < 1.5 ) {
3138			dat = (ac++) / (wg++) det;
3139			wg2++;
3140			} else {
3141			*dat = badval;
3142			ac++; wg++; wg2++;
3143			}
3144			dat += out->dimincs[ndims];
3145			} /* end of for loop */
3146			} /* end of flux flag set conditional */
3147			} /* end of convenience block */
3148
3149			/* End of code for normal pixels */
3150			} else {
3151			/* The pixel was ludicrously huge -- just set this pixel to nan */
3152			PDL_Long *dat = out->data;
3153			for(i=0;i
3154			dat += out->dimincs[i] * ovec[i];
3155			for(i=0;idims[ndims];i++) {
3156			dat = badval; / Should handle bad values too -- not yet */
3157			dat += out->dimincs[ndims];
3158			}
3159			}
3160
3161			/* Increment the pixel counter */
3162			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
3163			} while (1);
3164			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
3165			#line 3166 "lib/PDL/Transform-pp-map.c"
3166	5	50	}PDL_BROADCASTLOOP_END_map_readdata
		50
3167	5		} break;
3168	0		case PDL_UL: {
3169	0	0	PDL_DECLARE_PARAMS_map_1(PDL_ULong,K)
		0
		0
3170	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
3171			#line 611 "lib/PDL/Transform.pd"
3172
3173			/*
3174			* Pixel interpolation & averaging code
3175			*
3176			* Calls a common coordinate-transformation block (see following hdr)
3177			* that isn't dependent on the type of the input variable.
3178			*
3179			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
3180			* is handled internally. To simplify the broadcasting business, any
3181			* broadcast dimensions should all be collapsed to a single one by the
3182			* perl front-end.
3183			*
3184			*/
3185
3186			pdl *in = __params->in;
3187			pdl *out = __params->out;
3188			pdl *map = __params->map;
3189
3190			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
3191			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
3192			PDL_Indx ovec[ndims]; /* output pixel loop vector */
3193			PDL_Indx ivec[ndims]; /* input pixel loop vector */
3194			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
3195			PDL_Double dvec[ndims]; /* Residual vector for linearization */
3196			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
3197			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
3198			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
3199			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
3200			char bounds[ndims]; /* Boundary condition packed string */
3201			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
3202			char method; /* Method identifier (gets one of 'h','g') */
3203			PDL_Long big; /* Max size of input footprint for each pix */
3204			PDL_Double blur; /* Scaling of filter */
3205			PDL_Double sv_min; /* minimum singular value */
3206			char flux; /* Flag to indicate flux conservation */
3207			PDL_Indx i;
3208			PDL_ULong badval = SvNV(__params->bv);
3209			#define HANNING_LOOKUP_SIZE 2500
3210			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
3211			static int needs_hanning_calc = 1;
3212			PDL_Double zeta = 0;
3213			PDL_Double hanning_offset;
3214
3215			#define GAUSSIAN_LOOKUP_SIZE 4000
3216			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
3217			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
3218			static int needs_gaussian_calc = 1;
3219
3220			PDL_RETERROR(PDL_err, PDL->make_physical(in));
3221			PDL_RETERROR(PDL_err, PDL->make_physical(out));
3222			PDL_RETERROR(PDL_err, PDL->make_physical(map));
3223
3224			/***
3225			* Fill in the boundary condition array
3226			*/
3227			{
3228			char *bstr;
3229			STRLEN blen;
3230			bstr = SvPV(__params->boundary,blen);
3231
3232			if(blen == 0) {
3233			/* If no boundary is specified then every dim gets truncated */
3234			PDL_Indx i;
3235			for (i=0;i
3236			bounds[i] = 1;
3237			} else {
3238			PDL_Indx i;
3239			for(i=0;i
3240			switch(bstr[i < blen ? i : blen-1 ]) {
3241			case '0': case 'f': case 'F': /* forbid */
3242			bounds[i] = 0;
3243			break;
3244			case '1': case 't': case 'T': /* truncate */
3245			bounds[i] = 1;
3246			break;
3247			case '2': case 'e': case 'E': /* extend */
3248			bounds[i] = 2;
3249			break;
3250			case '3': case 'p': case 'P': /* periodic */
3251			bounds[i] = 3;
3252			break;
3253			case '4': case 'm': case 'M': /* mirror */
3254			bounds[i] = 4;
3255			break;
3256			default:
3257			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
3258			break;
3259			}
3260			}
3261			}
3262			}
3263
3264			/***
3265			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
3266			*/
3267			big = labs(__params->big);
3268			if(big <= 0)
3269			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
3270
3271			blur = fabs(__params->blur);
3272			if(blur < 0)
3273			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
3274
3275			sv_min = fabs(__params->sv_min);
3276			if(sv_min < 0)
3277			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
3278
3279			flux = __params->flux != 0;
3280
3281			{
3282			char *mstr;
3283			STRLEN mlen;
3284			mstr = SvPV(__params->method,mlen);
3285
3286			if(mlen==0)
3287			method = 'h';
3288			switch(*mstr) {
3289			case 'h': if( needs_hanning_calc ) {
3290			int i;
3291			for(i=0;i
3292			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
3293			}
3294			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
3295			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
3296			needs_hanning_calc = 0;
3297			}
3298			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
3299			case 'H': method = *mstr;
3300			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
3301			break;
3302
3303			case 'g': case 'j': method = 'g';
3304			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
3305
3306			if( needs_gaussian_calc ) {
3307			int i;
3308			for(i=0;i
3309			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
3310			}
3311			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
3312			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
3313			needs_gaussian_calc = 0;
3314			}
3315			break;
3316
3317			case 'G': case 'J': method = 'G'; break;
3318			default:
3319			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
3320			break;
3321			}
3322			}
3323
3324
3325
3326			/* End of initialization */
3327			/*************************************************************/
3328			/* Start of Real Work */
3329
3330			/* Initialize coordinate vector and map offset
3331			*/
3332			for(i=0;i
3333			ovec[i] = 0;
3334
3335			PDL_Double map_ptr = (PDL_Double )(map->data);
3336			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
3337			for (i=0; i < npdls; i++) offs[i] = 0;
3338			for (i=0; i < ndims; i++) {
3339			incs[i*npdls + 0] = map->dimincs[i+1];
3340			}
3341
3342			/* Main pixel loop (iterates over pixels in the output plane) */
3343			do {
3344			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
3345			/* Prefrobnicate the transformation matrix */
3346			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
3347
3348			#ifdef DEBUG_MAP
3349			{
3350			PDL_Indx k; PDL_Indx foo = 0;
3351			printf("ovec: [");
3352			for(k=0;k
3353			foo += ovec[k] * map->dimincs[k+1];
3354			printf(" %2td ",(PDL_Indx)(ovec[k]));
3355			}
3356			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
3357			for(k=0;k
3358			printf("%8.2f",(double)(((PDL_ULong )(map->data))[foo + kmap->dimincs[0]]));
3359			}
3360			printf("]\n");
3361			}
3362			#endif
3363
3364			/* Don't bother accumulating output if psize is too large */
3365			if(psize <= big) {
3366			/* Use the prefrobnicated matrix to generate a local linearization.
3367			* dvec gets the delta; ibvec gets the base.
3368			*/
3369			{
3370			PDL_Double *mp = map_ptr + offs[0];
3371			for (i=0;i
3372			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
3373			mp += map->dimincs[0];
3374			}
3375			}
3376
3377			/* Initialize input delta vector */
3378			for(i=0;i
3379			ivec[i] = -psize;
3380
3381			/* Initialize accumulators */
3382			{
3383			PDL_Double *ac = acc;
3384			for(i=0; i < in->dims[ndims]; i++)
3385			*(ac++) = 0.0;
3386
3387			}
3388			{
3389			PDL_Double *wg = wgt;
3390			for(i=0;i < in->dims[ndims]; i++)
3391			*(wg++) = 0.0;
3392			}
3393			{
3394			PDL_Double *wg = wgt2;
3395			for(i=0;i < in->dims[ndims]; i++)
3396			*(wg++) = 0.0;
3397			}
3398
3399
3400			/*
3401			* Calculate the original offset into the data array, to enable
3402			* delta calculations in the pixel loop
3403			*
3404			* i runs over dims; j holds the working integer index in the
3405			* current dim.
3406			*
3407			* This code matches the incrementation code at the bottom of the accumulation loop
3408			*/
3409
3410			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
3411			i_off = 0;
3412			for(i=0;i
3413			j = ivec[i] + ibvec[i];
3414			if(j<0 \|\| j >= in->dims[i]) {
3415			switch(bounds[i]) {
3416			case 0: /* no breakage allowed */
3417			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
3418			break;
3419			case 1: /* truncation */
3420			t_vio++;
3421			/* fall through */
3422			case 2: /* extension -- crop */
3423			if(j<0)
3424			j=0;
3425			else j = in->dims[i] - 1;
3426			break;
3427			case 3: /* periodic -- mod it */
3428			j %= in->dims[i];
3429			if(j<0)
3430			j += in->dims[i];
3431			break;
3432			case 4: /* mirror -- reflect off the edges */
3433			j += in->dims[i];
3434			j %= (in->dims[i]*2);
3435			if(j<0)
3436			j += in->dims[i]*2;
3437			j -= in->dims[i];
3438			if(j<0) {
3439			j *= -1;
3440			j -= 1;
3441			}
3442			break;
3443			default:
3444			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
3445			break;
3446			}
3447			}
3448			i_off += in->dimincs[i] * j;
3449			}
3450
3451			/* Initialize index stashes for later reference as we scan the footprint */
3452			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
3453			/* end of a dimensional scan. So we stash the index location at the */
3454			/* start of each dimensional scan here. When we finish incrementing */
3455			/* through a particular dim, we pull its value back out of the stash. */
3456			for(i=0;i
3457			index_stash[i] = i_off;
3458			}
3459
3460			/* The input accumulation loop is the hotspot for the entire operation. */
3461			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
3462			/* in the input array, use the linearized transform to bring each pixel center */
3463			/* forward to the output plane, and calculate a weighting based on the chosen */
3464			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
3465			/* table that is initialized the first time through the code. 'H' is the */
3466			/* same process, but explicitly calculated for each interation (~2x slower). */
3467			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
3468			/* into the input array fresh from the current input array vector each time, */
3469			/* we walk through the array using dimincs and the old offset. This saves */
3470			/* about half of the time spent on index calculation. */
3471
3472			do { /* Input accumulation loop */
3473			PDL_Double *cp;
3474			PDL_Double alpha;
3475			/* Calculate the weight of the current input point. Don't bother if we're
3476			* violating any truncation boundaries (in that case our value is zero, but
3477			* for the interpolation we also set the weight to zero).
3478			*/
3479			if( !t_vio ) {
3480
3481			PDL_Double *ap = tvec;
3482			PDL_Double *bp = dvec;
3483			PDL_Indx *ip = ivec;
3484			for(i=0; i
3485			(ap++) = (ip++) - *(bp++);
3486
3487			switch(method) {
3488			PDL_Double dd;
3489			case 'h':
3490			/* This is the Hanning window rolloff. It is a product of a simple */
3491			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
3492			/* is about 2x faster than using cos(theta) directly in each */
3493			/* weighting calculation, so we do. Using 2500 entries and linear */
3494			/* interpolation is accurate to about 10^-7, and should preserve */
3495			/* the contents of cache pretty well. */
3496			alpha = 1;
3497			cp = tmp;
3498			for (i=0; i
3499			PDL_Indx lodex;
3500			PDL_Indx hidex;
3501			PDL_Double beta;
3502			dd = 0;
3503			ap = tvec;
3504			/* Get the matrix-multiply element for this dimension */
3505			for (j=0;j
3506			dd += (ap++) *(cp++);
3507
3508			/* Do linear interpolation from the table */
3509			/* The table captures a hanning window centered 0.5 pixel from center. */
3510			/* We scale the filter by the blur parameter -- but if blur is less */
3511			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
3512			/* value on the pixel edge at 0.5. For blur greater than unity, we */
3513			/* scale simply. */
3514			beta = fabs(dd) - hanning_offset;
3515			if (beta > 0) {
3516			if (beta >= blur) {
3517			alpha = 0;
3518			i = ndims;
3519			} else {
3520			beta *= zeta;
3521			lodex = beta;
3522			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
3523			lodex = HANNING_LOOKUP_SIZE;
3524			hidex = lodex+1;
3525			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
3526			} /* end of interpolation branch */
3527			} /* end of beta > 0 branch */
3528			} /* end of dimension loop */
3529			break;
3530
3531			case 'H':
3532			/* This is the Hanning window rolloff with explicit calculation, preserved */
3533			/* in case someone actually wants the slower longer method. */
3534			alpha = 1;
3535			cp = tmp;
3536			for (i=0; i
3537			dd = 0;
3538			ap = tvec;
3539			for (j=0;j
3540			dd += (ap++) *(cp++);
3541			dd = (fabs(dd) - hanning_offset) / blur;
3542			if ( dd > 1 ) {
3543			alpha = 0;
3544			i = ndims;
3545			} else
3546			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
3547			}
3548			break;
3549
3550			case 'g':
3551			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
3552			{
3553			PDL_Double sum = 0;
3554			cp = tmp;
3555			for(i=0; i
3556			dd = 0;
3557			ap = tvec;
3558			for(j=0;j
3559			dd += (ap++) *(cp++);
3560			dd /= blur;
3561			sum += dd * dd;
3562			if(sum > GAUSSIAN_MAXVAL) {
3563			i = ndims; /* exit early if we're too far out */
3564			alpha = 0;
3565			}
3566			}
3567			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
3568			alpha = 0;
3569			} else {
3570			PDL_Indx lodex,hidex;
3571			PDL_Double beta = fabs(zeta * sum);
3572
3573			lodex = beta;
3574			beta -= lodex; hidex = lodex+1;
3575			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
3576
3577			}
3578			}
3579			break;
3580
3581			case 'G':
3582			/* This is the Gaussian rolloff with explicit calculation, preserved */
3583			/* in case someone actually wants the slower longer method. */
3584			{
3585			PDL_Double sum = 0;
3586			cp = tmp;
3587			for(i=0; i
3588			dd = 0;
3589			ap = tvec;
3590			for(j=0;j
3591			dd += (ap++) *(cp++);
3592			dd /= blur;
3593			sum += dd * dd;
3594			if(sum > 4) /* 2 pixels -- four half-widths */
3595			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
3596			}
3597
3598			if(sum > GAUSSIAN_MAXVAL)
3599			alpha = 0;
3600			else
3601			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
3602			}
3603			break;
3604			default:
3605			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
3606			}
3607
3608			{ /* convenience block -- accumulate the current point into the weighted sum. */
3609			/* This is more than simple assignment because we have our own explicit poor */
3610			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
3611			PDL_ULong dat = ((PDL_ULong )(in->data)) + i_off;
3612			PDL_Indx max = out->dims[ndims];
3613			for( i=0; i < max; i++ ) {
3614			if( (badval==0) \|\| (*dat != badval) ) {
3615			acc[i] += dat alpha;
3616			dat += in->dimincs[ndims];
3617			wgt[i] += alpha;
3618			}
3619			wgt2[i] += alpha; }
3620			}
3621			} /* end of t_vio check (i.e. of input accumulation) */
3622
3623
3624			/* Advance input accumulation loop. */
3625			/* We both increment the total vector and also advance the index. */
3626			carry = 1;
3627			for (i=0; i
3628			/* Advance the current element of the offset vector */
3629			ivec[i]++;
3630			j = ivec[i] + ibvec[i];
3631
3632			/* Advance the offset into the data array */
3633			if( j > 0 && j <= in->dims[i]-1 ) {
3634			/* Normal case -- just advance the input vector */
3635			i_off += in->dimincs[i];
3636			} else {
3637			/* Busted a boundary - either before or after. */
3638			switch(bounds[i]){
3639			case 0: /* no breakage allowed -- treat as truncation for interpolation */
3640			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
3641			if( j == 0 )
3642			t_vio--;
3643			else if( j == in->dims[i] )
3644			t_vio++;
3645			break;
3646			case 2: /* extension -- do nothing (so the same input point is re-used) */
3647			break;
3648			case 3: /* periodic -- advance and mod into the allowed range */
3649			if((j % in->dims[i]) == 0) {
3650			i_off -= in->dimincs[i] * (in->dims[i]-1);
3651			} else {
3652			i_off += in->dimincs[i];
3653			}
3654			break;
3655			case 4: /* mirror -- advance or retreat depending on phase */
3656			j += in->dims[i];
3657			j %= (in->dims[i]*2);
3658			j -= in->dims[i];
3659			if( j!=0 && j!= -in->dims[i] ) {
3660			if(j<0)
3661			i_off -= in->dimincs[i];
3662			else
3663			i_off += in->dimincs[i];
3664			}
3665			break;
3666			}
3667			}
3668
3669			/* Now check for carry */
3670			if (ivec[i] <= psize) {
3671			/* Normal case -- copy the current offset to the faster-running dim stashes */
3672			PDL_Indx k;
3673			for(k=0;k
3674			index_stash[k] = i_off;
3675			}
3676			carry = 0;
3677
3678			} else { /* End of this scan -- recover the last position, and mark carry */
3679			i_off = index_stash[i];
3680			if(bounds[i]==1) {
3681			j = ivec[i] + ibvec[i];
3682			if( j < 0 \|\| j >= in->dims[i] )
3683			t_vio--;
3684			ivec[i] = -psize;
3685			j = ivec[i] + ibvec[i];
3686			if( j < 0 \|\| j >= in->dims[i] )
3687			t_vio++;
3688			carry = 1;
3689			} else {
3690			ivec[i] = -psize;
3691			}
3692			}
3693			} /* End of counter-advance loop */
3694			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
3695
3696			{
3697			PDL_Double *ac = acc;
3698			PDL_Double *wg = wgt;
3699			PDL_Double *wg2 = wgt2;
3700			PDL_ULong *dat = out->data;
3701
3702			/* Calculate output vector offset */
3703			for(i=0;i
3704			dat += out->dimincs[i] * ovec[i];
3705
3706			if (!flux) {
3707			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
3708			for (i=0; i < out->dims[ndims]; i++) {
3709			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
3710			ac++; wg++; wg2++;
3711			dat += out->dimincs[ndims];
3712			}
3713			} else {
3714			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
3715			PDL_Double det = tmp[ndims*ndims];
3716			for(i=0; i < out->dims[ndims]; i++) {
3717			if(wg && (wg2 / *wg) < 1.5 ) {
3718			dat = (ac++) / (wg++) det;
3719			wg2++;
3720			} else {
3721			*dat = badval;
3722			ac++; wg++; wg2++;
3723			}
3724			dat += out->dimincs[ndims];
3725			} /* end of for loop */
3726			} /* end of flux flag set conditional */
3727			} /* end of convenience block */
3728
3729			/* End of code for normal pixels */
3730			} else {
3731			/* The pixel was ludicrously huge -- just set this pixel to nan */
3732			PDL_ULong *dat = out->data;
3733			for(i=0;i
3734			dat += out->dimincs[i] * ovec[i];
3735			for(i=0;idims[ndims];i++) {
3736			dat = badval; / Should handle bad values too -- not yet */
3737			dat += out->dimincs[ndims];
3738			}
3739			}
3740
3741			/* Increment the pixel counter */
3742			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
3743			} while (1);
3744			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
3745			#line 3746 "lib/PDL/Transform-pp-map.c"
3746	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
3747	0		} break;
3748	0		case PDL_IND: {
3749	0	0	PDL_DECLARE_PARAMS_map_1(PDL_Indx,N)
		0
		0
3750	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
3751			#line 611 "lib/PDL/Transform.pd"
3752
3753			/*
3754			* Pixel interpolation & averaging code
3755			*
3756			* Calls a common coordinate-transformation block (see following hdr)
3757			* that isn't dependent on the type of the input variable.
3758			*
3759			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
3760			* is handled internally. To simplify the broadcasting business, any
3761			* broadcast dimensions should all be collapsed to a single one by the
3762			* perl front-end.
3763			*
3764			*/
3765
3766			pdl *in = __params->in;
3767			pdl *out = __params->out;
3768			pdl *map = __params->map;
3769
3770			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
3771			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
3772			PDL_Indx ovec[ndims]; /* output pixel loop vector */
3773			PDL_Indx ivec[ndims]; /* input pixel loop vector */
3774			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
3775			PDL_Double dvec[ndims]; /* Residual vector for linearization */
3776			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
3777			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
3778			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
3779			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
3780			char bounds[ndims]; /* Boundary condition packed string */
3781			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
3782			char method; /* Method identifier (gets one of 'h','g') */
3783			PDL_Long big; /* Max size of input footprint for each pix */
3784			PDL_Double blur; /* Scaling of filter */
3785			PDL_Double sv_min; /* minimum singular value */
3786			char flux; /* Flag to indicate flux conservation */
3787			PDL_Indx i;
3788			PDL_Indx badval = SvNV(__params->bv);
3789			#define HANNING_LOOKUP_SIZE 2500
3790			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
3791			static int needs_hanning_calc = 1;
3792			PDL_Double zeta = 0;
3793			PDL_Double hanning_offset;
3794
3795			#define GAUSSIAN_LOOKUP_SIZE 4000
3796			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
3797			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
3798			static int needs_gaussian_calc = 1;
3799
3800			PDL_RETERROR(PDL_err, PDL->make_physical(in));
3801			PDL_RETERROR(PDL_err, PDL->make_physical(out));
3802			PDL_RETERROR(PDL_err, PDL->make_physical(map));
3803
3804			/***
3805			* Fill in the boundary condition array
3806			*/
3807			{
3808			char *bstr;
3809			STRLEN blen;
3810			bstr = SvPV(__params->boundary,blen);
3811
3812			if(blen == 0) {
3813			/* If no boundary is specified then every dim gets truncated */
3814			PDL_Indx i;
3815			for (i=0;i
3816			bounds[i] = 1;
3817			} else {
3818			PDL_Indx i;
3819			for(i=0;i
3820			switch(bstr[i < blen ? i : blen-1 ]) {
3821			case '0': case 'f': case 'F': /* forbid */
3822			bounds[i] = 0;
3823			break;
3824			case '1': case 't': case 'T': /* truncate */
3825			bounds[i] = 1;
3826			break;
3827			case '2': case 'e': case 'E': /* extend */
3828			bounds[i] = 2;
3829			break;
3830			case '3': case 'p': case 'P': /* periodic */
3831			bounds[i] = 3;
3832			break;
3833			case '4': case 'm': case 'M': /* mirror */
3834			bounds[i] = 4;
3835			break;
3836			default:
3837			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
3838			break;
3839			}
3840			}
3841			}
3842			}
3843
3844			/***
3845			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
3846			*/
3847			big = labs(__params->big);
3848			if(big <= 0)
3849			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
3850
3851			blur = fabs(__params->blur);
3852			if(blur < 0)
3853			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
3854
3855			sv_min = fabs(__params->sv_min);
3856			if(sv_min < 0)
3857			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
3858
3859			flux = __params->flux != 0;
3860
3861			{
3862			char *mstr;
3863			STRLEN mlen;
3864			mstr = SvPV(__params->method,mlen);
3865
3866			if(mlen==0)
3867			method = 'h';
3868			switch(*mstr) {
3869			case 'h': if( needs_hanning_calc ) {
3870			int i;
3871			for(i=0;i
3872			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
3873			}
3874			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
3875			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
3876			needs_hanning_calc = 0;
3877			}
3878			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
3879			case 'H': method = *mstr;
3880			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
3881			break;
3882
3883			case 'g': case 'j': method = 'g';
3884			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
3885
3886			if( needs_gaussian_calc ) {
3887			int i;
3888			for(i=0;i
3889			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
3890			}
3891			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
3892			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
3893			needs_gaussian_calc = 0;
3894			}
3895			break;
3896
3897			case 'G': case 'J': method = 'G'; break;
3898			default:
3899			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
3900			break;
3901			}
3902			}
3903
3904
3905
3906			/* End of initialization */
3907			/*************************************************************/
3908			/* Start of Real Work */
3909
3910			/* Initialize coordinate vector and map offset
3911			*/
3912			for(i=0;i
3913			ovec[i] = 0;
3914
3915			PDL_Double map_ptr = (PDL_Double )(map->data);
3916			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
3917			for (i=0; i < npdls; i++) offs[i] = 0;
3918			for (i=0; i < ndims; i++) {
3919			incs[i*npdls + 0] = map->dimincs[i+1];
3920			}
3921
3922			/* Main pixel loop (iterates over pixels in the output plane) */
3923			do {
3924			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
3925			/* Prefrobnicate the transformation matrix */
3926			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
3927
3928			#ifdef DEBUG_MAP
3929			{
3930			PDL_Indx k; PDL_Indx foo = 0;
3931			printf("ovec: [");
3932			for(k=0;k
3933			foo += ovec[k] * map->dimincs[k+1];
3934			printf(" %2td ",(PDL_Indx)(ovec[k]));
3935			}
3936			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
3937			for(k=0;k
3938			printf("%8.2f",(double)(((PDL_Indx )(map->data))[foo + kmap->dimincs[0]]));
3939			}
3940			printf("]\n");
3941			}
3942			#endif
3943
3944			/* Don't bother accumulating output if psize is too large */
3945			if(psize <= big) {
3946			/* Use the prefrobnicated matrix to generate a local linearization.
3947			* dvec gets the delta; ibvec gets the base.
3948			*/
3949			{
3950			PDL_Double *mp = map_ptr + offs[0];
3951			for (i=0;i
3952			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
3953			mp += map->dimincs[0];
3954			}
3955			}
3956
3957			/* Initialize input delta vector */
3958			for(i=0;i
3959			ivec[i] = -psize;
3960
3961			/* Initialize accumulators */
3962			{
3963			PDL_Double *ac = acc;
3964			for(i=0; i < in->dims[ndims]; i++)
3965			*(ac++) = 0.0;
3966
3967			}
3968			{
3969			PDL_Double *wg = wgt;
3970			for(i=0;i < in->dims[ndims]; i++)
3971			*(wg++) = 0.0;
3972			}
3973			{
3974			PDL_Double *wg = wgt2;
3975			for(i=0;i < in->dims[ndims]; i++)
3976			*(wg++) = 0.0;
3977			}
3978
3979
3980			/*
3981			* Calculate the original offset into the data array, to enable
3982			* delta calculations in the pixel loop
3983			*
3984			* i runs over dims; j holds the working integer index in the
3985			* current dim.
3986			*
3987			* This code matches the incrementation code at the bottom of the accumulation loop
3988			*/
3989
3990			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
3991			i_off = 0;
3992			for(i=0;i
3993			j = ivec[i] + ibvec[i];
3994			if(j<0 \|\| j >= in->dims[i]) {
3995			switch(bounds[i]) {
3996			case 0: /* no breakage allowed */
3997			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
3998			break;
3999			case 1: /* truncation */
4000			t_vio++;
4001			/* fall through */
4002			case 2: /* extension -- crop */
4003			if(j<0)
4004			j=0;
4005			else j = in->dims[i] - 1;
4006			break;
4007			case 3: /* periodic -- mod it */
4008			j %= in->dims[i];
4009			if(j<0)
4010			j += in->dims[i];
4011			break;
4012			case 4: /* mirror -- reflect off the edges */
4013			j += in->dims[i];
4014			j %= (in->dims[i]*2);
4015			if(j<0)
4016			j += in->dims[i]*2;
4017			j -= in->dims[i];
4018			if(j<0) {
4019			j *= -1;
4020			j -= 1;
4021			}
4022			break;
4023			default:
4024			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
4025			break;
4026			}
4027			}
4028			i_off += in->dimincs[i] * j;
4029			}
4030
4031			/* Initialize index stashes for later reference as we scan the footprint */
4032			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
4033			/* end of a dimensional scan. So we stash the index location at the */
4034			/* start of each dimensional scan here. When we finish incrementing */
4035			/* through a particular dim, we pull its value back out of the stash. */
4036			for(i=0;i
4037			index_stash[i] = i_off;
4038			}
4039
4040			/* The input accumulation loop is the hotspot for the entire operation. */
4041			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
4042			/* in the input array, use the linearized transform to bring each pixel center */
4043			/* forward to the output plane, and calculate a weighting based on the chosen */
4044			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
4045			/* table that is initialized the first time through the code. 'H' is the */
4046			/* same process, but explicitly calculated for each interation (~2x slower). */
4047			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
4048			/* into the input array fresh from the current input array vector each time, */
4049			/* we walk through the array using dimincs and the old offset. This saves */
4050			/* about half of the time spent on index calculation. */
4051
4052			do { /* Input accumulation loop */
4053			PDL_Double *cp;
4054			PDL_Double alpha;
4055			/* Calculate the weight of the current input point. Don't bother if we're
4056			* violating any truncation boundaries (in that case our value is zero, but
4057			* for the interpolation we also set the weight to zero).
4058			*/
4059			if( !t_vio ) {
4060
4061			PDL_Double *ap = tvec;
4062			PDL_Double *bp = dvec;
4063			PDL_Indx *ip = ivec;
4064			for(i=0; i
4065			(ap++) = (ip++) - *(bp++);
4066
4067			switch(method) {
4068			PDL_Double dd;
4069			case 'h':
4070			/* This is the Hanning window rolloff. It is a product of a simple */
4071			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
4072			/* is about 2x faster than using cos(theta) directly in each */
4073			/* weighting calculation, so we do. Using 2500 entries and linear */
4074			/* interpolation is accurate to about 10^-7, and should preserve */
4075			/* the contents of cache pretty well. */
4076			alpha = 1;
4077			cp = tmp;
4078			for (i=0; i
4079			PDL_Indx lodex;
4080			PDL_Indx hidex;
4081			PDL_Double beta;
4082			dd = 0;
4083			ap = tvec;
4084			/* Get the matrix-multiply element for this dimension */
4085			for (j=0;j
4086			dd += (ap++) *(cp++);
4087
4088			/* Do linear interpolation from the table */
4089			/* The table captures a hanning window centered 0.5 pixel from center. */
4090			/* We scale the filter by the blur parameter -- but if blur is less */
4091			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
4092			/* value on the pixel edge at 0.5. For blur greater than unity, we */
4093			/* scale simply. */
4094			beta = fabs(dd) - hanning_offset;
4095			if (beta > 0) {
4096			if (beta >= blur) {
4097			alpha = 0;
4098			i = ndims;
4099			} else {
4100			beta *= zeta;
4101			lodex = beta;
4102			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
4103			lodex = HANNING_LOOKUP_SIZE;
4104			hidex = lodex+1;
4105			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
4106			} /* end of interpolation branch */
4107			} /* end of beta > 0 branch */
4108			} /* end of dimension loop */
4109			break;
4110
4111			case 'H':
4112			/* This is the Hanning window rolloff with explicit calculation, preserved */
4113			/* in case someone actually wants the slower longer method. */
4114			alpha = 1;
4115			cp = tmp;
4116			for (i=0; i
4117			dd = 0;
4118			ap = tvec;
4119			for (j=0;j
4120			dd += (ap++) *(cp++);
4121			dd = (fabs(dd) - hanning_offset) / blur;
4122			if ( dd > 1 ) {
4123			alpha = 0;
4124			i = ndims;
4125			} else
4126			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
4127			}
4128			break;
4129
4130			case 'g':
4131			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
4132			{
4133			PDL_Double sum = 0;
4134			cp = tmp;
4135			for(i=0; i
4136			dd = 0;
4137			ap = tvec;
4138			for(j=0;j
4139			dd += (ap++) *(cp++);
4140			dd /= blur;
4141			sum += dd * dd;
4142			if(sum > GAUSSIAN_MAXVAL) {
4143			i = ndims; /* exit early if we're too far out */
4144			alpha = 0;
4145			}
4146			}
4147			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
4148			alpha = 0;
4149			} else {
4150			PDL_Indx lodex,hidex;
4151			PDL_Double beta = fabs(zeta * sum);
4152
4153			lodex = beta;
4154			beta -= lodex; hidex = lodex+1;
4155			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
4156
4157			}
4158			}
4159			break;
4160
4161			case 'G':
4162			/* This is the Gaussian rolloff with explicit calculation, preserved */
4163			/* in case someone actually wants the slower longer method. */
4164			{
4165			PDL_Double sum = 0;
4166			cp = tmp;
4167			for(i=0; i
4168			dd = 0;
4169			ap = tvec;
4170			for(j=0;j
4171			dd += (ap++) *(cp++);
4172			dd /= blur;
4173			sum += dd * dd;
4174			if(sum > 4) /* 2 pixels -- four half-widths */
4175			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
4176			}
4177
4178			if(sum > GAUSSIAN_MAXVAL)
4179			alpha = 0;
4180			else
4181			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
4182			}
4183			break;
4184			default:
4185			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
4186			}
4187
4188			{ /* convenience block -- accumulate the current point into the weighted sum. */
4189			/* This is more than simple assignment because we have our own explicit poor */
4190			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
4191			PDL_Indx dat = ((PDL_Indx )(in->data)) + i_off;
4192			PDL_Indx max = out->dims[ndims];
4193			for( i=0; i < max; i++ ) {
4194			if( (badval==0) \|\| (*dat != badval) ) {
4195			acc[i] += dat alpha;
4196			dat += in->dimincs[ndims];
4197			wgt[i] += alpha;
4198			}
4199			wgt2[i] += alpha; }
4200			}
4201			} /* end of t_vio check (i.e. of input accumulation) */
4202
4203
4204			/* Advance input accumulation loop. */
4205			/* We both increment the total vector and also advance the index. */
4206			carry = 1;
4207			for (i=0; i
4208			/* Advance the current element of the offset vector */
4209			ivec[i]++;
4210			j = ivec[i] + ibvec[i];
4211
4212			/* Advance the offset into the data array */
4213			if( j > 0 && j <= in->dims[i]-1 ) {
4214			/* Normal case -- just advance the input vector */
4215			i_off += in->dimincs[i];
4216			} else {
4217			/* Busted a boundary - either before or after. */
4218			switch(bounds[i]){
4219			case 0: /* no breakage allowed -- treat as truncation for interpolation */
4220			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
4221			if( j == 0 )
4222			t_vio--;
4223			else if( j == in->dims[i] )
4224			t_vio++;
4225			break;
4226			case 2: /* extension -- do nothing (so the same input point is re-used) */
4227			break;
4228			case 3: /* periodic -- advance and mod into the allowed range */
4229			if((j % in->dims[i]) == 0) {
4230			i_off -= in->dimincs[i] * (in->dims[i]-1);
4231			} else {
4232			i_off += in->dimincs[i];
4233			}
4234			break;
4235			case 4: /* mirror -- advance or retreat depending on phase */
4236			j += in->dims[i];
4237			j %= (in->dims[i]*2);
4238			j -= in->dims[i];
4239			if( j!=0 && j!= -in->dims[i] ) {
4240			if(j<0)
4241			i_off -= in->dimincs[i];
4242			else
4243			i_off += in->dimincs[i];
4244			}
4245			break;
4246			}
4247			}
4248
4249			/* Now check for carry */
4250			if (ivec[i] <= psize) {
4251			/* Normal case -- copy the current offset to the faster-running dim stashes */
4252			PDL_Indx k;
4253			for(k=0;k
4254			index_stash[k] = i_off;
4255			}
4256			carry = 0;
4257
4258			} else { /* End of this scan -- recover the last position, and mark carry */
4259			i_off = index_stash[i];
4260			if(bounds[i]==1) {
4261			j = ivec[i] + ibvec[i];
4262			if( j < 0 \|\| j >= in->dims[i] )
4263			t_vio--;
4264			ivec[i] = -psize;
4265			j = ivec[i] + ibvec[i];
4266			if( j < 0 \|\| j >= in->dims[i] )
4267			t_vio++;
4268			carry = 1;
4269			} else {
4270			ivec[i] = -psize;
4271			}
4272			}
4273			} /* End of counter-advance loop */
4274			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
4275
4276			{
4277			PDL_Double *ac = acc;
4278			PDL_Double *wg = wgt;
4279			PDL_Double *wg2 = wgt2;
4280			PDL_Indx *dat = out->data;
4281
4282			/* Calculate output vector offset */
4283			for(i=0;i
4284			dat += out->dimincs[i] * ovec[i];
4285
4286			if (!flux) {
4287			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
4288			for (i=0; i < out->dims[ndims]; i++) {
4289			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
4290			ac++; wg++; wg2++;
4291			dat += out->dimincs[ndims];
4292			}
4293			} else {
4294			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
4295			PDL_Double det = tmp[ndims*ndims];
4296			for(i=0; i < out->dims[ndims]; i++) {
4297			if(wg && (wg2 / *wg) < 1.5 ) {
4298			dat = (ac++) / (wg++) det;
4299			wg2++;
4300			} else {
4301			*dat = badval;
4302			ac++; wg++; wg2++;
4303			}
4304			dat += out->dimincs[ndims];
4305			} /* end of for loop */
4306			} /* end of flux flag set conditional */
4307			} /* end of convenience block */
4308
4309			/* End of code for normal pixels */
4310			} else {
4311			/* The pixel was ludicrously huge -- just set this pixel to nan */
4312			PDL_Indx *dat = out->data;
4313			for(i=0;i
4314			dat += out->dimincs[i] * ovec[i];
4315			for(i=0;idims[ndims];i++) {
4316			dat = badval; / Should handle bad values too -- not yet */
4317			dat += out->dimincs[ndims];
4318			}
4319			}
4320
4321			/* Increment the pixel counter */
4322			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
4323			} while (1);
4324			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
4325			#line 4326 "lib/PDL/Transform-pp-map.c"
4326	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
4327	0		} break;
4328	0		case PDL_ULL: {
4329	0	0	PDL_DECLARE_PARAMS_map_1(PDL_ULongLong,P)
		0
		0
4330	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
4331			#line 611 "lib/PDL/Transform.pd"
4332
4333			/*
4334			* Pixel interpolation & averaging code
4335			*
4336			* Calls a common coordinate-transformation block (see following hdr)
4337			* that isn't dependent on the type of the input variable.
4338			*
4339			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
4340			* is handled internally. To simplify the broadcasting business, any
4341			* broadcast dimensions should all be collapsed to a single one by the
4342			* perl front-end.
4343			*
4344			*/
4345
4346			pdl *in = __params->in;
4347			pdl *out = __params->out;
4348			pdl *map = __params->map;
4349
4350			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
4351			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
4352			PDL_Indx ovec[ndims]; /* output pixel loop vector */
4353			PDL_Indx ivec[ndims]; /* input pixel loop vector */
4354			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
4355			PDL_Double dvec[ndims]; /* Residual vector for linearization */
4356			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
4357			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
4358			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
4359			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
4360			char bounds[ndims]; /* Boundary condition packed string */
4361			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
4362			char method; /* Method identifier (gets one of 'h','g') */
4363			PDL_Long big; /* Max size of input footprint for each pix */
4364			PDL_Double blur; /* Scaling of filter */
4365			PDL_Double sv_min; /* minimum singular value */
4366			char flux; /* Flag to indicate flux conservation */
4367			PDL_Indx i;
4368			PDL_ULongLong badval = SvNV(__params->bv);
4369			#define HANNING_LOOKUP_SIZE 2500
4370			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
4371			static int needs_hanning_calc = 1;
4372			PDL_Double zeta = 0;
4373			PDL_Double hanning_offset;
4374
4375			#define GAUSSIAN_LOOKUP_SIZE 4000
4376			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
4377			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
4378			static int needs_gaussian_calc = 1;
4379
4380			PDL_RETERROR(PDL_err, PDL->make_physical(in));
4381			PDL_RETERROR(PDL_err, PDL->make_physical(out));
4382			PDL_RETERROR(PDL_err, PDL->make_physical(map));
4383
4384			/***
4385			* Fill in the boundary condition array
4386			*/
4387			{
4388			char *bstr;
4389			STRLEN blen;
4390			bstr = SvPV(__params->boundary,blen);
4391
4392			if(blen == 0) {
4393			/* If no boundary is specified then every dim gets truncated */
4394			PDL_Indx i;
4395			for (i=0;i
4396			bounds[i] = 1;
4397			} else {
4398			PDL_Indx i;
4399			for(i=0;i
4400			switch(bstr[i < blen ? i : blen-1 ]) {
4401			case '0': case 'f': case 'F': /* forbid */
4402			bounds[i] = 0;
4403			break;
4404			case '1': case 't': case 'T': /* truncate */
4405			bounds[i] = 1;
4406			break;
4407			case '2': case 'e': case 'E': /* extend */
4408			bounds[i] = 2;
4409			break;
4410			case '3': case 'p': case 'P': /* periodic */
4411			bounds[i] = 3;
4412			break;
4413			case '4': case 'm': case 'M': /* mirror */
4414			bounds[i] = 4;
4415			break;
4416			default:
4417			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
4418			break;
4419			}
4420			}
4421			}
4422			}
4423
4424			/***
4425			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
4426			*/
4427			big = labs(__params->big);
4428			if(big <= 0)
4429			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
4430
4431			blur = fabs(__params->blur);
4432			if(blur < 0)
4433			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
4434
4435			sv_min = fabs(__params->sv_min);
4436			if(sv_min < 0)
4437			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
4438
4439			flux = __params->flux != 0;
4440
4441			{
4442			char *mstr;
4443			STRLEN mlen;
4444			mstr = SvPV(__params->method,mlen);
4445
4446			if(mlen==0)
4447			method = 'h';
4448			switch(*mstr) {
4449			case 'h': if( needs_hanning_calc ) {
4450			int i;
4451			for(i=0;i
4452			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
4453			}
4454			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
4455			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
4456			needs_hanning_calc = 0;
4457			}
4458			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
4459			case 'H': method = *mstr;
4460			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
4461			break;
4462
4463			case 'g': case 'j': method = 'g';
4464			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
4465
4466			if( needs_gaussian_calc ) {
4467			int i;
4468			for(i=0;i
4469			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
4470			}
4471			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
4472			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
4473			needs_gaussian_calc = 0;
4474			}
4475			break;
4476
4477			case 'G': case 'J': method = 'G'; break;
4478			default:
4479			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
4480			break;
4481			}
4482			}
4483
4484
4485
4486			/* End of initialization */
4487			/*************************************************************/
4488			/* Start of Real Work */
4489
4490			/* Initialize coordinate vector and map offset
4491			*/
4492			for(i=0;i
4493			ovec[i] = 0;
4494
4495			PDL_Double map_ptr = (PDL_Double )(map->data);
4496			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
4497			for (i=0; i < npdls; i++) offs[i] = 0;
4498			for (i=0; i < ndims; i++) {
4499			incs[i*npdls + 0] = map->dimincs[i+1];
4500			}
4501
4502			/* Main pixel loop (iterates over pixels in the output plane) */
4503			do {
4504			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
4505			/* Prefrobnicate the transformation matrix */
4506			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
4507
4508			#ifdef DEBUG_MAP
4509			{
4510			PDL_Indx k; PDL_Indx foo = 0;
4511			printf("ovec: [");
4512			for(k=0;k
4513			foo += ovec[k] * map->dimincs[k+1];
4514			printf(" %2td ",(PDL_Indx)(ovec[k]));
4515			}
4516			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
4517			for(k=0;k
4518			printf("%8.2f",(double)(((PDL_ULongLong )(map->data))[foo + kmap->dimincs[0]]));
4519			}
4520			printf("]\n");
4521			}
4522			#endif
4523
4524			/* Don't bother accumulating output if psize is too large */
4525			if(psize <= big) {
4526			/* Use the prefrobnicated matrix to generate a local linearization.
4527			* dvec gets the delta; ibvec gets the base.
4528			*/
4529			{
4530			PDL_Double *mp = map_ptr + offs[0];
4531			for (i=0;i
4532			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
4533			mp += map->dimincs[0];
4534			}
4535			}
4536
4537			/* Initialize input delta vector */
4538			for(i=0;i
4539			ivec[i] = -psize;
4540
4541			/* Initialize accumulators */
4542			{
4543			PDL_Double *ac = acc;
4544			for(i=0; i < in->dims[ndims]; i++)
4545			*(ac++) = 0.0;
4546
4547			}
4548			{
4549			PDL_Double *wg = wgt;
4550			for(i=0;i < in->dims[ndims]; i++)
4551			*(wg++) = 0.0;
4552			}
4553			{
4554			PDL_Double *wg = wgt2;
4555			for(i=0;i < in->dims[ndims]; i++)
4556			*(wg++) = 0.0;
4557			}
4558
4559
4560			/*
4561			* Calculate the original offset into the data array, to enable
4562			* delta calculations in the pixel loop
4563			*
4564			* i runs over dims; j holds the working integer index in the
4565			* current dim.
4566			*
4567			* This code matches the incrementation code at the bottom of the accumulation loop
4568			*/
4569
4570			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
4571			i_off = 0;
4572			for(i=0;i
4573			j = ivec[i] + ibvec[i];
4574			if(j<0 \|\| j >= in->dims[i]) {
4575			switch(bounds[i]) {
4576			case 0: /* no breakage allowed */
4577			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
4578			break;
4579			case 1: /* truncation */
4580			t_vio++;
4581			/* fall through */
4582			case 2: /* extension -- crop */
4583			if(j<0)
4584			j=0;
4585			else j = in->dims[i] - 1;
4586			break;
4587			case 3: /* periodic -- mod it */
4588			j %= in->dims[i];
4589			if(j<0)
4590			j += in->dims[i];
4591			break;
4592			case 4: /* mirror -- reflect off the edges */
4593			j += in->dims[i];
4594			j %= (in->dims[i]*2);
4595			if(j<0)
4596			j += in->dims[i]*2;
4597			j -= in->dims[i];
4598			if(j<0) {
4599			j *= -1;
4600			j -= 1;
4601			}
4602			break;
4603			default:
4604			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
4605			break;
4606			}
4607			}
4608			i_off += in->dimincs[i] * j;
4609			}
4610
4611			/* Initialize index stashes for later reference as we scan the footprint */
4612			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
4613			/* end of a dimensional scan. So we stash the index location at the */
4614			/* start of each dimensional scan here. When we finish incrementing */
4615			/* through a particular dim, we pull its value back out of the stash. */
4616			for(i=0;i
4617			index_stash[i] = i_off;
4618			}
4619
4620			/* The input accumulation loop is the hotspot for the entire operation. */
4621			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
4622			/* in the input array, use the linearized transform to bring each pixel center */
4623			/* forward to the output plane, and calculate a weighting based on the chosen */
4624			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
4625			/* table that is initialized the first time through the code. 'H' is the */
4626			/* same process, but explicitly calculated for each interation (~2x slower). */
4627			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
4628			/* into the input array fresh from the current input array vector each time, */
4629			/* we walk through the array using dimincs and the old offset. This saves */
4630			/* about half of the time spent on index calculation. */
4631
4632			do { /* Input accumulation loop */
4633			PDL_Double *cp;
4634			PDL_Double alpha;
4635			/* Calculate the weight of the current input point. Don't bother if we're
4636			* violating any truncation boundaries (in that case our value is zero, but
4637			* for the interpolation we also set the weight to zero).
4638			*/
4639			if( !t_vio ) {
4640
4641			PDL_Double *ap = tvec;
4642			PDL_Double *bp = dvec;
4643			PDL_Indx *ip = ivec;
4644			for(i=0; i
4645			(ap++) = (ip++) - *(bp++);
4646
4647			switch(method) {
4648			PDL_Double dd;
4649			case 'h':
4650			/* This is the Hanning window rolloff. It is a product of a simple */
4651			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
4652			/* is about 2x faster than using cos(theta) directly in each */
4653			/* weighting calculation, so we do. Using 2500 entries and linear */
4654			/* interpolation is accurate to about 10^-7, and should preserve */
4655			/* the contents of cache pretty well. */
4656			alpha = 1;
4657			cp = tmp;
4658			for (i=0; i
4659			PDL_Indx lodex;
4660			PDL_Indx hidex;
4661			PDL_Double beta;
4662			dd = 0;
4663			ap = tvec;
4664			/* Get the matrix-multiply element for this dimension */
4665			for (j=0;j
4666			dd += (ap++) *(cp++);
4667
4668			/* Do linear interpolation from the table */
4669			/* The table captures a hanning window centered 0.5 pixel from center. */
4670			/* We scale the filter by the blur parameter -- but if blur is less */
4671			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
4672			/* value on the pixel edge at 0.5. For blur greater than unity, we */
4673			/* scale simply. */
4674			beta = fabs(dd) - hanning_offset;
4675			if (beta > 0) {
4676			if (beta >= blur) {
4677			alpha = 0;
4678			i = ndims;
4679			} else {
4680			beta *= zeta;
4681			lodex = beta;
4682			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
4683			lodex = HANNING_LOOKUP_SIZE;
4684			hidex = lodex+1;
4685			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
4686			} /* end of interpolation branch */
4687			} /* end of beta > 0 branch */
4688			} /* end of dimension loop */
4689			break;
4690
4691			case 'H':
4692			/* This is the Hanning window rolloff with explicit calculation, preserved */
4693			/* in case someone actually wants the slower longer method. */
4694			alpha = 1;
4695			cp = tmp;
4696			for (i=0; i
4697			dd = 0;
4698			ap = tvec;
4699			for (j=0;j
4700			dd += (ap++) *(cp++);
4701			dd = (fabs(dd) - hanning_offset) / blur;
4702			if ( dd > 1 ) {
4703			alpha = 0;
4704			i = ndims;
4705			} else
4706			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
4707			}
4708			break;
4709
4710			case 'g':
4711			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
4712			{
4713			PDL_Double sum = 0;
4714			cp = tmp;
4715			for(i=0; i
4716			dd = 0;
4717			ap = tvec;
4718			for(j=0;j
4719			dd += (ap++) *(cp++);
4720			dd /= blur;
4721			sum += dd * dd;
4722			if(sum > GAUSSIAN_MAXVAL) {
4723			i = ndims; /* exit early if we're too far out */
4724			alpha = 0;
4725			}
4726			}
4727			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
4728			alpha = 0;
4729			} else {
4730			PDL_Indx lodex,hidex;
4731			PDL_Double beta = fabs(zeta * sum);
4732
4733			lodex = beta;
4734			beta -= lodex; hidex = lodex+1;
4735			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
4736
4737			}
4738			}
4739			break;
4740
4741			case 'G':
4742			/* This is the Gaussian rolloff with explicit calculation, preserved */
4743			/* in case someone actually wants the slower longer method. */
4744			{
4745			PDL_Double sum = 0;
4746			cp = tmp;
4747			for(i=0; i
4748			dd = 0;
4749			ap = tvec;
4750			for(j=0;j
4751			dd += (ap++) *(cp++);
4752			dd /= blur;
4753			sum += dd * dd;
4754			if(sum > 4) /* 2 pixels -- four half-widths */
4755			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
4756			}
4757
4758			if(sum > GAUSSIAN_MAXVAL)
4759			alpha = 0;
4760			else
4761			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
4762			}
4763			break;
4764			default:
4765			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
4766			}
4767
4768			{ /* convenience block -- accumulate the current point into the weighted sum. */
4769			/* This is more than simple assignment because we have our own explicit poor */
4770			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
4771			PDL_ULongLong dat = ((PDL_ULongLong )(in->data)) + i_off;
4772			PDL_Indx max = out->dims[ndims];
4773			for( i=0; i < max; i++ ) {
4774			if( (badval==0) \|\| (*dat != badval) ) {
4775			acc[i] += dat alpha;
4776			dat += in->dimincs[ndims];
4777			wgt[i] += alpha;
4778			}
4779			wgt2[i] += alpha; }
4780			}
4781			} /* end of t_vio check (i.e. of input accumulation) */
4782
4783
4784			/* Advance input accumulation loop. */
4785			/* We both increment the total vector and also advance the index. */
4786			carry = 1;
4787			for (i=0; i
4788			/* Advance the current element of the offset vector */
4789			ivec[i]++;
4790			j = ivec[i] + ibvec[i];
4791
4792			/* Advance the offset into the data array */
4793			if( j > 0 && j <= in->dims[i]-1 ) {
4794			/* Normal case -- just advance the input vector */
4795			i_off += in->dimincs[i];
4796			} else {
4797			/* Busted a boundary - either before or after. */
4798			switch(bounds[i]){
4799			case 0: /* no breakage allowed -- treat as truncation for interpolation */
4800			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
4801			if( j == 0 )
4802			t_vio--;
4803			else if( j == in->dims[i] )
4804			t_vio++;
4805			break;
4806			case 2: /* extension -- do nothing (so the same input point is re-used) */
4807			break;
4808			case 3: /* periodic -- advance and mod into the allowed range */
4809			if((j % in->dims[i]) == 0) {
4810			i_off -= in->dimincs[i] * (in->dims[i]-1);
4811			} else {
4812			i_off += in->dimincs[i];
4813			}
4814			break;
4815			case 4: /* mirror -- advance or retreat depending on phase */
4816			j += in->dims[i];
4817			j %= (in->dims[i]*2);
4818			j -= in->dims[i];
4819			if( j!=0 && j!= -in->dims[i] ) {
4820			if(j<0)
4821			i_off -= in->dimincs[i];
4822			else
4823			i_off += in->dimincs[i];
4824			}
4825			break;
4826			}
4827			}
4828
4829			/* Now check for carry */
4830			if (ivec[i] <= psize) {
4831			/* Normal case -- copy the current offset to the faster-running dim stashes */
4832			PDL_Indx k;
4833			for(k=0;k
4834			index_stash[k] = i_off;
4835			}
4836			carry = 0;
4837
4838			} else { /* End of this scan -- recover the last position, and mark carry */
4839			i_off = index_stash[i];
4840			if(bounds[i]==1) {
4841			j = ivec[i] + ibvec[i];
4842			if( j < 0 \|\| j >= in->dims[i] )
4843			t_vio--;
4844			ivec[i] = -psize;
4845			j = ivec[i] + ibvec[i];
4846			if( j < 0 \|\| j >= in->dims[i] )
4847			t_vio++;
4848			carry = 1;
4849			} else {
4850			ivec[i] = -psize;
4851			}
4852			}
4853			} /* End of counter-advance loop */
4854			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
4855
4856			{
4857			PDL_Double *ac = acc;
4858			PDL_Double *wg = wgt;
4859			PDL_Double *wg2 = wgt2;
4860			PDL_ULongLong *dat = out->data;
4861
4862			/* Calculate output vector offset */
4863			for(i=0;i
4864			dat += out->dimincs[i] * ovec[i];
4865
4866			if (!flux) {
4867			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
4868			for (i=0; i < out->dims[ndims]; i++) {
4869			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
4870			ac++; wg++; wg2++;
4871			dat += out->dimincs[ndims];
4872			}
4873			} else {
4874			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
4875			PDL_Double det = tmp[ndims*ndims];
4876			for(i=0; i < out->dims[ndims]; i++) {
4877			if(wg && (wg2 / *wg) < 1.5 ) {
4878			dat = (ac++) / (wg++) det;
4879			wg2++;
4880			} else {
4881			*dat = badval;
4882			ac++; wg++; wg2++;
4883			}
4884			dat += out->dimincs[ndims];
4885			} /* end of for loop */
4886			} /* end of flux flag set conditional */
4887			} /* end of convenience block */
4888
4889			/* End of code for normal pixels */
4890			} else {
4891			/* The pixel was ludicrously huge -- just set this pixel to nan */
4892			PDL_ULongLong *dat = out->data;
4893			for(i=0;i
4894			dat += out->dimincs[i] * ovec[i];
4895			for(i=0;idims[ndims];i++) {
4896			dat = badval; / Should handle bad values too -- not yet */
4897			dat += out->dimincs[ndims];
4898			}
4899			}
4900
4901			/* Increment the pixel counter */
4902			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
4903			} while (1);
4904			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
4905			#line 4906 "lib/PDL/Transform-pp-map.c"
4906	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
4907	0		} break;
4908	0		case PDL_LL: {
4909	0	0	PDL_DECLARE_PARAMS_map_1(PDL_LongLong,Q)
		0
		0
4910	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
4911			#line 611 "lib/PDL/Transform.pd"
4912
4913			/*
4914			* Pixel interpolation & averaging code
4915			*
4916			* Calls a common coordinate-transformation block (see following hdr)
4917			* that isn't dependent on the type of the input variable.
4918			*
4919			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
4920			* is handled internally. To simplify the broadcasting business, any
4921			* broadcast dimensions should all be collapsed to a single one by the
4922			* perl front-end.
4923			*
4924			*/
4925
4926			pdl *in = __params->in;
4927			pdl *out = __params->out;
4928			pdl *map = __params->map;
4929
4930			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
4931			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
4932			PDL_Indx ovec[ndims]; /* output pixel loop vector */
4933			PDL_Indx ivec[ndims]; /* input pixel loop vector */
4934			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
4935			PDL_Double dvec[ndims]; /* Residual vector for linearization */
4936			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
4937			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
4938			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
4939			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
4940			char bounds[ndims]; /* Boundary condition packed string */
4941			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
4942			char method; /* Method identifier (gets one of 'h','g') */
4943			PDL_Long big; /* Max size of input footprint for each pix */
4944			PDL_Double blur; /* Scaling of filter */
4945			PDL_Double sv_min; /* minimum singular value */
4946			char flux; /* Flag to indicate flux conservation */
4947			PDL_Indx i;
4948			PDL_LongLong badval = SvNV(__params->bv);
4949			#define HANNING_LOOKUP_SIZE 2500
4950			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
4951			static int needs_hanning_calc = 1;
4952			PDL_Double zeta = 0;
4953			PDL_Double hanning_offset;
4954
4955			#define GAUSSIAN_LOOKUP_SIZE 4000
4956			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
4957			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
4958			static int needs_gaussian_calc = 1;
4959
4960			PDL_RETERROR(PDL_err, PDL->make_physical(in));
4961			PDL_RETERROR(PDL_err, PDL->make_physical(out));
4962			PDL_RETERROR(PDL_err, PDL->make_physical(map));
4963
4964			/***
4965			* Fill in the boundary condition array
4966			*/
4967			{
4968			char *bstr;
4969			STRLEN blen;
4970			bstr = SvPV(__params->boundary,blen);
4971
4972			if(blen == 0) {
4973			/* If no boundary is specified then every dim gets truncated */
4974			PDL_Indx i;
4975			for (i=0;i
4976			bounds[i] = 1;
4977			} else {
4978			PDL_Indx i;
4979			for(i=0;i
4980			switch(bstr[i < blen ? i : blen-1 ]) {
4981			case '0': case 'f': case 'F': /* forbid */
4982			bounds[i] = 0;
4983			break;
4984			case '1': case 't': case 'T': /* truncate */
4985			bounds[i] = 1;
4986			break;
4987			case '2': case 'e': case 'E': /* extend */
4988			bounds[i] = 2;
4989			break;
4990			case '3': case 'p': case 'P': /* periodic */
4991			bounds[i] = 3;
4992			break;
4993			case '4': case 'm': case 'M': /* mirror */
4994			bounds[i] = 4;
4995			break;
4996			default:
4997			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
4998			break;
4999			}
5000			}
5001			}
5002			}
5003
5004			/***
5005			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
5006			*/
5007			big = labs(__params->big);
5008			if(big <= 0)
5009			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
5010
5011			blur = fabs(__params->blur);
5012			if(blur < 0)
5013			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
5014
5015			sv_min = fabs(__params->sv_min);
5016			if(sv_min < 0)
5017			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
5018
5019			flux = __params->flux != 0;
5020
5021			{
5022			char *mstr;
5023			STRLEN mlen;
5024			mstr = SvPV(__params->method,mlen);
5025
5026			if(mlen==0)
5027			method = 'h';
5028			switch(*mstr) {
5029			case 'h': if( needs_hanning_calc ) {
5030			int i;
5031			for(i=0;i
5032			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
5033			}
5034			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
5035			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
5036			needs_hanning_calc = 0;
5037			}
5038			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
5039			case 'H': method = *mstr;
5040			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
5041			break;
5042
5043			case 'g': case 'j': method = 'g';
5044			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
5045
5046			if( needs_gaussian_calc ) {
5047			int i;
5048			for(i=0;i
5049			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
5050			}
5051			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
5052			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
5053			needs_gaussian_calc = 0;
5054			}
5055			break;
5056
5057			case 'G': case 'J': method = 'G'; break;
5058			default:
5059			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
5060			break;
5061			}
5062			}
5063
5064
5065
5066			/* End of initialization */
5067			/*************************************************************/
5068			/* Start of Real Work */
5069
5070			/* Initialize coordinate vector and map offset
5071			*/
5072			for(i=0;i
5073			ovec[i] = 0;
5074
5075			PDL_Double map_ptr = (PDL_Double )(map->data);
5076			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
5077			for (i=0; i < npdls; i++) offs[i] = 0;
5078			for (i=0; i < ndims; i++) {
5079			incs[i*npdls + 0] = map->dimincs[i+1];
5080			}
5081
5082			/* Main pixel loop (iterates over pixels in the output plane) */
5083			do {
5084			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
5085			/* Prefrobnicate the transformation matrix */
5086			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
5087
5088			#ifdef DEBUG_MAP
5089			{
5090			PDL_Indx k; PDL_Indx foo = 0;
5091			printf("ovec: [");
5092			for(k=0;k
5093			foo += ovec[k] * map->dimincs[k+1];
5094			printf(" %2td ",(PDL_Indx)(ovec[k]));
5095			}
5096			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
5097			for(k=0;k
5098			printf("%8.2f",(double)(((PDL_LongLong )(map->data))[foo + kmap->dimincs[0]]));
5099			}
5100			printf("]\n");
5101			}
5102			#endif
5103
5104			/* Don't bother accumulating output if psize is too large */
5105			if(psize <= big) {
5106			/* Use the prefrobnicated matrix to generate a local linearization.
5107			* dvec gets the delta; ibvec gets the base.
5108			*/
5109			{
5110			PDL_Double *mp = map_ptr + offs[0];
5111			for (i=0;i
5112			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
5113			mp += map->dimincs[0];
5114			}
5115			}
5116
5117			/* Initialize input delta vector */
5118			for(i=0;i
5119			ivec[i] = -psize;
5120
5121			/* Initialize accumulators */
5122			{
5123			PDL_Double *ac = acc;
5124			for(i=0; i < in->dims[ndims]; i++)
5125			*(ac++) = 0.0;
5126
5127			}
5128			{
5129			PDL_Double *wg = wgt;
5130			for(i=0;i < in->dims[ndims]; i++)
5131			*(wg++) = 0.0;
5132			}
5133			{
5134			PDL_Double *wg = wgt2;
5135			for(i=0;i < in->dims[ndims]; i++)
5136			*(wg++) = 0.0;
5137			}
5138
5139
5140			/*
5141			* Calculate the original offset into the data array, to enable
5142			* delta calculations in the pixel loop
5143			*
5144			* i runs over dims; j holds the working integer index in the
5145			* current dim.
5146			*
5147			* This code matches the incrementation code at the bottom of the accumulation loop
5148			*/
5149
5150			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
5151			i_off = 0;
5152			for(i=0;i
5153			j = ivec[i] + ibvec[i];
5154			if(j<0 \|\| j >= in->dims[i]) {
5155			switch(bounds[i]) {
5156			case 0: /* no breakage allowed */
5157			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
5158			break;
5159			case 1: /* truncation */
5160			t_vio++;
5161			/* fall through */
5162			case 2: /* extension -- crop */
5163			if(j<0)
5164			j=0;
5165			else j = in->dims[i] - 1;
5166			break;
5167			case 3: /* periodic -- mod it */
5168			j %= in->dims[i];
5169			if(j<0)
5170			j += in->dims[i];
5171			break;
5172			case 4: /* mirror -- reflect off the edges */
5173			j += in->dims[i];
5174			j %= (in->dims[i]*2);
5175			if(j<0)
5176			j += in->dims[i]*2;
5177			j -= in->dims[i];
5178			if(j<0) {
5179			j *= -1;
5180			j -= 1;
5181			}
5182			break;
5183			default:
5184			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
5185			break;
5186			}
5187			}
5188			i_off += in->dimincs[i] * j;
5189			}
5190
5191			/* Initialize index stashes for later reference as we scan the footprint */
5192			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
5193			/* end of a dimensional scan. So we stash the index location at the */
5194			/* start of each dimensional scan here. When we finish incrementing */
5195			/* through a particular dim, we pull its value back out of the stash. */
5196			for(i=0;i
5197			index_stash[i] = i_off;
5198			}
5199
5200			/* The input accumulation loop is the hotspot for the entire operation. */
5201			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
5202			/* in the input array, use the linearized transform to bring each pixel center */
5203			/* forward to the output plane, and calculate a weighting based on the chosen */
5204			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
5205			/* table that is initialized the first time through the code. 'H' is the */
5206			/* same process, but explicitly calculated for each interation (~2x slower). */
5207			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
5208			/* into the input array fresh from the current input array vector each time, */
5209			/* we walk through the array using dimincs and the old offset. This saves */
5210			/* about half of the time spent on index calculation. */
5211
5212			do { /* Input accumulation loop */
5213			PDL_Double *cp;
5214			PDL_Double alpha;
5215			/* Calculate the weight of the current input point. Don't bother if we're
5216			* violating any truncation boundaries (in that case our value is zero, but
5217			* for the interpolation we also set the weight to zero).
5218			*/
5219			if( !t_vio ) {
5220
5221			PDL_Double *ap = tvec;
5222			PDL_Double *bp = dvec;
5223			PDL_Indx *ip = ivec;
5224			for(i=0; i
5225			(ap++) = (ip++) - *(bp++);
5226
5227			switch(method) {
5228			PDL_Double dd;
5229			case 'h':
5230			/* This is the Hanning window rolloff. It is a product of a simple */
5231			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
5232			/* is about 2x faster than using cos(theta) directly in each */
5233			/* weighting calculation, so we do. Using 2500 entries and linear */
5234			/* interpolation is accurate to about 10^-7, and should preserve */
5235			/* the contents of cache pretty well. */
5236			alpha = 1;
5237			cp = tmp;
5238			for (i=0; i
5239			PDL_Indx lodex;
5240			PDL_Indx hidex;
5241			PDL_Double beta;
5242			dd = 0;
5243			ap = tvec;
5244			/* Get the matrix-multiply element for this dimension */
5245			for (j=0;j
5246			dd += (ap++) *(cp++);
5247
5248			/* Do linear interpolation from the table */
5249			/* The table captures a hanning window centered 0.5 pixel from center. */
5250			/* We scale the filter by the blur parameter -- but if blur is less */
5251			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
5252			/* value on the pixel edge at 0.5. For blur greater than unity, we */
5253			/* scale simply. */
5254			beta = fabs(dd) - hanning_offset;
5255			if (beta > 0) {
5256			if (beta >= blur) {
5257			alpha = 0;
5258			i = ndims;
5259			} else {
5260			beta *= zeta;
5261			lodex = beta;
5262			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
5263			lodex = HANNING_LOOKUP_SIZE;
5264			hidex = lodex+1;
5265			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
5266			} /* end of interpolation branch */
5267			} /* end of beta > 0 branch */
5268			} /* end of dimension loop */
5269			break;
5270
5271			case 'H':
5272			/* This is the Hanning window rolloff with explicit calculation, preserved */
5273			/* in case someone actually wants the slower longer method. */
5274			alpha = 1;
5275			cp = tmp;
5276			for (i=0; i
5277			dd = 0;
5278			ap = tvec;
5279			for (j=0;j
5280			dd += (ap++) *(cp++);
5281			dd = (fabs(dd) - hanning_offset) / blur;
5282			if ( dd > 1 ) {
5283			alpha = 0;
5284			i = ndims;
5285			} else
5286			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
5287			}
5288			break;
5289
5290			case 'g':
5291			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
5292			{
5293			PDL_Double sum = 0;
5294			cp = tmp;
5295			for(i=0; i
5296			dd = 0;
5297			ap = tvec;
5298			for(j=0;j
5299			dd += (ap++) *(cp++);
5300			dd /= blur;
5301			sum += dd * dd;
5302			if(sum > GAUSSIAN_MAXVAL) {
5303			i = ndims; /* exit early if we're too far out */
5304			alpha = 0;
5305			}
5306			}
5307			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
5308			alpha = 0;
5309			} else {
5310			PDL_Indx lodex,hidex;
5311			PDL_Double beta = fabs(zeta * sum);
5312
5313			lodex = beta;
5314			beta -= lodex; hidex = lodex+1;
5315			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
5316
5317			}
5318			}
5319			break;
5320
5321			case 'G':
5322			/* This is the Gaussian rolloff with explicit calculation, preserved */
5323			/* in case someone actually wants the slower longer method. */
5324			{
5325			PDL_Double sum = 0;
5326			cp = tmp;
5327			for(i=0; i
5328			dd = 0;
5329			ap = tvec;
5330			for(j=0;j
5331			dd += (ap++) *(cp++);
5332			dd /= blur;
5333			sum += dd * dd;
5334			if(sum > 4) /* 2 pixels -- four half-widths */
5335			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
5336			}
5337
5338			if(sum > GAUSSIAN_MAXVAL)
5339			alpha = 0;
5340			else
5341			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
5342			}
5343			break;
5344			default:
5345			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
5346			}
5347
5348			{ /* convenience block -- accumulate the current point into the weighted sum. */
5349			/* This is more than simple assignment because we have our own explicit poor */
5350			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
5351			PDL_LongLong dat = ((PDL_LongLong )(in->data)) + i_off;
5352			PDL_Indx max = out->dims[ndims];
5353			for( i=0; i < max; i++ ) {
5354			if( (badval==0) \|\| (*dat != badval) ) {
5355			acc[i] += dat alpha;
5356			dat += in->dimincs[ndims];
5357			wgt[i] += alpha;
5358			}
5359			wgt2[i] += alpha; }
5360			}
5361			} /* end of t_vio check (i.e. of input accumulation) */
5362
5363
5364			/* Advance input accumulation loop. */
5365			/* We both increment the total vector and also advance the index. */
5366			carry = 1;
5367			for (i=0; i
5368			/* Advance the current element of the offset vector */
5369			ivec[i]++;
5370			j = ivec[i] + ibvec[i];
5371
5372			/* Advance the offset into the data array */
5373			if( j > 0 && j <= in->dims[i]-1 ) {
5374			/* Normal case -- just advance the input vector */
5375			i_off += in->dimincs[i];
5376			} else {
5377			/* Busted a boundary - either before or after. */
5378			switch(bounds[i]){
5379			case 0: /* no breakage allowed -- treat as truncation for interpolation */
5380			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
5381			if( j == 0 )
5382			t_vio--;
5383			else if( j == in->dims[i] )
5384			t_vio++;
5385			break;
5386			case 2: /* extension -- do nothing (so the same input point is re-used) */
5387			break;
5388			case 3: /* periodic -- advance and mod into the allowed range */
5389			if((j % in->dims[i]) == 0) {
5390			i_off -= in->dimincs[i] * (in->dims[i]-1);
5391			} else {
5392			i_off += in->dimincs[i];
5393			}
5394			break;
5395			case 4: /* mirror -- advance or retreat depending on phase */
5396			j += in->dims[i];
5397			j %= (in->dims[i]*2);
5398			j -= in->dims[i];
5399			if( j!=0 && j!= -in->dims[i] ) {
5400			if(j<0)
5401			i_off -= in->dimincs[i];
5402			else
5403			i_off += in->dimincs[i];
5404			}
5405			break;
5406			}
5407			}
5408
5409			/* Now check for carry */
5410			if (ivec[i] <= psize) {
5411			/* Normal case -- copy the current offset to the faster-running dim stashes */
5412			PDL_Indx k;
5413			for(k=0;k
5414			index_stash[k] = i_off;
5415			}
5416			carry = 0;
5417
5418			} else { /* End of this scan -- recover the last position, and mark carry */
5419			i_off = index_stash[i];
5420			if(bounds[i]==1) {
5421			j = ivec[i] + ibvec[i];
5422			if( j < 0 \|\| j >= in->dims[i] )
5423			t_vio--;
5424			ivec[i] = -psize;
5425			j = ivec[i] + ibvec[i];
5426			if( j < 0 \|\| j >= in->dims[i] )
5427			t_vio++;
5428			carry = 1;
5429			} else {
5430			ivec[i] = -psize;
5431			}
5432			}
5433			} /* End of counter-advance loop */
5434			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
5435
5436			{
5437			PDL_Double *ac = acc;
5438			PDL_Double *wg = wgt;
5439			PDL_Double *wg2 = wgt2;
5440			PDL_LongLong *dat = out->data;
5441
5442			/* Calculate output vector offset */
5443			for(i=0;i
5444			dat += out->dimincs[i] * ovec[i];
5445
5446			if (!flux) {
5447			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
5448			for (i=0; i < out->dims[ndims]; i++) {
5449			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
5450			ac++; wg++; wg2++;
5451			dat += out->dimincs[ndims];
5452			}
5453			} else {
5454			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
5455			PDL_Double det = tmp[ndims*ndims];
5456			for(i=0; i < out->dims[ndims]; i++) {
5457			if(wg && (wg2 / *wg) < 1.5 ) {
5458			dat = (ac++) / (wg++) det;
5459			wg2++;
5460			} else {
5461			*dat = badval;
5462			ac++; wg++; wg2++;
5463			}
5464			dat += out->dimincs[ndims];
5465			} /* end of for loop */
5466			} /* end of flux flag set conditional */
5467			} /* end of convenience block */
5468
5469			/* End of code for normal pixels */
5470			} else {
5471			/* The pixel was ludicrously huge -- just set this pixel to nan */
5472			PDL_LongLong *dat = out->data;
5473			for(i=0;i
5474			dat += out->dimincs[i] * ovec[i];
5475			for(i=0;idims[ndims];i++) {
5476			dat = badval; / Should handle bad values too -- not yet */
5477			dat += out->dimincs[ndims];
5478			}
5479			}
5480
5481			/* Increment the pixel counter */
5482			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
5483			} while (1);
5484			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
5485			#line 5486 "lib/PDL/Transform-pp-map.c"
5486	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
5487	0		} break;
5488	0		case PDL_F: {
5489	0	0	PDL_DECLARE_PARAMS_map_1(PDL_Float,F)
		0
		0
5490	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
5491			#line 611 "lib/PDL/Transform.pd"
5492
5493			/*
5494			* Pixel interpolation & averaging code
5495			*
5496			* Calls a common coordinate-transformation block (see following hdr)
5497			* that isn't dependent on the type of the input variable.
5498			*
5499			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
5500			* is handled internally. To simplify the broadcasting business, any
5501			* broadcast dimensions should all be collapsed to a single one by the
5502			* perl front-end.
5503			*
5504			*/
5505
5506			pdl *in = __params->in;
5507			pdl *out = __params->out;
5508			pdl *map = __params->map;
5509
5510			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
5511			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
5512			PDL_Indx ovec[ndims]; /* output pixel loop vector */
5513			PDL_Indx ivec[ndims]; /* input pixel loop vector */
5514			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
5515			PDL_Double dvec[ndims]; /* Residual vector for linearization */
5516			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
5517			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
5518			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
5519			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
5520			char bounds[ndims]; /* Boundary condition packed string */
5521			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
5522			char method; /* Method identifier (gets one of 'h','g') */
5523			PDL_Long big; /* Max size of input footprint for each pix */
5524			PDL_Double blur; /* Scaling of filter */
5525			PDL_Double sv_min; /* minimum singular value */
5526			char flux; /* Flag to indicate flux conservation */
5527			PDL_Indx i;
5528			PDL_Float badval = SvNV(__params->bv);
5529			#define HANNING_LOOKUP_SIZE 2500
5530			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
5531			static int needs_hanning_calc = 1;
5532			PDL_Double zeta = 0;
5533			PDL_Double hanning_offset;
5534
5535			#define GAUSSIAN_LOOKUP_SIZE 4000
5536			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
5537			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
5538			static int needs_gaussian_calc = 1;
5539
5540			PDL_RETERROR(PDL_err, PDL->make_physical(in));
5541			PDL_RETERROR(PDL_err, PDL->make_physical(out));
5542			PDL_RETERROR(PDL_err, PDL->make_physical(map));
5543
5544			/***
5545			* Fill in the boundary condition array
5546			*/
5547			{
5548			char *bstr;
5549			STRLEN blen;
5550			bstr = SvPV(__params->boundary,blen);
5551
5552			if(blen == 0) {
5553			/* If no boundary is specified then every dim gets truncated */
5554			PDL_Indx i;
5555			for (i=0;i
5556			bounds[i] = 1;
5557			} else {
5558			PDL_Indx i;
5559			for(i=0;i
5560			switch(bstr[i < blen ? i : blen-1 ]) {
5561			case '0': case 'f': case 'F': /* forbid */
5562			bounds[i] = 0;
5563			break;
5564			case '1': case 't': case 'T': /* truncate */
5565			bounds[i] = 1;
5566			break;
5567			case '2': case 'e': case 'E': /* extend */
5568			bounds[i] = 2;
5569			break;
5570			case '3': case 'p': case 'P': /* periodic */
5571			bounds[i] = 3;
5572			break;
5573			case '4': case 'm': case 'M': /* mirror */
5574			bounds[i] = 4;
5575			break;
5576			default:
5577			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
5578			break;
5579			}
5580			}
5581			}
5582			}
5583
5584			/***
5585			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
5586			*/
5587			big = labs(__params->big);
5588			if(big <= 0)
5589			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
5590
5591			blur = fabs(__params->blur);
5592			if(blur < 0)
5593			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
5594
5595			sv_min = fabs(__params->sv_min);
5596			if(sv_min < 0)
5597			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
5598
5599			flux = __params->flux != 0;
5600
5601			{
5602			char *mstr;
5603			STRLEN mlen;
5604			mstr = SvPV(__params->method,mlen);
5605
5606			if(mlen==0)
5607			method = 'h';
5608			switch(*mstr) {
5609			case 'h': if( needs_hanning_calc ) {
5610			int i;
5611			for(i=0;i
5612			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
5613			}
5614			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
5615			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
5616			needs_hanning_calc = 0;
5617			}
5618			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
5619			case 'H': method = *mstr;
5620			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
5621			break;
5622
5623			case 'g': case 'j': method = 'g';
5624			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
5625
5626			if( needs_gaussian_calc ) {
5627			int i;
5628			for(i=0;i
5629			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
5630			}
5631			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
5632			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
5633			needs_gaussian_calc = 0;
5634			}
5635			break;
5636
5637			case 'G': case 'J': method = 'G'; break;
5638			default:
5639			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
5640			break;
5641			}
5642			}
5643
5644
5645
5646			/* End of initialization */
5647			/*************************************************************/
5648			/* Start of Real Work */
5649
5650			/* Initialize coordinate vector and map offset
5651			*/
5652			for(i=0;i
5653			ovec[i] = 0;
5654
5655			PDL_Double map_ptr = (PDL_Double )(map->data);
5656			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
5657			for (i=0; i < npdls; i++) offs[i] = 0;
5658			for (i=0; i < ndims; i++) {
5659			incs[i*npdls + 0] = map->dimincs[i+1];
5660			}
5661
5662			/* Main pixel loop (iterates over pixels in the output plane) */
5663			do {
5664			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
5665			/* Prefrobnicate the transformation matrix */
5666			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
5667
5668			#ifdef DEBUG_MAP
5669			{
5670			PDL_Indx k; PDL_Indx foo = 0;
5671			printf("ovec: [");
5672			for(k=0;k
5673			foo += ovec[k] * map->dimincs[k+1];
5674			printf(" %2td ",(PDL_Indx)(ovec[k]));
5675			}
5676			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
5677			for(k=0;k
5678			printf("%8.2f",(double)(((PDL_Float )(map->data))[foo + kmap->dimincs[0]]));
5679			}
5680			printf("]\n");
5681			}
5682			#endif
5683
5684			/* Don't bother accumulating output if psize is too large */
5685			if(psize <= big) {
5686			/* Use the prefrobnicated matrix to generate a local linearization.
5687			* dvec gets the delta; ibvec gets the base.
5688			*/
5689			{
5690			PDL_Double *mp = map_ptr + offs[0];
5691			for (i=0;i
5692			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
5693			mp += map->dimincs[0];
5694			}
5695			}
5696
5697			/* Initialize input delta vector */
5698			for(i=0;i
5699			ivec[i] = -psize;
5700
5701			/* Initialize accumulators */
5702			{
5703			PDL_Double *ac = acc;
5704			for(i=0; i < in->dims[ndims]; i++)
5705			*(ac++) = 0.0;
5706
5707			}
5708			{
5709			PDL_Double *wg = wgt;
5710			for(i=0;i < in->dims[ndims]; i++)
5711			*(wg++) = 0.0;
5712			}
5713			{
5714			PDL_Double *wg = wgt2;
5715			for(i=0;i < in->dims[ndims]; i++)
5716			*(wg++) = 0.0;
5717			}
5718
5719
5720			/*
5721			* Calculate the original offset into the data array, to enable
5722			* delta calculations in the pixel loop
5723			*
5724			* i runs over dims; j holds the working integer index in the
5725			* current dim.
5726			*
5727			* This code matches the incrementation code at the bottom of the accumulation loop
5728			*/
5729
5730			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
5731			i_off = 0;
5732			for(i=0;i
5733			j = ivec[i] + ibvec[i];
5734			if(j<0 \|\| j >= in->dims[i]) {
5735			switch(bounds[i]) {
5736			case 0: /* no breakage allowed */
5737			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
5738			break;
5739			case 1: /* truncation */
5740			t_vio++;
5741			/* fall through */
5742			case 2: /* extension -- crop */
5743			if(j<0)
5744			j=0;
5745			else j = in->dims[i] - 1;
5746			break;
5747			case 3: /* periodic -- mod it */
5748			j %= in->dims[i];
5749			if(j<0)
5750			j += in->dims[i];
5751			break;
5752			case 4: /* mirror -- reflect off the edges */
5753			j += in->dims[i];
5754			j %= (in->dims[i]*2);
5755			if(j<0)
5756			j += in->dims[i]*2;
5757			j -= in->dims[i];
5758			if(j<0) {
5759			j *= -1;
5760			j -= 1;
5761			}
5762			break;
5763			default:
5764			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
5765			break;
5766			}
5767			}
5768			i_off += in->dimincs[i] * j;
5769			}
5770
5771			/* Initialize index stashes for later reference as we scan the footprint */
5772			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
5773			/* end of a dimensional scan. So we stash the index location at the */
5774			/* start of each dimensional scan here. When we finish incrementing */
5775			/* through a particular dim, we pull its value back out of the stash. */
5776			for(i=0;i
5777			index_stash[i] = i_off;
5778			}
5779
5780			/* The input accumulation loop is the hotspot for the entire operation. */
5781			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
5782			/* in the input array, use the linearized transform to bring each pixel center */
5783			/* forward to the output plane, and calculate a weighting based on the chosen */
5784			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
5785			/* table that is initialized the first time through the code. 'H' is the */
5786			/* same process, but explicitly calculated for each interation (~2x slower). */
5787			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
5788			/* into the input array fresh from the current input array vector each time, */
5789			/* we walk through the array using dimincs and the old offset. This saves */
5790			/* about half of the time spent on index calculation. */
5791
5792			do { /* Input accumulation loop */
5793			PDL_Double *cp;
5794			PDL_Double alpha;
5795			/* Calculate the weight of the current input point. Don't bother if we're
5796			* violating any truncation boundaries (in that case our value is zero, but
5797			* for the interpolation we also set the weight to zero).
5798			*/
5799			if( !t_vio ) {
5800
5801			PDL_Double *ap = tvec;
5802			PDL_Double *bp = dvec;
5803			PDL_Indx *ip = ivec;
5804			for(i=0; i
5805			(ap++) = (ip++) - *(bp++);
5806
5807			switch(method) {
5808			PDL_Double dd;
5809			case 'h':
5810			/* This is the Hanning window rolloff. It is a product of a simple */
5811			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
5812			/* is about 2x faster than using cos(theta) directly in each */
5813			/* weighting calculation, so we do. Using 2500 entries and linear */
5814			/* interpolation is accurate to about 10^-7, and should preserve */
5815			/* the contents of cache pretty well. */
5816			alpha = 1;
5817			cp = tmp;
5818			for (i=0; i
5819			PDL_Indx lodex;
5820			PDL_Indx hidex;
5821			PDL_Double beta;
5822			dd = 0;
5823			ap = tvec;
5824			/* Get the matrix-multiply element for this dimension */
5825			for (j=0;j
5826			dd += (ap++) *(cp++);
5827
5828			/* Do linear interpolation from the table */
5829			/* The table captures a hanning window centered 0.5 pixel from center. */
5830			/* We scale the filter by the blur parameter -- but if blur is less */
5831			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
5832			/* value on the pixel edge at 0.5. For blur greater than unity, we */
5833			/* scale simply. */
5834			beta = fabs(dd) - hanning_offset;
5835			if (beta > 0) {
5836			if (beta >= blur) {
5837			alpha = 0;
5838			i = ndims;
5839			} else {
5840			beta *= zeta;
5841			lodex = beta;
5842			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
5843			lodex = HANNING_LOOKUP_SIZE;
5844			hidex = lodex+1;
5845			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
5846			} /* end of interpolation branch */
5847			} /* end of beta > 0 branch */
5848			} /* end of dimension loop */
5849			break;
5850
5851			case 'H':
5852			/* This is the Hanning window rolloff with explicit calculation, preserved */
5853			/* in case someone actually wants the slower longer method. */
5854			alpha = 1;
5855			cp = tmp;
5856			for (i=0; i
5857			dd = 0;
5858			ap = tvec;
5859			for (j=0;j
5860			dd += (ap++) *(cp++);
5861			dd = (fabs(dd) - hanning_offset) / blur;
5862			if ( dd > 1 ) {
5863			alpha = 0;
5864			i = ndims;
5865			} else
5866			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
5867			}
5868			break;
5869
5870			case 'g':
5871			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
5872			{
5873			PDL_Double sum = 0;
5874			cp = tmp;
5875			for(i=0; i
5876			dd = 0;
5877			ap = tvec;
5878			for(j=0;j
5879			dd += (ap++) *(cp++);
5880			dd /= blur;
5881			sum += dd * dd;
5882			if(sum > GAUSSIAN_MAXVAL) {
5883			i = ndims; /* exit early if we're too far out */
5884			alpha = 0;
5885			}
5886			}
5887			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
5888			alpha = 0;
5889			} else {
5890			PDL_Indx lodex,hidex;
5891			PDL_Double beta = fabs(zeta * sum);
5892
5893			lodex = beta;
5894			beta -= lodex; hidex = lodex+1;
5895			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
5896
5897			}
5898			}
5899			break;
5900
5901			case 'G':
5902			/* This is the Gaussian rolloff with explicit calculation, preserved */
5903			/* in case someone actually wants the slower longer method. */
5904			{
5905			PDL_Double sum = 0;
5906			cp = tmp;
5907			for(i=0; i
5908			dd = 0;
5909			ap = tvec;
5910			for(j=0;j
5911			dd += (ap++) *(cp++);
5912			dd /= blur;
5913			sum += dd * dd;
5914			if(sum > 4) /* 2 pixels -- four half-widths */
5915			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
5916			}
5917
5918			if(sum > GAUSSIAN_MAXVAL)
5919			alpha = 0;
5920			else
5921			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
5922			}
5923			break;
5924			default:
5925			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
5926			}
5927
5928			{ /* convenience block -- accumulate the current point into the weighted sum. */
5929			/* This is more than simple assignment because we have our own explicit poor */
5930			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
5931			PDL_Float dat = ((PDL_Float )(in->data)) + i_off;
5932			PDL_Indx max = out->dims[ndims];
5933			for( i=0; i < max; i++ ) {
5934			if( (badval==0) \|\| (*dat != badval) ) {
5935			acc[i] += dat alpha;
5936			dat += in->dimincs[ndims];
5937			wgt[i] += alpha;
5938			}
5939			wgt2[i] += alpha; }
5940			}
5941			} /* end of t_vio check (i.e. of input accumulation) */
5942
5943
5944			/* Advance input accumulation loop. */
5945			/* We both increment the total vector and also advance the index. */
5946			carry = 1;
5947			for (i=0; i
5948			/* Advance the current element of the offset vector */
5949			ivec[i]++;
5950			j = ivec[i] + ibvec[i];
5951
5952			/* Advance the offset into the data array */
5953			if( j > 0 && j <= in->dims[i]-1 ) {
5954			/* Normal case -- just advance the input vector */
5955			i_off += in->dimincs[i];
5956			} else {
5957			/* Busted a boundary - either before or after. */
5958			switch(bounds[i]){
5959			case 0: /* no breakage allowed -- treat as truncation for interpolation */
5960			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
5961			if( j == 0 )
5962			t_vio--;
5963			else if( j == in->dims[i] )
5964			t_vio++;
5965			break;
5966			case 2: /* extension -- do nothing (so the same input point is re-used) */
5967			break;
5968			case 3: /* periodic -- advance and mod into the allowed range */
5969			if((j % in->dims[i]) == 0) {
5970			i_off -= in->dimincs[i] * (in->dims[i]-1);
5971			} else {
5972			i_off += in->dimincs[i];
5973			}
5974			break;
5975			case 4: /* mirror -- advance or retreat depending on phase */
5976			j += in->dims[i];
5977			j %= (in->dims[i]*2);
5978			j -= in->dims[i];
5979			if( j!=0 && j!= -in->dims[i] ) {
5980			if(j<0)
5981			i_off -= in->dimincs[i];
5982			else
5983			i_off += in->dimincs[i];
5984			}
5985			break;
5986			}
5987			}
5988
5989			/* Now check for carry */
5990			if (ivec[i] <= psize) {
5991			/* Normal case -- copy the current offset to the faster-running dim stashes */
5992			PDL_Indx k;
5993			for(k=0;k
5994			index_stash[k] = i_off;
5995			}
5996			carry = 0;
5997
5998			} else { /* End of this scan -- recover the last position, and mark carry */
5999			i_off = index_stash[i];
6000			if(bounds[i]==1) {
6001			j = ivec[i] + ibvec[i];
6002			if( j < 0 \|\| j >= in->dims[i] )
6003			t_vio--;
6004			ivec[i] = -psize;
6005			j = ivec[i] + ibvec[i];
6006			if( j < 0 \|\| j >= in->dims[i] )
6007			t_vio++;
6008			carry = 1;
6009			} else {
6010			ivec[i] = -psize;
6011			}
6012			}
6013			} /* End of counter-advance loop */
6014			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
6015
6016			{
6017			PDL_Double *ac = acc;
6018			PDL_Double *wg = wgt;
6019			PDL_Double *wg2 = wgt2;
6020			PDL_Float *dat = out->data;
6021
6022			/* Calculate output vector offset */
6023			for(i=0;i
6024			dat += out->dimincs[i] * ovec[i];
6025
6026			if (!flux) {
6027			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
6028			for (i=0; i < out->dims[ndims]; i++) {
6029			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
6030			ac++; wg++; wg2++;
6031			dat += out->dimincs[ndims];
6032			}
6033			} else {
6034			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
6035			PDL_Double det = tmp[ndims*ndims];
6036			for(i=0; i < out->dims[ndims]; i++) {
6037			if(wg && (wg2 / *wg) < 1.5 ) {
6038			dat = (ac++) / (wg++) det;
6039			wg2++;
6040			} else {
6041			*dat = badval;
6042			ac++; wg++; wg2++;
6043			}
6044			dat += out->dimincs[ndims];
6045			} /* end of for loop */
6046			} /* end of flux flag set conditional */
6047			} /* end of convenience block */
6048
6049			/* End of code for normal pixels */
6050			} else {
6051			/* The pixel was ludicrously huge -- just set this pixel to nan */
6052			PDL_Float *dat = out->data;
6053			for(i=0;i
6054			dat += out->dimincs[i] * ovec[i];
6055			for(i=0;idims[ndims];i++) {
6056			dat = badval; / Should handle bad values too -- not yet */
6057			dat += out->dimincs[ndims];
6058			}
6059			}
6060
6061			/* Increment the pixel counter */
6062			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
6063			} while (1);
6064			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
6065			#line 6066 "lib/PDL/Transform-pp-map.c"
6066	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
6067	0		} break;
6068	7		case PDL_D: {
6069	7	50	PDL_DECLARE_PARAMS_map_1(PDL_Double,D)
		50
		50
6070	21	50	PDL_BROADCASTLOOP_START_map_readdata {
		50
		50
		50
		50
		100
		100
6071			#line 611 "lib/PDL/Transform.pd"
6072
6073			/*
6074			* Pixel interpolation & averaging code
6075			*
6076			* Calls a common coordinate-transformation block (see following hdr)
6077			* that isn't dependent on the type of the input variable.
6078			*
6079			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
6080			* is handled internally. To simplify the broadcasting business, any
6081			* broadcast dimensions should all be collapsed to a single one by the
6082			* perl front-end.
6083			*
6084			*/
6085
6086			pdl *in = __params->in;
6087			pdl *out = __params->out;
6088			pdl *map = __params->map;
6089
6090			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
6091			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
6092			PDL_Indx ovec[ndims]; /* output pixel loop vector */
6093			PDL_Indx ivec[ndims]; /* input pixel loop vector */
6094			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
6095			PDL_Double dvec[ndims]; /* Residual vector for linearization */
6096			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
6097			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
6098			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
6099			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
6100			char bounds[ndims]; /* Boundary condition packed string */
6101			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
6102			char method; /* Method identifier (gets one of 'h','g') */
6103			PDL_Long big; /* Max size of input footprint for each pix */
6104			PDL_Double blur; /* Scaling of filter */
6105			PDL_Double sv_min; /* minimum singular value */
6106			char flux; /* Flag to indicate flux conservation */
6107			PDL_Indx i;
6108			PDL_Double badval = SvNV(__params->bv);
6109			#define HANNING_LOOKUP_SIZE 2500
6110			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
6111			static int needs_hanning_calc = 1;
6112			PDL_Double zeta = 0;
6113			PDL_Double hanning_offset;
6114
6115			#define GAUSSIAN_LOOKUP_SIZE 4000
6116			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
6117			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
6118			static int needs_gaussian_calc = 1;
6119
6120			PDL_RETERROR(PDL_err, PDL->make_physical(in));
6121			PDL_RETERROR(PDL_err, PDL->make_physical(out));
6122			PDL_RETERROR(PDL_err, PDL->make_physical(map));
6123
6124			/***
6125			* Fill in the boundary condition array
6126			*/
6127			{
6128			char *bstr;
6129			STRLEN blen;
6130			bstr = SvPV(__params->boundary,blen);
6131
6132			if(blen == 0) {
6133			/* If no boundary is specified then every dim gets truncated */
6134			PDL_Indx i;
6135			for (i=0;i
6136			bounds[i] = 1;
6137			} else {
6138			PDL_Indx i;
6139			for(i=0;i
6140			switch(bstr[i < blen ? i : blen-1 ]) {
6141			case '0': case 'f': case 'F': /* forbid */
6142			bounds[i] = 0;
6143			break;
6144			case '1': case 't': case 'T': /* truncate */
6145			bounds[i] = 1;
6146			break;
6147			case '2': case 'e': case 'E': /* extend */
6148			bounds[i] = 2;
6149			break;
6150			case '3': case 'p': case 'P': /* periodic */
6151			bounds[i] = 3;
6152			break;
6153			case '4': case 'm': case 'M': /* mirror */
6154			bounds[i] = 4;
6155			break;
6156			default:
6157			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
6158			break;
6159			}
6160			}
6161			}
6162			}
6163
6164			/***
6165			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
6166			*/
6167			big = labs(__params->big);
6168			if(big <= 0)
6169			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
6170
6171			blur = fabs(__params->blur);
6172			if(blur < 0)
6173			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
6174
6175			sv_min = fabs(__params->sv_min);
6176			if(sv_min < 0)
6177			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
6178
6179			flux = __params->flux != 0;
6180
6181			{
6182			char *mstr;
6183			STRLEN mlen;
6184			mstr = SvPV(__params->method,mlen);
6185
6186			if(mlen==0)
6187			method = 'h';
6188			switch(*mstr) {
6189			case 'h': if( needs_hanning_calc ) {
6190			int i;
6191			for(i=0;i
6192			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
6193			}
6194			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
6195			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
6196			needs_hanning_calc = 0;
6197			}
6198			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
6199			case 'H': method = *mstr;
6200			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
6201			break;
6202
6203			case 'g': case 'j': method = 'g';
6204			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
6205
6206			if( needs_gaussian_calc ) {
6207			int i;
6208			for(i=0;i
6209			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
6210			}
6211			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
6212			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
6213			needs_gaussian_calc = 0;
6214			}
6215			break;
6216
6217			case 'G': case 'J': method = 'G'; break;
6218			default:
6219			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
6220			break;
6221			}
6222			}
6223
6224
6225
6226			/* End of initialization */
6227			/*************************************************************/
6228			/* Start of Real Work */
6229
6230			/* Initialize coordinate vector and map offset
6231			*/
6232			for(i=0;i
6233			ovec[i] = 0;
6234
6235			PDL_Double map_ptr = (PDL_Double )(map->data);
6236			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
6237			for (i=0; i < npdls; i++) offs[i] = 0;
6238			for (i=0; i < ndims; i++) {
6239			incs[i*npdls + 0] = map->dimincs[i+1];
6240			}
6241
6242			/* Main pixel loop (iterates over pixels in the output plane) */
6243			do {
6244			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
6245			/* Prefrobnicate the transformation matrix */
6246			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
6247
6248			#ifdef DEBUG_MAP
6249			{
6250			PDL_Indx k; PDL_Indx foo = 0;
6251			printf("ovec: [");
6252			for(k=0;k
6253			foo += ovec[k] * map->dimincs[k+1];
6254			printf(" %2td ",(PDL_Indx)(ovec[k]));
6255			}
6256			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
6257			for(k=0;k
6258			printf("%8.2f",(double)(((PDL_Double )(map->data))[foo + kmap->dimincs[0]]));
6259			}
6260			printf("]\n");
6261			}
6262			#endif
6263
6264			/* Don't bother accumulating output if psize is too large */
6265			if(psize <= big) {
6266			/* Use the prefrobnicated matrix to generate a local linearization.
6267			* dvec gets the delta; ibvec gets the base.
6268			*/
6269			{
6270			PDL_Double *mp = map_ptr + offs[0];
6271			for (i=0;i
6272			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
6273			mp += map->dimincs[0];
6274			}
6275			}
6276
6277			/* Initialize input delta vector */
6278			for(i=0;i
6279			ivec[i] = -psize;
6280
6281			/* Initialize accumulators */
6282			{
6283			PDL_Double *ac = acc;
6284			for(i=0; i < in->dims[ndims]; i++)
6285			*(ac++) = 0.0;
6286
6287			}
6288			{
6289			PDL_Double *wg = wgt;
6290			for(i=0;i < in->dims[ndims]; i++)
6291			*(wg++) = 0.0;
6292			}
6293			{
6294			PDL_Double *wg = wgt2;
6295			for(i=0;i < in->dims[ndims]; i++)
6296			*(wg++) = 0.0;
6297			}
6298
6299
6300			/*
6301			* Calculate the original offset into the data array, to enable
6302			* delta calculations in the pixel loop
6303			*
6304			* i runs over dims; j holds the working integer index in the
6305			* current dim.
6306			*
6307			* This code matches the incrementation code at the bottom of the accumulation loop
6308			*/
6309
6310			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
6311			i_off = 0;
6312			for(i=0;i
6313			j = ivec[i] + ibvec[i];
6314			if(j<0 \|\| j >= in->dims[i]) {
6315			switch(bounds[i]) {
6316			case 0: /* no breakage allowed */
6317			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
6318			break;
6319			case 1: /* truncation */
6320			t_vio++;
6321			/* fall through */
6322			case 2: /* extension -- crop */
6323			if(j<0)
6324			j=0;
6325			else j = in->dims[i] - 1;
6326			break;
6327			case 3: /* periodic -- mod it */
6328			j %= in->dims[i];
6329			if(j<0)
6330			j += in->dims[i];
6331			break;
6332			case 4: /* mirror -- reflect off the edges */
6333			j += in->dims[i];
6334			j %= (in->dims[i]*2);
6335			if(j<0)
6336			j += in->dims[i]*2;
6337			j -= in->dims[i];
6338			if(j<0) {
6339			j *= -1;
6340			j -= 1;
6341			}
6342			break;
6343			default:
6344			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
6345			break;
6346			}
6347			}
6348			i_off += in->dimincs[i] * j;
6349			}
6350
6351			/* Initialize index stashes for later reference as we scan the footprint */
6352			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
6353			/* end of a dimensional scan. So we stash the index location at the */
6354			/* start of each dimensional scan here. When we finish incrementing */
6355			/* through a particular dim, we pull its value back out of the stash. */
6356			for(i=0;i
6357			index_stash[i] = i_off;
6358			}
6359
6360			/* The input accumulation loop is the hotspot for the entire operation. */
6361			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
6362			/* in the input array, use the linearized transform to bring each pixel center */
6363			/* forward to the output plane, and calculate a weighting based on the chosen */
6364			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
6365			/* table that is initialized the first time through the code. 'H' is the */
6366			/* same process, but explicitly calculated for each interation (~2x slower). */
6367			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
6368			/* into the input array fresh from the current input array vector each time, */
6369			/* we walk through the array using dimincs and the old offset. This saves */
6370			/* about half of the time spent on index calculation. */
6371
6372			do { /* Input accumulation loop */
6373			PDL_Double *cp;
6374			PDL_Double alpha;
6375			/* Calculate the weight of the current input point. Don't bother if we're
6376			* violating any truncation boundaries (in that case our value is zero, but
6377			* for the interpolation we also set the weight to zero).
6378			*/
6379			if( !t_vio ) {
6380
6381			PDL_Double *ap = tvec;
6382			PDL_Double *bp = dvec;
6383			PDL_Indx *ip = ivec;
6384			for(i=0; i
6385			(ap++) = (ip++) - *(bp++);
6386
6387			switch(method) {
6388			PDL_Double dd;
6389			case 'h':
6390			/* This is the Hanning window rolloff. It is a product of a simple */
6391			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
6392			/* is about 2x faster than using cos(theta) directly in each */
6393			/* weighting calculation, so we do. Using 2500 entries and linear */
6394			/* interpolation is accurate to about 10^-7, and should preserve */
6395			/* the contents of cache pretty well. */
6396			alpha = 1;
6397			cp = tmp;
6398			for (i=0; i
6399			PDL_Indx lodex;
6400			PDL_Indx hidex;
6401			PDL_Double beta;
6402			dd = 0;
6403			ap = tvec;
6404			/* Get the matrix-multiply element for this dimension */
6405			for (j=0;j
6406			dd += (ap++) *(cp++);
6407
6408			/* Do linear interpolation from the table */
6409			/* The table captures a hanning window centered 0.5 pixel from center. */
6410			/* We scale the filter by the blur parameter -- but if blur is less */
6411			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
6412			/* value on the pixel edge at 0.5. For blur greater than unity, we */
6413			/* scale simply. */
6414			beta = fabs(dd) - hanning_offset;
6415			if (beta > 0) {
6416			if (beta >= blur) {
6417			alpha = 0;
6418			i = ndims;
6419			} else {
6420			beta *= zeta;
6421			lodex = beta;
6422			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
6423			lodex = HANNING_LOOKUP_SIZE;
6424			hidex = lodex+1;
6425			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
6426			} /* end of interpolation branch */
6427			} /* end of beta > 0 branch */
6428			} /* end of dimension loop */
6429			break;
6430
6431			case 'H':
6432			/* This is the Hanning window rolloff with explicit calculation, preserved */
6433			/* in case someone actually wants the slower longer method. */
6434			alpha = 1;
6435			cp = tmp;
6436			for (i=0; i
6437			dd = 0;
6438			ap = tvec;
6439			for (j=0;j
6440			dd += (ap++) *(cp++);
6441			dd = (fabs(dd) - hanning_offset) / blur;
6442			if ( dd > 1 ) {
6443			alpha = 0;
6444			i = ndims;
6445			} else
6446			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
6447			}
6448			break;
6449
6450			case 'g':
6451			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
6452			{
6453			PDL_Double sum = 0;
6454			cp = tmp;
6455			for(i=0; i
6456			dd = 0;
6457			ap = tvec;
6458			for(j=0;j
6459			dd += (ap++) *(cp++);
6460			dd /= blur;
6461			sum += dd * dd;
6462			if(sum > GAUSSIAN_MAXVAL) {
6463			i = ndims; /* exit early if we're too far out */
6464			alpha = 0;
6465			}
6466			}
6467			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
6468			alpha = 0;
6469			} else {
6470			PDL_Indx lodex,hidex;
6471			PDL_Double beta = fabs(zeta * sum);
6472
6473			lodex = beta;
6474			beta -= lodex; hidex = lodex+1;
6475			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
6476
6477			}
6478			}
6479			break;
6480
6481			case 'G':
6482			/* This is the Gaussian rolloff with explicit calculation, preserved */
6483			/* in case someone actually wants the slower longer method. */
6484			{
6485			PDL_Double sum = 0;
6486			cp = tmp;
6487			for(i=0; i
6488			dd = 0;
6489			ap = tvec;
6490			for(j=0;j
6491			dd += (ap++) *(cp++);
6492			dd /= blur;
6493			sum += dd * dd;
6494			if(sum > 4) /* 2 pixels -- four half-widths */
6495			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
6496			}
6497
6498			if(sum > GAUSSIAN_MAXVAL)
6499			alpha = 0;
6500			else
6501			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
6502			}
6503			break;
6504			default:
6505			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
6506			}
6507
6508			{ /* convenience block -- accumulate the current point into the weighted sum. */
6509			/* This is more than simple assignment because we have our own explicit poor */
6510			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
6511			PDL_Double dat = ((PDL_Double )(in->data)) + i_off;
6512			PDL_Indx max = out->dims[ndims];
6513			for( i=0; i < max; i++ ) {
6514			if( (badval==0) \|\| (*dat != badval) ) {
6515			acc[i] += dat alpha;
6516			dat += in->dimincs[ndims];
6517			wgt[i] += alpha;
6518			}
6519			wgt2[i] += alpha; }
6520			}
6521			} /* end of t_vio check (i.e. of input accumulation) */
6522
6523
6524			/* Advance input accumulation loop. */
6525			/* We both increment the total vector and also advance the index. */
6526			carry = 1;
6527			for (i=0; i
6528			/* Advance the current element of the offset vector */
6529			ivec[i]++;
6530			j = ivec[i] + ibvec[i];
6531
6532			/* Advance the offset into the data array */
6533			if( j > 0 && j <= in->dims[i]-1 ) {
6534			/* Normal case -- just advance the input vector */
6535			i_off += in->dimincs[i];
6536			} else {
6537			/* Busted a boundary - either before or after. */
6538			switch(bounds[i]){
6539			case 0: /* no breakage allowed -- treat as truncation for interpolation */
6540			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
6541			if( j == 0 )
6542			t_vio--;
6543			else if( j == in->dims[i] )
6544			t_vio++;
6545			break;
6546			case 2: /* extension -- do nothing (so the same input point is re-used) */
6547			break;
6548			case 3: /* periodic -- advance and mod into the allowed range */
6549			if((j % in->dims[i]) == 0) {
6550			i_off -= in->dimincs[i] * (in->dims[i]-1);
6551			} else {
6552			i_off += in->dimincs[i];
6553			}
6554			break;
6555			case 4: /* mirror -- advance or retreat depending on phase */
6556			j += in->dims[i];
6557			j %= (in->dims[i]*2);
6558			j -= in->dims[i];
6559			if( j!=0 && j!= -in->dims[i] ) {
6560			if(j<0)
6561			i_off -= in->dimincs[i];
6562			else
6563			i_off += in->dimincs[i];
6564			}
6565			break;
6566			}
6567			}
6568
6569			/* Now check for carry */
6570			if (ivec[i] <= psize) {
6571			/* Normal case -- copy the current offset to the faster-running dim stashes */
6572			PDL_Indx k;
6573			for(k=0;k
6574			index_stash[k] = i_off;
6575			}
6576			carry = 0;
6577
6578			} else { /* End of this scan -- recover the last position, and mark carry */
6579			i_off = index_stash[i];
6580			if(bounds[i]==1) {
6581			j = ivec[i] + ibvec[i];
6582			if( j < 0 \|\| j >= in->dims[i] )
6583			t_vio--;
6584			ivec[i] = -psize;
6585			j = ivec[i] + ibvec[i];
6586			if( j < 0 \|\| j >= in->dims[i] )
6587			t_vio++;
6588			carry = 1;
6589			} else {
6590			ivec[i] = -psize;
6591			}
6592			}
6593			} /* End of counter-advance loop */
6594			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
6595
6596			{
6597			PDL_Double *ac = acc;
6598			PDL_Double *wg = wgt;
6599			PDL_Double *wg2 = wgt2;
6600			PDL_Double *dat = out->data;
6601
6602			/* Calculate output vector offset */
6603			for(i=0;i
6604			dat += out->dimincs[i] * ovec[i];
6605
6606			if (!flux) {
6607			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
6608			for (i=0; i < out->dims[ndims]; i++) {
6609			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
6610			ac++; wg++; wg2++;
6611			dat += out->dimincs[ndims];
6612			}
6613			} else {
6614			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
6615			PDL_Double det = tmp[ndims*ndims];
6616			for(i=0; i < out->dims[ndims]; i++) {
6617			if(wg && (wg2 / *wg) < 1.5 ) {
6618			dat = (ac++) / (wg++) det;
6619			wg2++;
6620			} else {
6621			*dat = badval;
6622			ac++; wg++; wg2++;
6623			}
6624			dat += out->dimincs[ndims];
6625			} /* end of for loop */
6626			} /* end of flux flag set conditional */
6627			} /* end of convenience block */
6628
6629			/* End of code for normal pixels */
6630			} else {
6631			/* The pixel was ludicrously huge -- just set this pixel to nan */
6632			PDL_Double *dat = out->data;
6633			for(i=0;i
6634			dat += out->dimincs[i] * ovec[i];
6635			for(i=0;idims[ndims];i++) {
6636			dat = badval; / Should handle bad values too -- not yet */
6637			dat += out->dimincs[ndims];
6638			}
6639			}
6640
6641			/* Increment the pixel counter */
6642			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
6643			} while (1);
6644			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
6645			#line 6646 "lib/PDL/Transform-pp-map.c"
6646	7	50	}PDL_BROADCASTLOOP_END_map_readdata
		50
6647	7		} break;
6648	0		case PDL_LD: {
6649	0	0	PDL_DECLARE_PARAMS_map_1(PDL_LDouble,E)
		0
		0
6650	0	0	PDL_BROADCASTLOOP_START_map_readdata {
		0
		0
		0
		0
		0
		0
6651			#line 611 "lib/PDL/Transform.pd"
6652
6653			/*
6654			* Pixel interpolation & averaging code
6655			*
6656			* Calls a common coordinate-transformation block (see following hdr)
6657			* that isn't dependent on the type of the input variable.
6658			*
6659			* The inputs are SVs to avoid hassling with broadcastloops; broadcasting
6660			* is handled internally. To simplify the broadcasting business, any
6661			* broadcast dimensions should all be collapsed to a single one by the
6662			* perl front-end.
6663			*
6664			*/
6665
6666			pdl *in = __params->in;
6667			pdl *out = __params->out;
6668			pdl *map = __params->map;
6669
6670			PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
6671			PDL_Double tmp[3ndimsndims+ndims]; /* Workspace for prefrobnication */
6672			PDL_Indx ovec[ndims]; /* output pixel loop vector */
6673			PDL_Indx ivec[ndims]; /* input pixel loop vector */
6674			PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
6675			PDL_Double dvec[ndims]; /* Residual vector for linearization */
6676			PDL_Double tvec[ndims]; /* Temporary floating-point vector */
6677			PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
6678			PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
6679			PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval finding */
6680			char bounds[ndims]; /* Boundary condition packed string */
6681			PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
6682			char method; /* Method identifier (gets one of 'h','g') */
6683			PDL_Long big; /* Max size of input footprint for each pix */
6684			PDL_Double blur; /* Scaling of filter */
6685			PDL_Double sv_min; /* minimum singular value */
6686			char flux; /* Flag to indicate flux conservation */
6687			PDL_Indx i;
6688			PDL_LDouble badval = SvNV(__params->bv);
6689			#define HANNING_LOOKUP_SIZE 2500
6690			static PDL_Double hanning_lookup[HANNING_LOOKUP_SIZE + 2];
6691			static int needs_hanning_calc = 1;
6692			PDL_Double zeta = 0;
6693			PDL_Double hanning_offset;
6694
6695			#define GAUSSIAN_LOOKUP_SIZE 4000
6696			#define GAUSSIAN_MAXVAL 6.25 /* 2.5 HWHMs (square it) */
6697			static PDL_Double gaussian_lookup[GAUSSIAN_LOOKUP_SIZE + 2];
6698			static int needs_gaussian_calc = 1;
6699
6700			PDL_RETERROR(PDL_err, PDL->make_physical(in));
6701			PDL_RETERROR(PDL_err, PDL->make_physical(out));
6702			PDL_RETERROR(PDL_err, PDL->make_physical(map));
6703
6704			/***
6705			* Fill in the boundary condition array
6706			*/
6707			{
6708			char *bstr;
6709			STRLEN blen;
6710			bstr = SvPV(__params->boundary,blen);
6711
6712			if(blen == 0) {
6713			/* If no boundary is specified then every dim gets truncated */
6714			PDL_Indx i;
6715			for (i=0;i
6716			bounds[i] = 1;
6717			} else {
6718			PDL_Indx i;
6719			for(i=0;i
6720			switch(bstr[i < blen ? i : blen-1 ]) {
6721			case '0': case 'f': case 'F': /* forbid */
6722			bounds[i] = 0;
6723			break;
6724			case '1': case 't': case 'T': /* truncate */
6725			bounds[i] = 1;
6726			break;
6727			case '2': case 'e': case 'E': /* extend */
6728			bounds[i] = 2;
6729			break;
6730			case '3': case 'p': case 'P': /* periodic */
6731			bounds[i] = 3;
6732			break;
6733			case '4': case 'm': case 'M': /* mirror */
6734			bounds[i] = 4;
6735			break;
6736			default:
6737			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Error in map: Unknown boundary condition '%c'",bstr[i]);
6738			break;
6739			}
6740			}
6741			}
6742			}
6743
6744			/***
6745			* Parse out the 'method', 'big', 'blur', and 'sv_min' arguments
6746			*/
6747			big = labs(__params->big);
6748			if(big <= 0)
6749			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'big' parameter must be >0");
6750
6751			blur = fabs(__params->blur);
6752			if(blur < 0)
6753			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'blur' parameter must be >= 0");
6754
6755			sv_min = fabs(__params->sv_min);
6756			if(sv_min < 0)
6757			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "map: 'sv_min' parameter must be >= 0");
6758
6759			flux = __params->flux != 0;
6760
6761			{
6762			char *mstr;
6763			STRLEN mlen;
6764			mstr = SvPV(__params->method,mlen);
6765
6766			if(mlen==0)
6767			method = 'h';
6768			switch(*mstr) {
6769			case 'h': if( needs_hanning_calc ) {
6770			int i;
6771			for(i=0;i
6772			hanning_lookup[i] = 0.5 + 0.5 * cos(3.1415926536 / HANNING_LOOKUP_SIZE * i);
6773			}
6774			hanning_lookup[HANNING_LOOKUP_SIZE] = 0;
6775			hanning_lookup[HANNING_LOOKUP_SIZE+1] = 0;
6776			needs_hanning_calc = 0;
6777			}
6778			zeta = HANNING_LOOKUP_SIZE / blur; /* fall-through deliberate */
6779			case 'H': method = *mstr;
6780			hanning_offset = (blur >= 1) ? 0 : 0.5 * (1.0 - blur);
6781			break;
6782
6783			case 'g': case 'j': method = 'g';
6784			zeta = GAUSSIAN_LOOKUP_SIZE / GAUSSIAN_MAXVAL;
6785
6786			if( needs_gaussian_calc ) {
6787			int i;
6788			for(i=0;i
6789			gaussian_lookup[i] = exp( - i * 1.386294 / zeta );
6790			}
6791			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE] = 0;
6792			gaussian_lookup[GAUSSIAN_LOOKUP_SIZE+1] = 0;
6793			needs_gaussian_calc = 0;
6794			}
6795			break;
6796
6797			case 'G': case 'J': method = 'G'; break;
6798			default:
6799			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Bug in map: unknown method '%c'",*mstr);
6800			break;
6801			}
6802			}
6803
6804
6805
6806			/* End of initialization */
6807			/*************************************************************/
6808			/* Start of Real Work */
6809
6810			/* Initialize coordinate vector and map offset
6811			*/
6812			for(i=0;i
6813			ovec[i] = 0;
6814
6815			PDL_Double map_ptr = (PDL_Double )(map->data);
6816			PDL_Indx npdls = 1, incs[npdls*ndims], offs[npdls];
6817			for (i=0; i < npdls; i++) offs[i] = 0;
6818			for (i=0; i < ndims; i++) {
6819			incs[i*npdls + 0] = map->dimincs[i+1];
6820			}
6821
6822			/* Main pixel loop (iterates over pixels in the output plane) */
6823			do {
6824			PDL_Indx psize; PDL_Indx i_off; PDL_Indx j; char t_vio; char carry;
6825			/* Prefrobnicate the transformation matrix */
6826			psize = (PDL_Long)(blur * PDL_xform_aux(map, ovec, tmp, sv_min) + 0.5)+1; /* assignment */
6827
6828			#ifdef DEBUG_MAP
6829			{
6830			PDL_Indx k; PDL_Indx foo = 0;
6831			printf("ovec: [");
6832			for(k=0;k
6833			foo += ovec[k] * map->dimincs[k+1];
6834			printf(" %2td ",(PDL_Indx)(ovec[k]));
6835			}
6836			printf("]; psize is %ld; big is %d; blur is %8.2f; map is [",psize,big, blur);
6837			for(k=0;k
6838			printf("%8.2f",(double)(((PDL_LDouble )(map->data))[foo + kmap->dimincs[0]]));
6839			}
6840			printf("]\n");
6841			}
6842			#endif
6843
6844			/* Don't bother accumulating output if psize is too large */
6845			if(psize <= big) {
6846			/* Use the prefrobnicated matrix to generate a local linearization.
6847			* dvec gets the delta; ibvec gets the base.
6848			*/
6849			{
6850			PDL_Double *mp = map_ptr + offs[0];
6851			for (i=0;i
6852			dvec[i] = mp - ( ibvec[i] = (PDL_Long)(mp + 0.5)); /* assignment */
6853			mp += map->dimincs[0];
6854			}
6855			}
6856
6857			/* Initialize input delta vector */
6858			for(i=0;i
6859			ivec[i] = -psize;
6860
6861			/* Initialize accumulators */
6862			{
6863			PDL_Double *ac = acc;
6864			for(i=0; i < in->dims[ndims]; i++)
6865			*(ac++) = 0.0;
6866
6867			}
6868			{
6869			PDL_Double *wg = wgt;
6870			for(i=0;i < in->dims[ndims]; i++)
6871			*(wg++) = 0.0;
6872			}
6873			{
6874			PDL_Double *wg = wgt2;
6875			for(i=0;i < in->dims[ndims]; i++)
6876			*(wg++) = 0.0;
6877			}
6878
6879
6880			/*
6881			* Calculate the original offset into the data array, to enable
6882			* delta calculations in the pixel loop
6883			*
6884			* i runs over dims; j holds the working integer index in the
6885			* current dim.
6886			*
6887			* This code matches the incrementation code at the bottom of the accumulation loop
6888			*/
6889
6890			t_vio = 0; /* truncation-boundary violation count - don't bother if it is nonzero */
6891			i_off = 0;
6892			for(i=0;i
6893			j = ivec[i] + ibvec[i];
6894			if(j<0 \|\| j >= in->dims[i]) {
6895			switch(bounds[i]) {
6896			case 0: /* no breakage allowed */
6897			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "index out-of-bounds in map");
6898			break;
6899			case 1: /* truncation */
6900			t_vio++;
6901			/* fall through */
6902			case 2: /* extension -- crop */
6903			if(j<0)
6904			j=0;
6905			else j = in->dims[i] - 1;
6906			break;
6907			case 3: /* periodic -- mod it */
6908			j %= in->dims[i];
6909			if(j<0)
6910			j += in->dims[i];
6911			break;
6912			case 4: /* mirror -- reflect off the edges */
6913			j += in->dims[i];
6914			j %= (in->dims[i]*2);
6915			if(j<0)
6916			j += in->dims[i]*2;
6917			j -= in->dims[i];
6918			if(j<0) {
6919			j *= -1;
6920			j -= 1;
6921			}
6922			break;
6923			default:
6924			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "Unknown boundary condition %td in map -- bug alert!",(PDL_Indx)bounds[i]);
6925			break;
6926			}
6927			}
6928			i_off += in->dimincs[i] * j;
6929			}
6930
6931			/* Initialize index stashes for later reference as we scan the footprint */
6932			/* It's a pain in the ass to deal with boundaries, and doubly so at the */
6933			/* end of a dimensional scan. So we stash the index location at the */
6934			/* start of each dimensional scan here. When we finish incrementing */
6935			/* through a particular dim, we pull its value back out of the stash. */
6936			for(i=0;i
6937			index_stash[i] = i_off;
6938			}
6939
6940			/* The input accumulation loop is the hotspot for the entire operation. */
6941			/* We loop over pixels in the region of interest (+/- psize in each dimension) */
6942			/* in the input array, use the linearized transform to bring each pixel center */
6943			/* forward to the output plane, and calculate a weighting based on the chosen */
6944			/* filter function. 'h' is a fast Hanning window rolloff using a lookup */
6945			/* table that is initialized the first time through the code. 'H' is the */
6946			/* same process, but explicitly calculated for each interation (~2x slower). */
6947			/* 'g' uses a radial Gaussian filter. Rather than calculate the array offset */
6948			/* into the input array fresh from the current input array vector each time, */
6949			/* we walk through the array using dimincs and the old offset. This saves */
6950			/* about half of the time spent on index calculation. */
6951
6952			do { /* Input accumulation loop */
6953			PDL_Double *cp;
6954			PDL_Double alpha;
6955			/* Calculate the weight of the current input point. Don't bother if we're
6956			* violating any truncation boundaries (in that case our value is zero, but
6957			* for the interpolation we also set the weight to zero).
6958			*/
6959			if( !t_vio ) {
6960
6961			PDL_Double *ap = tvec;
6962			PDL_Double *bp = dvec;
6963			PDL_Indx *ip = ivec;
6964			for(i=0; i
6965			(ap++) = (ip++) - *(bp++);
6966
6967			switch(method) {
6968			PDL_Double dd;
6969			case 'h':
6970			/* This is the Hanning window rolloff. It is a product of a simple */
6971			/* cos^2(theta) rolloff in each dimension. Using a lookup table */
6972			/* is about 2x faster than using cos(theta) directly in each */
6973			/* weighting calculation, so we do. Using 2500 entries and linear */
6974			/* interpolation is accurate to about 10^-7, and should preserve */
6975			/* the contents of cache pretty well. */
6976			alpha = 1;
6977			cp = tmp;
6978			for (i=0; i
6979			PDL_Indx lodex;
6980			PDL_Indx hidex;
6981			PDL_Double beta;
6982			dd = 0;
6983			ap = tvec;
6984			/* Get the matrix-multiply element for this dimension */
6985			for (j=0;j
6986			dd += (ap++) *(cp++);
6987
6988			/* Do linear interpolation from the table */
6989			/* The table captures a hanning window centered 0.5 pixel from center. */
6990			/* We scale the filter by the blur parameter -- but if blur is less */
6991			/* than unity, we shrink the hanning blur window while keeping the 0.5 */
6992			/* value on the pixel edge at 0.5. For blur greater than unity, we */
6993			/* scale simply. */
6994			beta = fabs(dd) - hanning_offset;
6995			if (beta > 0) {
6996			if (beta >= blur) {
6997			alpha = 0;
6998			i = ndims;
6999			} else {
7000			beta *= zeta;
7001			lodex = beta;
7002			beta -= lodex; if (lodex > HANNING_LOOKUP_SIZE)
7003			lodex = HANNING_LOOKUP_SIZE;
7004			hidex = lodex+1;
7005			alpha = hanning_lookup[hidex]beta + hanning_lookup[lodex]*(1-beta);
7006			} /* end of interpolation branch */
7007			} /* end of beta > 0 branch */
7008			} /* end of dimension loop */
7009			break;
7010
7011			case 'H':
7012			/* This is the Hanning window rolloff with explicit calculation, preserved */
7013			/* in case someone actually wants the slower longer method. */
7014			alpha = 1;
7015			cp = tmp;
7016			for (i=0; i
7017			dd = 0;
7018			ap = tvec;
7019			for (j=0;j
7020			dd += (ap++) *(cp++);
7021			dd = (fabs(dd) - hanning_offset) / blur;
7022			if ( dd > 1 ) {
7023			alpha = 0;
7024			i = ndims;
7025			} else
7026			alpha = (0.5 + 0.5 cos( dd * 3.1415926536 ));
7027			}
7028			break;
7029
7030			case 'g':
7031			/* This is the Gaussian rolloff. It does lookup into a precalculated exponential. */
7032			{
7033			PDL_Double sum = 0;
7034			cp = tmp;
7035			for(i=0; i
7036			dd = 0;
7037			ap = tvec;
7038			for(j=0;j
7039			dd += (ap++) *(cp++);
7040			dd /= blur;
7041			sum += dd * dd;
7042			if(sum > GAUSSIAN_MAXVAL) {
7043			i = ndims; /* exit early if we're too far out */
7044			alpha = 0;
7045			}
7046			}
7047			if( sum > GAUSSIAN_MAXVAL \|\| !isfinite(sum) \|\| isnan(sum) ) {
7048			alpha = 0;
7049			} else {
7050			PDL_Indx lodex,hidex;
7051			PDL_Double beta = fabs(zeta * sum);
7052
7053			lodex = beta;
7054			beta -= lodex; hidex = lodex+1;
7055			alpha = gaussian_lookup[hidex]beta + gaussian_lookup[lodex](1 - beta);
7056
7057			}
7058			}
7059			break;
7060
7061			case 'G':
7062			/* This is the Gaussian rolloff with explicit calculation, preserved */
7063			/* in case someone actually wants the slower longer method. */
7064			{
7065			PDL_Double sum = 0;
7066			cp = tmp;
7067			for(i=0; i
7068			dd = 0;
7069			ap = tvec;
7070			for(j=0;j
7071			dd += (ap++) *(cp++);
7072			dd /= blur;
7073			sum += dd * dd;
7074			if(sum > 4) /* 2 pixels -- four half-widths */
7075			i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */
7076			}
7077
7078			if(sum > GAUSSIAN_MAXVAL)
7079			alpha = 0;
7080			else
7081			alpha = exp(-sum * 1.386294); /* Gaussian, rt(2)-pix HWHM */
7082			}
7083			break;
7084			default:
7085			return PDL->make_error(PDL_EUSERERROR, "Error in map:" "This can't happen: method='%c'",method);
7086			}
7087
7088			{ /* convenience block -- accumulate the current point into the weighted sum. */
7089			/* This is more than simple assignment because we have our own explicit poor */
7090			/* man's broadcastloop here, so we accumulate each broadcasted element separately. */
7091			PDL_LDouble dat = ((PDL_LDouble )(in->data)) + i_off;
7092			PDL_Indx max = out->dims[ndims];
7093			for( i=0; i < max; i++ ) {
7094			if( (badval==0) \|\| (*dat != badval) ) {
7095			acc[i] += dat alpha;
7096			dat += in->dimincs[ndims];
7097			wgt[i] += alpha;
7098			}
7099			wgt2[i] += alpha; }
7100			}
7101			} /* end of t_vio check (i.e. of input accumulation) */
7102
7103
7104			/* Advance input accumulation loop. */
7105			/* We both increment the total vector and also advance the index. */
7106			carry = 1;
7107			for (i=0; i
7108			/* Advance the current element of the offset vector */
7109			ivec[i]++;
7110			j = ivec[i] + ibvec[i];
7111
7112			/* Advance the offset into the data array */
7113			if( j > 0 && j <= in->dims[i]-1 ) {
7114			/* Normal case -- just advance the input vector */
7115			i_off += in->dimincs[i];
7116			} else {
7117			/* Busted a boundary - either before or after. */
7118			switch(bounds[i]){
7119			case 0: /* no breakage allowed -- treat as truncation for interpolation */
7120			case 1: /* truncation -- if we crossed the boundary mark ourselves out-of-bounds */
7121			if( j == 0 )
7122			t_vio--;
7123			else if( j == in->dims[i] )
7124			t_vio++;
7125			break;
7126			case 2: /* extension -- do nothing (so the same input point is re-used) */
7127			break;
7128			case 3: /* periodic -- advance and mod into the allowed range */
7129			if((j % in->dims[i]) == 0) {
7130			i_off -= in->dimincs[i] * (in->dims[i]-1);
7131			} else {
7132			i_off += in->dimincs[i];
7133			}
7134			break;
7135			case 4: /* mirror -- advance or retreat depending on phase */
7136			j += in->dims[i];
7137			j %= (in->dims[i]*2);
7138			j -= in->dims[i];
7139			if( j!=0 && j!= -in->dims[i] ) {
7140			if(j<0)
7141			i_off -= in->dimincs[i];
7142			else
7143			i_off += in->dimincs[i];
7144			}
7145			break;
7146			}
7147			}
7148
7149			/* Now check for carry */
7150			if (ivec[i] <= psize) {
7151			/* Normal case -- copy the current offset to the faster-running dim stashes */
7152			PDL_Indx k;
7153			for(k=0;k
7154			index_stash[k] = i_off;
7155			}
7156			carry = 0;
7157
7158			} else { /* End of this scan -- recover the last position, and mark carry */
7159			i_off = index_stash[i];
7160			if(bounds[i]==1) {
7161			j = ivec[i] + ibvec[i];
7162			if( j < 0 \|\| j >= in->dims[i] )
7163			t_vio--;
7164			ivec[i] = -psize;
7165			j = ivec[i] + ibvec[i];
7166			if( j < 0 \|\| j >= in->dims[i] )
7167			t_vio++;
7168			carry = 1;
7169			} else {
7170			ivec[i] = -psize;
7171			}
7172			}
7173			} /* End of counter-advance loop */
7174			} while (carry==0); /* end of total data accumulation loop (termination condition has carry on last dim) */
7175
7176			{
7177			PDL_Double *ac = acc;
7178			PDL_Double *wg = wgt;
7179			PDL_Double *wg2 = wgt2;
7180			PDL_LDouble *dat = out->data;
7181
7182			/* Calculate output vector offset */
7183			for(i=0;i
7184			dat += out->dimincs[i] * ovec[i];
7185
7186			if (!flux) {
7187			/* Flux flag is NOT set -- normal case. Copy the weighted accumulated data. */
7188			for (i=0; i < out->dims[ndims]; i++) {
7189			dat = (wg && (wg2 / wg) < 1.5) ? ac / wg : badval;
7190			ac++; wg++; wg2++;
7191			dat += out->dimincs[ndims];
7192			}
7193			} else {
7194			/* Flux flag is set - scale by the (unpadded) determinant of the Jacobian */
7195			PDL_Double det = tmp[ndims*ndims];
7196			for(i=0; i < out->dims[ndims]; i++) {
7197			if(wg && (wg2 / *wg) < 1.5 ) {
7198			dat = (ac++) / (wg++) det;
7199			wg2++;
7200			} else {
7201			*dat = badval;
7202			ac++; wg++; wg2++;
7203			}
7204			dat += out->dimincs[ndims];
7205			} /* end of for loop */
7206			} /* end of flux flag set conditional */
7207			} /* end of convenience block */
7208
7209			/* End of code for normal pixels */
7210			} else {
7211			/* The pixel was ludicrously huge -- just set this pixel to nan */
7212			PDL_LDouble *dat = out->data;
7213			for(i=0;i
7214			dat += out->dimincs[i] * ovec[i];
7215			for(i=0;idims[ndims];i++) {
7216			dat = badval; / Should handle bad values too -- not yet */
7217			dat += out->dimincs[ndims];
7218			}
7219			}
7220
7221			/* Increment the pixel counter */
7222			if (!pdl_broadcast_nd_step(npdls, offs, 0, ndims, incs, out->dims, ovec)) break;
7223			} while (1);
7224			PDL->changed(out, PDL_PARENTDATACHANGED, 0);
7225			#line 7226 "lib/PDL/Transform-pp-map.c"
7226	0	0	}PDL_BROADCASTLOOP_END_map_readdata
		0
7227	0		} break;
7228	0		default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in map: unhandled datatype(%d), only handles (ABSULKNPQFDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
7229			}
7230			#undef PDL_IF_BAD
7231	12		return PDL_err;
7232			}
7233
7234
7235			#line 1857 "lib/PDL/PP.pm"
7236			pdl_error pdl_map_free(pdl_trans *__privtrans, char destroy) {
7237			pdl_error PDL_err = {0, NULL, 0};
7238			#line 7239 "lib/PDL/Transform-pp-map.c"
7239	12		pdl_params_map *__params = __privtrans->params; (void)__params;
7240	12	50	PDL_FREE_CODE(__privtrans, destroy, SvREFCNT_dec(__params->boundary); /* CType.get_free */
7241			SvREFCNT_dec(__params->method); /* CType.get_free */
7242			SvREFCNT_dec(__params->bv); /* CType.get_free */
7243	12		, ) return PDL_err;
7244			}
7245
7246			static pdl_datatypes pdl_map_vtable_gentypes[] = { PDL_SB, PDL_B, PDL_S, PDL_US, PDL_L, PDL_UL, PDL_IND, PDL_ULL, PDL_LL, PDL_F, PDL_D, PDL_LD, -1 };
7247			static PDL_Indx pdl_map_vtable_realdims[] = { 0 };
7248			static char *pdl_map_vtable_parnames[] = { "k0" };
7249			static short pdl_map_vtable_parflags[] = {
7250			0
7251			};
7252			static pdl_datatypes pdl_map_vtable_partypes[] = { -1 };
7253			static PDL_Indx pdl_map_vtable_realdims_starts[] = { 0 };
7254			static PDL_Indx pdl_map_vtable_realdims_ind_ids[] = { 0 };
7255			static char *pdl_map_vtable_indnames[] = { "" };
7256			pdl_transvtable pdl_map_vtable = {
7257			PDL_TRANS_DO_BROADCAST, 0, pdl_map_vtable_gentypes, 1, 1, NULL /CORE21/,
7258			pdl_map_vtable_realdims, pdl_map_vtable_parnames,
7259			pdl_map_vtable_parflags, pdl_map_vtable_partypes,
7260			pdl_map_vtable_realdims_starts, pdl_map_vtable_realdims_ind_ids, 0,
7261			0, pdl_map_vtable_indnames,
7262			NULL, pdl_map_readdata, NULL,
7263			pdl_map_free,
7264			sizeof(pdl_params_map),"PDL::Transform::map"
7265			};
7266
7267
7268	12		pdl_error pdl_run_map(pdl k0,pdl in,pdl out,pdl map,SV boundary,SV method,long big,double blur,double sv_min,char flux,SV *bv) {
7269	12		pdl_error PDL_err = {0, NULL, 0};
7270	12	50	if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
7271	12		pdl_trans *__privtrans = PDL->create_trans(&pdl_map_vtable);
7272	12	50	if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
7273	12		pdl_params_map *__params = __privtrans->params;
7274	12		__privtrans->pdls[0] = k0;
7275	12	50	PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
7276	12		(__params->in) = (in); /* CType.get_copy */
7277	12		(__params->out) = (out); /* CType.get_copy */
7278	12		(__params->map) = (map); /* CType.get_copy */
7279	12		(__params->boundary) = newSVsv(boundary); /* CType.get_copy */
7280	12		(__params->method) = newSVsv(method); /* CType.get_copy */
7281	12		(__params->big) = (big); /* CType.get_copy */
7282	12		(__params->blur) = (blur); /* CType.get_copy */
7283	12		(__params->sv_min) = (sv_min); /* CType.get_copy */
7284	12		(__params->flux) = (flux); /* CType.get_copy */
7285	12		(__params->bv) = newSVsv(bv); /* CType.get_copy */
7286	12	50	PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
7287	12		return PDL_err;
7288			}