BiCubic Interpolation does not match any known filter!

[magick]

Thanks for the continuity advice. I use alpha / beta / gamma as temporary variables. In many cases I had already used alpha / beta so that left gamma for dealing with transparency.

[NicolasRobidoux]

Thanks for the continuity advice. I use alpha / beta / gamma as temporary variables. In many cases I had already used alpha / beta so that left gamma for dealing with transparency.

(Of course my second suggested fix was absurd given if gamma is transparency and transparency is nominally in [0.1].)

-----

Warning: When using a non-monotone method like Catmull-Rom, gamma (as the interpolated alpha=transparency channel is called in this code) can overshoot it's bounds, not only below but also above. This is not an issue with bilinear or cubic or quadratic spline smoothing, but it definitely is with Catmull-Rom (and my suggested Mitchell).

You may need to add (continuous) clamping at the "top" (which is always equal to 1?).

[anthony]

Sorry to stick my nose in here. I don't want to make a lot of work for you here, but just wondering if this is limited to catmull-rom or can be done for the whole cubic family with B,C control. Would each B,C combination need separate optimization code?

Fred

It would be nice, BUT, interpolation is not set-up like a full filter which has a setup and structure to hold settings while you are making use of it., and needs to be fast. As such when any setting for it is set, you have to decode the setting string and store the numbers immediately. Where as in a resize filter, you can do that on filter set-up for the operation.

Actually that is why -interpolation filter is so darn slow. The filter needs to be set-up for each pixel lookup!!!!
I am seriously considering removing that interpolation.

ASIDE: Using a resampling filter for kernel generation (2d-cylindrical or 1d-tensor kernels) however would be useful and 'cool', especially if you can specify a 'scale' for that specific kernel.

Anyway.... The point is you have to have interpolation to be prepared and ready to be used for a direct lookup call, and any extra values beyond the interpolation method need to be pre-decoded and stored in every image, and need to be reset for every new image, as well as after any settings 'pop' from the setting stack. As such storing 'expert cubic coefficients' for interpolations may not be a good idea.

I think we should just stick to pre-defined cubics with hard-coded coefficients. They can use the same code, with different coefficients selected by the current interpolation method.

[anthony]

gamma being used for the inverted alpha weights, is fairly standard throughout the IM code. Though its naming threw me when I first saw it.

In the past I have changed it, but it would generate go back to using "gamma" when Crisy re-factored the code later, which doing some global code change.



The tricky part is that often you need two weights when dealing with a weight function that may not be included in the 'neighbourhood' (because of edges effects). that is raw weight added up and inverted for alpha, and alpha multiplied weights for color pixels (unless the color pixels get pre-multiplied, which does nto always happen).

This will get worse, if we start implementing my suggestion of allowing Virtual pixels 'none' actually mean no virtual pixel component in the result. Even though we are still getting virtual pixels in the 'pixel lookup view array'. This suggestion basically meant that for some images such as Distort Resize, Bluring, Convolution, Morphology, etc. you don't need to 'double process' the images to remove virtual pixel contributions from the resulting image.

Remember resize is written very differently to almost every other 'convolving' image processing operator as it does NOT include virtual pixels in its results.

[NicolasRobidoux]

magick: My suggested fixes of

Code: Select all

gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);

do not make the result continuous, because "your" formula has two large jumps: one at MagickEpsilon (which my fixes took care of) and one at -MagickEpsilon (which I did not fix). The jumps are size abs(1/MagickEpsilon) - 1, that is, really big.

Let me briefly explain how I would fix the alpha issue. This being said, I may be missing some issues. GEGL does not fix this yet that I know of (I sidestepped the issue by programming a nonlinear monotone resampler) and VIPS ignores alpha issues, so my other code bases are no help here. MyPaint, I believe, carefully deals with alpha issues, but I don't, in detail, know how it does it. Maybe some other software packages could offer guidance?

In other words, I'm not really sure I know enough about how everything fits together in IM, or about what the standard solution is, to be sure that my fix is "the" fix.

Gist: Clamp alpha to its natural range (which I assume is [0,1], although I understand that it may make sense to use a range related to the Quantum) in whatever computation you do, and if you use its reciprocal, clamp it to [MagickEpsilon,1]. In this context, this gives, for non-monotone interpolators:

Code: Select all

      /* This gamma is the interpolated alpha. */
      gamma=cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*alpha[2]+
        cx[3]*alpha[3])+cy[1]*(cx[0]*alpha[4]+cx[1]*alpha[5]+
        cx[2]*alpha[6]+cx[3]*alpha[7])+cy[2]*(cx[0]*alpha[8]+
        cx[1]*alpha[9]+cx[2]*alpha[10]+cx[3]*alpha[11])+cy[3]*(
        cx[0]*alpha[12]+cx[1]*alpha[13]+cx[2]*alpha[14]+cx[3]*alpha[15]);
      /* Clamp overshoot */
      gamma=(gamma>(MagickRealType) 1. ? (MagickRealType) 1. : gamma);
      /* Make sure you are not negative or too close to 0 */
      gamma=(gamma<(MagickRealType) MagickEpsilon ? (MagickRealType) MagickEpsilon : gamma);
      /* Take the reciprocal of the interpolated transparency */
      gamma = (MagickRealType) 1./gamma;

