1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
LANCZOS_INTERP_FUNC(PyObject * py_ixi,PyObject * py_iyi,PyObject * py_dx,PyObject * py_dy,PyObject * loutputs,PyObject * linputs)5 static int LANCZOS_INTERP_FUNC(PyObject* py_ixi, PyObject* py_iyi,
6 PyObject* py_dx, PyObject* py_dy,
7 PyObject* loutputs, PyObject* linputs) {
8 npy_intp W,H, N;
9 npy_intp Nimages;
10 npy_intp i, j;
11 PyArray_Descr* dtype = PyArray_DescrFromType(NPY_FLOAT);
12 PyArray_Descr* itype = PyArray_DescrFromType(NPY_INT32);
13 int req = NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED |
14 NPY_ARRAY_NOTSWAPPED | NPY_ARRAY_ELEMENTSTRIDES;
15 int reqout = req | NPY_ARRAY_WRITEABLE | NPY_ARRAY_UPDATEIFCOPY;
16
17 const int32_t *ixi, *iyi;
18 const float *dx, *dy;
19
20 PyArrayObject *np_ixi, *np_iyi, *np_dx, *np_dy;
21
22 /*
23 dx,dy are in [-0.5, 0.5].
24
25 Lanczos-3 kernel is zero outside [-3,3].
26
27 We build a look-up table where [0] is L(-3.5).
28
29 And organized so that:
30 lut[0] = L(-3.5)
31 lut[1] = L(-2.5)
32 lut[2] = L(-1.5)
33 lut[3] = L(-0.5)
34 lut[4] = L( 0.5)
35 lut[5] = L( 1.5)
36 lut[6] = L( 2.5)
37 lut[7] is empty for padding
38 actually, = sum(lut[0:7])
39
40 lut[8] = L(-3.499)
41 lut[9] = L(-2.499)
42 lut[10] = L(-1.499)
43 ...
44 ...
45 lut[8184] = L(-2.501)
46 lut[8185] = L(-1.501)
47 lut[8186] = L(-0.501)
48 ...
49
50 This is annoying because [-3.5,3] and [3,3.5] are zero so we
51 have to sum 7 elements rather than 6. But the alternatives
52 seem worse.
53 */
54
55 //static const int L = 5;
56 // Nlutunit is number of bins per unit x
57 //static const int Nlutunit = 1024;
58 static const int Nlutunit = 2048;
59 static const double lut0 = -(L + 0.5); //-5.5; //-(L+0.5);
60 static const int Nunits = 2*(L+1); //12; //2*(L+1);
61 //static const int Nlut = Nunits * Nlutunit;
62 //static float lut[24576];
63 //static float lut[Nunits*Nlutunit];
64 // HACK -- 2048 here == Nlutunit... some versions of gcc don't believe Nunits*Nlutunit is constant
65 static float lut[2*(L+1)*2048];
66 static int initialized = 0;
67
68 if (!initialized) {
69 // this table has the elements you need to use together
70 // stored together: L(x[0]), L(x[0]+1), L(x[0]+2), ...;
71 // L(x[1]), L(x[1]+1), L(x[2]+2), ...
72 for (i=0; i<Nlutunit; i++) {
73 double x,f;
74 double acc = 0.;
75 x = (lut0 + ((i+0.5) / (double)Nlutunit));
76 for (j=0; j<Nunits; j++, x+=1.0) {
77 if (x <= -L || x >= L) {
78 f = 0.0;
79 } else if (x == 0) {
80 f = 1.0;
81 } else {
82 f = L * sin(M_PI * x) * sin(M_PI / L * x) /
83 (M_PI * M_PI * x * x);
84 }
85 lut[i * Nunits + j] = f;
86 acc += f;
87 }
88 lut[i*Nunits + Nunits-1] = acc;
89 //printf("acc: %f\n", acc);
90 }
91 /*
92 for (i=0; i<Nlut; i++) {
93 printf("lut[% 4li] = %f\n", i, lut[i]);
94 }
95 */
96 initialized = 1;
97 }
98
99 // CheckFromAny steals the dtype reference
100 Py_INCREF(itype);
101 np_ixi = (PyArrayObject*)PyArray_CheckFromAny(py_ixi, itype, 1, 1, req, NULL);
102 np_iyi = (PyArrayObject*)PyArray_CheckFromAny(py_iyi, itype, 1, 1, req, NULL);
103 // At this point, itype refcount = 0
104 Py_INCREF(dtype);
105 Py_INCREF(dtype);
106 np_dx = (PyArrayObject*)PyArray_CheckFromAny(py_dx, dtype, 1, 1, req, NULL);
107 np_dy = (PyArrayObject*)PyArray_CheckFromAny(py_dy, dtype, 1, 1, req, NULL);
108 // dtype refcount = 1 (we use it more below)
109 if (!np_ixi || !np_iyi) {
110 ERR("ixi,iyi arrays are wrong type / shape\n");
111 return -1;
112 }
113 if (!np_dx || !np_dy) {
114 ERR("dx,dy arrays are wrong type / shape\n");
115 return -1;
116 }
117 N = PyArray_DIM(np_ixi, 0);
118 if ((PyArray_DIM(np_iyi, 0) != N) ||
119 (PyArray_DIM(np_dx, 0) != N) ||
120 (PyArray_DIM(np_dy, 0) != N)) {
121 ERR("ixi,iyi,dx,dy arrays must be same size\n");
122 return -1;
123 }
124
125 if (!PyList_Check(loutputs) ||
126 !PyList_Check(linputs)) {
127 ERR("loutputs and linputs must be lists of np arrays\n");
128 return -1;
129 }
130 Nimages = PyList_Size(loutputs);
131 if (PyList_Size(linputs) != Nimages) {
132 ERR("loutputs and linputs must be same length\n");
133 return -1;
134 }
135
136 for (i=0; i<Nimages; i++) {
137 PyArrayObject* np_inimg;
138 PyArrayObject* np_outimg;
139 const float *inimg;
140 float *outimg;
141
142 ixi = PyArray_DATA(np_ixi);
143 iyi = PyArray_DATA(np_iyi);
144 dx = PyArray_DATA(np_dx);
145 dy = PyArray_DATA(np_dy);
146
147 Py_INCREF(dtype);
148 Py_INCREF(dtype);
149 np_inimg = (PyArrayObject*)PyArray_CheckFromAny(PyList_GetItem(linputs, i), dtype, 2, 2, req, NULL);
150 np_outimg = (PyArrayObject*)PyArray_CheckFromAny(PyList_GetItem(loutputs, i), dtype, 1, 1, reqout, NULL);
151 if (!np_inimg || !np_outimg) {
152 ERR("Failed to convert input and output images to right type/shape\n");
153 return -1;
154 }
155 if (PyArray_DIM(np_outimg, 0) != N) {
156 ERR("Output image must be same shape as ixo\n");
157 return -1;
158 }
159 H = PyArray_DIM(np_inimg, 0);
160 W = PyArray_DIM(np_inimg, 1);
161 inimg = PyArray_DATA(np_inimg);
162 outimg = PyArray_DATA(np_outimg);
163
164 for (j=0; j<N; j++, outimg++, ixi++, iyi++) {
165 // resample inimg[ iyi[j] + dy[j], ixi[j] + dx[j] ]
166 // to outimg[ j ]
167 npy_intp u,v;
168 int tx0, ty0;
169 float acc = 0.;
170 float nacc;
171 const float* ly;
172 int ix,iy;
173 tx0 = (int)((-(dx[j]+L) - lut0) * Nlutunit);
174 ty0 = (int)((-(dy[j]+L) - lut0) * Nlutunit);
175 // clip
176 tx0 = MAX(0, MIN(Nlutunit-1, tx0));
177 ty0 = MAX(0, MIN(Nlutunit-1, ty0));
178 tx0 *= Nunits;
179 ty0 *= Nunits;
180 ly = lut + ty0;
181 iy = *iyi;
182 ix = *ixi;
183 // special-case pixels near the image edges.
184 if (ix < L || ix >= (W-L) || iy < L || iy >= (H-L)) {
185 iy -= L;
186 // Lanczos kernel in y direction
187 for (v=0; v<2*L+1; v++, iy++, ly++) {
188 float accx = 0;
189 int clipiy = MAX(0, MIN((int)(H-1), iy));
190 int ix = *ixi - L;
191 const float* lx = lut + tx0;
192 const float* inpix = inimg + clipiy * W;
193 // Lanczos kernel in x direction
194 for (u=0; u<2*L+1; u++, ix++, lx++) {
195 int clipix = MAX(0, MIN((int)(W-1), ix));
196 accx += (*lx) * (inpix[clipix]);
197 }
198 acc += (*ly) * accx;
199 }
200 } else {
201 iy -= L;
202 // Lanczos kernel in y direction
203 for (v=0; v<2*L+1; v++,
204 iy++, ly++) {
205 float accx = 0;
206 int ix = *ixi - L;
207 const float* lx = lut + tx0;
208 const float* inpix = inimg + iy * W + ix;
209 // Lanczos kernel in x direction
210 for (u=0; u<2*L+1; u++,
211 lx++, inpix++) {
212 accx += (*lx) * (*inpix);
213 }
214 acc += (*ly) * accx;
215 }
216 }
217 nacc = lut[tx0 + Nunits-1] * lut[ty0 + Nunits-1];
218 *outimg = acc / nacc;
219 }
220 Py_DECREF(np_inimg);
221 Py_DECREF(np_outimg);
222 }
223 Py_DECREF(dtype);
224 Py_DECREF(np_ixi);
225 Py_DECREF(np_iyi);
226 Py_DECREF(np_dx);
227 Py_DECREF(np_dy);
228 return 0;
229 }
230