[SOLVED] making sense of sigmoidal-contrast

Discuss digital image processing techniques and algorithms. We encourage its application to ImageMagick but you can discuss any software solutions here.
User avatar
anthony
Posts: 8883
Joined: 2004-05-31T19:27:03-07:00
Authentication code: 8675308
Location: Brisbane, Australia

Re: making sense of sigmoidal-contrast

Post by anthony »

I always felt the LUT unnecessary when the math is not too complex. It is also less exact than 16bit (1024 entries or 10 bit resolution) and fails to handle HDRI values that range beyond the LUT bounds.


The 1/2 was either the 'mean' for handling the singularity (EG generate flat gray on zero argument), or it is a rounding constant (add 0.5 before converting to interger) in the LUT or range changes (normalized floating point value to quantum range integer).
Anthony Thyssen -- Webmaster for ImageMagick Example Pages
https://imagemagick.org/Usage/
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: making sense of sigmoidal-contrast

Post by NicolasRobidoux »

anthony wrote:I always felt the LUT unnecessary when the math is not too complex.
Agreed.
anthony wrote:It is also less exact than 16bit (1024 entries or 10 bit resolution) and fails to handle HDRI values that range beyond the LUT bounds.
Agreed.
anthony wrote:The 1/2 was either the 'mean' for handling the singularity (EG generate flat gray on zero argument), or it is a rounding constant (add 0.5 before converting to interger) in the LUT or range changes (normalized floating point value to quantum range integer).
That's what it was in the version I removed. In the version I cleaned up and "restored", it looked like rounding the result (not setting the input value to the middle of the relevant interval, which actually makes more sense to me) if rounding was done by clamping, but the result was stored in MagickRealType, and it did not occur to me to check whether MagickRealType could be an integer type.
-----
In any case, should I remove the LUT and organize things so that direct computation is as cheap as it can? I actually would prefer that. (I won't use my fast series approximations because we push values near singularities and I want almost perfect invertibility. Computing 1/(1+exp(a*(m-s)) is not expensive anyway, esp. if I rearrange things so as to minimize divisions (and hoist constants out, but the compiler should be doing that without assistance).
-----
P.S. MagickRealType can't be an integer type, so it could not have been "cheap C rounding of nonnegative values" (viewtopic.php?f=2&t=21196#p86192). Anyway, the question of what the "+0.5" was doing there is moot if I get rid of the LUT.
Last edited by NicolasRobidoux on 2012-08-06T04:22:18-07:00, edited 1 time in total.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: making sense of sigmoidal-contrast

Post by NicolasRobidoux »

Note: The way the scaled sigmoidal works, you can overrun the bounds by arbitrary amounts when using the contrast enhancing version (-sigmoidal-contrast), but when using the contrast-reducing version (+sigmoidal-contrast) you can only overrun them by quantities that depend on a (the contrast) and m (the midpoint), although when a is zero there is no bound (because the mapping is the identity). It should not take me long to figure out the exact range into which the input of +sigmoidal-contrast needs to be. (Calculus 1.) And then I could clamp things internally so as to extend the domain of the + version when we don't use LUTs. I'd document in the code comments. (Why publish articles when you can insert explanatory comments in published code?)
Also, I think that there is "post-mortem" clamping going on. I'd be enclined to remove it as well.
-----
Cristy: In another thread, I mentioned that "automatic -clamp" may be a good thing to have for the sigmoidal-contrast operation. It turns out I was wrong. One just needs to clamp "enough" internally when using the + version ("enough" depending on the value of contrast: more contrast requires more clamping). Otherwise, no need.
-----
(Blows my mind that something so apparently simple as sigmoidal-contrast took me so much time to fully figure out.)
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: making sense of sigmoidal-contrast

Post by NicolasRobidoux »

(Right now I'm only modifying IM7 so as to make comparison easy and give some time to testers to catch whatever bug I may be introducing. Will port to IM6 later.)
When I have a minute, I will change enhance.c so that sigmoidal-contrast not use a LUT in HDRI mode, and allow the maximum overshoots allowed (anything for -sigmoidal-contrast, bound to a range determined by the values of a (contrast) and possibly b (midpoint) when +sigmoidal-contrast is used, unless a is very close to zero in which case both - and +sigmoidal-contrast are the identity). I'll leave the LUT in for non-HDRI mode.
This will lead to better behavior when there are intermediate result overshoots worth preserving, for example when there is negative transparency that arises from the use of a resampling filter with negative lobes.
henrywho
Posts: 188
Joined: 2011-08-17T06:46:40-07:00
Authentication code: 8675308

Re: making sense of sigmoidal-contrast

Post by henrywho »

NicolasRobidoux wrote:When I have a minute, I will change enhance.c so that sigmoidal-contrast not use a LUT in HDRI mode
Hope we will have an official HDRI compile. :D
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: making sense of sigmoidal-contrast

Post by NicolasRobidoux »

henrywho wrote:...
Hope we will have an official HDRI compile. :D
viewtopic.php?f=2&t=21534&start=15#p88646
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: making sense of sigmoidal-contrast

Post by NicolasRobidoux »

I just ported the sigmoidal-contrast improvements from IM7 to IM6.
From svn rev 9287 both branches of ImageMagick have matching sigmoidal-contrast.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: [SOLVED] making sense of sigmoidal-contrast

Post by NicolasRobidoux »

I've branched sigmoidal-contrast again: As of IM7 svn rev 9292, it does not use a LUT; in IM6, it still does.
I want the more accurate results because I have it in my plans to find a good value of the sigmoidization contrast using a quantitative test suite. So, I don't want results to be corrupted by a 10 bit LUT when using exponentials.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: [SOLVED] making sense of sigmoidal-contrast

Post by NicolasRobidoux »

Well, well:
In HDRI mode, bleeding-edge IM6 with the LUT takes about .7 seconds to do the sigmoidized 800% enlargement of the dragon I show elsewhere, and bleeding-edge IM7 without the LUT takes about .8 seconds.
Of course this may have to do with other things than LUT VS direct in sigmoidal-contrast.
But at this point I sure can't claim that no LUT is faster (on a recent 2-core laptop).
P.S. With very large images, the LUT-free one is much slower.
Last edited by NicolasRobidoux on 2012-09-11T07:16:46-07:00, edited 1 time in total.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: [SOLVED] making sense of sigmoidal-contrast

Post by NicolasRobidoux »

I benchmarked: As I expected, the tanh/atanh clone of the logistic version of sigmoidal-contrast is faster.
User avatar
anthony
Posts: 8883
Joined: 2004-05-31T19:27:03-07:00
Authentication code: 8675308
Location: Brisbane, Australia

Re: making sense of sigmoidal-contrast

Post by anthony »

NicolasRobidoux wrote:Now I see that fudging things by using ( s(0) + s(1) ) / 2 computed from the actual m was a pretty nice and expedient solution.
Thank you. That was what I did to solve the alpha = 0 case so it produces uniform gray, rather than uniform black, when the function becomes a constant.

Of course, later on I discovered that when alpha is small btu not zero, the autoscaling, just expands the very slight gradient into a near linear (no-op) function in any case, so I suppose that limit case was not a great solution in any case.

As I tried to explained to Fred at the time, having an Alpha less than 1.0 is just not 'useful' in the sigmoidal function.
Anthony Thyssen -- Webmaster for ImageMagick Example Pages
https://imagemagick.org/Usage/
User avatar
fmw42
Posts: 25562
Joined: 2007-07-02T17:14:51-07:00
Authentication code: 1152
Location: Sunnyvale, California, USA

Re: making sense of sigmoidal-contrast

Post by fmw42 »

anthony wrote:
As I tried to explained to Fred at the time, having an Alpha less than 1.0 is just not 'useful' in the sigmoidal function.
But (in concept) it should approach a straight line as contrast goes to zero. The IM function does not. Thus as the contrast gets small enough (say .001 or so -- even 0.1 would be helpful) one can just substitute a straight line in the IM function. I believe that Nicholas proved/mentioned that earlier in one of the discussions.

Also I had at first created my sigmoidal script using tanh and atanh. But then I realized that that produced the same results as the IM sigmoidal function, but could not prove mathematically that they did. So I just used the IM function rather than create a lut from the tanh/atanh and apply that.

Just a few days ago I found a mathematical proof at the reference http://de.wikipedia.org/wiki/Sigmoidfunktion from a link on a reference to http://en.wikipedia.org/wiki/File:Gjl-t%28x%29.svg that Nicolas found.

If Nicolas has completed his development with the tanh and atanh and it replicates the IM sigmodal-contrast and if it does not need any special processing when the contrast gets too small and since it is faster, then perhaps we should switch -sigmoidal-contrast to use his version.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: [SOLVED] making sense of sigmoidal-contrast

Post by NicolasRobidoux »

I have already put in both IM6 and IM7 a sigmoidal that gracefully becomes the identity when contrast=0. That is: pixel values are unchanged when contrast=0, which makes complete sense.
I did not have to force it: This is what happens naturally to the scaled sigmoidal contrast when a -> 0. ("a" is the contrast value inside the C code. I did not change that.)
All I had to do is handle the removable singularity (inconsequential division by 0). See the source code comments:

Code: Select all

MagickExport MagickBooleanType SigmoidalContrastImage(Image *image,
  const MagickBooleanType sharpen,const double contrast,const double midpoint,
  ExceptionInfo *exception)
{
#define SigmoidalContrastImageTag  "SigmoidalContrast/Image"

  CacheView
    *image_view;

  MagickBooleanType
    status;

  MagickOffsetType
    progress;

  ssize_t
    y;

  /*
    Side effect: clamps values unless contrast<MagickEpsilon, in which
    case nothing is done.
  */
  if (contrast<MagickEpsilon)
    return(MagickTrue);

  /*
    Sigmoidal function 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 its 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 formula:
    ( 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 above 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, hence the branching.
    In addition, HDRI may have out of gamut values.
    InverseScaledSigmoidal is not a two-side inverse of
    ScaledSigmoidal: It is only a right inverse. This is unavoidable.
  */
#if defined(MAGICKCORE_HAVE_ATANH)
#define InverseScaledSigmoidal(a,b,x) ({                             \
  const double _argument =                                           \
    (Sigmoidal((a),(b),1.0)-Sigmoidal((a),(b),0.0)) * (x) +          \
    Sigmoidal((a),(b),0.0);                                          \
  const double _clamped_argument =                                   \
    ( _argument < -1+MagickEpsilon ? -1+MagickEpsilon :              \
    ( _argument > 1-MagickEpsilon ? 1-MagickEpsilon : _argument ) ); \
  (b) + (2.0/(a)) * atanh(_clamped_argument); })
#else
#define InverseScaledSigmoidal(a,b,x) ({                             \
  const double _argument =                                           \
    (Sigmoidal((a),(b),1.0)-Sigmoidal((a),(b),0.0)) * (x) +          \
    Sigmoidal((a),(b),0.0);                                          \
  const double _clamped_argument =                                   \
    ( _argument < MagickEpsilon ? MagickEpsilon :                    \
    ( _argument > 1-MagickEpsilon ? 1-MagickEpsilon : _argument ) ); \
  (b) + (-1.0/(a)) * log(1.0/_clamped_argument+-1.0); })
