search for a decent interpolatory EWA filter

[NicolasRobidoux]

I've been thinking (and tinkering) on and off about how to construct a high quality EWA filter that would be interpolatory (that is, which would leave pixel values alone under no-op).

Triangle (the "cone tent" filter in the EWA context) piles on the artifacts unless the image is smooth or one is downsampling in all directions. It's definitely worse than orthogonal Triangle (= bilinear when upsampling).

Hermite barely smooths out the rough edges of Triangle.

I've also tried an even (single piece) polynomial with support exactly sqrt(2) and the least reasonable degree given the usual smoothness constraints, and it did not work so well.

-----

The starting point of the current attempt is the BC-splines rescaled so that the support extends to sqrt(2), that is, with blur = sqrt(2)/2.

The constraint of interpolation then is quite simple: If outer(x) is the value of the BC-spline between 1 and 2, we want to enforce outer(sqrt(2)) = 0.

This gives B = 3 sqrt(2) C.

Given that the blur and B are set, this reduces the search for "best" to one parameter (C).

My first attempt will consist of minimizing the max error when upsampling an affine gradient aligned with one of the axes. I'll need to modify some of my earlier code for this.

[NicolasRobidoux]

Uhmmm! This one certainly is not going to be made "nice" without some kicking and screaming.

[NicolasRobidoux]

Not sure there is a member of the family that works well, but

Code: Select all

convert rose: -define filter:blur=.7071067811865475 -define filter:c=1. -define filter:b=4.2426406871929 -filter Cubic -distort Resize 3000% rose_c_c_is_1.png

is a neat effect.

[NicolasRobidoux]

It does not look like any of them is a high quality scheme for upsampling.

[NicolasRobidoux]

Doing a quick sweep, none of them is very good for upsampling (although it should be noted that at this point nobody knows how to make a cardinal EWA filter that works really well for upsampling).

Best so far:
Making the cardinal BC-spline a Keys cubic gives B = (9 - 3 sqrt(2))/7, C = (2 - 3 sqrt(2))/14:

Code: Select all

convert rose: -define filter:blur=.7071067811865475 -define filter:c=.16018862050852 -filter Cubic -distort Resize 3000% rose_ewa_cardinal_keys.png

which is representative of the best results that be obtained with the family of cardinal EWA filters (based on my really quick check) and is a lot better than Triangle and Hermite EWA.