The casts are to prevent time wasting conversions to and from float and double (if MagickRealType is float), although those involving MagickEpsilon may be a waste since it would make sense that this global be declared as a MagickRealType (and I realize that modern compilers/chips may be smart enough not to do that).

-----

Another side effect is that as things are right now the reciprocated alpha can "flip" halos, that is, it can replace "clipped wave ripples" by "skipping stones" (a bit like plotting the absolute value of the cosine function, clipping the result just above zero; think of this in the context of haloing). The reason for this is that gamma can be negative.

-----

Another comment: The fabs is unnecessary in the current code for bilinear because bilinear is a monotone method: gamma can't be negative.

[magick]

Ok, I'll get your bicubic in the Subversion trunk. You / Anthony are welcome to make changes once its operational.

[NicolasRobidoux]

I thought about it some more, and this is what I think would be a good overall solution:

The way I understand alpha in the context of resampling is a little bit like what physicists/engineers call an "envelope".
The solution I proposed for this issue for GEGL is that one not necessarily use the same resampling method for the alpha channel: Always use a monotone method (actually, monotone is not required: it is enough to be range bound to [0,1], and "local boundedness" does the job), trying to choose a good match. (There has to be a standard solution but I've not run across it.)

What I propose that I implement is this:

bilinear interpolation uses bilinear alpha (like now)

quadratic B-spline smoothing uses quadratic B-spline smoothing alpha (quadratic is a really fast scheme: 9 point stencil)

Catrom interpolation uses cubic B-spline smoothing alpha

If more than 3 schemes (besides nearest) are desired, they could be

Mitchell interpolation uses cubic B-spline smoothing alpha

cubic B-spline smoothing uses cubic B-spline smoothing alpha

Note: I'm totally happy implementing less as opposed to more. This is the full menu, including what I believe are good choices. I don't use these methods, so I don't particularly care what is on order.

[anthony]

magick: My suggested fixes of

Code: Select all

gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);

do not make the result continuous, because "your" formula has two large jumps: one at MagickEpsilon (which my fixes took care of) and one at -MagickEpsilon (which I did not fix). The jumps are size abs(1/MagickEpsilon) - 1, that is, really big.

Actually the input 'gamma' is just an addition of weights. what happens near zero is that you are dealing with a fully transparent pixel, and that the result does not matter. What happens is that the color does jump -- but color is goes to zero and produces fully-transparent black. That is the discontinuity. Basically fully-transparent pixels do not matter so just make it black.

The reason for this code is purely to prevent division by zero! Mind you 'MagickEpsilon' could probably be much smaller than 1e-10 for IM Q16.

Gist: Clamp alpha to its natural range (which I assume is [0,1], although I understand that it may make sense to use a range related to the Quantum) in whatever computation you do, and if you use its reciprocal, clamp it to [MagickEpsilon,1]. In this context, this gives, for non-monotone interpolators:

Code: Select all

      /* This gamma is the interpolated alpha. */
      gamma=cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*alpha[2]+
        cx[3]*alpha[3])+cy[1]*(cx[0]*alpha[4]+cx[1]*alpha[5]+
        cx[2]*alpha[6]+cx[3]*alpha[7])+cy[2]*(cx[0]*alpha[8]+
        cx[1]*alpha[9]+cx[2]*alpha[10]+cx[3]*alpha[11])+cy[3]*(
        cx[0]*alpha[12]+cx[1]*alpha[13]+cx[2]*alpha[14]+cx[3]*alpha[15]);
      /* Clamp overshoot */
      gamma=(gamma>(MagickRealType) 1. ? (MagickRealType) 1. : gamma);
      /* Make sure you are not negative or too close to 0 */
      gamma=(gamma<(MagickRealType) MagickEpsilon ? (MagickRealType) MagickEpsilon : gamma);
      /* Take the reciprocal of the interpolated transparency */
      gamma = (MagickRealType) 1./gamma;

The casts are to prevent time wasting conversions to and from float and double (if MagickRealType is float), although those involving MagickEpsilon may be a waste since it would make sense that this global be declared as a MagickRealType (and I realize that modern compilers/chips may be smart enough not to do that).

It is probably over complicated. Whatever it should be, it will be added to a very large number of operators!
from composition though to convolution, morphology, and much more.

Another comment: The fabs is unnecessary in the current code for bilinear because bilinear is a monotone method: gamma can't be negative.

