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);
}
