root / ImageMagick / trunk / magick / enhance.c

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%              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                %
%                                                                             %
%                                                                             %
%                    ImageMagick Image Enhancement Methods                    %
%                                                                             %
%                              Software Design                                %
%                                John Cristy                                  %
%                                 July 1992                                   %
%                                                                             %
%                                                                             %
%  Copyright 1999-2008 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.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
*/

/*
  Include declarations.
*/
#include "magick/studio.h"
#include "magick/cache-view.h"
#include "magick/color.h"
#include "magick/color-private.h"
#include "magick/colorspace.h"
#include "magick/enhance.h"
#include "magick/exception.h"
#include "magick/exception-private.h"
#include "magick/gem.h"
#include "magick/geometry.h"
#include "magick/image.h"
#include "magick/image-private.h"
#include "magick/memory_.h"
#include "magick/monitor.h"
#include "magick/monitor-private.h"
#include "magick/quantum.h"
#include "magick/quantum-private.h"
#include "magick/resample.h"
#include "magick/string_.h"

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     C l u t I m a g e                                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ClutImage() replaces colors in the image from a color lookup table.
%
%  The format of the ClutImage method is:
%
%      MagickBooleanType ClutImage(Image *image,Image *clut_image)
%      MagickBooleanType ClutImageChannel(Image *image,
%        const ChannelType channel,Image *clut_image)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o channel: the channel.
%
%    o clut_image: the color lookup.
%
*/

static ResampleFilter **DestroyResampleFilterThreadSet(ResampleFilter **filter)
{
  register long
    i;

  assert(filter != (ResampleFilter **) NULL);
  for (i=0; i < (long) GetCacheViewMaximumThreads(); i++)
    if (filter[i] != (ResampleFilter *) NULL)
      filter[i]=DestroyResampleFilter(filter[i]);
  return((ResampleFilter **) RelinquishMagickMemory(filter));
}

static ResampleFilter **AcquireResampleFilterThreadSet(const Image *image,
  ExceptionInfo *exception)
{
  register long
    i;

  ResampleFilter
    **filter;

  filter=(ResampleFilter **) AcquireQuantumMemory(GetCacheViewMaximumThreads(),
    sizeof(*filter));
  if (filter == (ResampleFilter **) NULL)
    return((ResampleFilter **) NULL);
  (void) ResetMagickMemory(filter,0,GetCacheViewMaximumThreads()*
    sizeof(*filter));
  for (i=0; i < (long) GetCacheViewMaximumThreads(); i++)
  {
    filter[i]=AcquireResampleFilter(image,exception);
    if (filter[i] == (ResampleFilter *) NULL)
      return(DestroyResampleFilterThreadSet(filter));
  }
  return(filter);
}

MagickExport MagickBooleanType ClutImage(Image *image,const Image *clut_image)
{
  return(ClutImageChannel(image,DefaultChannels,clut_image));
}

