fx.c

Go to the documentation of this file.

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%                                 FFFFF  X   X                                %
%                                 F       X X                                 %
%                                 FFF      X                                  %
%                                 F       X X                                 %
%                                 F      X   X                                %
%                                                                             %
%                                                                             %
%                   MagickCore Image Special Effects Methods                  %
%                                                                             %
%                               Software Design                               %
%                                 John Cristy                                 %
%                                 October 1996                                %
%                                                                             %
%                                                                             %
%  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/annotate.h"
#include "magick/cache.h"
#include "magick/cache-view.h"
#include "magick/color.h"
#include "magick/color-private.h"
#include "magick/composite.h"
#include "magick/decorate.h"
#include "magick/draw.h"
#include "magick/effect.h"
#include "magick/enhance.h"
#include "magick/exception.h"
#include "magick/exception-private.h"
#include "magick/fx.h"
#include "magick/fx-private.h"
#include "magick/gem.h"
#include "magick/geometry.h"
#include "magick/layer.h"
#include "magick/list.h"
#include "magick/log.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/option.h"
#include "magick/pixel-private.h"
#include "magick/property.h"
#include "magick/quantum.h"
#include "magick/random_.h"
#include "magick/resample.h"
#include "magick/resample-private.h"
#include "magick/resize.h"
#include "magick/shear.h"
#include "magick/splay-tree.h"
#include "magick/statistic.h"
#include "magick/string_.h"
#include "magick/transform.h"
#include "magick/utility.h"

/*
  Define declarations.
*/
#define LeftShiftOperator 0xf5
#define RightShiftOperator 0xf6
#define LessThanEqualOperator 0xf7
#define GreaterThanEqualOperator 0xf8
#define EqualOperator 0xf9
#define NotEqualOperator 0xfa
#define LogicalAndOperator 0xfb
#define LogicalOrOperator 0xfc

struct _FxInfo
{
  const Image
    *images;

  MagickBooleanType
    matte;

  char
    *expression;

  SplayTreeInfo
    *colors,
    *symbols;

  ResampleFilter
    **resample_filter;

