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