MagickExport MagickBooleanType ClutImageChannel(Image *image,
  const ChannelType channel,const Image *clut_image)
{
#define ClutImageTag  "Clut/Image"

  long
    adjust,
    progress,
    y;

  MagickBooleanType
    status;

  ResampleFilter
    **resample_filter;

  ViewInfo
    **image_view;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  assert(clut_image != (Image *) NULL);
  assert(clut_image->signature == MagickSignature);
  if (SetImageStorageClass(image,DirectClass) == MagickFalse)
    return(MagickFalse);
  /*
    Clut image.
  */
  status=MagickTrue;
  progress=0;
  adjust=clut_image->interpolate == IntegerInterpolatePixel ? 0 : 1;
  resample_filter=AcquireResampleFilterThreadSet(clut_image,&image->exception);
  image_view=AcquireCacheViewThreadSet(image);
  #pragma omp parallel for
  for (y=0; y < (long) image->rows; y++)
  {
    MagickPixelPacket
      pixel;

    register IndexPacket
      *indexes;

    register long
      id,
      x;

    register PixelPacket
      *q;

    if (status == MagickFalse)
      continue;
    id=GetCacheViewThreadId();
    q=GetCacheViewPixels(image_view[id],0,y,image->columns,1);
    if (q == (PixelPacket *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    indexes=GetCacheViewIndexes(image_view[id]);
    GetMagickPixelPacket(clut_image,&pixel);
    for (x=0; x < (long) image->columns; x++)
    {
      if ((channel & RedChannel) != 0)
        {
          pixel=ResamplePixelColor(resample_filter[id],QuantumScale*q->red*
            (clut_image->columns-adjust),QuantumScale*q->red*
            (clut_image->rows-adjust));
          q->red=RoundToQuantum(pixel.red);
        }
      if ((channel & GreenChannel) != 0)
        {
          pixel=ResamplePixelColor(resample_filter[id],QuantumScale*q->green*
            (clut_image->columns-adjust),QuantumScale*q->green*
            (clut_image->rows-adjust));
          q->green=RoundToQuantum(pixel.green);
        }
      if ((channel & BlueChannel) != 0)
        {
          pixel=ResamplePixelColor(resample_filter[id],QuantumScale*q->blue*
            (clut_image->columns-adjust),QuantumScale*q->blue*
            (clut_image->rows-adjust));
          q->blue=RoundToQuantum(pixel.blue);
        }
      if ((channel & OpacityChannel) != 0)
        {
          if (clut_image->matte == MagickFalse)
            {
              /*
                A gray-scale LUT replacement for an image alpha channel.
              */
              pixel=ResamplePixelColor(resample_filter[id],QuantumScale*
                (QuantumRange-q->opacity)*(clut_image->columns+adjust),
                QuantumScale*(QuantumRange-q->opacity)*(clut_image->rows+
                adjust));
              q->opacity=(Quantum) (QuantumRange-MagickPixelIntensityToQuantum(
                &pixel));
            }
          else
            if (image->matte == MagickFalse)
              {
                /*
                  A greyscale image being colored by a LUT with transparency.
                */
                pixel=ResamplePixelColor(resample_filter[id],QuantumScale*
                  PixelIntensity(q)*(clut_image->columns-adjust),QuantumScale*
                  PixelIntensity(q)*(clut_image->rows-adjust));
                q->opacity=RoundToQuantum(pixel.opacity);
              }
            else
              {
                /*
                  Direct alpha channel lookup.
                */
                pixel=ResamplePixelColor(resample_filter[id],QuantumScale*
                  q->opacity*(clut_image->columns-adjust),QuantumScale*
                  q->opacity* (clut_image->rows-adjust));
                q->opacity=RoundToQuantum(pixel.opacity);
              }
        }
      if (((channel & IndexChannel) != 0) &&
          (image->colorspace == CMYKColorspace))
        {
          pixel=ResamplePixelColor(resample_filter[id],QuantumScale*indexes[x]*
            (clut_image->columns-adjust),QuantumScale*indexes[x]*
            (clut_image->rows-adjust));
          indexes[x]=RoundToQuantum(pixel.index);
        }
      q++;
    }
    if (SyncCacheView(image_view[id]) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

        #pragma omp critical
        proceed=SetImageProgress(image,ClutImageTag,progress++,image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  image_view=DestroyCacheViewThreadSet(image_view);
  resample_filter=DestroyResampleFilterThreadSet(resample_filter);
  /*
    Enable alpha channel if CLUT image could enable it.
  */
  if ((clut_image->matte == MagickTrue) && ((channel & OpacityChannel) != 0))
    (void) SetImageAlphaChannel(image,ActivateAlphaChannel);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     C o n t r a s t I m a g e                                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ContrastImage() enhances the intensity differences between the lighter and
%  darker elements of the image.  Set sharpen to a MagickTrue to increase the
%  image contrast otherwise the contrast is reduced.
%
%  The format of the ContrastImage method is:
%
%      MagickBooleanType ContrastImage(Image *image,
%        const MagickBooleanType sharpen)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o sharpen: Increase or decrease image contrast.
%
*/

static void Contrast(const int sign,Quantum *red,Quantum *green,Quantum *blue)
{
  double
    brightness,
    hue,
    saturation;

  /*
    Enhance contrast: dark color become darker, light color become lighter.
  */
  assert(red != (Quantum *) NULL);
  assert(green != (Quantum *) NULL);
  assert(blue != (Quantum *) NULL);
  hue=0.0;
  saturation=0.0;
  brightness=0.0;
  ConvertRGBToHSB(*red,*green,*blue,&hue,&saturation,&brightness);
  brightness+=0.5*sign*(0.5*(sin(MagickPI*(brightness-0.5))+1.0)-brightness);
  if (brightness > 1.0)
    brightness=1.0;
  else
    if (brightness < 0.0)
      brightness=0.0;
  ConvertHSBToRGB(hue,saturation,brightness,red,green,blue);
}

MagickExport MagickBooleanType ContrastImage(Image *image,
  const MagickBooleanType sharpen)
{
#define ContrastImageTag  "Contrast/Image"

  int
    sign;

  long
    progress,
    y;

  MagickBooleanType
    status;

  register long
    i;

  ViewInfo
    **image_view;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  sign=sharpen != MagickFalse ? 1 : -1;
  if (image->storage_class == PseudoClass)
    {
      /*
        Contrast enhance colormap.
      */
      for (i=0; i < (long) image->colors; i++)
        Contrast(sign,&image->colormap[i].red,&image->colormap[i].green,
          &image->colormap[i].blue);
    }
  /*
    Contrast enhance image.
  */
  status=MagickTrue;
  progress=0;
  image_view=AcquireCacheViewThreadSet(image);
  #pragma omp parallel for
  for (y=0; y < (long) image->rows; y++)
  {
    register long
      id,
      x;

    register PixelPacket
      *q;

    if (status == MagickFalse)
      continue;
    id=GetCacheViewThreadId();
    q=GetCacheViewPixels(image_view[id],0,y,image->columns,1);
    if (q == (PixelPacket *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (long) image->columns; x++)
    {
      Contrast(sign,&q->red,&q->green,&q->blue);
      q++;
    }
    if (SyncCacheView(image_view[id]) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

        #pragma omp critical
        proceed=SetImageProgress(image,ContrastImageTag,progress++,image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  image_view=DestroyCacheViewThreadSet(image_view);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     C o n t r a s t S t r e t c h I m a g e                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  The ContrastStretchImage() is a simple image enhancement technique that
%  attempts to improve the contrast in an image by `stretching' the range of
%  intensity values it contains to span a desired range of values. It differs
%  from the more sophisticated histogram equalization in that it can only
%  apply %  a linear scaling function to the image pixel values.  As a result
%  the `enhancement' is less harsh.
%
%  The format of the ContrastStretchImage method is:
%
%      MagickBooleanType ContrastStretchImage(Image *image,
%        const char *levels)
%      MagickBooleanType ContrastStretchImageChannel(Image *image,
%        const unsigned long channel,const double black_point,
%        const double white_point)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o channel: the channel.
%
%    o black_point: the black point.
%
%    o white_point: the white point.
%
%    o levels: Specify the levels where the black and white points have the
%      range of 0 to number-of-pixels (e.g. 1%, 10x90%, etc.).
%
*/

MagickExport MagickBooleanType ContrastStretchImage(Image *image,
  const char *levels)
{
  double
    black_point,
    white_point;

  GeometryInfo
    geometry_info;

  MagickBooleanType
    status;

  MagickStatusType
    flags;

  /*
    Parse levels.
  */
  if (levels == (char *) NULL)
    return(MagickFalse);
  flags=ParseGeometry(levels,&geometry_info);
  black_point=geometry_info.rho;
  white_point=(double) image->columns*image->rows;
  if ((flags & SigmaValue) != 0)
    white_point=geometry_info.sigma;
  if ((flags & PercentValue) != 0)
    {
      black_point*=(double) QuantumRange/100.0;
      white_point*=(double) QuantumRange/100.0;
    }
  if ((flags & SigmaValue) == 0)
    white_point=(double) image->columns*image->rows-black_point;
  status=ContrastStretchImageChannel(image,DefaultChannels,black_point,
    white_point);
  return(status);
}

MagickExport MagickBooleanType ContrastStretchImageChannel(Image *image,
  const ChannelType channel,const double black_point,const double white_point)
{
#define MaxRange(color)  ((MagickRealType) ScaleQuantumToMap((Quantum) (color)))
#define ContrastStretchImageTag  "ContrastStretch/Image"

  double
    intensity;

  long
    progress,
    y;

  MagickBooleanType
    status;

  MagickPixelPacket
    black,
    *histogram,
    *stretch_map,
    white;

  register long
    i;

  ViewInfo
    **image_view;

  /*
    Allocate histogram and stretch map.
  */
  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  histogram=(MagickPixelPacket *) AcquireQuantumMemory(MaxMap+1UL,
    sizeof(*histogram));
  stretch_map=(MagickPixelPacket *) AcquireQuantumMemory(MaxMap+1UL,
    sizeof(*stretch_map));
  if ((histogram == (MagickPixelPacket *) NULL) ||
      (stretch_map == (MagickPixelPacket *) NULL))
    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
      image->filename);
  /*
    Form histogram.
  */
  status=MagickTrue;
  (void) ResetMagickMemory(histogram,0,(MaxMap+1)*sizeof(*histogram));
  image_view=AcquireCacheViewThreadSet(image);
  #pragma omp parallel for
  for (y=0; y < (long) image->rows; y++)
  {
    IndexPacket
      *indexes;

    register const PixelPacket
      *p;

    register long
      id,
      x;

    if (status == MagickFalse)
      continue;
    id=GetCacheViewThreadId();
    p=AcquireCacheViewPixels(image_view[id],0,y,image->columns,1,
      &image->exception);
    if (p == (const PixelPacket *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    indexes=GetCacheViewIndexes(image_view[id]);
    if (channel == DefaultChannels)
      #pragma omp parallel for
      for (x=0; x < (long) image->columns; x++)
      {
        Quantum
          intensity;

        intensity=PixelIntensityToQuantum(p+x);
        histogram[ScaleQuantumToMap(intensity)].red++;
        histogram[ScaleQuantumToMap(intensity)].green++;
        histogram[ScaleQuantumToMap(intensity)].blue++;
        histogram[ScaleQuantumToMap(intensity)].index++;
      }
    else
      #pragma omp parallel for
      for (x=0; x < (long) image->columns; x++)
      {
        if ((channel & RedChannel) != 0)
          histogram[ScaleQuantumToMap((p+x)->red)].red++;
        if ((channel & GreenChannel) != 0)
          histogram[ScaleQuantumToMap((p+x)->green)].green++;
        if ((channel & BlueChannel) != 0)
          histogram[ScaleQuantumToMap((p+x)->blue)].blue++;
        if ((channel & OpacityChannel) != 0)
          histogram[ScaleQuantumToMap((p+x)->opacity)].opacity++;
        if (((channel & IndexChannel) != 0) &&
            (image->colorspace == CMYKColorspace))
          histogram[ScaleQuantumToMap(indexes[x])].index++;
      }
  }
  /*
    Find the histogram boundaries by locating the black/white levels.
  */
  black.red=0.0;
  white.red=MaxRange(QuantumRange);
  if ((channel & RedChannel) != 0)
    {
      intensity=0.0;
      for (i=0; i <= (long) MaxMap; i++)
      {
        intensity+=histogram[i].red;
        if (intensity > black_point)
          break;
      }
      black.red=(MagickRealType) i;
      intensity=0.0;
      for (i=(long) MaxMap; i != 0; i--)
      {
        intensity+=histogram[i].red;
        if (intensity > ((double) image->columns*image->rows-white_point))
          break;
      }
      white.red=(MagickRealType) i;
    }
  black.green=0.0;
  white.green=MaxRange(QuantumRange);
  if ((channel & GreenChannel) != 0)
    {
      intensity=0.0;
      for (i=0; i <= (long) MaxMap; i++)
      {
        intensity+=histogram[i].green;
        if (intensity > black_point)
          break;
      }
      black.green=(MagickRealType) i;
      intensity=0.0;
      for (i=(long) MaxMap; i != 0; i--)
      {
        intensity+=histogram[i].green;
        if (intensity > ((double) image->columns*image->rows-white_point))
          break;
      }
      white.green=(MagickRealType) i;
    }
  black.blue=0.0;
  white.blue=MaxRange(QuantumRange);
  if ((channel & BlueChannel) != 0)
    {
      intensity=0.0;
      for (i=0; i <= (long) MaxMap; i++)
      {
        intensity+=histogram[i].blue;
        if (intensity > black_point)
          break;
      }
      black.blue=(MagickRealType) i;
      intensity=0.0;
      for (i=(long) MaxMap; i != 0; i--)
      {
        intensity+=histogram[i].blue;
        if (intensity > ((double) image->columns*image->rows-white_point))
          break;
      }
      white.blue=(MagickRealType) i;
    }
  black.opacity=0.0;
  white.opacity=MaxRange(QuantumRange);
  if ((channel & OpacityChannel) != 0)
    {
      intensity=0.0;
      for (i=0; i <= (long) MaxMap; i++)
      {
        intensity+=histogram[i].opacity;
        if (intensity > black_point)
          break;
      }
      black.opacity=(MagickRealType) i;
      intensity=0.0;
      for (i=(long) MaxMap; i != 0; i--)
      {
        intensity+=histogram[i].opacity;
        if (intensity > ((double) image->columns*image->rows-white_point))
          break;
      }
      white.opacity=(MagickRealType) i;
    }
  black.index=0.0;
  white.index=MaxRange(QuantumRange);
  if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
    {
      intensity=0.0;
      for (i=0; i <= (long) MaxMap; i++)
      {
        intensity+=histogram[i].index;
        if (intensity > black_point)
          break;
      }
      black.index=(MagickRealType) i;
      intensity=0.0;
      for (i=(long) MaxMap; i != 0; i--)
      {
        intensity+=histogram[i].index;
        if (intensity > ((double) image->columns*image->rows-white_point))
          break;
      }
      white.index=(MagickRealType) i;
    }
  histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
  /*
    Stretch the histogram to create the stretched image mapping.
  */
  (void) ResetMagickMemory(stretch_map,0,(MaxMap+1)*sizeof(*stretch_map));
  #pragma omp parallel for
  for (i=0; i <= (long) MaxMap; i++)
  {
    if ((channel & RedChannel) != 0)
      {
        if (i < (long) black.red)
          stretch_map[i].red=0.0;
        else
          if (i > (long) white.red)
            stretch_map[i].red=(MagickRealType) QuantumRange;
          else
            if (black.red != white.red)
              stretch_map[i].red=(MagickRealType) ScaleMapToQuantum(
                (MagickRealType) (MaxMap*(i-black.red)/(white.red-black.red)));
      }
    if ((channel & GreenChannel) != 0)
      {
        if (i < (long) black.green)
          stretch_map[i].green=0.0;
        else
          if (i > (long) white.green)
            stretch_map[i].green=(MagickRealType) QuantumRange;
          else
            if (black.green != white.green)
              stretch_map[i].green=(MagickRealType) ScaleMapToQuantum(
                (MagickRealType) (MaxMap*(i-black.green)/(white.green-
                black.green)));
      }
    if ((channel & BlueChannel) != 0)
      {
        if (i < (long) black.blue)
          stretch_map[i].blue=0.0;
        else
          if (i > (long) white.blue)
            stretch_map[i].blue=(MagickRealType) QuantumRange;
          else
            if (black.blue != white.blue)
              stretch_map[i].blue=(MagickRealType) ScaleMapToQuantum(
                (MagickRealType) (MaxMap*(i-black.blue)/(white.blue-
                black.blue)));
      }
    if ((channel & OpacityChannel) != 0)
      {
        if (i < (long) black.opacity)
          stretch_map[i].opacity=0.0;
        else
          if (i > (long) white.opacity)
            stretch_map[i].opacity=(MagickRealType) QuantumRange;
          else
            if (black.opacity != white.opacity)
              stretch_map[i].opacity=(MagickRealType) ScaleMapToQuantum(
                (MagickRealType) (MaxMap*(i-black.opacity)/(white.opacity-
                black.opacity)));
      }
    if (((channel & IndexChannel) != 0) &&
        (image->colorspace == CMYKColorspace))
      {
        if (i < (long) black.index)
          stretch_map[i].index=0.0;
        else
          if (i > (long) white.index)
            stretch_map[i].index=(MagickRealType) QuantumRange;
          else
            if (black.index != white.index)
              stretch_map[i].index=(MagickRealType) ScaleMapToQuantum(
                (MagickRealType) (MaxMap*(i-black.index)/(white.index-
                black.index)));
      }
  }
  /*
    Stretch the image.
  */
  if (((channel & OpacityChannel) != 0) || (((channel & IndexChannel) != 0) &&
      (image->colorspace == CMYKColorspace)))
    image->storage_class=DirectClass;
  if (image->storage_class == PseudoClass)
    {
      /*
        Stretch colormap.
      */
      #pragma omp parallel for
      for (i=0; i < (long) image->colors; i++)
      {
        if ((channel & RedChannel) != 0)
          {
            if (black.red != white.red)
              image->colormap[i].red=RoundToQuantum(stretch_map[
                ScaleQuantumToMap(image->colormap[i].red)].red);
          }
        if ((channel & GreenChannel) != 0)
          {
            if (black.green != white.green)
              image->colormap[i].green=RoundToQuantum(stretch_map[
                ScaleQuantumToMap(image->colormap[i].green)].green);
          }
        if ((channel & BlueChannel) != 0)
          {
            if (black.blue != white.blue)
              image->colormap[i].blue=RoundToQuantum(stretch_map[
                ScaleQuantumToMap(image->colormap[i].blue)].blue);
          }
        if ((channel & OpacityChannel) != 0)
          {
            if (black.opacity != white.opacity)
              image->colormap[i].opacity=RoundToQuantum(stretch_map[
                ScaleQuantumToMap(image->colormap[i].opacity)].opacity);
          }
      }
    }
  /*
    Stretch image.
  */
  status=MagickTrue;
  progress=0;
  #pragma omp parallel for
  for (y=0; y < (long) image->rows; y++)
  {
    IndexPacket
      *indexes;

    register long
      id,
      x;

    register PixelPacket
      *q;

    if (status == MagickFalse)
      continue;
    id=GetCacheViewThreadId();
    q=GetCacheViewPixels(image_view[id],0,y,image->columns,1);
    if (q == (PixelPacket *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    indexes=GetCacheViewIndexes(image_view[id]);
    for (x=0; x < (long) image->columns; x++)
    {
      if ((channel & RedChannel) != 0)
        {
          if (black.red != white.red)
            q->red=RoundToQuantum(stretch_map[ScaleQuantumToMap(
              q->red)].red);
        }
      if ((channel & GreenChannel) != 0)
        {
          if (black.green != white.green)
            q->green=RoundToQuantum(stretch_map[ScaleQuantumToMap(
              q->green)].green);
        }
      if ((channel & BlueChannel) != 0)
        {
          if (black.blue != white.blue)
            q->blue=RoundToQuantum(stretch_map[ScaleQuantumToMap(
              q->blue)].blue);
        }
      if ((channel & OpacityChannel) != 0)
        {
          if (black.opacity != white.opacity)
            q->opacity=RoundToQuantum(stretch_map[ScaleQuantumToMap(
              q->opacity)].opacity);
        }
      if (((channel & IndexChannel) != 0) &&
          (image->colorspace == CMYKColorspace))
        {
          if (black.index != white.index)
            indexes[x]=(IndexPacket) RoundToQuantum(stretch_map[
              ScaleQuantumToMap(indexes[x])].index);
        }
      q++;
    }
    if (SyncCacheView(image_view[id]) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

        #pragma omp critical
        proceed=SetImageProgress(image,ContrastStretchImageTag,progress++,
          image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  image_view=DestroyCacheViewThreadSet(image_view);
  stretch_map=(MagickPixelPacket *) RelinquishMagickMemory(stretch_map);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     E n h a n c e I m a g e                                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  EnhanceImage() applies a digital filter that improves the quality of a
%  noisy image.
%
%  The format of the EnhanceImage method is:
%
%      Image *EnhanceImage(const Image *image,ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o exception: Return any errors or warnings in this structure.
%
*/
MagickExport Image *EnhanceImage(const Image *image,ExceptionInfo *exception)
{
#define Enhance(weight) \
  mean=((MagickRealType) r->red+pixel.red)/2; \
  distance=(MagickRealType) r->red-(MagickRealType) pixel.red; \
  distance_squared=QuantumScale*(2.0*((MagickRealType) QuantumRange+1.0)+ \
     mean)*distance*distance; \
  mean=((MagickRealType) r->green+pixel.green)/2; \
  distance=(MagickRealType) r->green-(MagickRealType) pixel.green; \
  distance_squared+=4.0*distance*distance; \
  mean=((MagickRealType) r->blue+pixel.blue)/2; \
  distance=(MagickRealType) r->blue-(MagickRealType) pixel.blue; \
  distance_squared+=QuantumScale*(3.0*((MagickRealType) \
    QuantumRange+1.0)-1.0-mean)*distance*distance; \
  mean=((MagickRealType) r->opacity+pixel.opacity)/2; \
  distance=(MagickRealType) r->opacity-(MagickRealType) pixel.opacity; \
  distance_squared+=QuantumScale*(3.0*((MagickRealType) \
    QuantumRange+1.0)-1.0-mean)*distance*distance; \
  if (distance_squared < ((MagickRealType) QuantumRange*(MagickRealType) \
      QuantumRange/25.0f)) \
    { \
      aggregate.red+=(weight)*r->red; \
      aggregate.green+=(weight)*r->green; \
      aggregate.blue+=(weight)*r->blue; \
      aggregate.opacity+=(weight)*r->opacity; \
      total_weight+=(weight); \
    } \
  r++;
#define EnhanceImageTag  "Enhance/Image"

  Image
    *enhance_image;

  long
    progress,
    y;

  MagickBooleanType
    status;

  MagickPixelPacket
    zero;

  ViewInfo
    **enhance_view,
    **image_view;

  /*
    Initialize enhanced image attributes.
  */
  assert(image != (const Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  assert(exception != (ExceptionInfo *) NULL);
  assert(exception->signature == MagickSignature);
  if ((image->columns < 5) || (image->rows < 5))
    return((Image *) NULL);
  enhance_image=CloneImage(image,0,0,MagickTrue,exception);
  if (enhance_image == (Image *) NULL)
    return((Image *) NULL);
  if (SetImageStorageClass(enhance_image,DirectClass) == MagickFalse)
    {
      InheritException(exception,&enhance_image->exception);
      enhance_image=DestroyImage(enhance_image);
      return((Image *) NULL);
    }
  /*
    Enhance image.
  */
  status=MagickTrue;
  progress=0;
  (void) ResetMagickMemory(&zero,0,sizeof(zero));
  image_view=AcquireCacheViewThreadSet(image);
  enhance_view=AcquireCacheViewThreadSet(enhance_image);
  #pragma omp parallel for
  for (y=0; y < (long) image->rows; y++)
  {
    register const PixelPacket
      *p;

    register long
      id,
      x;

    register PixelPacket
      *q;

    /*
      Read another scan line.
    */
    if (status == MagickFalse)
      continue;
    id=GetCacheViewThreadId();
    p=AcquireCacheViewPixels(image_view[id],-2,y-2,image->columns+4,5,
      exception);
    q=GetCacheViewPixels(enhance_view[id],0,y,enhance_image->columns,1);
    if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (long) image->columns; x++)
    {
      MagickPixelPacket
        aggregate;

      MagickRealType
        distance,
        distance_squared,
        mean,
        total_weight;

      PixelPacket
        pixel;

      register const PixelPacket
        *r;

      /*
        Compute weighted average of target pixel color components.
      */
      aggregate=zero;
      total_weight=0.0;
      r=p+2*(image->columns+4)+2;
      pixel=(*r);
      r=p;
      Enhance(5.0); Enhance(8.0); Enhance(10.0); Enhance(8.0); Enhance(5.0);
      r=p+(image->columns+4);
      Enhance(8.0); Enhance(20.0); Enhance(40.0); Enhance(20.0); Enhance(8.0);