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