madVR - high quality video renderer (GPU assisted)

[NicolasRobidoux]

Mathias:
What is your "anchor" colorspace?
Actually, what is your "anchor" luma/luminance channel.
What I most care about is whether it is a linear light channel, or whether it belongs to a perceptual color space.
-----
It's cheaper to do something like a "colorspace transformation" just to luma/luminance, yes? Does it cost a lot to convert to linear light, or is it one of these things that are basically free? Are colorspace transformations implemented as LUTs or with computational code?
-----
This is the key code in ImageMagick (actually, the IM6 code computes a LUT; what I'm showing is the IM7 code). This is way too complicated for all sorts of reasons. But it gives the idea, hopefully.

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%              EEEEE  N   N  H   H   AAA   N   N   CCCC  EEEEE                %
%              E      NN  N  H   H  A   A  NN  N  C      E                    %
%              EEE    N N N  HHHHH  AAAAA  N N N  C      EEE                  %
%              E      N  NN  H   H  A   A  N  NN  C      E                    %
%              EEEEE  N   N  H   H  A   A  N   N   CCCC  EEEEE                %
%                                                                             %
%                                                                             %
%                    MagickCore Image Enhancement Methods                     %
%                                                                             %
%                              Software Design                                %
%                                John Cristy                                  %
%                                 July 1992                                   %
%                                                                             %
%                                                                             %
%  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
%  dedicated to making software imaging solutions freely available.           %
%                                                                             %
%  You may not use this file except in compliance with the License.  You may  %
%  obtain a copy of the License at                                            %
%                                                                             %
%    http://www.imagemagick.org/script/license.php                            %
%                                                                             %
%  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.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     S i g m o i d a l C o n t r a s t I m a g e                             %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  SigmoidalContrastImage() adjusts the contrast of an image with a non-linear
%  sigmoidal contrast algorithm.  Increase the contrast of the image using a
%  sigmoidal transfer function without saturating highlights or shadows.
%  Contrast indicates how much to increase the contrast (0 is none; 3 is
%  typical; 20 is pushing it); mid-point indicates where midtones fall in the
%  resultant image (0 is white; 50% is middle-gray; 100% is black).  Set
%  sharpen to MagickTrue to increase the image contrast otherwise the contrast
%  is reduced.
%
%  The format of the SigmoidalContrastImage method is:
%
%      MagickBooleanType SigmoidalContrastImage(Image *image,
%        const MagickBooleanType sharpen,const char *levels,
%        ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o sharpen: Increase or decrease image contrast.
%
%    o contrast: strength of the contrast, the larger the number the more
%      'threshold-like' it becomes.
%
%    o midpoint: midpoint of the function as a color value 0 to QuantumRange.
%
%    o exception: return any errors or warnings in this structure.
%
*/

/*
  ImageMagick 6 has a version of this function which uses LUTs.
*/

/*
  Sigmoidal function Sigmoidal with inflexion point moved to b and "slope
  constant" set to a.

  The first version, based on the hyperbolic tangent tanh, when combined with
  the scaling step, is an exact arithmetic clone of the the sigmoid function
  based on the logistic curve. The equivalence is based on the identity

    1/(1+exp(-t)) = (1+tanh(t/2))/2

  (http://de.wikipedia.org/wiki/Sigmoidfunktion) and the fact that the
  scaled sigmoidal derivation is invariant under affine transformations of
  the ordinate.

  The tanh version is almost certainly more accurate and cheaper.  The 0.5
  factor in the argument is to clone the legacy ImageMagick behavior. The
  reason for making the define depend on atanh even though it only uses tanh
  has to do with the construction of the inverse of the scaled sigmoidal.
*/
#if defined(MAGICKCORE_HAVE_ATANH)
#define Sigmoidal(a,b,x) ( tanh((0.5*(a))*((x)-(b))) )
#else
#define Sigmoidal(a,b,x) ( 1.0/(1.0+exp((a)*((b)-(x)))) )
#endif
/*
  Scaled sigmoidal function:

    ( Sigmoidal(a,b,x) - Sigmoidal(a,b,0) ) /
    ( Sigmoidal(a,b,1) - Sigmoidal(a,b,0) )

  See http://osdir.com/ml/video.image-magick.devel/2005-04/msg00006.html and
  http://www.cs.dartmouth.edu/farid/downloads/tutorials/fip.pdf.  The limit
  of ScaledSigmoidal as a->0 is the identity, but a=0 gives a division by
  zero. This is fixed below by exiting immediately when contrast is small,
  leaving the image (or colormap) unmodified. This appears to be safe because
  the series expansion of the logistic sigmoidal function around x=b is

  1/2-a*(b-x)/4+...

  so that the key denominator s(1)-s(0) is about a/4 (a/2 with tanh).
*/
#define ScaledSigmoidal(a,b,x) (                    \
  (Sigmoidal((a),(b),(x))-Sigmoidal((a),(b),0.0)) / \
  (Sigmoidal((a),(b),1.0)-Sigmoidal((a),(b),0.0)) )
/*
  Inverse of ScaledSigmoidal, used for +sigmoidal-contrast.  Because b
  may be 0 or 1, the argument of the hyperbolic tangent (resp. logistic
  sigmoidal) may be outside of the interval (-1,1) (resp. (0,1)), even
  when creating a LUT from in gamut values, hence the branching.  In
  addition, HDRI may have out of gamut values.
  InverseScaledSigmoidal is not a two-sided inverse of ScaledSigmoidal:
  It is only a right inverse. This is unavoidable.
*/
static inline double InverseScaledSigmoidal(const double a,const double b,
  const double x)
{
  const double sig0=Sigmoidal(a,b,0.0);
  const double sig1=Sigmoidal(a,b,1.0);
  const double argument=(sig1-sig0)*x+sig0;
  const double clamped=
    (
#if defined(MAGICKCORE_HAVE_ATANH)
      argument < -1+MagickEpsilon
      ?
      -1+MagickEpsilon
      :
      ( argument > 1-MagickEpsilon ? 1-MagickEpsilon : argument )
    );
  return(b+(2.0/a)*atanh(clamped));
#else
      argument < MagickEpsilon
      ?
      MagickEpsilon
      :
      ( argument > 1-MagickEpsilon ? 1-MagickEpsilon : argument )
    );
  return(b-log(1.0/clamped-1.0)/a);
#endif
}

