1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; -*- */
2 /* ***** BEGIN LICENSE BLOCK *****
3  * This file is part of openfx-supportext <https://github.com/devernay/openfx-supportext>,
4  * Copyright (C) 2013-2018 INRIA
5  *
6  * openfx-supportext is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * openfx-supportext is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with openfx-supportext.  If not, see <http://www.gnu.org/licenses/gpl-2.0.html>
18  * ***** END LICENSE BLOCK ***** */
19 
20 /*
21  * OFX Filter/Interpolation help functions
22  */
23 
24 #ifndef openfx_supportext_ofxsFilter_h
25 #define openfx_supportext_ofxsFilter_h
26 
27 #include <cmath>
28 #include <cassert>
29 #include <algorithm>
30 
31 #include "ofxsImageEffect.h"
32 
33 namespace OFX {
34 // GENERIC
35 #define kParamFilterType "filter"
36 #define kParamFilterTypeLabel "Filter"
37 #define kParamFilterTypeHint "Filtering algorithm - some filters may produce values outside of the initial range (*) or modify the values even if there is no movement (+)."
38 #define kParamFilterClamp "clamp"
39 #define kParamFilterClampLabel "Clamp"
40 #define kParamFilterClampHint "Clamp filter output within the original range - useful to avoid negative values in mattes"
41 #define kParamFilterBlackOutside "black_outside"
42 #define kParamFilterBlackOutsideLabel "Black outside"
43 #define kParamFilterBlackOutsideHint "Fill the area outside the source image with black"
44 
45 enum FilterEnum
46 {
47     eFilterImpulse,
48     eFilterBox,
49     eFilterBilinear,
50     eFilterCubic,
51     eFilterKeys,
52     eFilterSimon,
53     eFilterRifman,
54     eFilterMitchell,
55     eFilterParzen,
56     eFilterNotch,
57 };
58 
59 #define kFilterImpulse "Impulse", "(nearest neighbor / box) Use original values.", "impulse"
60 #define kFilterBox "Box", "Integrate the source image over the bounding box of the back-transformed pixel.", "box"
61 #define kFilterBilinear "Bilinear", "(tent / triangle) Bilinear interpolation between original values.", "bilinear"
62 #define kFilterCubic "Cubic", "(cubic spline) Some smoothing.", "cubic"
63 #define kFilterKeys "Keys", "(Catmull-Rom / Hermite spline) Some smoothing, plus minor sharpening (*).", "keys"
64 #define kFilterSimon "Simon", "Some smoothing, plus medium sharpening (*).", "simon"
65 #define kFilterRifman "Rifman", "Some smoothing, plus significant sharpening (*).", "rifman"
66 #define kFilterMitchell "Mitchell", "Some smoothing, plus blurring to hide pixelation (*+).", "mitchell"
67 #define kFilterParzen "Parzen", "(cubic B-spline) Greatest smoothing of all filters (+).", "parzen"
68 #define kFilterNotch "Notch", "Flat smoothing (which tends to hide moire' patterns) (+).", "notch"
69 
70 inline
71 void
72 ofxsFilterDescribeParamsInterpolate2D(OFX::ImageEffectDescriptor &desc,
73                                       OFX::PageParamDescriptor *page,
74                                       bool blackOutsideDefault = true)
75 {
76     // GENERIC PARAMETERS
77     //
78     {
79         OFX::ChoiceParamDescriptor* param = desc.defineChoiceParam(kParamFilterType);
80 
81         param->setLabel(kParamFilterTypeLabel);
82         param->setHint(kParamFilterTypeHint);
83         assert(param->getNOptions() == eFilterImpulse);
84         param->appendOption(kFilterImpulse);
85         assert(param->getNOptions() == eFilterBox);
86         param->appendOption(kFilterBox);
87         assert(param->getNOptions() == eFilterBilinear);
88         param->appendOption(kFilterBilinear);
89         assert(param->getNOptions() == eFilterCubic);
90         param->appendOption(kFilterCubic);
91         assert(param->getNOptions() == eFilterKeys);
92         param->appendOption(kFilterKeys);
93         assert(param->getNOptions() == eFilterSimon);
94         param->appendOption(kFilterSimon);
95         assert(param->getNOptions() == eFilterRifman);
96         param->appendOption(kFilterRifman);
97         assert(param->getNOptions() == eFilterMitchell);
98         param->appendOption(kFilterMitchell);
99         assert(param->getNOptions() == eFilterParzen);
100         param->appendOption(kFilterParzen);
101         assert(param->getNOptions() == eFilterNotch);
102         param->appendOption(kFilterNotch);
103         param->setDefault(eFilterCubic);
104         param->setAnimates(true);
105 #ifdef OFX_EXTENSIONS_NUKE
106         param->setLayoutHint(OFX::eLayoutHintNoNewLine, 1);
107 #endif
108         if (page) {
109             page->addChild(*param);
110         }
111     }
112 
113     // clamp
114     {
115         OFX::BooleanParamDescriptor* param = desc.defineBooleanParam(kParamFilterClamp);
116         param->setLabel(kParamFilterClampLabel);
117         param->setHint(kParamFilterClampHint);
118         param->setDefault(false);
119         param->setAnimates(true);
120 #ifdef OFX_EXTENSIONS_NUKE
121         param->setLayoutHint(OFX::eLayoutHintNoNewLine, 1);
122 #endif
123         if (page) {
124             page->addChild(*param);
125         }
126     }
127 
128     // blackOutside
129     {
130         OFX::BooleanParamDescriptor* param = desc.defineBooleanParam(kParamFilterBlackOutside);
131         param->setLabel(kParamFilterBlackOutsideLabel);
132         param->setHint(kParamFilterBlackOutsideHint);
133         param->setDefault(blackOutsideDefault);
134         param->setAnimates(true);
135         if (page) {
136             page->addChild(*param);
137         }
138     }
139 } // ofxsFilterDescribeParamsInterpolate2D
140 
141 /*
142    Maple code to compute the filters.
143 
144  # Mitchell, D. and A. Netravali, "Reconstruction Filters in Computer Graphics."
145  # http://www.cs.utexas.edu/users/fussell/courses/cs384g/lectures/mitchell/Mitchell.pdf
146  # Computer Graphics, Vol. 22, No. 4, pp. 221-228.
147  # (B, C)
148  # (1/3, 1/3) - Defaults recommended by Mitchell and Netravali
149  # (1, 0) - Equivalent to the Cubic B-Spline
150  # (0, 0.5) - Equivalent to the Catmull-Rom Spline
151  # (0, C) - The family of Cardinal Cubic Splines
152  # (B, 0) - Duff's tensioned B-Splines.
153    unassign('Ip'):unassign('Ic'):unassign('In'):unassign('Ia'):
154    unassign('Jp'):unassign('Jc'):unassign('Jn'):unassign('Ja'):
155    P:= x -> ((12-9*B-6*C)*x**3 + (-18+12*B+6*C)*x**2+(6-2*B))/6;
156    Q:= x -> ((-B-6*C)*x**3 + (6*B+30*C)*x**2 + (-12*B-48*C)*x + (8*B+24*C))/6;
157 
158    R := d -> Q(d+1)*Ip + P(d)*Ic + P(1-d) * In + Q(2-d)*Ia;
159 
160  # how does it perform on a linear function?
161    R0 :=  d -> Q(d+1)*(Ic-1) + P(d)*Ic + P(1-d) * (Ic+1) + Q(2-d)*(Ic+2);
162 
163  # Cubic (cubic splines - depends only on In and Ic, derivatives are 0 at the center of each sample)
164    collect(subs({B=0,C=0},R(d)),d);
165    collect(subs({B=0,C=0},R0(d)),d);
166 
167  # Catmull-Rom / Keys / Hermite spline - gives linear func if input is linear
168    collect(subs({B=0,C=0.5},R(d)),d);
169    collect(subs({B=0,C=0.5},R0(d)),d);
170 
171  # Simon
172    collect(subs({B=0,C=0.75},R(d)),d);
173    collect(subs({B=0,C=0.75},R0(d)),d);
174 
175  # Rifman
176    collect(subs({B=0,C=1.},R(d)),d);
177    collect(subs({B=0,C=1.},R0(d)),d);
178 
179  # Mitchell - gives linear func if input is linear
180    collect(subs({B=1/3, C=1/3},R(d)),d);
181    collect(subs({B=1/3, C=1/3},R0(d)),d);
182 
183  # Parzen (Cubic B-spline) - gives linear func if input is linear
184    collect(subs({B=1,C=0},R(d)),d);
185    collect(subs({B=1,C=0},R0(d)),d);
186 
187  # Notch - gives linear func if input is linear
188    collect(subs({B=3/2,C=-1/4},R(d)),d);
189    collect(subs({B=3/2,C=-1/4},R0(d)),d);
190  */
191 inline double
ofxsFilterClampVal(double I,double Ic,double In)192 ofxsFilterClampVal(double I,
193                    double Ic,
194                    double In)
195 {
196     double Imin = (std::min)(Ic, In);
197 
198     if (I < Imin) {
199         return Imin;
200     }
201     double Imax = (std::max)(Ic, In);
202     if (I > Imax) {
203         return Imax;
204     }
205 
206     return I;
207 }
208 
209 inline
210 double
ofxsFilterLinear(double Ic,double In,double d)211 ofxsFilterLinear(double Ic,
212                  double In,
213                  double d)
214 {
215     return Ic + d * (In - Ic);
216 }
217 
218 static inline
219 double
ofxsFilterCubic(double Ic,double In,double d,bool clamp)220 ofxsFilterCubic(double Ic,
221                 double In,
222                 double d,
223                 bool clamp)
224 {
225     double I = Ic + d * d * ( (-3 * Ic + 3 * In ) + d * (2 * Ic - 2 * In ) );
226 
227     if (clamp) {
228         I = ofxsFilterClampVal(I, Ic, In);
229     }
230 
231     return I;
232 }
233 
234 inline
235 double
ofxsFilterKeys(double Ip,double Ic,double In,double Ia,double d,bool clamp)236 ofxsFilterKeys(double Ip,
237                double Ic,
238                double In,
239                double Ia,
240                double d,
241                bool clamp)
242 {
243     double I = Ic  + d * ( (-Ip + In ) + d * ( (2 * Ip - 5 * Ic + 4 * In - Ia ) + d * (-Ip + 3 * Ic - 3 * In + Ia ) ) ) / 2;
244 
245     if (clamp) {
246         I = ofxsFilterClampVal(I, Ic, In);
247     }
248 
249     return I;
250 }
251 
252 inline
253 double
ofxsFilterSimon(double Ip,double Ic,double In,double Ia,double d,bool clamp)254 ofxsFilterSimon(double Ip,
255                 double Ic,
256                 double In,
257                 double Ia,
258                 double d,
259                 bool clamp)
260 {
261     double I = Ic  + d * ( (-3 * Ip + 3 * In ) + d * ( (6 * Ip - 9 * Ic + 6 * In - 3 * Ia ) + d * (-3 * Ip + 5 * Ic - 5 * In + 3 * Ia ) ) ) / 4;
262 
263     if (clamp) {
264         I = ofxsFilterClampVal(I, Ic, In);
265     }
266 
267     return I;
268 }
269 
270 inline
271 double
ofxsFilterRifman(double Ip,double Ic,double In,double Ia,double d,bool clamp)272 ofxsFilterRifman(double Ip,
273                  double Ic,
274                  double In,
275                  double Ia,
276                  double d,
277                  bool clamp)
278 {
279     double I = Ic  + d * ( (-Ip + In ) + d * ( (2 * Ip - 2 * Ic + In - Ia ) + d * (-Ip + Ic - In + Ia ) ) );
280 
281     if (clamp) {
282         I = ofxsFilterClampVal(I, Ic, In);
283     }
284 
285     return I;
286 }
287 
288 inline
289 double
ofxsFilterMitchell(double Ip,double Ic,double In,double Ia,double d,bool clamp)290 ofxsFilterMitchell(double Ip,
291                    double Ic,
292                    double In,
293                    double Ia,
294                    double d,
295                    bool clamp)
296 {
297     double I = ( Ip + 16 * Ic + In + d * ( (-9 * Ip + 9 * In ) + d * ( (15 * Ip - 36 * Ic + 27 * In - 6 * Ia ) + d * (-7 * Ip + 21 * Ic - 21 * In + 7 * Ia ) ) ) ) / 18;
298 
299     if (clamp) {
300         I = ofxsFilterClampVal(I, Ic, In);
301     }
302 
303     return I;
304 }
305 
306 inline
307 double
ofxsFilterParzen(double Ip,double Ic,double In,double Ia,double d,bool)308 ofxsFilterParzen(double Ip,
309                  double Ic,
310                  double In,
311                  double Ia,
312                  double d,
313                  bool /*clamp*/)
314 {
315     double I = ( Ip + 4 * Ic + In + d * ( (-3 * Ip + 3 * In ) + d * ( (3 * Ip - 6 * Ic + 3 * In ) + d * (-Ip + 3 * Ic - 3 * In + Ia ) ) ) ) / 6;
316 
317     // clamp is not necessary for Parzen
318     return I;
319 }
320 
321 inline
322 double
ofxsFilterNotch(double Ip,double Ic,double In,double Ia,double d,bool)323 ofxsFilterNotch(double Ip,
324                 double Ic,
325                 double In,
326                 double Ia,
327                 double d,
328                 bool /*clamp*/)
329 {
330     double I = ( Ip + 2 * Ic + In + d * ( (-2 * Ip + 2 * In ) + d * ( (Ip - Ic - In + Ia ) ) ) ) / 4;
331 
332     // clamp is not necessary for Notch
333     return I;
334 }
335 
336 
337 /////////////////////////////////////////////////
338 // BOX FILTER START
339 /////////////////////////////////////////////////
340 
341 /// @brief Add to vector v the integral of the signal contained in l, seen as piecewise constant, from x1 to x2.
342 /// x = 0 corresponds to the left of the first pixel, x = 1 corresponds to the right of the first pixel / left of the second pixel
343 ///
344 template <class PIX>
345 void
ofxsFilterIntegrate1d(const PIX * l,const size_t nsamples,const size_t stride,const size_t depth,const double x1,const double x2,const bool zeroOutside,float * v)346 ofxsFilterIntegrate1d(const PIX* l, // pointer to data start
347                       const size_t nsamples,  // number of samples in the line
348                       const size_t stride, // increment from one data point to the next
349                       const size_t depth, // dimension of each sample, also the dimension of result vector v
350                       const double x1,
351                       const double x2,
352                       const bool zeroOutside, // if true, outside of the data is zero. If false, use Neumann boundary conditions (outside is the closest data point)
353                       float *v) // vector of dimension depth containing the result
354 {
355     assert(x2 >= x1);
356     assert(stride >= depth);
357     size_t ifirst, ilast; // index of the first/last pixel
358     double fracfirst, fraclast; // fraction to remove from the first/last pixel
359     if (x1 < 0.) {
360         ifirst = 0;
361         fracfirst = 0.;
362     } else if (nsamples <= x1) {
363         ifirst = nsamples - 1;
364         fracfirst = 0.;
365     } else {
366         ifirst = (size_t)floor(x1);
367         fracfirst = x1 - ifirst;
368     }
369     if (x2 < 0.) {
370         ilast = 0;
371         fraclast = 0.;
372     } else if (nsamples <= x2) {
373         ilast = nsamples - 1;
374         fraclast = 0.;
375     } else {
376         ilast =  (size_t)floor(x2);
377         fraclast = ilast + 1 - x2;
378     }
379     // start border condition
380     if (x1 < ifirst && !zeroOutside) {
381         for (size_t j = 0; j < depth; ++j) {
382             v[j] += l[ifirst * stride + j] * (float)(ifirst-x1);
383         }
384     }
385     // pre-subtract partial first pixel
386     if (fracfirst > 0.) {
387         for (size_t j = 0; j < depth; ++j) {
388             v[j] -= l[ifirst * stride + j] * (float)fracfirst;
389         }
390     }
391     // sum all covered pixels
392     for (size_t i = ifirst; i <= ilast; ++i) {
393         for (size_t j = 0; j < depth; ++j) {
394             v[j] += l[i * stride + j];
395         }
396     }
397     // subtract partial last pixel
398     if (fraclast > 0.) {
399         for (size_t j = 0; j < depth; ++j) {
400             v[j] -= l[ilast * stride + j] * (float)fraclast;
401         }
402     }
403     // end border condition
404     if (x2 > nsamples && !zeroOutside) {
405         for (size_t j = 0; j < depth; ++j) {
406             v[j] += l[ilast * stride + j] * (float)(x2 - nsamples);
407         }
408     }
409 }
410 
411 /// @brief Compute the mean of the signal contained in l, seen as piecewise constant, in the rectangular area delimited by x1, x2, y1, y2.
412 /// x = 0 corresponds to the left of the first pixel, x = 1 corresponds to the right of the first pixel / left of the second pixel
413 ///
414 template <class PIX>
415 void
ofxsFilterIntegrate2d(const PIX * a,const size_t awidth,const size_t aheight,const size_t axstride,const size_t aystride,const size_t depth,const OfxRectD & area,const bool zeroOutside,float * p,float * v)416 ofxsFilterIntegrate2d(const PIX* a, // pointer to data start
417                       const size_t awidth,  // width of the array
418                       const size_t aheight,  // height of the array
419                       const size_t axstride, // increment from one data point to the next (must be >= depth)
420                       const size_t aystride, // increment from one data line to the next (usually awidth * axstride)
421                       const size_t depth, // dimension of each sample, also the dimension of result vector v
422                       const OfxRectD& area,
423                       const bool zeroOutside, // if true, outside of the data is zero. If false, use Neumann boundary conditions (outside is the closest data point)
424                       float *p, // temporary storage of size depth
425                       float *v) // vector of dimension depth containing the result
426 {
427     double x1 = area.x1;
428     double y1 = area.y1;
429     double x2 = area.x2;
430     double y2 = area.y2;
431     assert(y2 >= y1);
432     size_t ifirst, ilast; // index of the first/last line
433     double fracfirst, fraclast; // fraction to remove from the first/last line
434     if (y1 < 0.) {
435         ifirst = 0;
436         fracfirst = 0.;
437     } else if (aheight <= y1) {
438         ifirst = aheight - 1;
439         fracfirst = 0.;
440     } else {
441         ifirst = (size_t)floor(y1);
442         fracfirst = y1 - ifirst;
443     }
444     if (y2 < 0.) {
445         ilast = 0;
446         fraclast = 0.;
447     } else if (aheight <= y2) {
448         ilast = aheight - 1;
449         fraclast = 0.;
450     } else {
451         ilast =  (size_t)floor(y2);
452         fraclast = ilast + 1 - y2;
453     }
454 
455     // compute result for first line
456     for (size_t j = 0; j < depth; ++j) {
457         p[j] = 0.;
458     }
459     ofxsFilterIntegrate1d(&a[ifirst * aystride], awidth, axstride, depth, x1, x2, zeroOutside, p);
460     // start border condition
461     if (y1 < ifirst && !zeroOutside) {
462         for (size_t j = 0; j < depth; ++j) {
463             v[j] += p[j] * (float)(ifirst - y1);
464         }
465     }
466     // subtract partial first line
467     if (fracfirst > 0.) {
468         for (size_t j = 0; j < depth; ++j) {
469             v[j] -= p[j] * (float)fracfirst;
470         }
471     }
472     // sum all covered lines
473     // first line
474     for (size_t j = 0; j < depth; ++j) {
475         v[j] += p[j];
476     }
477     // all lines except first and last
478     for (size_t i = ifirst + 1; i < ilast; ++i) {
479         // (results accumulates in v)
480         ofxsFilterIntegrate1d(&a[i * aystride], awidth, axstride, depth, x1, x2, zeroOutside, v);
481     }
482     // last line
483     if (ilast > ifirst) { // (if equal, the result is already in p)
484         for (size_t j = 0; j < depth; ++j) {
485             p[j] = 0.;
486         }
487         ofxsFilterIntegrate1d(&a[ilast * aystride], awidth, axstride, depth, x1, x2, zeroOutside, p);
488         for (size_t j = 0; j < depth; ++j) {
489             v[j] += p[j];
490         }
491     }
492     // subtract partial last pixel
493     if (fraclast > 0.) {
494         for (size_t j = 0; j < depth; ++j) {
495             v[j] -= p[j] * (float)fraclast;
496         }
497     }
498     // end border condition
499     if (y2 > aheight && !zeroOutside) {
500         for (size_t j = 0; j < depth; ++j) {
501             v[j] += p[j] * (float)(y2 - aheight);
502         }
503     }
504 }
505 
506 /// @brief resize the area from image a indicated by from and put it in image b at to.
507 /// If @param from is partially outside of a, pixels are considered to be black and transparent if zeroOutside is true,
508 /// else they take the value of the closest pixel in a.
509 /// The @param to may be partially outside of b.
510 template <class PIX>
511 void
ofxsFilterResize2d(const PIX * a,const size_t awidth,const size_t aheight,const size_t axstride,const size_t aystride,const size_t depth,const OfxRectD & from,const bool zeroOutside,float * b,const size_t bwidth,const size_t bheight,const size_t bxstride,const size_t bystride,const OfxRectI & to)512 ofxsFilterResize2d(const PIX* a, // pointer to data start
513                    const size_t awidth,  // number of samples in the line
514                    const size_t aheight,  // number of samples in the line
515                    const size_t axstride, // increment from one data point to the next (must be >= depth)
516                    const size_t aystride, // increment from one data line to the next (usually awidth * axstride)
517                    const size_t depth, // dimension of each sample, also the dimension of result vector v
518                    const OfxRectD& from,
519                    const bool zeroOutside, // if true, outside of the data is zero (Dirichlet boundary conditions). If false, outside is the closest data point (Neumann boundary conditions).
520                    float* b, // pointer to output start
521                    const size_t bwidth,  // number of samples in the line
522                    const size_t bheight,  // number of samples in the line
523                    const size_t bxstride, // inscrement from one data point to the next (must be >= depth)
524                    const size_t bystride,
525                    const OfxRectI& to)
526 
527 {
528     assert(awidth > 0 && aheight > 0 && axstride > 0 && aystride > 0 && depth > 0);
529     assert(bwidth > 0 && bheight > 0 && bxstride > 0 && bystride > 0);
530     double x1 = from.x1;
531     double y1 = from.y1;
532     double x2 = from.x2;
533     double y2 = from.y2;
534     assert(x2 >= x1);
535     assert(y2 >= y1);
536     int ox1 = to.x1;
537     int oy1 = to.y1;
538     int ox2 = to.x2;
539     int oy2 = to.y2;
540     assert(ox2 > ox1);
541     assert(oy2 > oy1);
542     // pixel factor
543     double vwidth = (x2 - x1) / (ox2 - ox1);
544     double vheight = (y2 - y1) / (oy2 - oy1);
545 
546     // adjust output to valid areas of b
547     if (ox1 < 0) {
548         x1 -= vwidth * ox1;
549         ox1 = 0;
550     }
551     if (ox2 > (int)bwidth) {
552         x2 -= vwidth * ((int)bwidth - ox2);
553         ox2 = (int)bwidth;
554     }
555     assert(x2 >= x1);
556     assert(ox2 >= ox1);
557     if (ox2 <= ox1) {
558         // nothing to draw
559         return;
560     }
561     if (oy1 < 0) {
562         y1 -= vheight * oy1;
563         oy1 = 0;
564     }
565     if (oy2 > (int)bheight) {
566         y2 -= vheight * ((int)bheight - oy2);
567         oy2 = (int)bheight;
568     }
569     assert(y2 >= y1);
570     assert(oy2 >= oy1);
571     if (oy2 <= oy1) {
572         // nothing to draw
573         return;
574     }
575 
576     float *p = new float[depth];
577     // #pragma parallel for
578     for (int j = oy1; j < oy2; ++j) {
579         OfxRectD area;
580         area.y1 = y1 + (j - oy1) * vheight;
581         area.y2 = area.y1 + vheight;
582         for (int i = ox1; i < ox2; ++i) {
583             area.x1 = x1 + (i - ox1) * vwidth;
584             area.x2 = area.x1 + vwidth;
585             // compute one pixel of the resized image
586             float *v = &b[j * bystride + i * bxstride];
587             // zero the result, since integrate_2d accumulates
588             for (size_t k = 0; k < depth; ++k) {
589                 v[k] = 0.;
590             }
591             ofxsFilterIntegrate2d(a, awidth, aheight, axstride, aystride, depth,
592                                   area,
593                                   zeroOutside,
594                                   p,
595                                   v);
596             // normalize by the surface of the pixel
597             for (size_t k = 0; k < depth; ++k) {
598                 v[k] /= vwidth * vheight;
599             }
600         }
601     }
602     delete [] p;
603 }
604 
605 /////////////////////////////////////////////////
606 // BOX FILTER END
607 /////////////////////////////////////////////////
608 
609 
610 #define OFXS_APPLY4(f, j) double I ## j = f(Ip ## j, Ic ## j, In ## j, Ia ## j, dx, clamp)
611 
612 #define OFXS_CUBIC2D(f)                                      \
613     inline                                           \
614     double                                                  \
615     f ## 2D (double Ipp, double Icp, double Inp, double Iap, \
616              double Ipc, double Icc, double Inc, double Iac, \
617              double Ipn, double Icn, double Inn, double Ian, \
618              double Ipa, double Ica, double Ina, double Iaa, \
619              double dx, double dy, bool clamp)               \
620     {                                                       \
621         OFXS_APPLY4(f, p); OFXS_APPLY4(f, c); OFXS_APPLY4(f, n); OFXS_APPLY4(f, a); \
622         return f(Ip, Ic, In, Ia, dy, clamp);            \
623     }
624 
625 OFXS_CUBIC2D(ofxsFilterKeys);
626 OFXS_CUBIC2D(ofxsFilterSimon);
627 OFXS_CUBIC2D(ofxsFilterRifman);
628 OFXS_CUBIC2D(ofxsFilterMitchell);
629 OFXS_CUBIC2D(ofxsFilterParzen);
630 OFXS_CUBIC2D(ofxsFilterNotch);
631 
632 #undef OFXS_CUBIC2D
633 #undef OFXS_APPLY4
634 
635 template <class PIX>
636 PIX
ofxsGetPixComp(const PIX * p,int c)637 ofxsGetPixComp(const PIX* p,
638                int c)
639 {
640     return p ? p[c] : PIX();
641 }
642 
643 // Macros used in ofxsFilterInterpolate2D
644 #define OFXS_CLAMPXY(m) \
645     m ## x = (std::max)( srcImg->getBounds().x1, (std::min)(m ## x, srcImg->getBounds().x2 - 1) ); \
646     m ## y = (std::max)( srcImg->getBounds().y1, (std::min)(m ## y, srcImg->getBounds().y2 - 1) )
647 
648 #define OFXS_GETPIX(i, j) PIX * P ## i ## j = (PIX *)srcImg->getPixelAddress(i ## x, j ## y)
649 
650 #define OFXS_GETI(i, j)   const double I ## i ## j = ofxsGetPixComp(P ## i ## j, c)
651 
652 #define OFXS_GETPIX4(i)  OFXS_GETPIX(i, p); OFXS_GETPIX(i, c); OFXS_GETPIX(i, n); OFXS_GETPIX(i, a);
653 
654 #define OFXS_GETI4(i)    OFXS_GETI(i, p); OFXS_GETI(i, c); OFXS_GETI(i, n); OFXS_GETI(i, a);
655 
656 
657 #define OFXS_I44         Ipp, Icp, Inp, Iap, \
658     Ipc, Icc, Inc, Iac, \
659     Ipn, Icn, Inn, Ian, \
660     Ipa, Ica, Ina, Iaa
661 
662 // note that the center of pixel (0,0) has pixel coordinates (0.5,0.5)
663 template <class PIX, int nComponents, FilterEnum filter, bool clamp>
664 bool
ofxsFilterInterpolate2D(double fx,double fy,const OFX::Image * srcImg,bool blackOutside,float * tmpPix)665 ofxsFilterInterpolate2D(double fx,
666                         double fy,            //!< coordinates of the pixel to be interpolated in srcImg in pixel coordinates
667                         const OFX::Image *srcImg, //!< image to be transformed
668                         bool blackOutside,
669                         float *tmpPix) //!< destination pixel in float format
670 {
671     if (!srcImg) {
672         for (int c = 0; c < nComponents; ++c) {
673             tmpPix[c] = 0;
674         }
675 
676         return false;
677     }
678     bool inside = true; // return true, except if outside and black
679     // GENERIC TRANSFORM
680     // from here on, everything is generic, and should be moved to a generic transform class
681     // Important: (0,0) is the *corner*, not the *center* of the first pixel (see OpenFX specs)
682     switch (filter) {
683     case eFilterImpulse:
684     case eFilterBox: {
685         ///nearest neighboor
686         // the center of pixel (0,0) has coordinates (0.5,0.5)
687         int mx = (int)std::floor(fx);     // don't add 0.5
688         int my = (int)std::floor(fy);     // don't add 0.5
689 
690         if (!blackOutside) {
691             OFXS_CLAMPXY(m);
692         }
693         OFXS_GETPIX(m, m);
694         if (Pmm) {
695             for (int c = 0; c < nComponents; ++c) {
696                 tmpPix[c] = Pmm[c];
697             }
698         } else {
699             for (int c = 0; c < nComponents; ++c) {
700                 tmpPix[c] = 0;
701             }
702             inside = false;
703         }
704         break;
705     }
706     case eFilterBilinear:
707     case eFilterCubic: {
708         // bilinear or cubic
709         // the center of pixel (0,0) has coordinates (0.5,0.5)
710         int cx = (int)std::floor(fx - 0.5);
711         int cy = (int)std::floor(fy - 0.5);
712         int nx = cx + 1;
713         int ny = cy + 1;
714         if (!blackOutside) {
715             OFXS_CLAMPXY(c);
716             OFXS_CLAMPXY(n);
717         }
718 
719         const double dx = (std::max)( 0., (std::min)(fx - 0.5 - cx, 1.) );
720         const double dy = (std::max)( 0., (std::min)(fy - 0.5 - cy, 1.) );
721 
722         OFXS_GETPIX(c, c); OFXS_GETPIX(n, c); OFXS_GETPIX(c, n); OFXS_GETPIX(n, n);
723         if (Pcc || Pnc || Pcn || Pnn) {
724             for (int c = 0; c < nComponents; ++c) {
725                 OFXS_GETI(c, c); OFXS_GETI(n, c); OFXS_GETI(c, n); OFXS_GETI(n, n);
726                 if (filter == eFilterBilinear) {
727                     double Ic = ofxsFilterLinear(Icc, Inc, dx);
728                     double In = ofxsFilterLinear(Icn, Inn, dx);
729                     tmpPix[c] = (float)ofxsFilterLinear(Ic, In, dy);
730                 } else if (filter == eFilterCubic) {
731                     double Ic = ofxsFilterCubic(Icc, Inc, dx, clamp);
732                     double In = ofxsFilterCubic(Icn, Inn, dx, clamp);
733                     tmpPix[c] = (float)ofxsFilterCubic(Ic, In, dy, clamp);
734                 } else {
735                     assert(0);
736                 }
737             }
738         } else {
739             for (int c = 0; c < nComponents; ++c) {
740                 tmpPix[c] = 0;
741             }
742             inside = false;
743         }
744         break;
745     }
746 
747     // (B,C) cubic filters
748     case eFilterKeys:
749     case eFilterSimon:
750     case eFilterRifman:
751     case eFilterMitchell:
752     case eFilterParzen:
753     case eFilterNotch: {
754         // the center of pixel (0,0) has coordinates (0.5,0.5)
755         int cx = (int)std::floor(fx - 0.5);
756         int cy = (int)std::floor(fy - 0.5);
757         int px = cx - 1;
758         int py = cy - 1;
759         int nx = cx + 1;
760         int ny = cy + 1;
761         int ax = cx + 2;
762         int ay = cy + 2;
763         if (!blackOutside) {
764             OFXS_CLAMPXY(c);
765             OFXS_CLAMPXY(p);
766             OFXS_CLAMPXY(n);
767             OFXS_CLAMPXY(a);
768         }
769         const double dx = (std::max)( 0., (std::min)(fx - 0.5 - cx, 1.) );
770         const double dy = (std::max)( 0., (std::min)(fy - 0.5 - cy, 1.) );
771 
772         OFXS_GETPIX4(p); OFXS_GETPIX4(c); OFXS_GETPIX4(n); OFXS_GETPIX4(a);
773         if (Ppp || Pcp || Pnp || Pap || Ppc || Pcc || Pnc || Pac || Ppn || Pcn || Pnn || Pan || Ppa || Pca || Pna || Paa) {
774             for (int c = 0; c < nComponents; ++c) {
775                 //double Ipp = get(Ppp,c);, etc.
776                 OFXS_GETI4(p); OFXS_GETI4(c); OFXS_GETI4(n); OFXS_GETI4(a);
777                 double I = 0.;
778                 switch (filter) {
779                 case eFilterKeys:
780                     I = ofxsFilterKeys2D(OFXS_I44, dx, dy, clamp);
781                     break;
782                 case eFilterSimon:
783                     I = ofxsFilterSimon2D(OFXS_I44, dx, dy, clamp);
784                     break;
785                 case eFilterRifman:
786                     I = ofxsFilterRifman2D(OFXS_I44, dx, dy, clamp);
787                     break;
788                 case eFilterMitchell:
789                     I = ofxsFilterMitchell2D(OFXS_I44, dx, dy, clamp);
790                     break;
791                 case eFilterParzen:
792                     I = ofxsFilterParzen2D(OFXS_I44, dx, dy, false);
793                     break;
794                 case eFilterNotch:
795                     I = ofxsFilterNotch2D(OFXS_I44, dx, dy, false);
796                     break;
797                 default:
798                     assert(0);
799                 }
800                 tmpPix[c] = (float)I;
801             }
802         } else {
803             for (int c = 0; c < nComponents; ++c) {
804                 tmpPix[c] = 0;
805             }
806             inside = false;
807         }
808         break;
809     }
810 
811     default:
812         assert(0);
813         break;
814     } // switch
815 
816     return inside;
817 } // ofxsFilterInterpolate2D
818 
819 /*
820  * Interpolation with SuperSampling, to avoid moire artifacts when minimizing.
821  *
822 
823    ofxsFilterInterpolate2D() does not take into account scaling or distortion effects.
824    A consequence is that the transform nodes may produce aliasing artifacts when downscaling by a factor of 2 or more.
825 
826    There are several solutions to this problem is the case where the same texture has to be mapped *several times*:
827 
828  * Trilinear mipmapping (as implemented by OpenGL) still produces artifacts when scaling is anisotropic (i.e. the scaling factor is different along two directions)
829  * [Feline (McCormack, 1999)](http://www.hpl.hp.com/techreports/Compaq-DEC/WRL-99-1.pdf), which is close to what is proposed in [OpenGL's anisotropic texture filter](http://www.opengl.org/registry/specs/EXT/texture_filter_anisotropic.txt) is probably 4-5 times slower than mipmapping, but produces less artifacts
830  * [EWA (Heckbert 1989)](https://www.cs.cmu.edu/~ph/texfund/texfund.pdf) would give the highest quality, but is probably 20 times slower than mipmapping.
831 
832    A sample implementation of the three methods is given in [Mesa 3D](http://mesa3d.org/)'s [software rasterizer, src/mesa/swrast/s_texfilter.c](http://cgit.freedesktop.org/mesa/mesa/tree/src/mesa/swrast/s_texfilter.c).
833 
834  * However*, in our case, the texture usually has to be mapped only once. Thus it is more appropriate to use one of the techniques described in this document: <http://people.cs.clemson.edu/~dhouse/courses/405/notes/antialiasing.pdf>.
835 
836  # Our implementation:
837 
838    We chose to use a standard supersampling method without jitter (to avoid randomness in rendered images).
839 
840    The first implementation was interpolating accross scales between supersampled pixels (see OFX_FILTER_SUPERSAMPLING_TRILINEAR below), but since we noticed that using the highest scale produces less moire, and it even costs a bit less (less tests in the processing).
841 
842    We also noticed that supersampled pixels don't need to use anything better than bilinear filter. The impulse filter still produces moire, and other filters are overkill or may even produce more moire.
843 
844  */
845 
846 #ifdef DEBUG
847 #define _DBG_COUNT(x) (x)
848 #else
849 #define _DBG_COUNT(x) ( (void)0 )
850 #endif
851 
852 // Internal function for supersampling (should never be called by the user)
853 // note that the center of pixel (0,0) has pixel coordinates (0.5,0.5)
854 template <class PIX, int nComponents, FilterEnum filter, int subx, int suby>
855 void
ofxsFilterInterpolate2DSuperInternal(double fx,double fy,double Jxx,double Jxy,double Jyx,double Jyy,double sx,double sy,int isx,int isy,const OFX::Image * srcImg,bool blackOutside,float * tmpPix)856 ofxsFilterInterpolate2DSuperInternal(double fx,
857                                      double fy,            //!< coordinates of the pixel to be interpolated in srcImg in pixel coordinates
858                                      double Jxx, //!< derivative of fx over x
859                                      double Jxy, //!< derivative of fx over y
860                                      double Jyx, //!< derivative of fy over x
861                                      double Jyy, //!< derivative of fy over y
862                                      double sx, //!< scale over x as a power of 3
863                                      double sy, //!< scale over y as a power of 3
864                                      int isx, //!< floor(sx)
865                                      int isy,  //!< floor(sy)
866                                      const OFX::Image *srcImg, //!< image to be transformed
867                                      bool blackOutside,
868                                      float *tmpPix) //!< input: interpolated center filter. output: destination pixel in float format
869 {
870     // do supersampling.
871     // All values are computed using nearest neighbor interpolation, except for the center value
872 
873     // compute number of samples over each dimension, i.e. pow(nis*,3)
874     // see http://stackoverflow.com/a/101613/2607517
875     int nisx;
876     {
877         int base = 3;
878         int exp = isx;
879         int result = 1;
880         while (exp) {
881             if (exp & 1) {
882                 result *= base;
883             }
884             exp >>= 1;
885             base *= base;
886         }
887         nisx = result;
888     }
889     /// linear version:
890     //nisx = 1;
891     //for (int p = 0; p < isx; ++p) {
892     //    nisx *= 3;
893     //}
894     int nisy;
895     {
896         int base = 3;
897         int exp = isy;
898         int result = 1;
899         while (exp) {
900             if (exp & 1) {
901                 result *= base;
902             }
903             exp >>= 1;
904             base *= base;
905         }
906         nisy = result;
907     }
908 
909     /// linear version:
910     //nisy = 1;
911     //for (int p = 0; p < isy; ++p) {
912     //    nisy *= 3;
913     //}
914     assert( nisx == std::pow( (double)3, (double)isx ) && nisy == std::pow( (double)3, (double)isy ) );
915 
916     // compute the pixel value at scales (isx,isy) (nsx,isy) (isx,nsy) (nsx,nsy), and interpolate bilinearly using sx,sy
917     float *pii = tmpPix;
918     float pni[nComponents];
919     float pin[nComponents];
920     float pnn[nComponents];
921 #ifdef DEBUG
922     int piicount = 1;
923     int pnicount = 0;
924     int pincount = 0;
925     int pnncount = 0;
926 #endif
927 
928     // initialize to value of center pixel
929     if (subx) {
930         std::copy(tmpPix, tmpPix + nComponents, pni);
931         _DBG_COUNT(pnicount = 1);
932         if (suby) {
933             std::copy(tmpPix, tmpPix + nComponents, pnn);
934             _DBG_COUNT(pnncount = 1);
935         }
936     }
937     if (suby) {
938         std::copy(tmpPix, tmpPix + nComponents, pin);
939         _DBG_COUNT(pincount = 1);
940     }
941 
942     // accumulate
943     for (int y = -nisy / 2; y <= nisy / 2; ++y) {
944         for (int x = -nisx / 2; x <= nisx / 2; ++x) {
945             // subsample position
946             double sfx = fx + (Jxx * x) / nisx + (Jxy * y) / nisy;
947             double sfy = fy + (Jyx * x) / nisx + (Jyy * y) / nisy;
948             if ( (x != 0) || (y != 0) ) { // center pixel was already accumulated
949                 float tmp[nComponents];
950                 ofxsFilterInterpolate2D<PIX, nComponents, filter, false>(sfx, sfy, srcImg, blackOutside, tmp);
951                 for (int c = 0; c < nComponents; ++c) {
952                     pii[c] += tmp[c];
953                     _DBG_COUNT( piicount += (c == 0) );
954                     // other scales
955                     if (subx) {
956                         pni[c] += tmp[c];
957                         _DBG_COUNT( pnicount += (c == 0) );
958                         if (suby) {
959                             pnn[c] += tmp[c];
960                             _DBG_COUNT( pnncount += (c == 0) );
961                         }
962                     }
963                     if (suby) {
964                         pin[c] += tmp[c];
965                         _DBG_COUNT( pincount += (c == 0) );
966                     }
967                 }
968             }
969             // subsamples from next scales
970             for (int j = -suby; j <= suby; ++j) {
971                 for (int i = -subx; i <= subx; ++i) {
972                     if ( (i != 0) || (j != 0) ) { // center subsample was already accumulated
973                         double subfx = sfx + (Jxx * i) / (nisx * 3) + (Jxy * j) / (nisy * 3);
974                         double subfy = sfy + (Jyx * i) / (nisx * 3) + (Jyy * j) / (nisy * 3);
975                         {
976                             float tmp[nComponents];
977                             ofxsFilterInterpolate2D<PIX, nComponents, filter, false>(subfx, subfy, srcImg, blackOutside, tmp);
978                             for (int c = 0; c < nComponents; ++c) {
979                                 // other scales
980                                 if (subx) {
981                                     if (j == 0) {
982                                         pni[c] += tmp[c];
983                                         _DBG_COUNT( pnicount += (c == 0) );
984                                     }
985                                     if (suby) {
986                                         pnn[c] += tmp[c];
987                                         _DBG_COUNT( pnncount += (c == 0) );
988                                     }
989                                 }
990                                 if (suby) {
991                                     if (i == 0) {
992                                         pin[c] += tmp[c];
993                                         _DBG_COUNT( pincount += (c == 0) );
994                                     }
995                                 }
996                             }
997                         }
998                     }
999                 }
1000             }
1001         }
1002     }
1003 
1004     // divide by number of values
1005     int insamples = nisx * nisy;
1006 
1007 #ifdef DEBUG
1008     assert(piicount == insamples);
1009     if (subx) {
1010         assert(pnicount == insamples * 3);
1011         if (suby) {
1012             assert(pnncount == insamples * 9);
1013         }
1014     }
1015     if (suby) {
1016         assert(pincount == insamples * 3);
1017     }
1018 #endif
1019 
1020     for (int c = 0; c < nComponents; ++c) {
1021         pii[c] /= insamples;
1022         if (subx) {
1023             pni[c] /= insamples * 3;
1024             if (suby) {
1025                 pnn[c] /= insamples * 9;
1026             }
1027         }
1028         if (suby) {
1029             pin[c] /= insamples * 3;
1030         }
1031     }
1032     if (subx) {
1033         // interpolate linearly over sx
1034         float alpha = (float)(sx - isx);
1035         for (int c = 0; c < nComponents; ++c) {
1036             pii[c] = pii[c] + alpha * (pni[c] - pii[c]);
1037         }
1038         if (suby) {
1039             for (int c = 0; c < nComponents; ++c) {
1040                 pin[c] = pin[c] + alpha * (pnn[c] - pin[c]);
1041             }
1042         }
1043     }
1044     if (suby) {
1045         // interpolate linearly over sy
1046         float alpha = (float)(sy - isy);
1047         for (int c = 0; c < nComponents; ++c) {
1048             pii[c] = pii[c] + alpha * (pin[c] - pii[c]);
1049         }
1050     }
1051 
1052     // pii is actually an alias to tmpPix, so everything is done
1053 } // ofxsFilterInterpolate2DSuperInternal
1054 
1055 #undef _DBG_COUNT
1056 
1057 inline bool
ofxsFilterOutside(double x,double y,const OfxRectI & bounds)1058 ofxsFilterOutside(double x,
1059                   double y,
1060                   const OfxRectI &bounds)
1061 {
1062     return x < bounds.x1 || bounds.x2 <= x || y < bounds.y1 || bounds.y2 <= y;
1063 }
1064 
1065 // Interpolation using the given filter and supersampling for minification
1066 // note that the center of pixel (0,0) has pixel coordinates (0.5,0.5)
1067 template <class PIX, int nComponents, FilterEnum filter, bool clamp>
1068 void
ofxsFilterInterpolate2DSuper(double fx,double fy,double Jxx,double Jxy,double Jyx,double Jyy,const OFX::Image * srcImg,bool blackOutside,float * tmpPix)1069 ofxsFilterInterpolate2DSuper(double fx,
1070                              double fy,            //!< coordinates of the pixel to be interpolated in srcImg in pixel coordinates
1071                              double Jxx, //!< derivative of fx over x
1072                              double Jxy, //!< derivative of fx over y
1073                              double Jyx, //!< derivative of fy over x
1074                              double Jyy, //!< derivative of fy over y
1075                              const OFX::Image *srcImg, //!< image to be transformed
1076                              bool blackOutside,
1077                              float *tmpPix) //!< destination pixel in float format
1078 {
1079     if ( !srcImg || !srcImg->getPixelData() ) {
1080         for (int c = 0; c < nComponents; ++c) {
1081             tmpPix[c] = 0.;
1082         }
1083         return;
1084     }
1085     if (Jxx == 0. && Jxy == 0. && Jyx == 0. && Jyy == 0.) {
1086         ofxsFilterInterpolate2D<PIX,nComponents,filter,clamp>(fx, fy, srcImg, blackOutside, tmpPix);
1087 
1088         return;
1089     }
1090     if (filter == eFilterBox) {
1091         for (int c = 0; c < nComponents; ++c) {
1092             tmpPix[c] = 0.;
1093         }
1094         // Box filter is a special case:
1095         // 1- compute the bounding box of the backtransformed pixel
1096         // 2- integrate the input image over this bounding box
1097         //
1098         //
1099         double x, y;
1100         double x1, x2, y1, y2;
1101         x1 = x2 = fx - Jxx * 0.5 - Jxy * 0.5;
1102         y1 = y2 = fy - Jyx * 0.5 - Jyy * 0.5;
1103         x = fx + Jxx * 0.5 - Jxy * 0.5;
1104         y = fy + Jyx * 0.5 - Jyy * 0.5;
1105         x1 = (std::min)(x1, x);
1106         y1 = (std::min)(y1, y);
1107         x2 = (std::max)(x2, x);
1108         y2 = (std::max)(y2, y);
1109         x = fx - Jxx * 0.5 + Jxy * 0.5;
1110         y = fy - Jyx * 0.5 + Jyy * 0.5;
1111         x1 = (std::min)(x1, x);
1112         y1 = (std::min)(y1, y);
1113         x2 = (std::max)(x2, x);
1114         y2 = (std::max)(y2, y);
1115         x = fx + Jxx * 0.5 + Jxy * 0.5;
1116         y = fy + Jyx * 0.5 + Jyy * 0.5;
1117         x1 = (std::min)(x1, x);
1118         y1 = (std::min)(y1, y);
1119         x2 = (std::max)(x2, x);
1120         y2 = (std::max)(y2, y);
1121         if (x2 <= x1 || y2 <= y1) {
1122             // empty pixel
1123             ofxsFilterInterpolate2D<PIX,nComponents,filter,clamp>(fx, fy, srcImg, blackOutside, tmpPix);
1124 
1125             return;
1126         }
1127 
1128         PIX* a = (PIX*)srcImg->getPixelData();
1129         const OfxRectI& srcBounds = srcImg->getBounds();
1130         const size_t awidth = srcBounds.x2 - srcBounds.x1;
1131         const size_t aheight = srcBounds.y2 - srcBounds.y1;
1132         x1 -= srcBounds.x1;
1133         y1 -= srcBounds.y1;
1134         x2 -= srcBounds.x1;
1135         y2 -= srcBounds.y1;
1136         const size_t axstride = srcImg->getPixelBytes() / sizeof(PIX);
1137         const size_t aystride = srcImg->getRowBytes() / sizeof(PIX);
1138         float p[nComponents];
1139         OfxRectD area = { x1, y1, x2, y2 };
1140         ofxsFilterIntegrate2d(a, awidth, aheight, axstride, aystride, nComponents,
1141                               area,
1142                               blackOutside,
1143                               p,
1144                               tmpPix);
1145         // normalize by the surface of the pixel
1146         float s = (float)((x2 - x1) * (y2 - y1));
1147         if (s != 0.f) {
1148             for (int c = 0; c < nComponents; ++c) {
1149                 tmpPix[c] /= s;
1150             }
1151         }
1152 
1153         return;
1154     }
1155     // first, compute the center value
1156     bool inside = ofxsFilterInterpolate2D<PIX, nComponents, filter, clamp>(fx, fy, srcImg, blackOutside, tmpPix);
1157 
1158     if (!inside) {
1159         // Center of the pixel is outside.
1160         // no supersampling if we're outside (we don't want to supersample black and transparent areas)
1161         // ... but we still have to check wether the entire pixel is outside
1162         const OfxRectI &bounds = srcImg->getBounds();
1163         // we check the four corners of the pixel
1164         if ( ofxsFilterOutside(fx - Jxx * 0.5 - Jxy * 0.5, fy - Jyx * 0.5 - Jyy * 0.5, bounds) &&
1165              ofxsFilterOutside(fx + Jxx * 0.5 - Jxy * 0.5, fy + Jyx * 0.5 - Jyy * 0.5, bounds) &&
1166              ofxsFilterOutside(fx - Jxx * 0.5 + Jxy * 0.5, fy - Jyx * 0.5 + Jyy * 0.5, bounds) &&
1167              ofxsFilterOutside(fx + Jxx * 0.5 + Jxy * 0.5, fy + Jyx * 0.5 + Jyy * 0.5, bounds) ) {
1168             return;
1169         }
1170     }
1171 
1172     double dx = Jxx * Jxx + Jyx * Jyx; // squared norm of the derivative over x
1173     double dy = Jxy * Jxy + Jyy * Jyy; // squared norm of the derivative over x
1174 
1175     if ( (dx <= 1.) && (dy <= 1.) ) {
1176         // no minificationin either direction, means no supersampling
1177         return;
1178     }
1179 
1180     // maximum scale is 4, which is 81x81 pixels for a scale factor < 1/81
1181     // rather than taking sqrt(dx), we divide its log by 2
1182     double sx = (dx <= 1.) ? 0. : (std::min)(std::log(dx) / ( 2 * std::log(3.) ), 4.); // scale over x as a power of 3
1183     double sy = (dy <= 1.) ? 0. : (std::min)(std::log(dy) / ( 2 * std::log(3.) ), 4.); // scale over y as a power of 3
1184 //#define OFX_FILTER_SUPERSAMPLING_TRILINEAR
1185 #ifdef OFX_FILTER_SUPERSAMPLING_TRILINEAR
1186     // produces artifacts
1187     int isx = std::floor(sx);
1188     int isy = std::floor(sy);
1189     int subx = (sx > isx);
1190     int suby = (sy > isy);
1191 
1192     // we use bilinear filtering for the supersamples (except for the center point).
1193     if (subx) {
1194         if (suby) {
1195             return ofxsFilterInterpolate2DSuperInternal<PIX, nComponents, eFilterBilinear, true, true>(fx, fy, Jxx, Jxy, Jyx, Jyy, sx, sy, isx, isy, srcImg, blackOutside, tmpPix);
1196         } else {
1197             return ofxsFilterInterpolate2DSuperInternal<PIX, nComponents, eFilterBilinear, true, false>(fx, fy, Jxx, Jxy, Jyx, Jyy, sx, sy, isx, isy, srcImg, blackOutside, tmpPix);
1198         }
1199     } else {
1200         if (suby) {
1201             return ofxsFilterInterpolate2DSuperInternal<PIX, nComponents, eFilterBilinear, false, true>(fx, fy, Jxx, Jxy, Jyx, Jyy, sx, sy, isx, isy, srcImg, blackOutside, tmpPix);
1202         } else {
1203             return ofxsFilterInterpolate2DSuperInternal<PIX, nComponents, eFilterBilinear, false, false>(fx, fy, Jxx, Jxy, Jyx, Jyy, sx, sy, isx, isy, srcImg, blackOutside, tmpPix);
1204         }
1205     }
1206 #else
1207     // always use the supersampled data
1208     // produces less artifacts, costs less
1209     // the problem is that sx = 1.0001 is supersampled, which gives a result very different from sx=1
1210     //int isx = std::ceil(sx);
1211     //int isy = std::ceil(sy);
1212     // This is why we prefer rounding. The jump will be at sx=sqrt(3)=1.732.
1213     // This produces quicker renders too, since we supersample less.
1214     int isx = (int)std::ceil(sx-0.5);
1215     int isy = (int)std::ceil(sy-0.5);
1216 
1217     return ofxsFilterInterpolate2DSuperInternal<PIX, nComponents, eFilterBilinear, false, false>(fx, fy, Jxx, Jxy, Jyx, Jyy, isx, isy, isx, isy, srcImg, blackOutside, tmpPix);
1218 #endif
1219 } // ofxsFilterInterpolate2DSuper
1220 
1221 #undef OFXS_CLAMPXY
1222 #undef OFXS_GETPIX
1223 #undef OFXS_GETI
1224 #undef OFXS_GETPIX4
1225 #undef OFXS_GETI
1226 #undef OFXS_I44
1227 
1228 
1229 inline void
ofxsFilterExpandRoD(OFX::ImageEffect *,double pixelAspectRatio,const OfxPointD & renderScale,bool blackOutside,OfxRectD * rod)1230 ofxsFilterExpandRoD(OFX::ImageEffect* /*effect*/,
1231                     double pixelAspectRatio,
1232                     const OfxPointD & renderScale,
1233                     bool blackOutside,
1234                     OfxRectD *rod)
1235 {
1236     if ( (rod->x2 <= rod->x1) || (rod->y2 <= rod->y1) ) {
1237         // empty rod
1238 
1239         return;
1240     }
1241 
1242     // No need to round things up here, we must give the *actual* RoD
1243 
1244     if (!blackOutside) {
1245         // if it's not black outside, the RoD should contain the project (we can't rely on the host to fill it).
1246         // [FD] 2014/09/01: disabled this. The transformed object may have a size which is very different from the project size
1247         /*
1248            OfxPointD size = effect->getProjectSize();
1249            OfxPointD offset = effect->getProjectOffset();
1250 
1251            rod->x1 = (std::min)(rod->x1, offset.x);
1252            rod->x2 = (std::max)(rod->x2, offset.x + size.x);
1253            rod->y1 = (std::min)(rod->y1, offset.y);
1254            rod->y2 = (std::max)(rod->y2, offset.y + size.y);
1255          */
1256     } else {
1257         // expand the RoD to get at least one black pixel
1258         double pixelSizeX = pixelAspectRatio / renderScale.x;
1259         double pixelSizeY = 1. / renderScale.y;
1260         if (rod->x1 > kOfxFlagInfiniteMin) {
1261             rod->x1 = rod->x1 - pixelSizeX;
1262         }
1263         if (rod->x2 < kOfxFlagInfiniteMax) {
1264             rod->x2 = rod->x2 + pixelSizeX;
1265         }
1266         if (rod->y1 > kOfxFlagInfiniteMin) {
1267             rod->y1 = rod->y1 - pixelSizeY;
1268         }
1269         if (rod->y2 < kOfxFlagInfiniteMax) {
1270             rod->y2 = rod->y2 + pixelSizeY;
1271         }
1272     }
1273 #if 0
1274     // The following code may be needed for hosts which do not
1275     // round correctly the RoD/RoI
1276     // This should correspond to pixel boundaries at the given renderScale,
1277     // which is why we have to round things up/down.
1278     if (rod->x1 > kOfxFlagInfiniteMin) {
1279         rod->x1 = ( std::floor(rod->x1 / pixelSizeX) ) * pixelSizeX;
1280     }
1281     if (rod->x2 < kOfxFlagInfiniteMax) {
1282         rod->x2 = ( std::ceil(rod->x2 / pixelSizeX) )  * pixelSizeX;
1283     }
1284     if (rod->y1 > kOfxFlagInfiniteMin) {
1285         rod->y1 = ( std::floor(rod->y1 / pixelSizeY) ) * pixelSizeY;
1286     }
1287     if (rod->y2 < kOfxFlagInfiniteMax) {
1288         rod->y2 = ( std::ceil(rod->y2 / pixelSizeY) )  * pixelSizeY;
1289     }
1290 #endif
1291     assert(rod->x1 <= rod->x2 && rod->y1 <= rod->y2);
1292 }
1293 
1294 inline void
ofxsFilterExpandRoI(const OfxRectD & roi,double pixelAspectRatio,const OfxPointD & renderScale,FilterEnum filter,bool doMasking,double mix,OfxRectD * srcRoI)1295 ofxsFilterExpandRoI(const OfxRectD &roi,
1296                     double pixelAspectRatio,
1297                     const OfxPointD & renderScale,
1298                     FilterEnum filter,
1299                     bool doMasking,
1300                     double mix,
1301                     OfxRectD *srcRoI)
1302 {
1303     // No need to round things up here, we must give the *actual* RoI,
1304     // the host should compute the right image region from it (by rounding it up/down).
1305 
1306     if ( (roi.x2 <= roi.x1) || (roi.y2 <= roi.y1) ) {
1307         *srcRoI = roi;
1308 
1309         return;
1310     }
1311 
1312     double pixelSizeX = pixelAspectRatio / renderScale.x;
1313     double pixelSizeY = 1. / renderScale.y;
1314 
1315     switch (filter) {
1316     case eFilterImpulse:
1317         // nearest neighbor, the exact region is OK
1318         break;
1319     case eFilterBox:
1320         // box filter, the exact region is OK
1321         break;
1322     case eFilterBilinear:
1323     case eFilterCubic:
1324         // bilinear or cubic, expand by 0.5 pixels
1325         if (srcRoI->x1 > kOfxFlagInfiniteMin) {
1326             srcRoI->x1 -= 0.5 * pixelSizeX;
1327         }
1328         if (srcRoI->x2 < kOfxFlagInfiniteMax) {
1329             srcRoI->x2 += 0.5 * pixelSizeX;
1330         }
1331         if (srcRoI->y1 > kOfxFlagInfiniteMin) {
1332             srcRoI->y1 -= 0.5 * pixelSizeY;
1333         }
1334         if (srcRoI->y2 < kOfxFlagInfiniteMax) {
1335             srcRoI->y2 += 0.5 * pixelSizeY;
1336         }
1337         break;
1338     case eFilterKeys:
1339     case eFilterSimon:
1340     case eFilterRifman:
1341     case eFilterMitchell:
1342     case eFilterParzen:
1343     case eFilterNotch:
1344         // bicubic, expand by 1.5 pixels
1345         if (srcRoI->x1 > kOfxFlagInfiniteMin) {
1346             srcRoI->x1 -= 1.5 * pixelSizeX;
1347         }
1348         if (srcRoI->x2 < kOfxFlagInfiniteMax) {
1349             srcRoI->x2 += 1.5 * pixelSizeX;
1350         }
1351         if (srcRoI->y1 > kOfxFlagInfiniteMin) {
1352             srcRoI->y1 -= 1.5 * pixelSizeY;
1353         }
1354         if (srcRoI->y2 < kOfxFlagInfiniteMax) {
1355             srcRoI->y2 += 1.5 * pixelSizeY;
1356         }
1357         break;
1358     }
1359     if ( doMasking || (mix != 1.) ) {
1360         // for masking or mixing, we also need the source image for that same roi.
1361         // compute the union of both ROIs
1362         srcRoI->x1 = (std::min)(srcRoI->x1, roi.x1);
1363         srcRoI->x2 = (std::max)(srcRoI->x2, roi.x2);
1364         srcRoI->y1 = (std::min)(srcRoI->y1, roi.y1);
1365         srcRoI->y2 = (std::max)(srcRoI->y2, roi.y2);
1366     }
1367 #if 0
1368     // The following code may be needed for hosts which do not
1369     // round correctly the RoD/RoI
1370     // This should correspond to pixel boundaries at the given renderScale,
1371     // which is why we have to round things up/down.
1372     if (srcRoI->x1 > kOfxFlagInfiniteMin) {
1373         srcRoI->x1 = std::floor(srcRoI->x1 / pixelSizeX) * pixelSizeX;
1374     }
1375     if (srcRoI->x2 < kOfxFlagInfiniteMax) {
1376         srcRoI->x2 = std::ceil(srcRoI->x2 / pixelSizeX)  * pixelSizeX;
1377     }
1378     if (srcRoI->y1 > kOfxFlagInfiniteMin) {
1379         srcRoI->y1 = std::floor(srcRoI->y1 / pixelSizeY) * pixelSizeY;
1380     }
1381     if (srcRoI->y2 < kOfxFlagInfiniteMax) {
1382         srcRoI->y2 = std::ceil(srcRoI->y2 / pixelSizeY)  * pixelSizeY;
1383     }
1384 #endif
1385     assert(srcRoI->x1 < srcRoI->x2 && srcRoI->y1 < srcRoI->y2);
1386 } // ofxsFilterExpandRoI
1387 } // OFX
1388 
1389 #endif // ifndef openfx_supportext_ofxsFilter_h
1390