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