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