I suspect that the "optimized for minimal error on affine gradients..." will be slightly better, but not by much. (And I don't have time to compute it right now. I believe that the optimized one will be very close, maybe right at, the above Keys one.)

This being said, the RobidouxSharp discussed in another thread is definitely better if it is the "look and feel" of an interpolatory method that you want (without it actually being interpolatory).

[NicolasRobidoux]

Anthony:

TTBOMK, the above Keys Cubic (with blur = sqrt(2) / 2) is the best interpolatory EWA filter currently available EWA. Not that I or someone else may not manage to find a better one.

If you are not upsampling a lot in any direction (no more than something like 3x) and you want the default EWA to be the identity under no-op (which I imagine you don't since your orthogonal default is Mitchell which is not the identity under no-op), this would be my first choice at this point.

[NicolasRobidoux]

An interesting one. Again, nothing to write home about if upsampling a lot.

Code: Select all

convert rose: -define filter:blur=.7071067811865475 -define filter:b=2 -define filter:c=.471404520791 -filter Cubic -distort Resize 3000% rose_c_p471404520791.png

[imaggie]

Hi , just chipping in on this , catmull-rom spline is used for interpolation in gimp. By definition it passes though all four reference points so would seem to fit your no-op criterion.

It can tend to overshoot on steep changes but that is more or less a result of containing it to pass through all four points, so I guess you can't have it all.

It's also very fast to do being a relatively simple calculation, unlike the [strike]beziers [/strike] Bessels in jinc et al.

PS I've seen the other thread you are all well aware of the properties on catmul-rom. sorry for the noise.

[fmw42]

unlike the beziers in jinc et al.

To be precise, jinc uses Bessel functions.

[NicolasRobidoux]

Hi , just chipping in on this , catmull-rom spline is used for interpolation in gimp. By definition it passes though all four reference points so would seem to fit your no-op criterion.

It can tend to overshoot on steep changes but that is more or less a result of containing it to pass through all four points, so I guess you can't have it all.

It's also very fast to do being a relatively simple calculation, unlike the [strike]beziers [/strike] Bessels in jinc et al.

PS I've seen the other thread you are all well aware of the properties on catmul-rom. sorry for the noise.

Quick reply:

• Indeed, Catmull-Rom is interpolatory, and consequently satisfies the "no-op" criterion. It is, however, an "orthogonal" method. This thread has to do with Elliptical Weighted Averaging (EWA) methods (what you get with distort e.g. -distort resize). The Catmull-Rom cubic is not interpolatory when used as an EWA filter.
• Indeed, Catmull-Rom is fast.
• FYI, I've written the code that does Catmull-Rom for the VIPS library, and if I remember correctly, touched the GEGL code that does the same. One of my masters students just defended her thesis which, among other things, aims at reducing Catmull-Rom overshoots. Not programmed yet. Let's say that I am a little more than familiar with Catmull-Rom.
• EWA is not as cheap as an orthogonal method implemented so as to take advantage of the orthogonality. And Jincs indeed are not cheap (unless one uses tricks discussed in the thesis of my student Chantal Racette, or uses LUTs, which is standard).
• Noise is OK: It's the risk one takes trying to be helpful. Thank you for trying

[anthony]

Note sure there is a member of the family that works well, but

Code: Select all

convert rose: -define filter:blur=.7071067811865475 -define filter:c=1. -define filter:b=4.2426406871929 -filter Cubic -distort Resize 3000% rose_c_c_is_1.png

is a neat effect.

The problem here is your B-C curve has stronge negative weights in it!
The result is much like what you get due to the EWA resampling filter missing all pixels!
See the experts note in...
http://www.imagemagick.org/Usage/resize ... terpolated
ASIDE: There is a code switch in 'resample.c' called "DEBUG NO_PIXEL_HIT" that will replace the default action of using a 'interpolated color' with a 'red' color for debugging. Though in this case you can never miss pixels, but get negative pixel colors instead!

Graph the filter...

Code: Select all

   convert null:  -define filter:blur=.7071067811865475 \
              -define filter:c=1. -define filter:b=4.2426406871929 -filter Cubic \
              -define filter:verbose=1 -distort SRT 0 null: > data.dat
  gnuplot
   set grid
   plot "data.dat" with lines
    set term png
    set output "filter_wierd.png"
    replot

Much of the weirdness is caused by the the negative weight weightings at 0 distance from the sample point
As such when you get a pixel in that negative area, you get weird effects. I do not have a check for a negative
weighted result! As such I think we are seeing an programmed 'integer modulus' effect

Actually I rather like what you get with a blur of 0.84 Its a very funny effect!

[NicolasRobidoux]

Indeed, I explored quickly without thinking too hard. And I should have taken the time to look at the plots.

Math:

If one enforces that the value of the B-spline at 0 be > 0, then one must have that C<1/sqrt(2)=.7071... which certainly eliminates the "neat effect" ones right off the bat.

[NicolasRobidoux]

I think I'm not going to be able to escape programming the "minimize the affine gradient error" code. Manana.

And it's not going to be a stellar upsampling scheme no matter what. Robidoux and RobidouxSharp (and Lanczos and LanczosSharp) are almost certainly (99% certain) better than any interpolatory EWA built with BC-splines (unless you absolutely need the interpolatory property, obviously).

[NicolasRobidoux]

Finally got around to programming the "compute the max error reconstructing a linear gradient with unit slope parallel to one of the axes (the x-axis, actually)" C program.

With blur = 1/sqrt(2), the cardinal cubic which is also a Keys cubic gives a max error of about .174. The Keys one is not the minimizer within the family of cardinal ones.

[NicolasRobidoux]

Finally, here is a pretty good approximation to the optimal pair. Although this may be the best interpolatory EWA filter ever devised, I can't help to think that it must be possible to put one together that is even better. (Actually, I now think that the BC-splines rescaled by 1/sqrt(2) is the wrong family to optimize over.) In particular, if you are fine with letting go of being interpolatory, there are much better EWA methods: Robidoux, RobidouxSharp, Lanczos and LanczosSharp, for example.

The (approximate) optimal pair (with respect to reconstructing affine gradients with gradient parallel to one of the axes): C = 0.49257366 and B = 3 sqrt(2) C = 2.089813051319261.

When reconstructing z(x,y) = a + x (or z(x,y) = a + y) with a a constant, the maximum error over the entire plane is approximately .141622596233478. The C code which verifies this is at the end of this post. Compare with a maximum error of .5 when using Nearest Neighbour, and of 0 for (orthogonal) bilinear or Catmull-Rom.

The results of enlarging a lot with this method are qualitatively different to those obtained with any standard method I know of. Unless the image is very smooth, the look and feel is that of slightly blurred half-size square tiles which are laid out diagonally, with corners at the original pixel positions and "rounded" just enough to get the interpolation property. The good news is that this scheme is almost halo free. Example:

Code: Select all

convert rose: -define filter:blur=.7071067811865475 -define filter:c=.49257366 -define filter:b=2.089813051319261 -filter Cubic -distort Resize 3000% rose_c_optimal.png

Compare, say, with nearest neighbour:

Code: Select all

convert rose: -filter Point -resize 3000% rose_o_point.png

If you don't believe me that this is the best interpolatory EWA method so far, try your interpolatory EWA challenger with the above 3000% enlargement task. (Please let me know if I happen to be wrong.)

If there is interest in making this a mainstream method (and honestly I don't understand why there should be), I'll fully converge the coefficients (when I've installed a decent OS of my new, faster, computers). My opinion, however, is that this is not something that should be used unless one absolutely wants an interpolatory EWA method. It works well for mild enlargements and very mild reductions, but I don't recommend it for anything else.

The other thing that needs to be double checked (requires a trivial mod of the following C code) is that the sum of the weights is positive everywhere. I'm pretty sure it is but I've not checked this directly (yet).

The optimization code:

Code: Select all

#include <stdio.h>
#include <math.h>

/*
 * ewaCardinal
 *
 * Program estimating the maximum EWA (Elliptical Weighted Averaging)
 * resampling error for raster image data which is equal to the column index,
 * when the filter kernel is a BC-spline with B = 3 sqrt(2) C and the
 * extend is shrunk to sqrt(2) from the standard 2, which corresponds to
 * an ImageMagick "blur" set to .5 sqrt(2).
 *
 * The unscaled spline kernel is written in the following form:
 *
 * when 0 <= x < 1, it is
 *
 *   1 - sqrt(2) C +
 *   x^2 ( -3 + ( 6 sqrt(2) + 1 ) C + x ( 2 - ( 2 + 9 sqrt(2) ) / 2 C ) )
 *
 * when 1 <= x < 2, it is
 *
 *    ( ( - sqrt(2) - 2 ) C / 2 ) ( x - 2 )^2
 *    ( x - ( 2 sqrt(2) + 2 ) / ( sqrt(2) + 2 ) )
 */

/*
 * Copyright 2011 Nicolas Robidoux
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

/*
 * Compile with
 *
 * gcc -O2 -march=native ewaCardinal.c -lm -o ewaCardinal
 */

/*
 * MIN_C is the smallest value of C to be tested.
 */               
#define MIN_C (0.49257365)

/*
 * MAX_C is the largest value of C to be tested.
 */
#define MAX_C (0.49257367)

/*
 * C_SAMPLING is one less than the number of evenly spaced values of C
 * to be tested.
 */
#define C_SAMPLING 2

/*
 * X_SAMPLING is one less than the number of samples of the error, in
 * both the horizontal and vertical direction, to be taken, for each
 * tested value of C, in the square [0,1/2]x[0,1/2]. The max error on
 * this square is the max error over the entire plane.
 */
// #define X_SAMPLING 65536
#define X_SAMPLING 16384
// #define X_SAMPLING 5048

static inline double fin(const double c,
                         const double x)
{
     const double result =
          1.
          -
          sqrt(2.) * c
          +
          x * x * ( -3.
                    + ( 6. * sqrt(2.) + 1. ) * c
                    + x * ( 2. - ( 2. + 9. * sqrt(2.) ) * .5 * c ) );
     return(result);
}

static inline double fout(const double c,
                          const double x)
{
     const double x_minus_2 = x - 2;
     const double result =
          ( -0.5 * ( 2. + sqrt(2.) ) )
          * c
          * ( x - ( 2. * sqrt(2.) + 2. ) / ( sqrt(2.) + 2. ) )
          * x_minus_2 * x_minus_2;
     return(result);
}
 
static inline double pt(const int i)
{
     const double result = i + .5;
     return(result);
}

static inline double bicubic(const double c,
                             const double x)
{
     double result = 0;
     if ( x < 1. ) {
          result = fin(c,x);
     }
     else if ( x < 2. ) {
          result = fout(c,x);
     }
     return(result);
}

static inline double radial(const double c,
                            const double x,
                            const double y)
{
     const double rsq = x*x+y*y;
     const double result =
          bicubic(c,sqrt(rsq+rsq));
     return(result);
}

int main()
{
     const double ch = ( MAX_C - MIN_C) / C_SAMPLING;
     const double h = .5 / X_SAMPLING;

     double best = 1;
     double best_c = MIN_C;

     int k = 0;
     do {
          double worst = 0.;
          const double c = MIN_C + k * ch;

          int i = 0;
          do {
               const double x = i * h;

               int j = 0;
               do {
                    const double y = j * h;
                    double numerator = 0;
                    double denominator = 0;
                    int l = -1;
                    do {
                         const double xpt = pt(l);
                         int m = -1;
                         do {
                              const double ypt = pt(m);
                              const double weight = radial(c,x-xpt,y-ypt);
                              numerator += xpt * weight;
                              denominator += weight;
                         } while ( ++m < 2 );
                    } while ( ++l < 2 );
                    {
                         const double error = numerator / denominator - x;
                         const double error_sq = error * error;
                         if ( error_sq > worst ) {
                              worst = error_sq;
                         }
                    }
               } while ( ++j <= X_SAMPLING );

          } while ( ++i <= X_SAMPLING );

          if ( worst < best ) {
               best = worst;
               best_c = c;
          }

     } while ( ++k <= C_SAMPLING );

     printf("least absolute error = %-17.15f\n", sqrt(best));
     printf("best C = %-17.15f\n", best_c);

     return(0);
}