  ExceptionInfo
    *exception;
};

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   A c q u i r e F x I n f o                                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AcquireFxInfo() allocates the FxInfo structure.
%
%  The format of the AcquireFxInfo method is:
%
%      FxInfo *AcquireFxInfo(Image *image,const char *expression)
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o expression: the expression.
%
*/
MagickExport FxInfo *AcquireFxInfo(const Image *image,const char *expression)
{
  char
    fx_op[2];

  FxInfo
    *fx_info;

  register long
    i;

  fx_info=(FxInfo *) AcquireMagickMemory(sizeof(*fx_info));
  if (fx_info == (FxInfo *) NULL)
    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
  (void) ResetMagickMemory(fx_info,0,sizeof(*fx_info));
  fx_info->exception=AcquireExceptionInfo();
  fx_info->images=image;
  fx_info->matte=image->matte;
  fx_info->colors=NewSplayTree(CompareSplayTreeString,RelinquishMagickMemory,
    RelinquishMagickMemory);
  fx_info->symbols=NewSplayTree(CompareSplayTreeString,RelinquishMagickMemory,
    RelinquishMagickMemory);
  fx_info->resample_filter=(ResampleFilter **) AcquireQuantumMemory(
    GetImageListLength(fx_info->images),sizeof(*fx_info->resample_filter));
  if (fx_info->resample_filter == (ResampleFilter **) NULL)
    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
  for (i=0; i < (long) GetImageListLength(fx_info->images); i++)
    fx_info->resample_filter[i]=AcquireResampleFilter(GetImageFromList(
      fx_info->images,i),fx_info->exception);
  fx_info->expression=ConstantString(expression);
  (void) SubstituteString(&fx_info->expression," ","");  /* compact string */
  if ((strstr(fx_info->expression,"e+") != (char *) NULL) ||
      (strstr(fx_info->expression,"e-") != (char *) NULL))
    {
      /*
        Convert scientific notation.
      */
      (void) SubstituteString(&fx_info->expression,"0e+","0*10^");
      (void) SubstituteString(&fx_info->expression,"1e+","1*10^");
      (void) SubstituteString(&fx_info->expression,"2e+","2*10^");
      (void) SubstituteString(&fx_info->expression,"3e+","3*10^");
      (void) SubstituteString(&fx_info->expression,"4e+","4*10^");
      (void) SubstituteString(&fx_info->expression,"5e+","5*10^");
      (void) SubstituteString(&fx_info->expression,"6e+","6*10^");
      (void) SubstituteString(&fx_info->expression,"7e+","7*10^");
      (void) SubstituteString(&fx_info->expression,"8e+","8*10^");
      (void) SubstituteString(&fx_info->expression,"9e+","9*10^");
      (void) SubstituteString(&fx_info->expression,"0e-","0*10^-");
      (void) SubstituteString(&fx_info->expression,"1e-","1*10^-");
      (void) SubstituteString(&fx_info->expression,"2e-","2*10^-");
      (void) SubstituteString(&fx_info->expression,"3e-","3*10^-");
      (void) SubstituteString(&fx_info->expression,"4e-","4*10^-");
      (void) SubstituteString(&fx_info->expression,"5e-","5*10^-");
      (void) SubstituteString(&fx_info->expression,"6e-","6*10^-");
      (void) SubstituteString(&fx_info->expression,"7e-","7*10^-");
      (void) SubstituteString(&fx_info->expression,"8e-","8*10^-");
      (void) SubstituteString(&fx_info->expression,"9e-","9*10^-");
    }
  /*
    Convert complex to simple operators.
  */
  fx_op[1]='\0';
  *fx_op=(char) LeftShiftOperator;
  (void) SubstituteString(&fx_info->expression,"<<",fx_op);
  *fx_op=(char) RightShiftOperator;
  (void) SubstituteString(&fx_info->expression,">>",fx_op);
  *fx_op=(char) LessThanEqualOperator;
  (void) SubstituteString(&fx_info->expression,"<=",fx_op);
  *fx_op=(char) GreaterThanEqualOperator;
  (void) SubstituteString(&fx_info->expression,">=",fx_op);
  *fx_op=(char) EqualOperator;
  (void) SubstituteString(&fx_info->expression,"==",fx_op);
  *fx_op=(char) NotEqualOperator;
  (void) SubstituteString(&fx_info->expression,"!=",fx_op);
  *fx_op=(char) LogicalAndOperator;
  (void) SubstituteString(&fx_info->expression,"&&",fx_op);
  *fx_op=(char) LogicalOrOperator;
  (void) SubstituteString(&fx_info->expression,"||",fx_op);
  return(fx_info);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     A d d N o i s e I m a g e                                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AddNoiseImage() adds random noise to the image.
%
%  The format of the AddNoiseImage method is:
%
%      Image *AddNoiseImage(const Image *image,const NoiseType noise_type,
%        ExceptionInfo *exception)
%      Image *AddNoiseImageChannel(const Image *image,const ChannelType channel,
%        const NoiseType noise_type,ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o channel: the channel type.
%
%    o noise_type:  The type of noise: Uniform, Gaussian, Multiplicative,
%      Impulse, Laplacian, or Poisson.
%
%    o exception: return any errors or warnings in this structure.
%
*/

static Quantum GenerateNoise(const Quantum pixel,const NoiseType noise_type,
  const MagickRealType attenuate)
{
#define NoiseEpsilon  (attenuate*1.0e-5)
#define SigmaUniform  ScaleCharToQuantum((unsigned char) (attenuate*4.0+0.5))
#define SigmaGaussian  ScaleCharToQuantum((unsigned char) (attenuate*4.0+0.5))
#define SigmaImpulse  (attenuate*0.10)
#define SigmaLaplacian ScaleCharToQuantum((unsigned char) (attenuate*10.0+0.5))
#define SigmaMultiplicativeGaussian \
  ScaleCharToQuantum((unsigned char) (attenuate*1.0+0.5))
#define SigmaPoisson  (attenuate*0.05)
#define TauGaussian  ScaleCharToQuantum((unsigned char) (attenuate*20.0+0.5))

  MagickRealType
    alpha,
    beta,
    noise,
    sigma;

  alpha=GetPseudoRandomValue();
  if (alpha == 0.0)
    alpha=1.0;
  switch (noise_type)
  {
    case UniformNoise:
    default:
    {
      noise=(MagickRealType) pixel+SigmaUniform*(alpha-0.5);
      break;
    }
    case GaussianNoise:
    {
      MagickRealType
        tau;

      beta=GetPseudoRandomValue();
      sigma=sqrt(-2.0*log((double) alpha))*cos((double) (2.0*MagickPI*beta));
      tau=sqrt(-2.0*log((double) alpha))*sin((double) (2.0*MagickPI*beta));
      noise=(MagickRealType) pixel+sqrt((double) pixel)*SigmaGaussian*sigma+
        TauGaussian*tau;
      break;
    }
    case MultiplicativeGaussianNoise:
    {
      if (alpha <= NoiseEpsilon)
        sigma=(MagickRealType) QuantumRange;
      else
        sigma=sqrt(-2.0*log((double) alpha));
      beta=GetPseudoRandomValue();
      noise=(MagickRealType) pixel+pixel*SigmaMultiplicativeGaussian*sigma/2.0*
        cos((double) (2.0*MagickPI*beta));
      break;
    }
    case ImpulseNoise:
    {
      if (alpha < (SigmaImpulse/2.0))
        noise=0.0;
       else
         if (alpha >= (1.0-(SigmaImpulse/2.0)))
           noise=(MagickRealType) QuantumRange;
         else
           noise=(MagickRealType) pixel;
      break;
    }
    case LaplacianNoise:
    {
      if (alpha <= 0.5)
        {
          if (alpha <= NoiseEpsilon)
            noise=(MagickRealType) pixel-(MagickRealType) QuantumRange;
          else
            noise=(MagickRealType) pixel+ScaleCharToQuantum((unsigned char)
              (SigmaLaplacian*log((double) (2.0*alpha))+0.5));
          break;
        }
      beta=1.0-alpha;
      if (beta <= (0.5*NoiseEpsilon))
        noise=(MagickRealType) (pixel+QuantumRange);
      else
        noise=(MagickRealType) pixel-ScaleCharToQuantum((unsigned char)
          (SigmaLaplacian*log((double) (2.0*beta))+0.5));
      break;
    }
    case PoissonNoise:
    {
      MagickRealType
        poisson;

      register long
        i;

      poisson=exp(-SigmaPoisson*(double) ScaleQuantumToChar(pixel));
      for (i=0; alpha > poisson; i++)
      {
        beta=GetPseudoRandomValue();
        alpha=alpha*beta;
      }
      noise=(MagickRealType) ScaleCharToQuantum((unsigned char)
        (i/SigmaPoisson));
      break;
    }
    case RandomNoise:
    {
      noise=(MagickRealType) QuantumRange*GetPseudoRandomValue();
      break;
    }
  }
  return(RoundToQuantum(noise));
}

MagickExport Image *AddNoiseImage(const Image *image,const NoiseType noise_type,
  ExceptionInfo *exception)
{
  Image
    *noise_image;

  noise_image=AddNoiseImageChannel(image,DefaultChannels,noise_type,exception);
  return(noise_image);
}

MagickExport Image *AddNoiseImageChannel(const Image *image,
  const ChannelType channel,const NoiseType noise_type,ExceptionInfo *exception)
{
#define AddNoiseImageTag  "AddNoise/Image"

  const char
    *option;

  Image
    *noise_image;

  long
    progress,
    y;

  MagickBooleanType
    status;

  MagickRealType
    attenuate;

  ViewInfo
    *image_view,
    *noise_view;

  /*
    Initialize noise 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);
  noise_image=CloneImage(image,0,0,MagickTrue,exception);
  if (noise_image == (Image *) NULL)
    return((Image *) NULL);
  if (SetImageStorageClass(noise_image,DirectClass) == MagickFalse)
    {
      InheritException(exception,&noise_image->exception);
      noise_image=DestroyImage(noise_image);
      return((Image *) NULL);
    }
  /*
    Add noise in each row.
  */
  attenuate=1.0;
  option=GetImageProperty(image,"attenuate");
  if (option != (char *) NULL)
    attenuate=atof(option);
  status=MagickTrue;
  progress=0;
  image_view=AcquireCacheView(image);
  noise_view=AcquireCacheView(noise_image);
  for (y=0; y < (long) image->rows; y++)
  {
    register const IndexPacket
      *indexes;

    register const PixelPacket
      *p;

    register IndexPacket
      *noise_indexes;

    register long
      x;

    register PixelPacket
      *q;

    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
    q=GetCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
      exception);
    if ((p == (PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
      {
        status=MagickFalse;
        break;
      }
    indexes=GetCacheViewVirtualIndexes(image_view);
    noise_indexes=GetCacheViewAuthenticIndexes(noise_view);
    for (x=0; x < (long) image->columns; x++)
    {
      if ((channel & RedChannel) != 0)
        q->red=GenerateNoise(p->red,noise_type,attenuate);
      if ((channel & GreenChannel) != 0)
        q->green=GenerateNoise(p->green,noise_type,attenuate);
      if ((channel & BlueChannel) != 0)
        q->blue=GenerateNoise(p->blue,noise_type,attenuate);
      if ((channel & OpacityChannel) != 0)
        q->opacity=GenerateNoise(p->opacity,noise_type,attenuate);
      if (((channel & IndexChannel) != 0) &&
          (image->colorspace == CMYKColorspace))
        noise_indexes[x]=(IndexPacket) GenerateNoise(indexes[x],noise_type,
          attenuate);
      p++;
      q++;
    }
    if (SyncCacheViewAuthenticPixels(noise_view,exception) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

        proceed=SetImageProgress(image,AddNoiseImageTag,progress++,
          image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
    if (status == MagickFalse)
      break;
  }
  noise_view=DestroyCacheView(noise_view);
  image_view=DestroyCacheView(image_view);
  if (status == MagickFalse)
    noise_image=DestroyImage(noise_image);
  return(noise_image);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     C h a r c o a l I m a g e                                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  CharcoalImage() creates a new image that is a copy of an existing one with
%  the edge highlighted.  It allocates the memory necessary for the new Image
%  structure and returns a pointer to the new image.
%
%  The format of the CharcoalImage method is:
%
%      Image *CharcoalImage(const Image *image,const double radius,
%        const double sigma,ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o radius: the radius of the pixel neighborhood.
%
%    o sigma: the standard deviation of the Gaussian, in pixels.
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport Image *CharcoalImage(const Image *image,const double radius,
  const double sigma,ExceptionInfo *exception)
{
  Image
    *charcoal_image,
    *clone_image,
    *edge_image;

  assert(image != (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);
  clone_image=CloneImage(image,0,0,MagickTrue,exception);
  if (clone_image == (Image *) NULL)
    return((Image *) NULL);
  (void) SetImageType(clone_image,GrayscaleType);
  edge_image=EdgeImage(clone_image,radius,exception);
  clone_image=DestroyImage(clone_image);
  if (edge_image == (Image *) NULL)
    return((Image *) NULL);
  charcoal_image=BlurImage(edge_image,radius,sigma,exception);
  edge_image=DestroyImage(edge_image);
  if (charcoal_image == (Image *) NULL)
    return((Image *) NULL);
  (void) NormalizeImage(charcoal_image);
  (void) NegateImage(charcoal_image,MagickFalse);
  (void) SetImageType(charcoal_image,GrayscaleType);
  return(charcoal_image);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     C o l o r i z e I m a g e                                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ColorizeImage() blends the fill color with each pixel in the image.
%  A percentage blend is specified with opacity.  Control the application
%  of different color components by specifying a different percentage for
%  each component (e.g. 90/100/10 is 90% red, 100% green, and 10% blue).
%
%  The format of the ColorizeImage method is:
%
%      Image *ColorizeImage(const Image *image,const char *opacity,
%        const PixelPacket colorize,ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o opacity:  A character string indicating the level of opacity as a
%      percentage.
%
%    o colorize: A color value.
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport Image *ColorizeImage(const Image *image,const char *opacity,
  const PixelPacket colorize,ExceptionInfo *exception)
{
#define ColorizeImageTag  "Colorize/Image"

  GeometryInfo
    geometry_info;

  Image
    *colorize_image;

  long
    progress,
    y;

  MagickBooleanType
    status;

  MagickPixelPacket
    pixel;

  MagickStatusType
    flags;

  ViewInfo
    **colorize_view,
    **image_view;

  /*
    Allocate colorized image.
  */
  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);
  colorize_image=CloneImage(image,image->columns,image->rows,MagickTrue,
    exception);
  if (colorize_image == (Image *) NULL)
    return((Image *) NULL);
  if (SetImageStorageClass(colorize_image,DirectClass) == MagickFalse)
    {
      InheritException(exception,&colorize_image->exception);
      colorize_image=DestroyImage(colorize_image);
      return((Image *) NULL);
    }
  if (opacity == (const char *) NULL)
    return(colorize_image);
  /*
    Determine RGB values of the pen color.
  */
  flags=ParseGeometry(opacity,&geometry_info);
  pixel.red=geometry_info.rho;
  pixel.green=geometry_info.rho;
  pixel.blue=geometry_info.rho;
  pixel.opacity=(MagickRealType) OpaqueOpacity;
  if ((flags & SigmaValue) != 0)
    pixel.green=geometry_info.sigma;
  if ((flags & XiValue) != 0)
    pixel.blue=geometry_info.xi;
  if ((flags & PsiValue) != 0)
    pixel.opacity=geometry_info.psi;
  /*
    Colorize DirectClass image.
  */
  status=MagickTrue;
  progress=0;
  image_view=AcquireCacheViewThreadSet(image);
  colorize_view=AcquireCacheViewThreadSet(colorize_image);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(dynamic,8) shared(progress,status)
#endif
  for (y=0; y < (long) image->rows; y++)
  {
    MagickBooleanType
      sync;

    register const PixelPacket
      *p;

    register long
      id,
      x;

    register PixelPacket
      *q;

    if (status == MagickFalse)
      continue;
    id=GetCacheViewThreadId();
    p=GetCacheViewVirtualPixels(image_view[id],0,y,image->columns,1,exception);
    q=QueueCacheViewAuthenticPixels(colorize_view[id],0,y,
      colorize_image->columns,1,exception);
    if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (long) image->columns; x++)
    {
      q->red=(Quantum) ((p->red*(100.0-pixel.red)+
        colorize.red*pixel.red)/100.0);
      q->green=(Quantum) ((p->green*(100.0-pixel.green)+
        colorize.green*pixel.green)/100.0);
      q->blue=(Quantum) ((p->blue*(100.0-pixel.blue)+
        colorize.blue*pixel.blue)/100.0);
      q->opacity=(Quantum) ((p->opacity*(100.0-pixel.opacity)+
        colorize.opacity*pixel.opacity)/100.0);
      p++;
      q++;
    }
    sync=SyncCacheViewAuthenticPixels(colorize_view[id],exception);
    if (sync == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp critical
#endif
        proceed=SetImageProgress(image,ColorizeImageTag,progress++,image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  image_view=DestroyCacheViewThreadSet(image_view);
  colorize_view=DestroyCacheViewThreadSet(colorize_view);
  if (status == MagickFalse)
    colorize_image=DestroyImage(colorize_image);
  return(colorize_image);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     C o n v o l v e I m a g e                                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ConvolveImage() applies a custom convolution kernel to the image.
%
%  The format of the ConvolveImage method is:
%
%      Image *ConvolveImage(const Image *image,const unsigned long order,
%        const double *kernel,ExceptionInfo *exception)
%      Image *ConvolveImageChannel(const Image *image,const ChannelType channel,
%        const unsigned long order,const double *kernel,
%        ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o channel: the channel type.
%
%    o order: the number of columns and rows in the filter kernel.
%
%    o kernel: An array of double representing the convolution kernel.
%
%    o exception: return any errors or warnings in this structure.
%
*/

MagickExport Image *ConvolveImage(const Image *image,const unsigned long order,
  const double *kernel,ExceptionInfo *exception)
{
  Image
    *convolve_image;

  convolve_image=ConvolveImageChannel(image,DefaultChannels,order,kernel,
    exception);
  return(convolve_image);
}

MagickExport Image *ConvolveImageChannel(const Image *image,
  const ChannelType channel,const unsigned long order,const double *kernel,
  ExceptionInfo *exception)
{
#define ConvolveImageTag  "Convolve/Image"

  double
    *normal_kernel;

  Image
    *convolve_image;

  long
    progress,
    y;

  MagickBooleanType
    status;

  MagickPixelPacket
    zero;

  MagickRealType
    bias,
    gamma;

  register long
    i;

  unsigned long
    width;

  /*
    Initialize convolve image attributes.
  */
  assert(image != (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);
  width=order;
  if ((width % 2) == 0)
    ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
  convolve_image=CloneImage(image,0,0,MagickTrue,exception);
  if (convolve_image == (Image *) NULL)
    return((Image *) NULL);
  if (SetImageStorageClass(convolve_image,DirectClass) == MagickFalse)
    {
      InheritException(exception,&convolve_image->exception);
      convolve_image=DestroyImage(convolve_image);
      return((Image *) NULL);
    }
  if (image->debug != MagickFalse)
    {
      char
        format[MaxTextExtent],
        *message;

      long
        u,
        v;

      register const double
        *k;

      (void) LogMagickEvent(TransformEvent,GetMagickModule(),
        "  ConvolveImage with %ldx%ld kernel:",width,width);
      message=AcquireString("");
      k=kernel;
      for (v=0; v < (long) width; v++)
      {
        *message='\0';
        (void) FormatMagickString(format,MaxTextExtent,"%ld: ",v);
        (void) ConcatenateString(&message,format);
        for (u=0; u < (long) width; u++)
        {
          (void) FormatMagickString(format,MaxTextExtent,"%+f ",*k++);
          (void) ConcatenateString(&message,format);
        }
        (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
      }
      message=DestroyString(message);
    }
  /*
    Normalize kernel.
  */
  normal_kernel=(double *) AcquireQuantumMemory(width*width,
    sizeof(*normal_kernel));
  if (normal_kernel == (double *) NULL)
    {
      convolve_image=DestroyImage(convolve_image);
      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
    }
  gamma=0.0;
  for (i=0; i < (long) (width*width); i++)
    gamma+=kernel[i];
  gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
  for (i=0; i < (long) (width*width); i++)
    normal_kernel[i]=gamma*kernel[i];
  /*
    Convolve image.
  */
  status=MagickTrue;
  progress=0;
  GetMagickPixelPacket(image,&zero);
  bias=image->bias;
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(dynamic,8) shared(progress,status)
#endif
  for (y=0; y < (long) image->rows; y++)
  {
    MagickBooleanType
      sync;

    register const IndexPacket
      *indexes;

    register const PixelPacket
      *p;

    register IndexPacket
      *convolve_indexes;

    register long
      x;

    register PixelPacket
      *q;

    if (status == MagickFalse)
      continue;
    p=GetVirtualPixels(image,-((long) width/2L),y-(long) (width/2L),
      image->columns+width,width,exception);
    q=GetAuthenticPixels(convolve_image,0,y,convolve_image->columns,1,
      exception);
    if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    indexes=GetVirtualIndexQueue(image);
    convolve_indexes=GetAuthenticIndexQueue(convolve_image);
    for (x=0; x < (long) image->columns; x++)
    {
      long
        j,
        v;

      MagickPixelPacket
        pixel;

      register const double
        *k;

      register long
        u;

      pixel=zero;
      k=normal_kernel;
      j=0;
      if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
        {
          for (v=0; v < (long) width; v++)
          {
            for (u=0; u < (long) width; u++)
            {
              pixel.red+=(*k)*(p+u+j)->red;
              pixel.green+=(*k)*(p+u+j)->green;
              pixel.blue+=(*k)*(p+u+j)->blue;
              pixel.opacity+=(*k)*(p+u+j)->opacity;
              k++;
            }
            j+=image->columns+width;
          }
          if ((channel & RedChannel) != 0)
            q->red=RoundToQuantum(pixel.red+bias);
          if ((channel & GreenChannel) != 0)
            q->green=RoundToQuantum(pixel.green+bias);
          if ((channel & BlueChannel) != 0)
            q->blue=RoundToQuantum(pixel.blue+bias);
          if ((channel & OpacityChannel) != 0)
            q->opacity=RoundToQuantum(pixel.opacity+bias);
          if (((channel & IndexChannel) != 0) &&
              (image->colorspace == CMYKColorspace))
            {
              k=normal_kernel;
              j=0;
              for (v=0; v < (long) width; v++)
              {
                for (u=0; u < (long) width; u++)
                {
                  pixel.index+=(*k)*indexes[x+u+j];
                  k++;
                }
                j+=image->columns+width;
              }
              convolve_indexes[x]=RoundToQuantum(pixel.index+bias);
            }
        }
      else
        {
          MagickRealType
            alpha,
            gamma;

          gamma=0.0;
          for (v=0; v < (long) width; v++)
          {
            for (u=0; u < (long) width; u++)
            {
              alpha=(MagickRealType) (QuantumScale*(QuantumRange-
                (p+u+j)->opacity));
              pixel.red+=(*k)*alpha*(p+u+j)->red;
              pixel.green+=(*k)*alpha*(p+u+j)->green;
              pixel.blue+=(*k)*alpha*(p+u+j)->blue;
              pixel.opacity+=(*k)*(p+u+j)->opacity;
              gamma+=alpha;
              k++;
            }
            j+=image->columns+width;
          }
          gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
          if ((channel & RedChannel) != 0)
            q->red=RoundToQuantum(gamma*pixel.red+bias);
          if ((channel & GreenChannel) != 0)
            q->green=RoundToQuantum(gamma*pixel.green+bias);
          if ((channel & BlueChannel) != 0)
            q->blue=RoundToQuantum(gamma*pixel.blue+bias);
          if ((channel & OpacityChannel) != 0)
            q->opacity=RoundToQuantum(pixel.opacity+bias);
          if (((channel & IndexChannel) != 0) &&
              (image->colorspace == CMYKColorspace))
            {
              k=normal_kernel;
              j=0;
              for (v=0; v < (long) width; v++)
              {