search for a decent interpolatory EWA filter

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
search for a decent interpolatory EWA filter
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 noop).
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 BCsplines 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 BCspline 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.
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 BCsplines 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 BCspline 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.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
Uhmmm! This one certainly is not going to be made "nice" without some kicking and screaming.
Last edited by NicolasRobidoux on 20111217T09:17:1707:00, edited 1 time in total.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
Not sure there is a member of the family that works well, but
is a neat effect.
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
Last edited by NicolasRobidoux on 20111223T10:19:3607:00, edited 3 times in total.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
It does not look like any of them is a high quality scheme for upsampling.
Last edited by NicolasRobidoux on 20111217T09:18:0907:00, edited 1 time in total.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
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 BCspline a Keys cubic gives B = (9  3 sqrt(2))/7, C = (2  3 sqrt(2))/14:
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).
Best so far:
Making the cardinal BCspline 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
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).

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
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 noop (which I imagine you don't since your orthogonal default is Mitchell which is not the identity under noop), this would be my first choice at this point.
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 noop (which I imagine you don't since your orthogonal default is Mitchell which is not the identity under noop), this would be my first choice at this point.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
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
Re: search for a decent interpolatory EWA filter
Hi , just chipping in on this , catmullrom spline is used for interpolation in gimp. By definition it passes though all four reference points so would seem to fit your noop 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 catmulrom. sorry for the noise.
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 catmulrom. sorry for the noise.
Last edited by imaggie on 20111219T14:00:5007:00, edited 2 times in total.
 fmw42
 Posts: 23675
 Joined: 20070702T17:14:5107:00
 Authentication code: 1152
 Location: Sunnyvale, California, USA
Re: search for a decent interpolatory EWA filter
To be precise, jinc uses Bessel functions.unlike the beziers in jinc et al.
Last edited by fmw42 on 20111220T11:36:1307:00, edited 1 time in total.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
Quick reply:imaggie wrote:Hi , just chipping in on this , catmullrom spline is used for interpolation in gimp. By definition it passes though all four reference points so would seem to fit your noop 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 catmulrom. sorry for the noise.
 Indeed, CatmullRom is interpolatory, and consequently satisfies the "noop" 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 CatmullRom cubic is not interpolatory when used as an EWA filter.
 Indeed, CatmullRom is fast.
 FYI, I've written the code that does CatmullRom 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 CatmullRom overshoots. Not programmed yet. Let's say that I am a little more than familiar with CatmullRom.
 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
 Posts: 8870
 Joined: 20040531T19:27:0307:00
 Authentication code: 8675308
 Location: Brisbane, Australia
Re: search for a decent interpolatory EWA filter
The problem here is your BC curve has stronge negative weights in it!NicolasRobidoux wrote:Note sure there is a member of the family that works well, butis a neat effect.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
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!
Anthony Thyssen  Webmaster for ImageMagick Example Pages
http://www.imagemagick.org/Usage/
http://www.imagemagick.org/Usage/

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
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 Bspline 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.
Math:
If one enforces that the value of the Bspline 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.
Last edited by NicolasRobidoux on 20111225T15:42:1907:00, edited 3 times in total.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
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 BCsplines (unless you absolutely need the interpolatory property, obviously).
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 BCsplines (unless you absolutely need the interpolatory property, obviously).
Last edited by NicolasRobidoux on 20111225T15:43:0007:00, edited 1 time in total.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
Finally got around to programming the "compute the max error reconstructing a linear gradient with unit slope parallel to one of the axes (the xaxis, 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.
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.

 Posts: 1944
 Joined: 20100828T11:16:0007:00
 Authentication code: 8675308
 Location: Montreal, Canada
Re: search for a decent interpolatory EWA filter
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 BCsplines 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 CatmullRom.
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 halfsize 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:Compare, say, with nearest neighbour:
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:
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 CatmullRom.
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 halfsize 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
Code: Select all
convert rose: filter Point resize 3000% rose_o_point.png
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 BCspline 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/LICENSE2.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,xxpt,yypt);
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);
}
Last edited by NicolasRobidoux on 20111227T08:33:4207:00, edited 3 times in total.