File Coverage

lib/PDL/Image2D-pp-warp2d.c
Criterion Covered Total %
statement 116 254 45.6
branch 53 230 23.0
condition n/a
subroutine n/a
pod n/a
total 169 484 34.9


line stmt bran cond sub pod time code
1              
2             #line 453 "lib/PDL/PP.pm"
3             /*
4             * THIS FILE WAS GENERATED BY PDL::PP from lib/PDL/Image2D.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_Image2D
21             extern Core* PDL; /* Structure hold core C functions */
22             #line 23 "lib/PDL/Image2D-pp-warp2d.c"
23              
24             /* Fast Modulus with proper negative behaviour */
25              
26             #define REALMOD(a,b) {while ((a)>=(b)) (a) -= (b); while ((a)<0) (a) += (b);}
27              
28             #define X(symbol, ctype, ppsym, ...) \
29             ctype quick_select_ ## ppsym(ctype arr[], int n);
30             PDL_TYPELIST_REAL(X)
31             #undef X
32             #define EZ(x) ez ? 0 : (x)
33             int getnewsize(PDL_Indx cols, PDL_Indx rows, float fangle, PDL_Indx *newcols, PDL_Indx *newrows);
34             typedef unsigned char imT; /* image type */
35             int rotate(imT *im, imT *out, int cols, int rows, int nc, int nr,
36             float fangle, imT bgval, int antialias);
37             #include "resample.h"
38              
39             #line 1846 "lib/PDL/PP.pm"
40             typedef struct pdl_params_warp2d {
41             #line 42 "lib/PDL/Image2D-pp-warp2d.c"
42             char *kernel_type;
43             double noval;
44             IV nsamples;
45             } pdl_params_warp2d;
46              
47              
48             #line 1857 "lib/PDL/PP.pm"
49             pdl_error pdl_warp2d_redodims(pdl_trans *__privtrans) {
50             pdl_error PDL_err = {0, NULL, 0};
51             #line 52 "lib/PDL/Image2D-pp-warp2d.c"
52 6           pdl_params_warp2d *__params = __privtrans->params; (void)__params;
53 6           __privtrans->ind_sizes[3] = __params->nsamples;
54 6 50         PDL_RETERROR(PDL_err, PDL->redodims_default(__privtrans));
55 6           return PDL_err;
56             }
57              
58              
59             #line 1857 "lib/PDL/PP.pm"
60             pdl_error pdl_warp2d_readdata(pdl_trans *__privtrans) {
61             pdl_error PDL_err = {0, NULL, 0};
62             #line 63 "lib/PDL/Image2D-pp-warp2d.c"
63 6           pdl_params_warp2d *__params = __privtrans->params; (void)__params;
64 6           register PDL_Indx __m_size = __privtrans->ind_sizes[0];
65 6           register PDL_Indx __n_size = __privtrans->ind_sizes[1];
66 6           register PDL_Indx __np_size = __privtrans->ind_sizes[2];
67 6 50         if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in warp2d:" "broadcast.incs NULL");
68             /* broadcastloop declarations */
69             int __brcloopval;
70             register PDL_Indx __tind0,__tind1; /* counters along dim */
71 6           register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
72             /* dims here are how many steps along those dims */
73 6           register PDL_Indx __tinc0_img = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
74 6           register PDL_Indx __tinc0_px = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
75 6           register PDL_Indx __tinc0_py = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
76 6           register PDL_Indx __tinc0_warp = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,0);
77 6           register PDL_Indx __tinc0_poly = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,0);
78 6           register PDL_Indx __tinc0_kernel = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,0);
79 6           register PDL_Indx __tinc1_img = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
80 6           register PDL_Indx __tinc1_px = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
81 6           register PDL_Indx __tinc1_py = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
82 6           register PDL_Indx __tinc1_warp = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,1);
83 6           register PDL_Indx __tinc1_poly = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,1);
84 6           register PDL_Indx __tinc1_kernel = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,1);
85             #define PDL_BROADCASTLOOP_START_warp2d_readdata PDL_BROADCASTLOOP_START( \
86             readdata, \
87             __privtrans->broadcast, \
88             __privtrans->vtable, \
89             img_datap += __offsp[0]; \
90             px_datap += __offsp[1]; \
91             py_datap += __offsp[2]; \
92             warp_datap += __offsp[3]; \
93             poly_datap += __offsp[4]; \
94             kernel_datap += __offsp[5]; \
95             , \
96             ( ,img_datap += __tinc1_img - __tinc0_img * __tdims0 \
97             ,px_datap += __tinc1_px - __tinc0_px * __tdims0 \
98             ,py_datap += __tinc1_py - __tinc0_py * __tdims0 \
99             ,warp_datap += __tinc1_warp - __tinc0_warp * __tdims0 \
100             ,poly_datap += __tinc1_poly - __tinc0_poly * __tdims0 \
101             ,kernel_datap += __tinc1_kernel - __tinc0_kernel * __tdims0 \
102             ), \
103             ( ,img_datap += __tinc0_img \
104             ,px_datap += __tinc0_px \
105             ,py_datap += __tinc0_py \
106             ,warp_datap += __tinc0_warp \
107             ,poly_datap += __tinc0_poly \
108             ,kernel_datap += __tinc0_kernel \
109             ) \
110             )
111             #define PDL_BROADCASTLOOP_END_warp2d_readdata PDL_BROADCASTLOOP_END( \
112             __privtrans->broadcast, \
113             img_datap -= __tinc1_img * __tdims1 + __offsp[0]; \
114             px_datap -= __tinc1_px * __tdims1 + __offsp[1]; \
115             py_datap -= __tinc1_py * __tdims1 + __offsp[2]; \
116             warp_datap -= __tinc1_warp * __tdims1 + __offsp[3]; \
117             poly_datap -= __tinc1_poly * __tdims1 + __offsp[4]; \
118             kernel_datap -= __tinc1_kernel * __tdims1 + __offsp[5]; \
119             )
120 6           register PDL_Indx __inc_img_m = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,0)]; (void)__inc_img_m;register PDL_Indx __inc_img_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,1)]; (void)__inc_img_n;
121 6           register PDL_Indx __inc_kernel_ns = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,5,0)]; (void)__inc_kernel_ns;
122 6           register PDL_Indx __inc_poly_np = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,4,0)]; (void)__inc_poly_np;
123 6           register PDL_Indx __inc_px_np0 = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,0)]; (void)__inc_px_np0;register PDL_Indx __inc_px_np1 = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,1)]; (void)__inc_px_np1;
124 6           register PDL_Indx __inc_py_np0 = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,2,0)]; (void)__inc_py_np0;register PDL_Indx __inc_py_np1 = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,2,1)]; (void)__inc_py_np1;
125 6           register PDL_Indx __inc_warp_m = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,3,0)]; (void)__inc_warp_m;register PDL_Indx __inc_warp_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,3,1)]; (void)__inc_warp_n;
126             #ifndef PDL_DECLARE_PARAMS_warp2d_1
127             #define PDL_DECLARE_PARAMS_warp2d_1(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_px,PDL_PPSYM_PARAM_px,PDL_TYPE_PARAM_py,PDL_PPSYM_PARAM_py,PDL_TYPE_PARAM_poly,PDL_PPSYM_PARAM_poly,PDL_TYPE_PARAM_kernel,PDL_PPSYM_PARAM_kernel) \
128             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, img, (__privtrans->pdls[0]), 1, PDL_PPSYM_OP) \
129             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_px, px, (__privtrans->pdls[1]), 1, PDL_PPSYM_PARAM_px) \
130             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_py, py, (__privtrans->pdls[2]), 1, PDL_PPSYM_PARAM_py) \
131             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, warp, (__privtrans->pdls[3]), 1, PDL_PPSYM_OP) \
132             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_poly, poly, (__privtrans->pdls[4]), 1, PDL_PPSYM_PARAM_poly) \
133             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_kernel, kernel, (__privtrans->pdls[5]), 1, PDL_PPSYM_PARAM_kernel)
134             #endif
135             #define PDL_IF_BAD(t,f) f
136 6           switch (__privtrans->__datatype) { /* Start generic switch */
137 0           case PDL_F: {
138 0 0         PDL_DECLARE_PARAMS_warp2d_1(PDL_Float,F,PDL_LDouble,E,PDL_LDouble,E,PDL_LDouble,E,PDL_LDouble,E)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
139             { PDL_Indx k, lx_3, ly_3, px, py, tabx, taby;
140             PDL_Indx da[16], db[16] ;
141 0           PDL_LDouble cur, neighbors[16], rsc[8], sumrs, x, y, *kernel = kernel_datap, *poly = poly_datap;
142              
143             /* Generate interpolation kernel */
144 0 0         if (!generate_interpolation_kernel(__params->kernel_type, __privtrans->ind_sizes[3], kernel))
145 0           return PDL->make_error(PDL_EUSERERROR, "Error in warp2d:" "Invalid kernel type '%s'",__params->kernel_type);
146              
147             /* Compute sizes */
148 0           lx_3 = __privtrans->ind_sizes[0] - 3;
149 0           ly_3 = __privtrans->ind_sizes[1] - 3;
150              
151             /* Pre compute leaps for 16 closest neighbors positions */
152 0           da[0] = -1; db[0] = -1;
153 0           da[1] = 0; db[1] = -1;
154 0           da[2] = 1; db[2] = -1;
155 0           da[3] = 2; db[3] = -1;
156              
157 0           da[4] = -1; db[4] = 0;
158 0           da[5] = 0; db[5] = 0;
159 0           da[6] = 1; db[6] = 0;
160 0           da[7] = 2; db[7] = 0;
161              
162 0           da[8] = -1; db[8] = 1;
163 0           da[9] = 0; db[9] = 1;
164 0           da[10] = 1; db[10] = 1;
165 0           da[11] = 2; db[11] = 1;
166              
167 0           da[12] = -1; db[12] = 2;
168 0           da[13] = 0; db[13] = 2;
169 0           da[14] = 1; db[14] = 2;
170 0           da[15] = 2; db[15] = 2;
171              
172 0           poly[0] = 1.0;
173              
174             /* Loop over the output image */
175 0 0         PDL_BROADCASTLOOP_START_warp2d_readdata
    0          
    0          
    0          
    0          
    0          
    0          
176 0 0         {/* Open n */ PDL_EXPAND2(register PDL_Indx n=0, __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
177              
178             /* fill in poly array */
179 0 0         {/* Open np=1 */ PDL_EXPAND2(register PDL_Indx np=PDLMAX((1),0), __np_stop=(__np_size)); for(; np<__np_stop; np+=1) {
180 0           poly[np] = (PDL_LDouble) n * poly[np-1];
181             }} /* Close np=1 */
182              
183 0 0         {/* Open m */ PDL_EXPAND2(register PDL_Indx m=0, __m_stop=(__m_size)); for(; m<__m_stop; m+=1) {
184              
185             /* Compute the original source for this pixel */
186 0           x = poly2d_compute( __privtrans->ind_sizes[2], px_datap, m, poly );
187 0           y = poly2d_compute( __privtrans->ind_sizes[2], py_datap, m, poly );
188              
189             /* Which is the closest integer positioned neighbor? */
190 0           px = (int)x ;
191 0           py = (int)y ;
192              
193 0 0         if ((px < 1) || (px > lx_3) || (py < 1) || (py > ly_3))
    0          
    0          
    0          
194 0           (warp_datap)[0+(__inc_warp_m*(m))+(__inc_warp_n*(n))] = (PDL_Float) __params->noval;
195             else {
196              
197             /* Now feed the positions for the closest 16 neighbors */
198 0 0         for (k=0 ; k<16 ; k++) {
199 0           neighbors[k] = (PDL_LDouble) (img_datap)[0+(__inc_img_m*(px+da[k]))+(__inc_img_n*(py+db[k]))];
200             }
201              
202             /* Which tabulated value index shall we use? */
203 0           tabx = (x - (PDL_LDouble)px) * (PDL_LDouble)(TABSPERPIX) ;
204 0           taby = (y - (PDL_LDouble)py) * (PDL_LDouble)(TABSPERPIX) ;
205              
206             /* Compute resampling coefficients */
207             /* rsc[0..3] in x, rsc[4..7] in y */
208              
209 0           rsc[0] = kernel[TABSPERPIX + tabx] ;
210 0           rsc[1] = kernel[tabx] ;
211 0           rsc[2] = kernel[TABSPERPIX - tabx] ;
212 0           rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
213 0           rsc[4] = kernel[TABSPERPIX + taby] ;
214 0           rsc[5] = kernel[taby] ;
215 0           rsc[6] = kernel[TABSPERPIX - taby] ;
216 0           rsc[7] = kernel[2 * TABSPERPIX - taby] ;
217              
218 0           sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
219 0           (rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;
220              
221             /* Compute interpolated pixel now */
222 0           cur = rsc[4] * ( rsc[0]*neighbors[0] +
223 0           rsc[1]*neighbors[1] +
224 0           rsc[2]*neighbors[2] +
225 0           rsc[3]*neighbors[3] ) +
226 0           rsc[5] * ( rsc[0]*neighbors[4] +
227 0           rsc[1]*neighbors[5] +
228 0           rsc[2]*neighbors[6] +
229 0           rsc[3]*neighbors[7] ) +
230 0           rsc[6] * ( rsc[0]*neighbors[8] +
231 0           rsc[1]*neighbors[9] +
232 0           rsc[2]*neighbors[10] +
233 0           rsc[3]*neighbors[11] ) +
234 0           rsc[7] * ( rsc[0]*neighbors[12] +
235 0           rsc[1]*neighbors[13] +
236 0           rsc[2]*neighbors[14] +
237 0           rsc[3]*neighbors[15] ) ;
238              
239             /* Copy the value to the output image */
240 0           (warp_datap)[0+(__inc_warp_m*(m))+(__inc_warp_n*(n))] = (PDL_Float) (cur/sumrs);
241              
242             } /* if: edge or interior */
243              
244             }} /* Close m */ /* loop(m) */
245             }} /* Close n */ /* loop(n) */
246 0 0         PDL_BROADCASTLOOP_END_warp2d_readdata /* broadcastloop */
    0          
247             }
248 0           } break;
249 6           case PDL_D: {
250 6 50         PDL_DECLARE_PARAMS_warp2d_1(PDL_Double,D,PDL_LDouble,E,PDL_LDouble,E,PDL_LDouble,E,PDL_LDouble,E)
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
251             { PDL_Indx k, lx_3, ly_3, px, py, tabx, taby;
252             PDL_Indx da[16], db[16] ;
253 6           PDL_LDouble cur, neighbors[16], rsc[8], sumrs, x, y, *kernel = kernel_datap, *poly = poly_datap;
254              
255             /* Generate interpolation kernel */
256 6 50         if (!generate_interpolation_kernel(__params->kernel_type, __privtrans->ind_sizes[3], kernel))
257 0           return PDL->make_error(PDL_EUSERERROR, "Error in warp2d:" "Invalid kernel type '%s'",__params->kernel_type);
258              
259             /* Compute sizes */
260 6           lx_3 = __privtrans->ind_sizes[0] - 3;
261 6           ly_3 = __privtrans->ind_sizes[1] - 3;
262              
263             /* Pre compute leaps for 16 closest neighbors positions */
264 6           da[0] = -1; db[0] = -1;
265 6           da[1] = 0; db[1] = -1;
266 6           da[2] = 1; db[2] = -1;
267 6           da[3] = 2; db[3] = -1;
268              
269 6           da[4] = -1; db[4] = 0;
270 6           da[5] = 0; db[5] = 0;
271 6           da[6] = 1; db[6] = 0;
272 6           da[7] = 2; db[7] = 0;
273              
274 6           da[8] = -1; db[8] = 1;
275 6           da[9] = 0; db[9] = 1;
276 6           da[10] = 1; db[10] = 1;
277 6           da[11] = 2; db[11] = 1;
278              
279 6           da[12] = -1; db[12] = 2;
280 6           da[13] = 0; db[13] = 2;
281 6           da[14] = 1; db[14] = 2;
282 6           da[15] = 2; db[15] = 2;
283              
284 6           poly[0] = 1.0;
285              
286             /* Loop over the output image */
287 24 50         PDL_BROADCASTLOOP_START_warp2d_readdata
    50          
    50          
    50          
    50          
    100          
    100          
288 306 100         {/* Open n */ PDL_EXPAND2(register PDL_Indx n=0, __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
289              
290             /* fill in poly array */
291 900 100         {/* Open np=1 */ PDL_EXPAND2(register PDL_Indx np=PDLMAX((1),0), __np_stop=(__np_size)); for(; np<__np_stop; np+=1) {
292 600           poly[np] = (PDL_LDouble) n * poly[np-1];
293             }} /* Close np=1 */
294              
295 15300 100         {/* Open m */ PDL_EXPAND2(register PDL_Indx m=0, __m_stop=(__m_size)); for(; m<__m_stop; m+=1) {
296              
297             /* Compute the original source for this pixel */
298 15000           x = poly2d_compute( __privtrans->ind_sizes[2], px_datap, m, poly );
299 15000           y = poly2d_compute( __privtrans->ind_sizes[2], py_datap, m, poly );
300              
301             /* Which is the closest integer positioned neighbor? */
302 15000           px = (int)x ;
303 15000           py = (int)y ;
304              
305 15000 100         if ((px < 1) || (px > lx_3) || (py < 1) || (py > ly_3))
    100          
    100          
    100          
306 1882           (warp_datap)[0+(__inc_warp_m*(m))+(__inc_warp_n*(n))] = (PDL_Double) __params->noval;
307             else {
308              
309             /* Now feed the positions for the closest 16 neighbors */
310 223006 100         for (k=0 ; k<16 ; k++) {
311 209888           neighbors[k] = (PDL_LDouble) (img_datap)[0+(__inc_img_m*(px+da[k]))+(__inc_img_n*(py+db[k]))];
312             }
313              
314             /* Which tabulated value index shall we use? */
315 13118           tabx = (x - (PDL_LDouble)px) * (PDL_LDouble)(TABSPERPIX) ;
316 13118           taby = (y - (PDL_LDouble)py) * (PDL_LDouble)(TABSPERPIX) ;
317              
318             /* Compute resampling coefficients */
319             /* rsc[0..3] in x, rsc[4..7] in y */
320              
321 13118           rsc[0] = kernel[TABSPERPIX + tabx] ;
322 13118           rsc[1] = kernel[tabx] ;
323 13118           rsc[2] = kernel[TABSPERPIX - tabx] ;
324 13118           rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
325 13118           rsc[4] = kernel[TABSPERPIX + taby] ;
326 13118           rsc[5] = kernel[taby] ;
327 13118           rsc[6] = kernel[TABSPERPIX - taby] ;
328 13118           rsc[7] = kernel[2 * TABSPERPIX - taby] ;
329              
330 13118           sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
331 13118           (rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;
332              
333             /* Compute interpolated pixel now */
334 13118           cur = rsc[4] * ( rsc[0]*neighbors[0] +
335 13118           rsc[1]*neighbors[1] +
336 13118           rsc[2]*neighbors[2] +
337 13118           rsc[3]*neighbors[3] ) +
338 13118           rsc[5] * ( rsc[0]*neighbors[4] +
339 13118           rsc[1]*neighbors[5] +
340 13118           rsc[2]*neighbors[6] +
341 13118           rsc[3]*neighbors[7] ) +
342 13118           rsc[6] * ( rsc[0]*neighbors[8] +
343 13118           rsc[1]*neighbors[9] +
344 13118           rsc[2]*neighbors[10] +
345 13118           rsc[3]*neighbors[11] ) +
346 13118           rsc[7] * ( rsc[0]*neighbors[12] +
347 13118           rsc[1]*neighbors[13] +
348 13118           rsc[2]*neighbors[14] +
349 13118           rsc[3]*neighbors[15] ) ;
350              
351             /* Copy the value to the output image */
352 13118           (warp_datap)[0+(__inc_warp_m*(m))+(__inc_warp_n*(n))] = (PDL_Double) (cur/sumrs);
353              
354             } /* if: edge or interior */
355              
356             }} /* Close m */ /* loop(m) */
357             }} /* Close n */ /* loop(n) */
358 6 50         PDL_BROADCASTLOOP_END_warp2d_readdata /* broadcastloop */
    50          
359             }
360 6           } break;
361 0           case PDL_LD: {
362 0 0         PDL_DECLARE_PARAMS_warp2d_1(PDL_LDouble,E,PDL_LDouble,E,PDL_LDouble,E,PDL_LDouble,E,PDL_LDouble,E)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
363             { PDL_Indx k, lx_3, ly_3, px, py, tabx, taby;
364             PDL_Indx da[16], db[16] ;
365 0           PDL_LDouble cur, neighbors[16], rsc[8], sumrs, x, y, *kernel = kernel_datap, *poly = poly_datap;
366              
367             /* Generate interpolation kernel */
368 0 0         if (!generate_interpolation_kernel(__params->kernel_type, __privtrans->ind_sizes[3], kernel))
369 0           return PDL->make_error(PDL_EUSERERROR, "Error in warp2d:" "Invalid kernel type '%s'",__params->kernel_type);
370              
371             /* Compute sizes */
372 0           lx_3 = __privtrans->ind_sizes[0] - 3;
373 0           ly_3 = __privtrans->ind_sizes[1] - 3;
374              
375             /* Pre compute leaps for 16 closest neighbors positions */
376 0           da[0] = -1; db[0] = -1;
377 0           da[1] = 0; db[1] = -1;
378 0           da[2] = 1; db[2] = -1;
379 0           da[3] = 2; db[3] = -1;
380              
381 0           da[4] = -1; db[4] = 0;
382 0           da[5] = 0; db[5] = 0;
383 0           da[6] = 1; db[6] = 0;
384 0           da[7] = 2; db[7] = 0;
385              
386 0           da[8] = -1; db[8] = 1;
387 0           da[9] = 0; db[9] = 1;
388 0           da[10] = 1; db[10] = 1;
389 0           da[11] = 2; db[11] = 1;
390              
391 0           da[12] = -1; db[12] = 2;
392 0           da[13] = 0; db[13] = 2;
393 0           da[14] = 1; db[14] = 2;
394 0           da[15] = 2; db[15] = 2;
395              
396 0           poly[0] = 1.0;
397              
398             /* Loop over the output image */
399 0 0         PDL_BROADCASTLOOP_START_warp2d_readdata
    0          
    0          
    0          
    0          
    0          
    0          
400 0 0         {/* Open n */ PDL_EXPAND2(register PDL_Indx n=0, __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
401              
402             /* fill in poly array */
403 0 0         {/* Open np=1 */ PDL_EXPAND2(register PDL_Indx np=PDLMAX((1),0), __np_stop=(__np_size)); for(; np<__np_stop; np+=1) {
404 0           poly[np] = (PDL_LDouble) n * poly[np-1];
405             }} /* Close np=1 */
406              
407 0 0         {/* Open m */ PDL_EXPAND2(register PDL_Indx m=0, __m_stop=(__m_size)); for(; m<__m_stop; m+=1) {
408              
409             /* Compute the original source for this pixel */
410 0           x = poly2d_compute( __privtrans->ind_sizes[2], px_datap, m, poly );
411 0           y = poly2d_compute( __privtrans->ind_sizes[2], py_datap, m, poly );
412              
413             /* Which is the closest integer positioned neighbor? */
414 0           px = (int)x ;
415 0           py = (int)y ;
416              
417 0 0         if ((px < 1) || (px > lx_3) || (py < 1) || (py > ly_3))
    0          
    0          
    0          
418 0           (warp_datap)[0+(__inc_warp_m*(m))+(__inc_warp_n*(n))] = (PDL_LDouble) __params->noval;
419             else {
420              
421             /* Now feed the positions for the closest 16 neighbors */
422 0 0         for (k=0 ; k<16 ; k++) {
423 0           neighbors[k] = (PDL_LDouble) (img_datap)[0+(__inc_img_m*(px+da[k]))+(__inc_img_n*(py+db[k]))];
424             }
425              
426             /* Which tabulated value index shall we use? */
427 0           tabx = (x - (PDL_LDouble)px) * (PDL_LDouble)(TABSPERPIX) ;
428 0           taby = (y - (PDL_LDouble)py) * (PDL_LDouble)(TABSPERPIX) ;
429              
430             /* Compute resampling coefficients */
431             /* rsc[0..3] in x, rsc[4..7] in y */
432              
433 0           rsc[0] = kernel[TABSPERPIX + tabx] ;
434 0           rsc[1] = kernel[tabx] ;
435 0           rsc[2] = kernel[TABSPERPIX - tabx] ;
436 0           rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
437 0           rsc[4] = kernel[TABSPERPIX + taby] ;
438 0           rsc[5] = kernel[taby] ;
439 0           rsc[6] = kernel[TABSPERPIX - taby] ;
440 0           rsc[7] = kernel[2 * TABSPERPIX - taby] ;
441              
442 0           sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
443 0           (rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;
444              
445             /* Compute interpolated pixel now */
446 0           cur = rsc[4] * ( rsc[0]*neighbors[0] +
447 0           rsc[1]*neighbors[1] +
448 0           rsc[2]*neighbors[2] +
449 0           rsc[3]*neighbors[3] ) +
450 0           rsc[5] * ( rsc[0]*neighbors[4] +
451 0           rsc[1]*neighbors[5] +
452 0           rsc[2]*neighbors[6] +
453 0           rsc[3]*neighbors[7] ) +
454 0           rsc[6] * ( rsc[0]*neighbors[8] +
455 0           rsc[1]*neighbors[9] +
456 0           rsc[2]*neighbors[10] +
457 0           rsc[3]*neighbors[11] ) +
458 0           rsc[7] * ( rsc[0]*neighbors[12] +
459 0           rsc[1]*neighbors[13] +
460 0           rsc[2]*neighbors[14] +
461 0           rsc[3]*neighbors[15] ) ;
462              
463             /* Copy the value to the output image */
464 0           (warp_datap)[0+(__inc_warp_m*(m))+(__inc_warp_n*(n))] = (PDL_LDouble) (cur/sumrs);
465              
466             } /* if: edge or interior */
467              
468             }} /* Close m */ /* loop(m) */
469             }} /* Close n */ /* loop(n) */
470 0 0         PDL_BROADCASTLOOP_END_warp2d_readdata /* broadcastloop */
    0          
471             }
472 0           } break;
473 0           default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in warp2d: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
474             }
475             #undef PDL_IF_BAD
476 6           return PDL_err;
477             }
478              
479              
480             #line 1857 "lib/PDL/PP.pm"
481             pdl_error pdl_warp2d_free(pdl_trans *__privtrans, char destroy) {
482             pdl_error PDL_err = {0, NULL, 0};
483             #line 484 "lib/PDL/Image2D-pp-warp2d.c"
484 6           pdl_params_warp2d *__params = __privtrans->params; (void)__params;
485 6 50         PDL_FREE_CODE(__privtrans, destroy, free(__params->kernel_type); /* CType.get_free */
486 6           , ) return PDL_err;
487             }
488              
489             static pdl_datatypes pdl_warp2d_vtable_gentypes[] = { PDL_F, PDL_D, PDL_LD, -1 };
490             static PDL_Indx pdl_warp2d_vtable_realdims[] = { 2, 2, 2, 2, 1, 1 };
491             static char *pdl_warp2d_vtable_parnames[] = { "img","px","py","warp","poly","kernel" };
492             static short pdl_warp2d_vtable_parflags[] = {
493             0,
494             PDL_PARAM_ISPHYS|PDL_PARAM_ISTYPED,
495             PDL_PARAM_ISPHYS|PDL_PARAM_ISTYPED,
496             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISWRITE,
497             PDL_PARAM_ISCREAT|PDL_PARAM_ISPHYS|PDL_PARAM_ISTEMP|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE,
498             PDL_PARAM_ISCREAT|PDL_PARAM_ISPHYS|PDL_PARAM_ISTEMP|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE
499             };
500             static pdl_datatypes pdl_warp2d_vtable_partypes[] = { -1, PDL_LD, PDL_LD, -1, PDL_LD, PDL_LD };
501             static PDL_Indx pdl_warp2d_vtable_realdims_starts[] = { 0, 2, 4, 6, 8, 9 };
502             static PDL_Indx pdl_warp2d_vtable_realdims_ind_ids[] = { 0, 1, 2, 2, 2, 2, 0, 1, 2, 3 };
503             static char *pdl_warp2d_vtable_indnames[] = { "m","n","np","ns" };
504             pdl_transvtable pdl_warp2d_vtable = {
505             PDL_TRANS_DO_BROADCAST|PDL_TRANS_BADIGNORE, 0, pdl_warp2d_vtable_gentypes, 3, 6, NULL /*CORE21*/,
506             pdl_warp2d_vtable_realdims, pdl_warp2d_vtable_parnames,
507             pdl_warp2d_vtable_parflags, pdl_warp2d_vtable_partypes,
508             pdl_warp2d_vtable_realdims_starts, pdl_warp2d_vtable_realdims_ind_ids, 10,
509             4, pdl_warp2d_vtable_indnames,
510             pdl_warp2d_redodims, pdl_warp2d_readdata, NULL,
511             pdl_warp2d_free,
512             sizeof(pdl_params_warp2d),"PDL::Image2D::warp2d"
513             };
514              
515              
516 6           pdl_error pdl_run_warp2d(pdl *img,pdl *px,pdl *py,pdl *warp,char *kernel_type,double noval,IV nsamples) {
517 6           pdl_error PDL_err = {0, NULL, 0};
518 6 50         if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
519 6           pdl_trans *__privtrans = PDL->create_trans(&pdl_warp2d_vtable);
520 6 50         if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
521 6           pdl_params_warp2d *__params = __privtrans->params;
522 6           __privtrans->pdls[0] = img;
523 6           __privtrans->pdls[1] = px;
524 6           __privtrans->pdls[2] = py;
525 6           __privtrans->pdls[3] = warp;
526 6 50         PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
527 6           (__params->kernel_type) = malloc(strlen(kernel_type)+1); strcpy(__params->kernel_type,kernel_type); /* CType.get_copy */
528 6           (__params->noval) = (noval); /* CType.get_copy */
529 6           (__params->nsamples) = (nsamples); /* CType.get_copy */
530 6 50         PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
531 6           return PDL_err;
532             }