Gamma errors in resizing filters?

[Dark Shikari]

An interesting read (http://www.4p8.com/eric.brasseur/gamma.html).

Any chance we could get a properly gamma-corrected Avisynth filter?

[MfA]

Not just resizing, anything that averages pixels ... so denoisers, deinterlacers, the works.

Hell, even video codec subpixel interpolation introduces errors (needing extra bits to encode) because of the nonlinear color space (although it probably won't make a huge difference on average).

Another nice site with some examples and conversion source code :
http://mysite.verizon.net/spitzak/conversion/index.html

[Wilbert]

They claim that the filters need to operate at a higher bit-depth. Luckely 16 bit will be available in v2.6, so that should be could enough.

[IanB]

Don't be mistaken, don't be mislead ...

2.6.0 will only have API infrastructure for future 16 and 32 bit formats. There will be no 16 bit code for any filters.

In the interim...
Levels(16, 1/2.2, 235, 0, 255, False) # Use PC range to avoid truncation
...Resize(640, 480)
Levels(0, 2.2, 255, 16, 235, False)This can cause banding with some images, a compromise is to lower the Gamma from 2.2 to trade quantization error for gamma error.

[Wilbert]

Don't be mistaken, don't be mislead ...

2.6.0 will only have API infrastructure for future 16 and 32 bit formats. There will be no 16 bit code for any filters.
I know there is no 16/32 bit code yet. But once the API infrastructure is ready for it, people can start coding. I can't wait to start with some :)

[2Bdecided]

Idiots question: can't the code work in 16-bits or whatever internally, as long as it read and writes 8-bit data from/to what's around it?

I.e. couldn't you make a gamma correct resize filter which slotted into AVIsynth?

(I know nothing).


btw, those pages are fascinating. With a lot of the features in photoshop-like programs which don't work quite as you expect, or are difficult to set properly - it's because of this error! I love the torture images in the first link to show just how bad it can get.

Cheers,
David.

[IanB]

@2Bdecided,

Yes it would be possible. The current resizers do the intermediate calculations with 16 bit arithmetic already.

The real problem is doing the gamma exponetiation quickly. The fastest method I have in my kit bag is to use precalculated lookup tables, while this is fairly quick as things go, it is significantly slower than the current resizer code.

And yes the plan for 16bit and 32bit per pixel channel colour spaces is for them to be gamma=1.0

[MfA]

You can do really fast exponentiation with int/float casting tricks (and log on the way back). In a pinch squaring is pretty close too :)

[IanB]

@MfA,

Yes, X*X is fast in MMX/SSE, which is great for gamma=2.0, but no use in the general case.

However for Sqrt(X) I don't have a MMX/SSE that is as fast as lookup table.

Any and all code fragment donations accepted :D

Damn Intel, I'd kill for PXLAT[BW] instructions or Movq mm0, [ebx+mm0] type of addessing.

[MfA]

Do you have IEEExplore access?

http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=595279

Lack of SIMD LUT is painful, but understandable ... multiple very narrow ports to the cache? Forget it.

PS. GPUs can do all this a lot faster of course ...

PPS. LUT is faster for square root than SSE RSQRTPS + RCPPS?

[Dark Shikari]

@MfA,

Yes, X*X is fast in MMX/SSE, which is great for gamma=2.0, but no use in the general case.

However for Sqrt(X) I don't have a MMX/SSE that is as fast as lookup table.

Any and all code fragment donations accepted :D

Damn Intel, I'd kill for PXLAT[BW] instructions or Movq mm0, [ebx+mm0] type of addessing.Sqrt(x)? I remember a fast bit of code Pengvado tossed me a while back that did it in about 10 clock cycles.

[IanB]

IEEExplore access?No.
LUT is faster for square root than SSE RSQRTPS + RCPPS?It's border line but generally 4 LUT's beat the SIMD Sqrt, and it only for gamma=2.0 not general values.
Sqrt(x)? I remember a fast bit of code Pengvado tossed me a while back that did it in about 10 clock cycles.Well don't horde it, share it.

[Dark Shikari]

Ah, nevermind, looked it up--turned out it was EXP2, not square root.

inline float fast_exp2_sse2(float x)
{
static const float ss_bias = 12582912.0; // 3<<22
static const float ss_ln2 = 0.693147180559945;
static const float ss_1 = 1.0;
static const float ss_0_5 = 0.5;
float t, u;
asm volatile (
"movaps %0, %1 \n\t"
"addss %3, %0 \n\t"
"movaps %0, %2 \n\t"
"subss %3, %0 \n\t" // round(x)
"pslld $23, %2 \n\t" // round(x) in the exponent
"subss %0, %1 \n\t" // frac(x)
"mulss %4, %1 \n\t" // frac(x)*ln2 = y
"movaps %1, %0 \n\t"
"mulss %1, %1 \n\t" // y*y
"addss %5, %0 \n\t" // 1+y
"mulss %6, %1 \n\t" // y*y*.5
"addss %1, %0 \n\t" // 1+y+y*y*.5
"paddd %2, %0 \n\t" // (1+y+y*y*.5)<<round(x)
:"+x"(x), "=x"(t), "=x"(u)
: "m"(ss_bias), "m"(ss_ln2), "m"(ss_1), "m"(ss_0_5)
);
return x;
}

Couldn't one use similar code to the "fast inverse square root" used in the Quake 3 engine?

[Jawed]

Presumably any filter, such as fft3dfilter, can work in linear space provided that it contains code at the start to linearise and code at the end to encode as gamma-2.2.

I'm also wondering if algorithms that use sum of absolute differences or sum or squares, would stand to gain.

Jawed

[MfA]

Code written by Chris Lomont for the casting tricks I mentioned can be found here :

http://www.lomont.org/Software/Misc/FloatTricks/FloatHack.zip

[sunitram]

Hi, I have an approximation to exp that might be useful here. I use this in a Java program which works quite well. Here is a working version in C (I dont have much experience in C so I am sure this can be written in a better way):

double fast_exp2(double a) {
int tmp = (int)(1048576.0*a + 1072632447);
double p = 0.0;
*(1 + (int *)&p) = tmp;
return p;
}

This approximates 2^x with an average precision of about 2.5% or so, and is extremely fast at it. It is basically the same as using a lookup table with 2048 values. The line where tmp is calculated is basically the same as

int tmp = (int)((1<<20) * a + (1023* (1<<20) - 60801));

For more info see this paper: http://citeseer.ist.psu.edu/schraudolph98fast.html
I have a compileable example here: http://martin.ankerl.com/wp-content/uploads/2008/02/exp.cpp
(compiles with g++ -O3 -fno-strict-aliasing exp.cpp)

Actually this technique can be easily used to do all kinds of approximations, e^x, log(x), a^b and so on.

[MfA]

Which is exactly what Blinn did in his paper a year before that one ;)

