Legacy ImageMagick Discussions Archive

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