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