Partial revert '#if 0' cleanup
[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                         copy_vn_fl(float_output, components, 0.0f);
116                 }
117                 if (byte_output) {
118                         copy_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, bool wrap_x, bool wrap_y)
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                 /* pixel value must be already wrapped, however values at boundaries may flip */
283                 if (wrap_x) {
284                         if (x1 < 0) x1 = width  - 1;
285                         if (x2 >= width) x2 = 0;
286                 }
287                 else if (x2 < 0 || x1 >= width) {
288                         copy_vn_fl(float_output, components, 0.0f);
289                         return;
290                 }
291
292                 if (wrap_y) {
293                         if (y1 < 0) y1 = height - 1;
294                         if (y2 >= height) y2 = 0;
295                 }
296                 else if (y2 < 0 || y1 >= height) {
297                         copy_vn_fl(float_output, components, 0.0f);
298                         return;
299                 }
300
301                 /* sample including outside of edges of image */
302                 if (x1 < 0 || y1 < 0) row1 = empty;
303                 else row1 = float_buffer + width * y1 * components + components * x1;
304
305                 if (x1 < 0 || y2 > height - 1) row2 = empty;
306                 else row2 = float_buffer + width * y2 * components + components * x1;
307
308                 if (x2 > width - 1 || y1 < 0) row3 = empty;
309                 else row3 = float_buffer + width * y1 * components + components * x2;
310
311                 if (x2 > width - 1 || y2 > height - 1) row4 = empty;
312                 else row4 = float_buffer + width * y2 * components + components * x2;
313
314                 a = u - floorf(u);
315                 b = v - floorf(v);
316                 a_b = a * b; ma_b = (1.0f - a) * b; a_mb = a * (1.0f - b); ma_mb = (1.0f - a) * (1.0f - b);
317
318                 if (components == 1) {
319                         float_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
320                 }
321                 else if (components == 3) {
322                         float_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
323                         float_output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1];
324                         float_output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2];
325                 }
326                 else {
327                         float_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
328                         float_output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1];
329                         float_output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2];
330                         float_output[3] = ma_mb * row1[3] + a_mb * row3[3] + ma_b * row2[3] + a_b * row4[3];
331                 }
332         }
333         else {
334                 const unsigned char *row1, *row2, *row3, *row4;
335                 unsigned char empty[4] = {0, 0, 0, 0};
336
337                 /* pixel value must be already wrapped, however values at boundaries may flip */
338                 if (wrap_x) {
339                         if (x1 < 0) x1 = width  - 1;
340                         if (x2 >= width) x2 = 0;
341                 }
342                 else if (x2 < 0 || x1 >= width) {
343                         copy_vn_uchar(byte_output, components, 0);
344                         return;
345                 }
346
347                 if (wrap_y) {
348                         if (y1 < 0) y1 = height - 1;
349                         if (y2 >= height) y2 = 0;
350                 }
351                 else if (y2 < 0 || y1 >= height) {
352                         copy_vn_uchar(byte_output, components, 0);
353                         return;
354                 }
355
356                 /* sample including outside of edges of image */
357                 if (x1 < 0 || y1 < 0) row1 = empty;
358                 else row1 = byte_buffer + width * y1 * components + components * x1;
359
360                 if (x1 < 0 || y2 > height - 1) row2 = empty;
361                 else row2 = byte_buffer + width * y2 * components + components * x1;
362
363                 if (x2 > width - 1 || y1 < 0) row3 = empty;
364                 else row3 = byte_buffer + width * y1 * components + components * x2;
365
366                 if (x2 > width - 1 || y2 > height - 1) row4 = empty;
367                 else row4 = byte_buffer + width * y2 * components + components * x2;
368
369                 a = u - floorf(u);
370                 b = v - floorf(v);
371                 a_b = a * b; ma_b = (1.0f - a) * b; a_mb = a * (1.0f - b); ma_mb = (1.0f - a) * (1.0f - b);
372
373                 if (components == 1) {
374                         byte_output[0] = (unsigned char)(ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0] + 0.5f);
375                 }
376                 else if (components == 3) {
377                         byte_output[0] = (unsigned char)(ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0] + 0.5f);
378                         byte_output[1] = (unsigned char)(ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1] + 0.5f);
379                         byte_output[2] = (unsigned char)(ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2] + 0.5f);
380                 }
381                 else {
382                         byte_output[0] = (unsigned char)(ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0] + 0.5f);
383                         byte_output[1] = (unsigned char)(ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1] + 0.5f);
384                         byte_output[2] = (unsigned char)(ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2] + 0.5f);
385                         byte_output[3] = (unsigned char)(ma_mb * row1[3] + a_mb * row3[3] + ma_b * row2[3] + a_b * row4[3] + 0.5f);
386                 }
387         }
388 }
389
390 void BLI_bilinear_interpolation_fl(const float *buffer, float *output, int width, int height,
391                                    int components, float u, float v)
392 {
393         bilinear_interpolation(NULL, buffer, NULL, output, width, height, components, u, v, false, false);
394 }
395
396 void BLI_bilinear_interpolation_char(const unsigned char *buffer, unsigned char *output, int width, int height,
397                                      int components, float u, float v)
398 {
399         bilinear_interpolation(buffer, NULL, output, NULL, width, height, components, u, v, false, false);
400 }
401
402 void BLI_bilinear_interpolation_wrap_fl(const float *buffer, float *output, int width, int height,
403                                         int components, float u, float v, bool wrap_x, bool wrap_y)
404 {
405         bilinear_interpolation(NULL, buffer, NULL, output, width, height, components, u, v, wrap_x, wrap_y);
406 }
407
408 void BLI_bilinear_interpolation_wrap_char(const unsigned char *buffer, unsigned char *output, int width, int height,
409                                           int components, float u, float v, bool wrap_x, bool wrap_y)
410 {
411         bilinear_interpolation(buffer, NULL, output, NULL, width, height, components, u, v, wrap_x, wrap_y);
412 }
413
414
415 /**************************************************************************
416  * Filtering method based on
417  * "Creating raster omnimax images from multiple perspective views using the elliptical weighted average filter"
418  * by Ned Greene and Paul S. Heckbert (1986)
419  ***************************************************************************/
420
421 /* table of (exp(ar) - exp(a)) / (1 - exp(a)) for r in range [0, 1] and a = -2
422  * used instead of actual gaussian, otherwise at high texture magnifications circular artifacts are visible */
423 #define EWA_MAXIDX 255
424 const float EWA_WTS[EWA_MAXIDX + 1] = {
425         1.f, 0.990965f, 0.982f, 0.973105f, 0.96428f, 0.955524f, 0.946836f, 0.938216f, 0.929664f,
426         0.921178f, 0.912759f, 0.904405f, 0.896117f, 0.887893f, 0.879734f, 0.871638f, 0.863605f,
427         0.855636f, 0.847728f, 0.839883f, 0.832098f, 0.824375f, 0.816712f, 0.809108f, 0.801564f,
428         0.794079f, 0.786653f, 0.779284f, 0.771974f, 0.76472f, 0.757523f, 0.750382f, 0.743297f,
429         0.736267f, 0.729292f, 0.722372f, 0.715505f, 0.708693f, 0.701933f, 0.695227f, 0.688572f,
430         0.68197f, 0.67542f, 0.66892f, 0.662471f, 0.656073f, 0.649725f, 0.643426f, 0.637176f,
431         0.630976f, 0.624824f, 0.618719f, 0.612663f, 0.606654f, 0.600691f, 0.594776f, 0.588906f,
432         0.583083f, 0.577305f, 0.571572f, 0.565883f, 0.56024f, 0.55464f, 0.549084f, 0.543572f,
433         0.538102f, 0.532676f, 0.527291f, 0.521949f, 0.516649f, 0.511389f, 0.506171f, 0.500994f,
434         0.495857f, 0.490761f, 0.485704f, 0.480687f, 0.475709f, 0.470769f, 0.465869f, 0.461006f,
435         0.456182f, 0.451395f, 0.446646f, 0.441934f, 0.437258f, 0.432619f, 0.428017f, 0.42345f,
436         0.418919f, 0.414424f, 0.409963f, 0.405538f, 0.401147f, 0.39679f, 0.392467f, 0.388178f,
437         0.383923f, 0.379701f, 0.375511f, 0.371355f, 0.367231f, 0.363139f, 0.359079f, 0.355051f,
438         0.351055f, 0.347089f, 0.343155f, 0.339251f, 0.335378f, 0.331535f, 0.327722f, 0.323939f,
439         0.320186f, 0.316461f, 0.312766f, 0.3091f, 0.305462f, 0.301853f, 0.298272f, 0.294719f,
440         0.291194f, 0.287696f, 0.284226f, 0.280782f, 0.277366f, 0.273976f, 0.270613f, 0.267276f,
441         0.263965f, 0.26068f, 0.257421f, 0.254187f, 0.250979f, 0.247795f, 0.244636f, 0.241502f,
442         0.238393f, 0.235308f, 0.232246f, 0.229209f, 0.226196f, 0.223206f, 0.220239f, 0.217296f,
443         0.214375f, 0.211478f, 0.208603f, 0.20575f, 0.20292f, 0.200112f, 0.197326f, 0.194562f,
444         0.191819f, 0.189097f, 0.186397f, 0.183718f, 0.18106f, 0.178423f, 0.175806f, 0.17321f,
445         0.170634f, 0.168078f, 0.165542f, 0.163026f, 0.16053f, 0.158053f, 0.155595f, 0.153157f,
446         0.150738f, 0.148337f, 0.145955f, 0.143592f, 0.141248f, 0.138921f, 0.136613f, 0.134323f,
447         0.132051f, 0.129797f, 0.12756f, 0.125341f, 0.123139f, 0.120954f, 0.118786f, 0.116635f,
448         0.114501f, 0.112384f, 0.110283f, 0.108199f, 0.106131f, 0.104079f, 0.102043f, 0.100023f,
449         0.0980186f, 0.09603f, 0.094057f, 0.0920994f, 0.0901571f, 0.08823f, 0.0863179f, 0.0844208f,
450         0.0825384f, 0.0806708f, 0.0788178f, 0.0769792f, 0.0751551f, 0.0733451f, 0.0715493f, 0.0697676f,
451         0.0679997f, 0.0662457f, 0.0645054f, 0.0627786f, 0.0610654f, 0.0593655f, 0.0576789f, 0.0560055f,
452         0.0543452f, 0.0526979f, 0.0510634f, 0.0494416f, 0.0478326f, 0.0462361f, 0.0446521f, 0.0430805f,
453         0.0415211f, 0.039974f, 0.0384389f, 0.0369158f, 0.0354046f, 0.0339052f, 0.0324175f, 0.0309415f,
454         0.029477f, 0.0280239f, 0.0265822f, 0.0251517f, 0.0237324f, 0.0223242f, 0.020927f, 0.0195408f,
455         0.0181653f, 0.0168006f, 0.0154466f, 0.0141031f, 0.0127701f, 0.0114476f, 0.0101354f, 0.00883339f,
456         0.00754159f, 0.00625989f, 0.00498819f, 0.00372644f, 0.00247454f, 0.00123242f, 0.f
457 };
458
459 static void radangle2imp(float a2, float b2, float th, float *A, float *B, float *C, float *F)
460 {
461         float ct2 = cosf(th);
462         const float st2 = 1.0f - ct2 * ct2;     /* <- sin(th)^2 */
463         ct2 *= ct2;
464         *A = a2 * st2 + b2 * ct2;
465         *B = (b2 - a2) * sinf(2.0f * th);
466         *C = a2 * ct2 + b2 * st2;
467         *F = a2 * b2;
468 }
469
470 /* all tests here are done to make sure possible overflows are hopefully minimized */
471 void BLI_ewa_imp2radangle(float A, float B, float C, float F, float *a, float *b, float *th, float *ecc)
472 {
473         if (F <= 1e-5f) {       /* use arbitrary major radius, zero minor, infinite eccentricity */
474                 *a = sqrtf(A > C ? A : C);
475                 *b = 0.0f;
476                 *ecc = 1e10f;
477                 *th = 0.5f * (atan2f(B, A - C) + (float)M_PI);
478         }
479         else {
480                 const float AmC = A - C, ApC = A + C, F2 = F * 2.0f;
481                 const float r = sqrtf(AmC * AmC + B * B);
482                 float d = ApC - r;
483                 *a = (d <= 0.0f) ? sqrtf(A > C ? A : C) : sqrtf(F2 / d);
484                 d = ApC + r;
485                 if (d <= 0.0f) {
486                         *b = 0.0f;
487                         *ecc = 1e10f;
488                 }
489                 else {
490                         *b = sqrtf(F2 / d);
491                         *ecc = *a / *b;
492                 }
493                 /* incr theta by 0.5*pi (angle of major axis) */
494                 *th = 0.5f * (atan2f(B, AmC) + (float)M_PI);
495         }
496 }
497
498 void BLI_ewa_filter(const int width, const int height,
499                     const bool intpol,
500                     const bool use_alpha,
501                     const float uv[2],
502                     const float du[2],
503                     const float dv[2],
504                     ewa_filter_read_pixel_cb read_pixel_cb,
505                     void *userdata,
506                     float result[4])
507 {
508         /* scaling dxt/dyt by full resolution can cause overflow because of huge A/B/C and esp. F values,
509          * scaling by aspect ratio alone does the opposite, so try something in between instead... */
510         const float ff2 = (float)width, ff = sqrtf(ff2), q = (float)height / ff;
511         const float Ux = du[0] * ff, Vx = du[1] * q, Uy = dv[0] * ff, Vy = dv[1] * q;
512         float A = Vx * Vx + Vy * Vy;
513         float B = -2.0f * (Ux * Vx + Uy * Vy);
514         float C = Ux * Ux + Uy * Uy;
515         float F = A * C - B * B * 0.25f;
516         float a, b, th, ecc, a2, b2, ue, ve, U0, V0, DDQ, U, ac1, ac2, BU, d;
517         int u, v, u1, u2, v1, v2;
518
519         /* The so-called 'high' quality ewa method simply adds a constant of 1 to both A & C,
520          * so the ellipse always covers at least some texels. But since the filter is now always larger,
521          * it also means that everywhere else it's also more blurry then ideally should be the case.
522          * So instead here the ellipse radii are modified instead whenever either is too low.
523          * Use a different radius based on interpolation switch, just enough to anti-alias when interpolation is off,
524          * and slightly larger to make result a bit smoother than bilinear interpolation when interpolation is on
525          * (minimum values: const float rmin = intpol ? 1.f : 0.5f;) */
526         const float rmin = (intpol ? 1.5625f : 0.765625f) / ff2;
527         BLI_ewa_imp2radangle(A, B, C, F, &a, &b, &th, &ecc);
528         if ((b2 = b * b) < rmin) {
529                 if ((a2 = a * a) < rmin) {
530                         B = 0.0f;
531                         A = C = rmin;
532                         F = A * C;
533                 }
534                 else {
535                         b2 = rmin;
536                         radangle2imp(a2, b2, th, &A, &B, &C, &F);
537                 }
538         }
539
540         ue = ff * sqrtf(C);
541         ve = ff * sqrtf(A);
542         d = (float)(EWA_MAXIDX + 1) / (F * ff2);
543         A *= d;
544         B *= d;
545         C *= d;
546
547         U0 = uv[0] * (float)width;
548         V0 = uv[1] * (float)height;
549         u1 = (int)(floorf(U0 - ue));
550         u2 = (int)(ceilf(U0 + ue));
551         v1 = (int)(floorf(V0 - ve));
552         v2 = (int)(ceilf(V0 + ve));
553
554         /* sane clamping to avoid unnecessarily huge loops */
555         /* note: if eccentricity gets clamped (see above),
556          * the ue/ve limits can also be lowered accordingly
557          */
558         if (U0 - (float)u1 > EWA_MAXIDX) u1 = (int)U0 - EWA_MAXIDX;
559         if ((float)u2 - U0 > EWA_MAXIDX) u2 = (int)U0 + EWA_MAXIDX;
560         if (V0 - (float)v1 > EWA_MAXIDX) v1 = (int)V0 - EWA_MAXIDX;
561         if ((float)v2 - V0 > EWA_MAXIDX) v2 = (int)V0 + EWA_MAXIDX;
562
563         /* Early output check for cases the whole region is outside of the buffer. */
564         if ((u2 < 0 || u1 >= width) ||  (v2 < 0 || v1 >= height)) {
565                 zero_v4(result);
566                 return;
567         }
568
569         U0 -= 0.5f;
570         V0 -= 0.5f;
571         DDQ = 2.0f * A;
572         U = (float)u1 - U0;
573         ac1 = A * (2.0f * U + 1.0f);
574         ac2 = A * U * U;
575         BU = B * U;
576
577         d = 0.0f;
578         zero_v4(result);
579         for (v = v1; v <= v2; ++v) {
580                 const float V = (float)v - V0;
581                 float DQ = ac1 + B * V;
582                 float Q = (C * V + BU) * V + ac2;
583                 for (u = u1; u <= u2; ++u) {
584                         if (Q < (float)(EWA_MAXIDX + 1)) {
585                                 float tc[4];
586                                 const float wt = EWA_WTS[(Q < 0.0f) ? 0 : (unsigned int)Q];
587                                 read_pixel_cb(userdata, u, v, tc);
588                                 madd_v3_v3fl(result, tc, wt);
589                                 result[3] += use_alpha ? tc[3] * wt : 0.0f;
590                                 d += wt;
591                         }
592                         Q += DQ;
593                         DQ += DDQ;
594                 }
595         }
596
597         /* d should hopefully never be zero anymore */
598         d = 1.0f / d;
599         mul_v3_fl(result, d);
600         /* clipping can be ignored if alpha used, texr->ta already includes filtered edge */
601         result[3] = use_alpha ? result[3] * d : 1.0f;
602 }