MagickExport MagickBooleanType SigmoidalContrastImage(Image *image,
  const MagickBooleanType sharpen,const double contrast,const double midpoint,
  ExceptionInfo *exception)
{
#define SigmoidalContrastImageTag  "SigmoidalContrast/Image"
#define ScaledSig(x) ( ClampToQuantum(QuantumRange* \
  ScaledSigmoidal(contrast,QuantumScale*midpoint,QuantumScale*(x))) )
#define InverseScaledSig(x) ( ClampToQuantum(QuantumRange* \
  InverseScaledSigmoidal(contrast,QuantumScale*midpoint,QuantumScale*(x))) )

  CacheView
    *image_view;

  MagickBooleanType
    status;

  MagickOffsetType
    progress;

  ssize_t
    y;

  /*
    Convenience macros.
  */
  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  /*
    Side effect: may clamp values unless contrast<MagickEpsilon, in which
    case nothing is done.
  */
  if (contrast < MagickEpsilon)
    return(MagickTrue);
  /*
    Sigmoidal-contrast enhance colormap.
  */
  if (image->storage_class == PseudoClass)
    {
      register ssize_t
        i;

      if (sharpen != MagickFalse)
        for (i=0; i < (ssize_t) image->colors; i++)
        {
          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
            image->colormap[i].red=ScaledSig(image->colormap[i].red);
          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
            image->colormap[i].green=ScaledSig(image->colormap[i].green);
          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
            image->colormap[i].blue=ScaledSig(image->colormap[i].blue);
          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
            image->colormap[i].alpha=ScaledSig(image->colormap[i].alpha);
        }
      else
        for (i=0; i < (ssize_t) image->colors; i++)
        {
          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
            image->colormap[i].red=InverseScaledSig(image->colormap[i].red);
          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
            image->colormap[i].green=InverseScaledSig(image->colormap[i].green);
          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
            image->colormap[i].blue=InverseScaledSig(image->colormap[i].blue);
          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
            image->colormap[i].alpha=InverseScaledSig(image->colormap[i].alpha);
        }
    }
  /*
    Sigmoidal-contrast enhance image.
  */
  status=MagickTrue;
  progress=0;
  image_view=AcquireAuthenticCacheView(image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static,4) shared(progress,status) \
    dynamic_number_threads(image,image->columns,image->rows,1)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    register Quantum
      *restrict q;

    register ssize_t
      x;

    if (status == MagickFalse)
      continue;
    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
    if (q == (Quantum *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      register ssize_t
        i;

      if (GetPixelMask(image,q) != 0)
        {
          q+=GetPixelChannels(image);
          continue;
        }
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        PixelChannel
          channel;

        PixelTrait
          traits;

        channel=GetPixelChannelChannel(image,i);
        traits=GetPixelChannelTraits(image,channel);
        if ((traits & UpdatePixelTrait) == 0)
          continue;
        if (sharpen != MagickFalse)
          q[i]=ScaledSig(q[i]);
        else
          q[i]=InverseScaledSig(q[i]);
      }
      q+=GetPixelChannels(image);
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
        #pragma omp critical (MagickCore_SigmoidalContrastImage)
#endif
        proceed=SetImageProgress(image,SigmoidalContrastImageTag,progress++,
          image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  image_view=DestroyCacheView(image_view);
  return(status);
}

[TheLion]

Originally Posted by NicolasRobidoux

I've added (tensor) Ginseng, Lanczos, Lanczos4, and EWA QuadraticJinc, produced similarly to the earlier ones, to http://web.cs.laurentian.ca/nrobidou.../madVR/glyphs/
-----
Mathias:
Ginseng is very often within JND of Lanczos, so much so that you probably could substitute it for Lanczos, and nobody would know better. (Sorry: I should have mentioned this earlier.) This being said, it does have measurably less halo, and it is slightly less jaggy, although the difference is only noticeable with sharp high contrast features. Although the difference is small, it is (on my brash Samsung monitor), noticeable on the "glyph" enlargements.
Tensor Lanczos is definitely too sharpening when enlarging. Ginseng backs the dial down from 11, just a nudge.

In my opinion the basic issue is that people use "sharp" resampling methods like Lanczos4 to actually sharpen the image. But this shouldn't be the prime function of a resampling algorithm. When Mathias implements a high quality "sharpening filter" I expect people will start to prefer resampling algorithms with the least artifacts (ringing, aliasing) and not so much the one with the most "perceived sharpness".

Availability of a special purpose sharpness filter will definitely shift resampling preferences.

[NicolasRobidoux]

Mathias:
Sorry. I was not clear.
What I invented was using sigmoidal-contrast to do better upsampling. Another way to put it is that I invented a parameterized family of color spaces which help get better results when enlarging (upsampling, actually: does not need to be orthogonal).
The code I show above was mostly written by me. It's based on earlier version written by, I believe, Anthony Thyssen and Cristy with some help from others who I can't all list without digging (Fred Weinhaus is one of them, and there is a textbook involved).
To give you the most optimal code, I need to know a few things about how you do things. As I mentioned, the ImageMagick code is way too complicated because it is mostly used to do other things than what I call "sigmoidizing".
-----
Makes more sense?

[Xaurus]

Originally Posted by nevcairiel

With madVR, the smartest thing to do is keep everything checked. Defaults, use them.

Haha thanks nev, I can't even remember what the defaults are. I guess there's no "reset to defaults" button within LAV Video.

So everything checked, including 10-bit and 16-bit formats as well as RGB 32 & 24?

[madshi]

Yes, that makes a lot of sense, thanks!

Originally Posted by NicolasRobidoux

What is your "anchor" colorspace?
Actually, what is your "anchor" luma/luminance channel.
What I most care about is whether it is a linear light channel, or whether it belongs to a perceptual color space.
-----
It's cheaper to do something like a "colorspace transformation" just to luma/luminance, yes? Does it cost a lot to convert to linear light, or is it one of these things that are basically free? Are colorspace transformations implemented as LUTs or with computational code?

I usually get Y'CbCr 4:2:0 from the decoder. I convert that to R'G'B', meaning gamma corrected RGB, usually with BT.601 or BT.709 primaries. I can convert that to linear light. Conversion to linear light is not free, but it's fast enough to not have to worry about it too much. I usually scale in R'G'B', usually with BT.709 primaries. Colorspace transformations? I don't convert to XYZ or things like that. I only convert from Y'CbCr to R'G'B' which is done by using simple matrix math.

Originally Posted by NicolasRobidoux

This is the key code in ImageMagick (actually, the IM6 code computes a LUT; what I'm showing is the IM7 code). This is way too complicated for all sorts of reasons. But it gives the idea, hopefully.

Thanks. So this seems to be the key formula?

#define Sigmoidal(a,b,x) ( tanh((0.5*(a))*((x)-(b))) )
#define Sigmoidal(a,b,x) ( 1.0/(1.0+exp((a)*((b)-(x)))) )

I suppose "x" is the image data, probably in the range 0.0 (black) to 1.0 (white), correct? But which parameters are used for "a" and "b" when you tell ImageMagick to do e.g. "+sigmoidal-contrast 7.5"? Basically all I need is the exact formula to use for your sigmoidization trick. Thanks a lot!!

[madshi]

P.S: Here's the parking meter image in full/original size:

http://madshi.net/parking-big.png

I usually upscale by 400% to check upscaling algorithms.

[NicolasRobidoux]

Mathias:

==================
Sigmoidization HOW TO
==================

Let me explain things thus, and we can sort out the nitty gritty later.

I've always done sigmoidization separately on each channel of a
physically linear color space which in ImageMagick is linear RGB
with sRGB primaries.

So, let's suppose that we have one channel of image data x which is
strongly constrained to the interval [0,1] and which can reasonably be
interpreted as one component of linear light. I understand that in
general there may be mild over- or undershoots, resulting from color space
conversion (it has to do with negative coefficients in transformation
matrices when things are properly normalized). If there may be over
and undershoots, let me know and I'll tell you how to fix things.

First, choose a contrast. Like in the ImageMagick code, let's call the
contrast "a". a should be a float > 0. As a approaches 0,
sigmoidization has less and less effect, so you may as well skip it.
I generally like values between 7 and 12.

If tanh is in your math library, it should compute faster than the
formula that uses exponentials, because tanh is a very easy function
to approximate. I can also provide fast approximations (given a bit of
time). Or I can give you the formulation in terms of the exp function.
Ideally, though, you have access to both tanh and atanh (the
hyperbolic tangent and its inverse).

Have the following constant set:

s = tanh(0.25*a)

(This is equal to Sigmoidal(a,1) = -Sigmoidal(a,0) because tanh is odd.)

Step 1: Convert your input to the "sigmoidal color space", that is,
"replace" every x by the value y computed as follows:

y = 0.5 + (2.0/a) * atanh(-s+(2.0*s)*x) <- friendly to multiply-add

Now, we are in the "sigmoidal color space". You don't need to clamp:
If your original x was in [0,1], y is in [0,1] too.

(Again: I know how to fix things if there may be over and undershoots
in the x data.)

Step 2:

Upsample as usual, using the y values. Let's say that we get a pixel
value of w.

Step 3: Convert w back to "linear color space":

z = tanh((0.5*a)*(w-0.5)) = tanh((-0.25*a)+(0.5*a)*w) <- friendly to multiply-add

final value = (z+s)/2s = 0.5+(0.5/s)*z

Only clamp if you want to. Step 3 does not care about what interval w is in.

-----

That's it!

[madshi]

Thank you, that looks easy enough!!

Two problems:

(1) Pixel shaders can do "atan", "tanh" and "exp", but not "atanh", unfortunately.
(2) Black is 0.0 and white is 1.0, correct? In that case: Yes, there may be values outside of [0,1]. Not only from YCbCr -> RGB conversion, but some movies also contain BTB (black than black) and WTW data. How can I solve that?

Thanks!!

[NicolasRobidoux]

Originally Posted by madshi

(2) Black is 0.0 and white is 1.0, correct? In that case: Yes, there may be values outside of [0,1]. Not only from YCbCr -> RGB conversion, but some movies also contain BTB (black than black) and WTW data. How can I solve that?

Yes, 0 is black and 1 is white.
Can you determine ahead of time an (ideally reasonably tight) interval which each channel cannot overshoot (after conversion to linear light or what resembles it)?
(I'm sure computing the min and max in each channel for the whole movie is not an option.)
-----
The atanh issue is not a big deal. Just need to sit down and do algebra. Atanh is often not included because it's so easy to reconstruct.

[NicolasRobidoux]

Mathias:
No need to do algebra: http://mathworld.wolfram.com/Inverse...icTangent.html.
atanh(x) = 0.5 * log((1+x)/(1-x))

[NicolasRobidoux]

Originally Posted by madshi

P.S: Here's the parking meter image in full/original size:

http://madshi.net/parking-big.png

I usually upscale by 400% to check upscaling algorithms.

It's actually the small version I'd like, so I can re-enlarge by 400% like you do.

[madshi]

Originally Posted by NicolasRobidoux

Can you determine ahead of time an (ideally reasonably tight) interval which each channel cannot overshoot (after conversion to linear light or what resembles it)?

If my math is correct, the possible value range for all 3 channels in *linear light* should be:

[-0.0029842504073062995363224069370219,1.2143440435408780674380747705496]

Originally Posted by NicolasRobidoux

Mathias:
No need to do algebra: http://mathworld.wolfram.com/Inverse...icTangent.html.
atanh(x) = 0.5 * log((1+x)/(1-x))

Cool. One problem solved.

[madshi]

Originally Posted by NicolasRobidoux

It's actually the small version I'd like, so I can re-enlarge by 400% like you do.

small parking meter

[Luv]

XP32
Catalyst 12.4,HD 4770
ZP 5.02

Hello,Madshi

everything is fine here but I just can't select the Intel MPEG2 decoder (While the libav can be without any pb).ZP uses the Cyberlink decoder (Wich ouputs 4:2:2...).
Is it a question of priority ?
Thank you.

[madshi]

Ok, I've implemented the Sigmoidal contrast stuff and it seems to work pretty well. I like the results better compared to the "pretend linear" hack. Here they are for comparison:

jinc3 normal -|- jinc3 hack -|- jinc3 sigmoidal

To my eyes the "sigmoidal" image clearly looks best. I've still left the anti-ringing filter enabled, though. Although the "sigmoidal" stuff does reduce the ringing artifacts, it doesn't completely remove them. So image quality is still noticeably better with the added anti-ringing filter on top.

So, only 2 questions remaining:

(1) Which is the "best" sigmoidal contrast? I've seen 7.5 (named with "safe") and 11.5. The image from above is with 11.5, which I prefer with this test image. Should I always use 11.5? Or would that be a bad idea?
(2) Still need to handle those out of [0,1] values somehow.

Thanks a lot, Nicolas, with your help upscaling quality just made another step forward in madVR...

[madshi]

@Luv, you need the Intel decoder DLL, see first page of this thread for a download link. Personally, I prefer libav, though.

[sunnah]

Originally Posted by digitech

Sunnah im interested about your HTPC, can you give me details about processor type, speeed, amount of ram etc, etc, pls? im about to build a new htpc just for the main use of playing movies with madvr, any info about your htpc details will be very apprecciated, thanks!

Well I bought a game laptop, mainly because of the high-end specs overall.
I watch my Korean tv shows always 1080p quality, on the laptop's full hd screen or by hdmi on a full hd led tv.
The laptop also has an overclock button on it, so that if it ever stutters(never happened), then I just push it then the processor speed will increase plus the graphic processing speed.

This is the laptop, http://nl.msi.com/product/nb/GT683.html

Cheers.

Don't forget to enable the "Hardware acceleration: LAV CUVID"

[Luv]

Originally Posted by madshi

@Luv, you need the Intel decoder DLL, see first page of this thread for a download link. Personally, I prefer libav, though.

I forgot (again) to add it.Thanks,Madshi.
I'm going to test the libav this evening.Could you say why it is to be preferred?

[madshi]

It's much faster, decodes 4:2:2 etc etc...

[NicolasRobidoux]

Mathias:

Here's the fix (or, more precisely, here is "a" fix, since I am not
enforcing symmetry w.r.t black (0) and white (1) and I'm not 100% sure
this is the best choice): The constant s of my HOW TO is strictly
between 0 and 1. (In the interval (0,1), which I think some Europeans
write ]0,1[.)

Let's assume that the interval in which every single channel belongs
to is contained in (A,B). Note that I don't want the bounds to be
attained. In the case of
[-0.0029842504073062995363224069370219,1.2143440435408780674380747705496],
I'd use something like A=-0.00298426 and B=1.214345. That is, I'd
enlarge the interval by some epsilon on each side to make sure that
the bounds cannot be reached. If I was in double precision, I'd use
something tighter, of course.

What you are going to do is that you are going to fix things so that
(A,B) is mapped to (-1,1) by the transformation

x -> -s+(2.0*s)*(b+m*x)

A little bit of algebra says that this is accomplished by choosing

m = (1/L)/s

and

b = 0.5 + (D/L)/s

where L = B-A and D = -0.5*(B+A)

-----

So, all you have to do is add a preliminary step that replaces x by
b+m*x, and add a final step that replaces "final value" by
(-b/m)+(final value)/m.

So, it adds two multiply-adds: one at the beginning, and one at the
end.