[Archimedes]

An interesting read (http://www.4p8.com/eric.brasseur/gamma.html).

Any chance we could get a properly gamma-corrected Avisynth filter?
As a workaround, you can try this script. Input is a RGB32 image. Requires MaskTools 2.

function PhotoResize(clip input, int Width, int Height, bool GammaCorrection) {
Width = (Width * input.Height / input.Width <= Height) ? Width : Height * input.Width / input.Height
Height = (Width * input.Height / input.Width <= Height) ? Width * input.Height / input.Width : Height
clip1 = input.Spline36Resize(Width, Height)
clip2 = (GammaCorrection == true) ? input.Levels(0, 1/2.2, 255, 0, 255).Spline36Resize(Width, Height).Levels(0, 2.2, 255, 0, 255) : clip1
difference = (GammaCorrection == true) ? Subtract(input.Levels(0, 1/2.2, 255, 0, 255).Levels(0, 2.2, 255, 0, 255), input).Spline36Resize(Width, Height) : clip1
GammaCorrection == true ? Subtract(clip2, difference) : clip1
clip1 = (GammaCorrection == true) ? clip1.RGBtoYV12() : last
clip2 = (GammaCorrection == true) ? last.RGBtoYV12() : last
GammaCorrection == true ? mt_clamp(clip2, clip1, clip1, 255, 0, u=3, v=3) : last
isYV12() ? YV12toRGB() : last

function RGBtoYV12(clip input) {
input.Crop(0, 0, -(input.Width % 2), -(input.Height % 2))
PointResize(last.Width * 2, last.Height * 2)
ConvertToYV12(matrix="pc.601")
}
function YV12toRGB(clip input) {
input.ConvertToRGB32(matrix="pc.601")
PointResize(last.Width / 2, last.Height / 2)
}
}

LanczosResize, incorrect scaling (GammaCorrection=false):
http://img89.imageshack.us/img89/8840/img19850512x0384lanczospi4.png

ImageMagick, correct scaling (with gamma correction), 16 bit:
http://img246.imageshack.us/img246/2780/img19850512x0384imagemasq0.png

LanczosResize, correct scaling (GammaCorrection=true) without error correction, 8 bit:
http://img99.imageshack.us/img99/8138/img19850512x0384lanczoszn8.png

LanczosResize, correct scaling (GammaCorrection=true) with error correction, 8 bit:
http://img263.imageshack.us/img263/9395/img19850512x0384lanczosyr5.png

[IanB]

@Archimedes,

As I pointed out above range scaling YUV images to PC levels as part of the Gamma conversion is very worthwhile, you should try it in addition to your error correction.Levels(16, 1/2.2, 235, 0, 255, False) # Use PC range to avoid truncation
...Resize(640, 480)
Levels(0, 2.2, 255, 16, 235, False)

[Didée]

@ IanB,

Archimedes is dealing with RGB input, so that kind of range expansion is not appliable.


@ Archimedes:

Oh, that's not the kind of error protection that I had thought of. But thinking about it, this is probably the better way, and ... then you can simplify it.

> mt_clamp(clip2, clip1, clip1, 255, 0, u=3, v=3)

is identical to, but slower as

> mt_lutxy(clip2, clip1, "x y > x y ?", u=3,v=3)

which in turn is identical to, but slower as

> mt_logic(clip2, clip1, "max", u=3,v=3)

Seeing this, it should be possible to completely ditch out masktools together with the 2*supersampled RGB->YY12 conversion, and just use

> Overlay(clip2, clip1, mode="Lighten", pc_range=true)

[Archimedes]

There was an error in my function. clip2 variable has to be set as follow.

clip2 = (GammaCorrection == true) ? last.RGBtoYV12() : last

Regarding Overlay. I first used it, but it seems, Overlay doesn't work well with RGB input. First, it needs RGB24 input (not really a problem), second, images looses color informations (really a problem) when using it. So i used the MaskTools.

[Wilbert]

Regarding Overlay. I first used it, but it seems, Overlay doesn't work well with RGB input. First, it needs RGB24 input (not really a problem), second, images looses color informations (really a problem) when using it.
Yes we know. It's fixed in the CVS (wait for the next v2.58).

[IanB]

Also, Overlay() works in YUV colour space, so there will be a conversion (fixed in next 2.5.8) from and to RGB colour space with the associated small loss of data.

Many of the functions of Overlay() are available in the older internal Layer() filter, which works directly with YUY2 and RGB32.

[Archimedes]

Also, Overlay() works in YUV colour space, so there will be a conversion (fixed in next 2.5.8) from and to RGB colour space with the associated small loss of data.
Exactly that was my sneaking suspicion.