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