#endif
  /*
    Convenience macros.
  */
#define ScaledSig(x) ( ClampToQuantum(QuantumRange* \
  ScaledSigmoidal(contrast,QuantumScale*midpoint,QuantumScale*(x))) )
#define InverseScaledSig(x) ( ClampToQuantum(QuantumRange* \
  InverseScaledSigmoidal(contrast,QuantumScale*midpoint,QuantumScale*(x))) )

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  /*
    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);
}
This is the IM7 version, which does not use a LUT. The IM6 version still uses a LUT. IM7 is slower than IM6 when the images are large, and vice versa. And IM7 is more precise (which should be: default is HDRI).
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: [SOLVED] making sense of sigmoidal-contrast

Post by NicolasRobidoux »

There is a significant difference between the IM7 (no LUT) and the IM6 (LUT) sigmoidal-contrast as I implemented them: In IM7, out of nominal gamut pixel values will only be clamped if they risk overflowing or underflowing InverseScaledSigmoidal, and this only happens with +sigmoidal-contrast, not -sigmoidal-contrast. Basically, IM7 leaves alone pixel values that are out of the nominal gamut, as much as possible. Of course, this only really matters in HDRI.
Summary: In HDRI, the IM7 version of sigmoidal contrast will process negative and large pixel values without clamping if possible. Not so in IM6, because of the use of a LUT: values are always clamped by sigmoidal-contrast (even in HDRI), unless contrast is extremely close to 0 (<MagickEpsilon, actually), in which case sigmoidal-contrast is a no-op pass through.
Last edited by NicolasRobidoux on 2012-09-12T08:57:24-07:00, edited 1 time in total.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: [SOLVED] making sense of sigmoidal-contrast

Post by NicolasRobidoux »

Cristy just informed me that the above code does not compile under Windows. I used gcc compatible macro tricks :(
Cristy moved things to static inline functions, which of course is better.
Post Reply