that may be the case, but as I said it is common code in a LOT of places, that often gets searched/replaced for updates and changes. Making it different for purely positive cases makes it harder to maintain. Typically by Cristy.

[NicolasRobidoux]

Opinion:
Ideally,

Code: Select all

gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);

should be replaced by

Code: Select all

gamma=(gamma>(MagickRealType) 1. ? (MagickRealType) 1. : gamma);
gamma=(MagickRealType) 1./(gamma<(MagickRealType) MagickEpsilon ? (MagickRealType) MagickEpsilon : gamma);

with the first line omitted when the method used to interpolate transparency is monotone (which I am suggesting should be always, with the consequence that alpha would be interpolated with cubic or quadratic B-spline smoothing or even bilinear when interpolating RGB with CatRom).
In other words, there never is a need for a fabs within the replacement line of code, which is the same for monotone and non-monotone methods. In addition, fabs most likely will slow down the code, esp. when MagickRealType is float (if this is compiled against the C math library, not the C++ math library, because its use requires an additional cast back to float... which was already there because 1. was not cast to MagickRealType, but is not in my version).

[magick]

Based on these discussions, I changed the following defines:

• #define MagickEpsilon ((MagickRealType) 1.0e-16)

and added a new inline method:

Code: Select all

static inline MagickRealType MagickReciprocal(const MagickRealType x)
{
  return((MagickRealType) 1.0/(((x) > (MagickRealType) 0.0 ? (x) : -(x)) < MagickEpsilon ? MagickEpsilon : x));
}

MagickReciprocal() is designed to avoid divide by zero and to use multiply rather than divide which is typically slower on modern computers. Individual algorithms will need to clamp to [MagickEpsilon,1] since the general case stuffs up a number of existing algorithms in IM (e.g. adaptive sharpen).

[NicolasRobidoux]

Suggestion: Call this MagickAbsReciprocal or MagickPositiveReciprocal. The name is then fully descriptive.

[NicolasRobidoux]

Another suggestion on the same theme: I feel rather strongly that spurious negative values of a quantity, like alpha (transparency) that is used nonlinearly and that is, as far as the eventual product goes, clamped to nonnegative values, should be assigned, with any operation, what 0 is assigned. In the case of taking reciprocals, this gives

Code: Select all

static inline MagickRealType MagickPositiveReciprocal(const MagickRealType x)
{
 return(x < MagickEpsilon ? (MagickRealType) 1.0/MagickEpsilon : (MagickRealType) 1.0/x);
}

(I used a name that I find specific, but of course this is MagickReciprocal.) Besides being faster (no nested branching; the first quotient will be computed at compile time hence costs nothing), this does exactly that: The "magick reciprocal" of -.1 will be exactly the same as the "magick reciprocal" of 0, namely 10^16. This should prevent unexpected oscillations to some extent.
I am not sure that I agree with the principle that says that it's OK to put inferior values in the colour channels when the transparency is, for example, 0. Looks to me like something that should be avoided if possible, otherwise some users will get a bad surprise. The above is safer (within the constraints that I have not tracked down everywhere this is used).
Generally, unless negative values are meaningful, anything that takes a reciprocal should not only stay "away" from zero, but actually negative values as well.
(Ignore if I'm missing something.)

[NicolasRobidoux]

Trying to figure out how to deal with the alpha channel, I am starting to think that the best solution is to only use monotone methods, e.g. bilinear with bilinear, Hermite with Hermite, quadratic B-spline with quadratic B-spline, but cubic B-spline with cubic B-spline, Mitchell and Catmull-Rom.
Overshoots and oscillations should be more deadly when something is used nonlinearly.

[anthony]

Trying to figure out how to deal with the alpha channel, I am starting to think that the best solution is to only use monotone methods, e.g. bilinear with bilinear, Hermite with Hermite, quadratic B-spline with quadratic B-spline, but cubic B-spline with cubic B-spline, Mitchell and Catmull-Rom.
Overshoots and oscillations should be more deadly when something is used nonlinearly.

I don't follow this!

[NicolasRobidoux]

Trying to figure out how to deal with the alpha channel, I am starting to think that the best solution is to only use monotone methods, e.g. bilinear with bilinear, Hermite with Hermite, quadratic B-spline with quadratic B-spline, but cubic B-spline with cubic B-spline, Mitchell and Catmull-Rom.
Overshoots and oscillations should be more deadly when something is used nonlinearly.

The code which interpolates R, G and B also interpolates alpha because it multiplies the interpolated R, G or B by gamma=1/alpha.
It is my opinion that the method used to interpolate alpha should always be a monotone (no over/undershoots) method, even when the method used to R, G or B is a method with overshoots.
(Maybe "deadly" is an overstated overstatement. Let's say "damaging".)