Cleanup: Remove more #if 0 blocks
[blender.git] / source / blender / blenlib / intern / math_interp.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2012 by Blender Foundation.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Sergey Sharybin
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  *
27  */
28
29 /** \file blender/blenlib/intern/math_interp.c
30  *  \ingroup bli
31  */
32
33 #include <math.h>
34
35 #include "BLI_math.h"
36
37 #include "BLI_strict_flags.h"
38
39 /**************************************************************************
40  *                            INTERPOLATIONS
41  *
42  * Reference and docs:
43  * http://wiki.blender.org/index.php/User:Damiles#Interpolations_Algorithms
44  ***************************************************************************/
45
46 /* BICUBIC Interpolation functions
47  * More info: http://wiki.blender.org/index.php/User:Damiles#Bicubic_pixel_interpolation
48  * function assumes out to be zero'ed, only does RGBA */
49
50 static float P(float k)
51 {
52         float p1, p2, p3, p4;
53         p1 = max_ff(k + 2.0f, 0.0f);
54         p2 = max_ff(k + 1.0f, 0.0f);
55         p3 = max_ff(k, 0.0f);
56         p4 = max_ff(k - 1.0f, 0.0f);
57         return (float)(1.0f / 6.0f) * (p1 * p1 * p1 - 4.0f * p2 * p2 * p2 + 6.0f * p3 * p3 * p3 - 4.0f * p4 * p4 * p4);
58 }
59
60 static void vector_from_float(const float *data, float vector[4], int components)
61 {
62         if (components == 1) {
63                 vector[0] = data[0];
64         }
65         else if (components == 3) {
66                 copy_v3_v3(vector, data);
67         }
68         else {
69                 copy_v4_v4(vector, data);
70         }
71 }
72
73 static void vector_from_byte(const unsigned char *data, float vector[4], int components)
74 {
75         if (components == 1) {
76                 vector[0] = data[0];
77         }
78         else if (components == 3) {
79                 vector[0] = data[0];
80                 vector[1] = data[1];
81                 vector[2] = data[2];
82         }
83         else {
84                 vector[0] = data[0];
85                 vector[1] = data[1];
86                 vector[2] = data[2];
87                 vector[3] = data[3];
88         }
89 }
90
91 /* BICUBIC INTERPOLATION */
92 BLI_INLINE void bicubic_interpolation(const unsigned char *byte_buffer, const float *float_buffer,
93                                       unsigned char *byte_output, float *float_output, int width, int height,
94                                       int components, float u, float v)
95 {
96         int i, j, n, m, x1, y1;
97         float a, b, w, wx, wy[4], out[4];
98
99         /* sample area entirely outside image? */
100         if (ceil(u) < 0 || floor(u) > width - 1 || ceil(v) < 0 || floor(v) > height - 1) {
101                 if (float_output) {
102                         copy_vn_fl(float_output, components, 0.0f);
103                 }
104                 if (byte_output) {
105                         copy_vn_uchar(byte_output, components, 0);
106                 }
107                 return;
108         }
109
110         i = (int)floor(u);
111         j = (int)floor(v);
112         a = u - (float)i;
113         b = v - (float)j;
114
115         zero_v4(out);
116
117 /* Optimized and not so easy to read */
118
119         /* avoid calling multiple times */
120         wy[0] = P(b - (-1));
121         wy[1] = P(b -  0);
122         wy[2] = P(b -  1);
123         wy[3] = P(b -  2);
124
125         for (n = -1; n <= 2; n++) {
126                 x1 = i + n;
127                 CLAMP(x1, 0, width - 1);
128                 wx = P((float)n - a);
129                 for (m = -1; m <= 2; m++) {
130                         float data[4];
131
132                         y1 = j + m;
133                         CLAMP(y1, 0, height - 1);
134                         /* normally we could do this */
135                         /* w = P(n-a) * P(b-m); */
136                         /* except that would call P() 16 times per pixel therefor pow() 64 times, better precalc these */
137                         w = wx * wy[m + 1];
138
139                         if (float_output) {
140                                 const float *float_data = float_buffer + width * y1 * components + components * x1;
141
142                                 vector_from_float(float_data, data, components);
143                         }
144                         else {
145                                 const unsigned char *byte_data = byte_buffer + width * y1 * components + components * x1;
146
147                                 vector_from_byte(byte_data, data, components);
148                         }
149
150                         if (components == 1) {
151                                 out[0] += data[0] * w;
152                         }
153                         else if (components == 3) {
154                                 out[0] += data[0] * w;
155                                 out[1] += data[1] * w;
156                                 out[2] += data[2] * w;
157                         }
158                         else {
159                                 out[0] += data[0] * w;
160                                 out[1] += data[1] * w;
161                                 out[2] += data[2] * w;
162                                 out[3] += data[3] * w;
163                         }
164                 }
165         }
166
167 /* Done with optimized part */
168
169         if (float_output) {
170                 if (components == 1) {
171                         float_output[0] = out[0];
172                 }
173                 else if (components == 3) {
174                         copy_v3_v3(float_output, out);
175                 }
176                 else {
177                         copy_v4_v4(float_output, out);
178                 }
179         }
180         else {
181                 if (components == 1) {
182                         byte_output[0] = (unsigned char)(out[0] + 0.5f);
183                 }
184                 else if (components == 3) {
185                         byte_output[0] = (unsigned char)(out[0] + 0.5f);
186                         byte_output[1] = (unsigned char)(out[1] + 0.5f);
187                         byte_output[2] = (unsigned char)(out[2] + 0.5f);
188                 }
189                 else {
190                         byte_output[0] = (unsigned char)(out[0] + 0.5f);
191                         byte_output[1] = (unsigned char)(out[1] + 0.5f);
192                         byte_output[2] = (unsigned char)(out[2] + 0.5f);
193                         byte_output[3] = (unsigned char)(out[3] + 0.5f);
194                 }
195         }
196 }
197
198 void BLI_bicubic_interpolation_fl(const float *buffer, float *output, int width, int height,
199                                   int components, float u, float v)
200 {
201         bicubic_interpolation(NULL, buffer, NULL, output, width, height, components, u, v);
202 }
203
204 void BLI_bicubic_interpolation_char(const unsigned char *buffer, unsigned char *output, int width, int height,
205                                     int components, float u, float v)
206 {
207         bicubic_interpolation(buffer, NULL, output, NULL, width, height, components, u, v);
208 }
209
210 /* BILINEAR INTERPOLATION */
211 BLI_INLINE void bilinear_interpolation(const unsigned char *byte_buffer, const float *float_buffer,
212                                        unsigned char *byte_output, float *float_output, int width, int height,
213                                        int components, float u, float v, bool wrap_x, bool wrap_y)
214 {
215         float a, b;
216         float a_b, ma_b, a_mb, ma_mb;
217         int y1, y2, x1, x2;
218
219         /* ImBuf in must have a valid rect or rect_float, assume this is already checked */
220
221         x1 = (int)floor(u);
222         x2 = (int)ceil(u);
223         y1 = (int)floor(v);
224         y2 = (int)ceil(v);
225
226         if (float_output) {
227                 const float *row1, *row2, *row3, *row4;
228                 float empty[4] = {0.0f, 0.0f, 0.0f, 0.0f};
229
230                 /* pixel value must be already wrapped, however values at boundaries may flip */
231                 if (wrap_x) {
232                         if (x1 < 0) x1 = width  - 1;
233                         if (x2 >= width) x2 = 0;
234                 }
235                 else if (x2 < 0 || x1 >= width) {
236                         copy_vn_fl(float_output, components, 0.0f);
237                         return;
238                 }
239
240                 if (wrap_y) {
241                         if (y1 < 0) y1 = height - 1;
242                         if (y2 >= height) y2 = 0;
243                 }
244                 else if (y2 < 0 || y1 >= height) {
245                         copy_vn_fl(float_output, components, 0.0f);
246                         return;
247                 }
248
249                 /* sample including outside of edges of image */
250                 if (x1 < 0 || y1 < 0) row1 = empty;
251                 else row1 = float_buffer + width * y1 * components + components * x1;
252
253                 if (x1 < 0 || y2 > height - 1) row2 = empty;
254                 else row2 = float_buffer + width * y2 * components + components * x1;
255
256                 if (x2 > width - 1 || y1 < 0) row3 = empty;
257                 else row3 = float_buffer + width * y1 * components + components * x2;
258
259                 if (x2 > width - 1 || y2 > height - 1) row4 = empty;
260                 else row4 = float_buffer + width * y2 * components + components * x2;
261
262                 a = u - floorf(u);
263                 b = v - floorf(v);
264                 a_b = a * b; ma_b = (1.0f - a) * b; a_mb = a * (1.0f - b); ma_mb = (1.0f - a) * (1.0f - b);
265
266                 if (components == 1) {
267                         float_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
268                 }
269                 else if (components == 3) {
270                         float_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
271                         float_output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1];
272                         float_output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2];
273                 }
274                 else {
275                         float_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
276                         float_output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1];
277                         float_output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2];
278                         float_output[3] = ma_mb * row1[3] + a_mb * row3[3] + ma_b * row2[3] + a_b * row4[3];
279                 }
280         }
281         else {
282                 const unsigned char *row1, *row2, *row3, *row4;
283                 unsigned char empty[4] = {0, 0, 0, 0};
284
285                 /* pixel value must be already wrapped, however values at boundaries may flip */
286                 if (wrap_x) {
287                         if (x1 < 0) x1 = width  - 1;
288                         if (x2 >= width) x2 = 0;
289                 }
290                 else if (x2 < 0 || x1 >= width) {
291                         copy_vn_uchar(byte_output, components, 0);
292                         return;
293                 }
294
295                 if (wrap_y) {
296                         if (y1 < 0) y1 = height - 1;
297                         if (y2 >= height) y2 = 0;
298                 }
299                 else if (y2 < 0 || y1 >= height) {
300                         copy_vn_uchar(byte_output, components, 0);
301                         return;
302                 }
303
304                 /* sample including outside of edges of image */
305                 if (x1 < 0 || y1 < 0) row1 = empty;
306                 else row1 = byte_buffer + width * y1 * components + components * x1;
307
308                 if (x1 < 0 || y2 > height - 1) row2 = empty;
309                 else row2 = byte_buffer + width * y2 * components + components * x1;
310
311                 if (x2 > width - 1 || y1 < 0) row3 = empty;
312                 else row3 = byte_buffer + width * y1 * components + components * x2;
313
314                 if (x2 > width - 1 || y2 > height - 1) row4 = empty;
315                 else row4 = byte_buffer + width * y2 * components + components * x2;
316
317                 a = u - floorf(u);
318                 b = v - floorf(v);
319                 a_b = a * b; ma_b = (1.0f - a) * b; a_mb = a * (1.0f - b); ma_mb = (1.0f - a) * (1.0f - b);
320
321                 if (components == 1) {
322                         byte_output[0] = (unsigned char)(ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0] + 0.5f);
323                 }
324                 else if (components == 3) {
325                         byte_output[0] = (unsigned char)(ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0] + 0.5f);
326                         byte_output[1] = (unsigned char)(ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1] + 0.5f);
327                         byte_output[2] = (unsigned char)(ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2] + 0.5f);
328                 }
329                 else {
330                         byte_output[0] = (unsigned char)(ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0] + 0.5f);
331                         byte_output[1] = (unsigned char)(ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1] + 0.5f);
332                         byte_output[2] = (unsigned char)(ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2] + 0.5f);
333                         byte_output[3] = (unsigned char)(ma_mb * row1[3] + a_mb * row3[3] + ma_b * row2[3] + a_b * row4[3] + 0.5f);
334                 }
335         }
336 }
337
338 void BLI_bilinear_interpolation_fl(const float *buffer, float *output, int width, int height,
339                                    int components, float u, float v)
340 {
341         bilinear_interpolation(NULL, buffer, NULL, output, width, height, components, u, v, false, false);
342 }
343
344 void BLI_bilinear_interpolation_char(const unsigned char *buffer, unsigned char *output, int width, int height,
345                                      int components, float u, float v)
346 {
347         bilinear_interpolation(buffer, NULL, output, NULL, width, height, components, u, v, false, false);
348 }
349
350 void BLI_bilinear_interpolation_wrap_fl(const float *buffer, float *output, int width, int height,
351                                         int components, float u, float v, bool wrap_x, bool wrap_y)
352 {
353         bilinear_interpolation(NULL, buffer, NULL, output, width, height, components, u, v, wrap_x, wrap_y);
354 }
355
356 void BLI_bilinear_interpolation_wrap_char(const unsigned char *buffer, unsigned char *output, int width, int height,
357                                           int components, float u, float v, bool wrap_x, bool wrap_y)
358 {
359         bilinear_interpolation(buffer, NULL, output, NULL, width, height, components, u, v, wrap_x, wrap_y);
360 }
361
362
363 /**************************************************************************
364  * Filtering method based on
365  * "Creating raster omnimax images from multiple perspective views using the elliptical weighted average filter"
366  * by Ned Greene and Paul S. Heckbert (1986)
367  ***************************************************************************/
368
369 /* table of (exp(ar) - exp(a)) / (1 - exp(a)) for r in range [0, 1] and a = -2
370  * used instead of actual gaussian, otherwise at high texture magnifications circular artifacts are visible */
371 #define EWA_MAXIDX 255
372 const float EWA_WTS[EWA_MAXIDX + 1] = {
373         1.f, 0.990965f, 0.982f, 0.973105f, 0.96428f, 0.955524f, 0.946836f, 0.938216f, 0.929664f,
374         0.921178f, 0.912759f, 0.904405f, 0.896117f, 0.887893f, 0.879734f, 0.871638f, 0.863605f,
375         0.855636f, 0.847728f, 0.839883f, 0.832098f, 0.824375f, 0.816712f, 0.809108f, 0.801564f,
376         0.794079f, 0.786653f, 0.779284f, 0.771974f, 0.76472f, 0.757523f, 0.750382f, 0.743297f,
377         0.736267f, 0.729292f, 0.722372f, 0.715505f, 0.708693f, 0.701933f, 0.695227f, 0.688572f,
378         0.68197f, 0.67542f, 0.66892f, 0.662471f, 0.656073f, 0.649725f, 0.643426f, 0.637176f,
379         0.630976f, 0.624824f, 0.618719f, 0.612663f, 0.606654f, 0.600691f, 0.594776f, 0.588906f,
380         0.583083f, 0.577305f, 0.571572f, 0.565883f, 0.56024f, 0.55464f, 0.549084f, 0.543572f,
381         0.538102f, 0.532676f, 0.527291f, 0.521949f, 0.516649f, 0.511389f, 0.506171f, 0.500994f,
382         0.495857f, 0.490761f, 0.485704f, 0.480687f, 0.475709f, 0.470769f, 0.465869f, 0.461006f,
383         0.456182f, 0.451395f, 0.446646f, 0.441934f, 0.437258f, 0.432619f, 0.428017f, 0.42345f,
384         0.418919f, 0.414424f, 0.409963f, 0.405538f, 0.401147f, 0.39679f, 0.392467f, 0.388178f,
385         0.383923f, 0.379701f, 0.375511f, 0.371355f, 0.367231f, 0.363139f, 0.359079f, 0.355051f,
386         0.351055f, 0.347089f, 0.343155f, 0.339251f, 0.335378f, 0.331535f, 0.327722f, 0.323939f,
387         0.320186f, 0.316461f, 0.312766f, 0.3091f, 0.305462f, 0.301853f, 0.298272f, 0.294719f,
388         0.291194f, 0.287696f, 0.284226f, 0.280782f, 0.277366f, 0.273976f, 0.270613f, 0.267276f,
389         0.263965f, 0.26068f, 0.257421f, 0.254187f, 0.250979f, 0.247795f, 0.244636f, 0.241502f,
390         0.238393f, 0.235308f, 0.232246f, 0.229209f, 0.226196f, 0.223206f, 0.220239f, 0.217296f,
391         0.214375f, 0.211478f, 0.208603f, 0.20575f, 0.20292f, 0.200112f, 0.197326f, 0.194562f,
392         0.191819f, 0.189097f, 0.186397f, 0.183718f, 0.18106f, 0.178423f, 0.175806f, 0.17321f,
393         0.170634f, 0.168078f, 0.165542f, 0.163026f, 0.16053f, 0.158053f, 0.155595f, 0.153157f,
394         0.150738f, 0.148337f, 0.145955f, 0.143592f, 0.141248f, 0.138921f, 0.136613f, 0.134323f,
395         0.132051f, 0.129797f, 0.12756f, 0.125341f, 0.123139f, 0.120954f, 0.118786f, 0.116635f,
396         0.114501f, 0.112384f, 0.110283f, 0.108199f, 0.106131f, 0.104079f, 0.102043f, 0.100023f,
397         0.0980186f, 0.09603f, 0.094057f, 0.0920994f, 0.0901571f, 0.08823f, 0.0863179f, 0.0844208f,
398         0.0825384f, 0.0806708f, 0.0788178f, 0.0769792f, 0.0751551f, 0.0733451f, 0.0715493f, 0.0697676f,
399         0.0679997f, 0.0662457f, 0.0645054f, 0.0627786f, 0.0610654f, 0.0593655f, 0.0576789f, 0.0560055f,
400         0.0543452f, 0.0526979f, 0.0510634f, 0.0494416f, 0.0478326f, 0.0462361f, 0.0446521f, 0.0430805f,
401         0.0415211f, 0.039974f, 0.0384389f, 0.0369158f, 0.0354046f, 0.0339052f, 0.0324175f, 0.0309415f,
402         0.029477f, 0.0280239f, 0.0265822f, 0.0251517f, 0.0237324f, 0.0223242f, 0.020927f, 0.0195408f,
403         0.0181653f, 0.0168006f, 0.0154466f, 0.0141031f, 0.0127701f, 0.0114476f, 0.0101354f, 0.00883339f,
404         0.00754159f, 0.00625989f, 0.00498819f, 0.00372644f, 0.00247454f, 0.00123242f, 0.f
405 };
406
407 static void radangle2imp(float a2, float b2, float th, float *A, float *B, float *C, float *F)
408 {
409         float ct2 = cosf(th);
410         const float st2 = 1.0f - ct2 * ct2;     /* <- sin(th)^2 */
411         ct2 *= ct2;
412         *A = a2 * st2 + b2 * ct2;
413         *B = (b2 - a2) * sinf(2.0f * th);
414         *C = a2 * ct2 + b2 * st2;
415         *F = a2 * b2;
416 }
417
418 /* all tests here are done to make sure possible overflows are hopefully minimized */
419 void BLI_ewa_imp2radangle(float A, float B, float C, float F, float *a, float *b, float *th, float *ecc)
420 {
421         if (F <= 1e-5f) {       /* use arbitrary major radius, zero minor, infinite eccentricity */
422                 *a = sqrtf(A > C ? A : C);
423                 *b = 0.0f;
424                 *ecc = 1e10f;
425                 *th = 0.5f * (atan2f(B, A - C) + (float)M_PI);
426         }
427         else {
428                 const float AmC = A - C, ApC = A + C, F2 = F * 2.0f;
429                 const float r = sqrtf(AmC * AmC + B * B);
430                 float d = ApC - r;
431                 *a = (d <= 0.0f) ? sqrtf(A > C ? A : C) : sqrtf(F2 / d);
432                 d = ApC + r;
433                 if (d <= 0.0f) {
434                         *b = 0.0f;
435                         *ecc = 1e10f;
436                 }
437                 else {
438                         *b = sqrtf(F2 / d);
439                         *ecc = *a / *b;
440                 }
441                 /* incr theta by 0.5*pi (angle of major axis) */
442                 *th = 0.5f * (atan2f(B, AmC) + (float)M_PI);
443         }
444 }
445
446 void BLI_ewa_filter(const int width, const int height,
447                     const bool intpol,
448                     const bool use_alpha,
449                     const float uv[2],
450                     const float du[2],
451                     const float dv[2],
452                     ewa_filter_read_pixel_cb read_pixel_cb,
453                     void *userdata,
454                     float result[4])
455 {
456         /* scaling dxt/dyt by full resolution can cause overflow because of huge A/B/C and esp. F values,
457          * scaling by aspect ratio alone does the opposite, so try something in between instead... */
458         const float ff2 = (float)width, ff = sqrtf(ff2), q = (float)height / ff;
459         const float Ux = du[0] * ff, Vx = du[1] * q, Uy = dv[0] * ff, Vy = dv[1] * q;
460         float A = Vx * Vx + Vy * Vy;
461         float B = -2.0f * (Ux * Vx + Uy * Vy);
462         float C = Ux * Ux + Uy * Uy;
463         float F = A * C - B * B * 0.25f;
464         float a, b, th, ecc, a2, b2, ue, ve, U0, V0, DDQ, U, ac1, ac2, BU, d;
465         int u, v, u1, u2, v1, v2;
466
467         /* The so-called 'high' quality ewa method simply adds a constant of 1 to both A & C,
468          * so the ellipse always covers at least some texels. But since the filter is now always larger,
469          * it also means that everywhere else it's also more blurry then ideally should be the case.
470          * So instead here the ellipse radii are modified instead whenever either is too low.
471          * Use a different radius based on interpolation switch, just enough to anti-alias when interpolation is off,
472          * and slightly larger to make result a bit smoother than bilinear interpolation when interpolation is on
473          * (minimum values: const float rmin = intpol ? 1.f : 0.5f;) */
474         const float rmin = (intpol ? 1.5625f : 0.765625f) / ff2;
475         BLI_ewa_imp2radangle(A, B, C, F, &a, &b, &th, &ecc);
476         if ((b2 = b * b) < rmin) {
477                 if ((a2 = a * a) < rmin) {
478                         B = 0.0f;
479                         A = C = rmin;
480                         F = A * C;
481                 }
482                 else {
483                         b2 = rmin;
484                         radangle2imp(a2, b2, th, &A, &B, &C, &F);
485                 }
486         }
487
488         ue = ff * sqrtf(C);
489         ve = ff * sqrtf(A);
490         d = (float)(EWA_MAXIDX + 1) / (F * ff2);
491         A *= d;
492         B *= d;
493         C *= d;
494
495         U0 = uv[0] * (float)width;
496         V0 = uv[1] * (float)height;
497         u1 = (int)(floorf(U0 - ue));
498         u2 = (int)(ceilf(U0 + ue));
499         v1 = (int)(floorf(V0 - ve));
500         v2 = (int)(ceilf(V0 + ve));
501
502         /* sane clamping to avoid unnecessarily huge loops */
503         /* note: if eccentricity gets clamped (see above),
504          * the ue/ve limits can also be lowered accordingly
505          */
506         if (U0 - (float)u1 > EWA_MAXIDX) u1 = (int)U0 - EWA_MAXIDX;
507         if ((float)u2 - U0 > EWA_MAXIDX) u2 = (int)U0 + EWA_MAXIDX;
508         if (V0 - (float)v1 > EWA_MAXIDX) v1 = (int)V0 - EWA_MAXIDX;
509         if ((float)v2 - V0 > EWA_MAXIDX) v2 = (int)V0 + EWA_MAXIDX;
510
511         /* Early output check for cases the whole region is outside of the buffer. */
512         if ((u2 < 0 || u1 >= width) ||  (v2 < 0 || v1 >= height)) {
513                 zero_v4(result);
514                 return;
515         }
516
517         U0 -= 0.5f;
518         V0 -= 0.5f;
519         DDQ = 2.0f * A;
520         U = (float)u1 - U0;
521         ac1 = A * (2.0f * U + 1.0f);
522         ac2 = A * U * U;
523         BU = B * U;
524
525         d = 0.0f;
526         zero_v4(result);
527         for (v = v1; v <= v2; ++v) {
528                 const float V = (float)v - V0;
529                 float DQ = ac1 + B * V;
530                 float Q = (C * V + BU) * V + ac2;
531                 for (u = u1; u <= u2; ++u) {
532                         if (Q < (float)(EWA_MAXIDX + 1)) {
533                                 float tc[4];
534                                 const float wt = EWA_WTS[(Q < 0.0f) ? 0 : (unsigned int)Q];
535                                 read_pixel_cb(userdata, u, v, tc);
536                                 madd_v3_v3fl(result, tc, wt);
537                                 result[3] += use_alpha ? tc[3] * wt : 0.0f;
538                                 d += wt;
539                         }
540                         Q += DQ;
541                         DQ += DDQ;
542                 }
543         }
544
545         /* d should hopefully never be zero anymore */
546         d = 1.0f / d;
547         mul_v3_fl(result, d);
548         /* clipping can be ignored if alpha used, texr->ta already includes filtered edge */
549         result[3] = use_alpha ? result[3] * d : 1.0f;
550 }