Bicubic magnification - how?

[Rolf N]

It's frustrating. Hitting "bicubic filter" into Google returns tons of related conversation, even a link to these very forums, where, as usual, everyone agrees that bicubic is a great -- if not the best practical -- approach to magnification. But I can't find any straight information, just lots of talk.

I've found a link to an NVIDIA demo that does some mutated variant in a pixel shader.
I've found discussions about HP printers.
I've found millions of discussions about image manipulation software.

What I'm looking for is a method, or an algorithm if you will, that can scale a source image to quadruple its original size, and fits the description "bicubic magnification filter".

This is about textures, so the source images will always have power-of-two dimensions, and let's assume texture repeats (=>unconstrained coords), if that helps.

I'm not asking for code donations
It would be perfectly enough if someone could just describe the process for, say, a single component image. Something along these lines maybe:
Bilinear 2x2 upsampling - crummy conjecture by me, that would help me tremendously if it were for bicubic

dest(0;0):=  0.75*(0.75*src(0;0)+0.25*src(0;1))+
             0.25*(0.75*src(-1;0)+0.25*src(-1;1));
dest(1;0):=  0.75*(0.75*src(0;0)+0.25*(src(0;1))+
             0.25*(0.75*src(1;0)+0.25*src(1;1));
dest(1;-1):= 0.75*(0.75*src(0;0)+0.25*(src(0;-1))+
             0.25*(0.75*src(1;0)+0.25*src(1;-1));
dest(0;-1):= 0.75*(0.75*src(0;0)+0.25*(src(0;-1))+
             0.25*(0.75*src(-1;0)+0.25*src(-1;-1);

Please

 

[Humus]

float cerp(const float u0, const float u1, const float u2, const float u3, float x){
	float p = (u3 - u2) - (u0 - u1);
	float q = (u0 - u1) - p;
	float r = u2 - u0;
	return x * (x * (x * p + q) + r) + u1;
}

Cubic interpolation in one direction:

float result = cerp(pic[y][x-1], pic[y][x], pic[y][x+1], pic[y][x+2], t);

t is the interpolation factor. 0 returns pic[y][x], 1 returns pic[y][x+1].

Bicubic interpolation is simply that interpolation in two directions:

float res0 = cerp(pic[y-1][x-1], pic[y-1][x], pic[y-1][x+1], pic[y-1][x+2], tx);
float res1 = cerp(pic[y  ][x-1], pic[y  ][x], pic[y  ][x+1], pic[y  ][x+2], tx);
float res2 = cerp(pic[y+1][x-1], pic[y+1][x], pic[y+1][x+1], pic[y+1][x+2], tx);
float res3 = cerp(pic[y+2][x-1], pic[y+2][x], pic[y+2][x+1], pic[y+2][x+2], tx);

float result = cerp(res0, res1, res2, res3, ty);

 

[DemoCoder]

I love Horner's rule.

 

[Rolf N]

Thanks, Humus!
*big hug*

 

[KimB]

What sort of cubic interpolation does that do, Humus? There are two ways to do it, of course.

One is to match the cubic function to all four sample points.

Another is to match the cubic function only to the nearest two sample points, and also match the derivative of the function at the nearest two points, with some measure of slope (easiest is just to simply get the slope from two other sample points, i.e. if we are looking at point 2, take the slope of the line from point 1 to 3 as the slope at point 2).

Additionally, you can do the cubic interpolation via matrix multiplication relatively easily. It's been a while since I did the math, but I think it worked out to fewer operations than doing the interpolation in the way Humus showed it....

 

[Simon F]

Zeckensack: Do you want "interpolating" or "approximating" curves? The former means that the upscaled signal contains the original points (eg Bilinear will do this) while "approximating" just means that the upscaled signal will only "tend to follow" the original signal.

The problem with "interpolating" curves is that with anything >Bilinear, you'll get over/undershoot which means you'll have to clamp your results.

If you want "interpolating" look for "Catmull-Rom Splines" otherwise try "Uniform BSplines" if you are happy with an approximating curve.

 

[Simon F]

Additionally, you can do the cubic interpolation via matrix multiplication relatively easily.

You can express it mathematically that way but don't implement it that way. It's both inefficient and mathematically unstable.

 

[Simon F]

Zeck',
Just re-read your post and you're just after, say, 2x2 scaling? As Humus said, Bicubic is a separable filter which means you can do the X and Y passes independently which makes it much cheaper.


Let's say you want to do the X pass and you have a row of pixels:

 .... P[i]  P[i+1]  P[i+2]  P[i+3] ....

and want to manufacture twice as many pixels...

 .... P'[i]   Q[i]   P'[i+1]   Q[i+1] ... etc ...

Then for a Catmull-Rom spline (/me has a rummage through some old code) you have

P'[i] = P[i];
and 
Q[i] = (9 * (P[i] + P[i+1]) - (P[i-1] + P[i+2]))/16; //Repeat for RG&B channels

You'll have to check for over/underflow with the Qs and you'll probably need to repeat the start and end pixels at the begining and end of each row/column.

 

[MfA]

If you've got ACM access read this paper. Mitchell & Netravali is pretty much what most people think about when talking about cubic spline interpolation in image processing and computer graphics (their method is also commonly referred to as bc-splines).

The continuous B-spline methods can do better in a MSE sense, but you can forget about implementing them in a shader.

 

[Rolf N]

Zeckensack: Do you want "interpolating" or "approximating" curves? The former means that the upscaled signal contains the original points (eg Bilinear will do this) while "approximating" just means that the upscaled signal will only "tend to follow" the original signal.

I suppose I want an "approximating" fit. This is what even (my idea of) bilinear upscaling does, and I believe that's the way to go. Let me explain.

The undesirable thing about "my" bilinear upscaling is that it's not reversible. If I take the 4 texels generated by upscaling and boxfilter them back, the result may be different from the single source texel I originally had. But upscaling by just inserting new texels has another, and IMO much more severe drawback: the new, larger mipmap would have a "drift" towards the origin.

I'm a big fan of the theory that says that the bottom left texel in a texture map is actually a point at coords (0.5;0.5). If I overlay an insert-upscaled version, and scale it to fit, that point (exact same color) would have moved to (0.25;0.25). To the point at (0.75;0.25) and to its two other friends, which are just as close as the one at (0.25;0.25), the "original" point would only contribute 50 per cent. I think that's just wrong. The original signal texel should IMO contribute the same amount to its four "new neighbours". They are the same distance away.

The problem with "interpolating" curves is that with anything >Bilinear, you'll get over/undershoot which means you'll have to clamp your results.

Thanks for pointing that out, but it wouldn't be much of a problem. I'll be using MMX anyway, and that makes saturation completely free.

If you want "interpolating" look for "Catmull-Rom Splines" otherwise try "Uniform BSplines" if you are happy with an approximating curve.

But that's not the same thing as bicubic, or is it? One step at a time, please

Zeck',
Just re-read your post and you're just after, say, 2x2 scaling? As Humus said, Bicubic is a separable filter which means you can do the X and Y passes independently which makes it much cheaper.

Yes, it seems obvious that there are a few shared terms that can be reused. I'm actually going to do the exact same weights as in my bilinear conjecture.

But I won't just insert. Sorry

You'll have to check for over/underflow with the Qs and you'll probably need to repeat the start and end pixels at the begining and end of each row/column.

I will try to do the right thing and either "clamp to edge" or "repeat", depending on texture state

Many thanks for the input.
*bows*

 

[KimB]

Additionally, you can do the cubic interpolation via matrix multiplication relatively easily.

You can express it mathematically that way but don't implement it that way. It's both inefficient and mathematically unstable.

Well, looking at the number of operations, it looks like while the matrix multiplication method had a few more operations, it was more balanced in terms of multiplication and addition ops, and so, I would think, could end up being higher in performance. At least, for the case where you take the interpolating function to go through all four sample points....it didn't look as easy to align the ops in the other interpolation scheme I mentioned.

Edit:
Btw, when I looked at the math, I figured out that you really don't need to do a full matrix multiplication. Since you only want one of the final four values you'd get from the multiplication (the interpolated value), you can instead get by with a four-vector dot product for producing each sample point. The thing that causes this technique to require more ops is the generation of the four-vector to dot the four samples with. This does become less of a problem if the four-vector is reused (which it is when you do bicubic, since you do four bicubic samples in one direction with the same x (or y, of course) value). So, I suppose, you could, if you wanted to, create a lookup table for these sampling weights for even more efficient bicubic filtering....

 

[MfA]

I'm a big fan of the theory that says that the bottom left texel in a texture map is actually a point at coords (0.5;0.5).

Well .5/N anyway, that isnt so much a theory as it is exactly what DirectX does AFAIR. From a resizing point of view a bad idea, you will always get less distortion by keeping samples colocated between scales during resizing no matter how well you do it.

 

[Rolf N]

I'm a big fan of the theory that says that the bottom left texel in a texture map is actually a point at coords (0.5;0.5).

Well .5/N anyway, that isnt so much a theory as it is exactly what DirectX does AFAIR. From a resizing point of view a bad idea, you will always get less distortion by keeping samples colocated between scales during resizing no matter how well you do it.

The intention is to play along as nicely as possible with fixed function mipmap filters, and those are boxes around texel centers, not boxes around texel corners. I really don't want surface detail moving around depending on distance.

If OTOH I were aiming for a Photoshop plugin, I'd agree instantly.

 

[Simon F]

Zeckensack: Do you want "interpolating" or "approximating" curves? The former means that the upscaled signal contains the original points (eg Bilinear will do this) while "approximating" just means that the upscaled signal will only "tend to follow" the original signal.

I suppose I want an "approximating" fit. This is what even (my idea of) bilinear upscaling does, and I believe that's the way to go. Let me explain.

The undesirable thing about "my" bilinear upscaling is that it's not reversible. If I take the 4 texels generated by upscaling and boxfilter them back, the result may be different from the single source texel I originally had. But upscaling by just inserting new texels has another, and IMO much more severe drawback: the new, larger mipmap would have a "drift" towards the origin.

Ahhh! You want to go look at wavelets then, in particular, the "lifting" process. I was experimenting using Catmull-Rom as the basis for the filtering in such a wavelet and it all worked fine.

I'm a big fan of the theory that says that the bottom left texel in a texture map is actually a point at coords (0.5;0.5). ....

Yes that's fine and is where HW will assume it lies as well. It's all because Williams' MIP mapping paper used box filtering.

BTW, you don't have to use a box filter to do the down sampling - It's cheap and cheerful, but it's by no means the best method.

The problem with "interpolating" curves is that with anything >Bilinear, you'll get over/undershoot which means you'll have to clamp your results.

Thanks for pointing that out, but it wouldn't be much of a problem. I'll be using MMX anyway, and that makes saturation completely free.

Fair enough - it was only meant as something to watch out for.

If you want "interpolating" look for "Catmull-Rom Splines" otherwise try "Uniform BSplines" if you are happy with an approximating curve.

But that's not the same thing as bicubic, or is it? One step at a time, please

Catmull-Rom Splines are cubic and BSplines can be made to be any order - I was assuming cubic.

 

[Simon F]

If you've got ACM access read this paper. Mitchell & Netravali is pretty much what most people think about when talking about cubic spline interpolation in image processing and computer graphics (their method is also commonly referred to as bc-splines).

Interesting to see that Catmull-Rom falls inside their "satisfactory" region.

 

[Rolf N]

The problem with "interpolating" curves is that with anything >Bilinear, you'll get over/undershoot which means you'll have to clamp your results.

Thanks for pointing that out, but it wouldn't be much of a problem. I'll be using MMX anyway, and that makes saturation completely free.

Fair enough - it was only meant as something to watch out for.

And I will. That was worded a little strangely, what I was actually trying to say was that clamping won't be expensive or hard to do. But I wasn't aware of that issue at all, so it was indeed useful info.