root / ImageMagick / branches / ImageMagick-6.3.5 / magick / effect.c

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%                   EEEEE  FFFFF  FFFFF  EEEEE  CCCC  TTTTT                   %
%                   E      F      F      E     C        T                     %
%                   EEE    FFF    FFF    EEE   C        T                     %
%                   E      F      F      E     C        T                     %
%                   EEEEE  F      F      EEEEE  CCCC    T                     %
%                                                                             %
%                                                                             %
%                      ImageMagick Image Effects Methods                      %
%                                                                             %
%                               Software Design                               %
%                                 John Cristy                                 %
%                                 October 1996                                %
%                                                                             %
%                                                                             %
%  Copyright 1999-2007 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/property.h"
#include "magick/blob.h"
#include "magick/cache-view.h"
#include "magick/color.h"
#include "magick/color-private.h"
#include "magick/colorspace.h"
#include "magick/constitute.h"
#include "magick/decorate.h"
#include "magick/draw.h"
#include "magick/enhance.h"
#include "magick/exception.h"
#include "magick/exception-private.h"
#include "magick/effect.h"
#include "magick/fx.h"
#include "magick/gem.h"
#include "magick/geometry.h"
#include "magick/image-private.h"
#include "magick/list.h"
#include "magick/log.h"
#include "magick/memory_.h"
#include "magick/monitor.h"
#include "magick/montage.h"
#include "magick/pixel-private.h"
#include "magick/property.h"
#include "magick/quantize.h"
#include "magick/quantum.h"
#include "magick/random_.h"
#include "magick/resize.h"
#include "magick/resource_.h"
#include "magick/segment.h"
#include "magick/shear.h"
#include "magick/signature.h"
#include "magick/string_.h"
#include "magick/transform.h"
#include "magick/threshold.h"

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     A d a p t i v e B l u r I m a g e                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AdaptiveBlurImage() adaptively blurs the image by blurring less
%  intensely near image edges and more intensely far from edges.  We blur the
%  image with a Gaussian operator of the given radius and standard deviation
%  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
%  radius of 0 and AdaptiveBlurImage() selects a suitable radius for you.
%
%  The format of the AdaptiveBlurImage method is:
%
%      Image *AdaptiveBlurImage(const Image *image,const double radius,
%        const double sigma,ExceptionInfo *exception)
%      Image *AdaptiveBlurImageChannel(const Image *image,
%        const ChannelType channel,double radius,const double sigma,
%        ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: The image.
%
%    o channel: The channel type.
%
%    o radius: The radius of the Gaussian, in pixels, not counting the center
%      pixel.
%
%    o sigma: The standard deviation of the Laplacian, in pixels.
%
%    o exception: Return any errors or warnings in this structure.
%
*/

MagickExport Image *AdaptiveBlurImage(const Image *image,const double radius,
  const double sigma,ExceptionInfo *exception)
{
  Image
    *blur_image;

  blur_image=AdaptiveBlurImageChannel(image,DefaultChannels,radius,sigma,
    exception);
  return(blur_image);
}

MagickExport Image *AdaptiveBlurImageChannel(const Image *image,
  const ChannelType channel,const double radius,const double sigma,
  ExceptionInfo *exception)
{
#define AdaptiveBlurImageTag  "Convolve/Image"

  double
    **kernel;

  Image
    *blur_image,
    *edge_image,
    *gaussian_image;

  IndexPacket
    *indexes,
    *blur_indexes;

  long
    j,
    u,
    v,
    y;

  MagickBooleanType
    status;

  MagickPixelPacket
    pixel;

  MagickRealType
    alpha,
    gamma,
    normalize;

  register const double
    *k;

  register const PixelPacket
    *p,
    *r;

  register long
    i,
    x;

  register PixelPacket
    *q;

  unsigned long
    width;

  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);
  blur_image=CloneImage(image,0,0,MagickTrue,exception);
  if (blur_image == (Image *) NULL)
    return((Image *) NULL);
  if (fabs(sigma) <= MagickEpsilon)
    return(blur_image);
  if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
    {
      InheritException(exception,&blur_image->exception);
      blur_image=DestroyImage(blur_image);
      return((Image *) NULL);
    }
  /*
    Edge detect the image brighness channel, level, blur, and level again.
  */
  edge_image=EdgeImage(image,radius,exception);
  if (edge_image == (Image *) NULL)
    {
      blur_image=DestroyImage(blur_image);
      return((Image *) NULL);
    }
  (void) LevelImage(edge_image,"20%,95%");
  gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
  if (gaussian_image != (Image *) NULL)
    {
      edge_image=DestroyImage(edge_image);
      edge_image=gaussian_image;
    }
  (void) LevelImage(edge_image,"10%,95%");
  /*
    Create a set of kernels from maximum (radius,sigma) to minimum.
  */
  width=GetOptimalKernelWidth2D(radius,sigma);
  kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
  if (kernel == (double **) NULL)
    {
      edge_image=DestroyImage(edge_image);
      blur_image=DestroyImage(blur_image);
      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
    }
  (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
  for (i=0; i < (long) width; i+=2)
  {
    kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
      sizeof(**kernel));
    if (kernel[i] == (double *) NULL)
      break;
    j=0;
    normalize=0.0;
    for (v=(-((long) (width-i)/2)); v <= (long) ((width-i)/2); v++)
    {
      for (u=(-((long) (width-i)/2)); u <= (long) ((width-i)/2); u++)
      {
        alpha=exp(-((double) u*u+v*v)/(2.0*sigma*sigma));
        kernel[i][j]=(double) (alpha/(2.0*MagickPI*sigma*sigma));
        if (((width-i) < 3) || (u != 0) || (v != 0))
          normalize+=kernel[i][j];
        j++;
      }
    }
    kernel[i][j/2]=(double) ((-2.0)*normalize);
    normalize=0.0;
    for (j=0; j < (long) ((width-i)*(width-i)); j++)
      normalize+=kernel[i][j];
    if (fabs(normalize) <= MagickEpsilon)
      normalize=1.0;
    normalize=1.0/normalize;
    for (j=0; j < (long) ((width-i)*(width-i)); j++)
      kernel[i][j]=(double) (normalize*kernel[i][j]);
  }
  if (i < (long) width)
    {
      for (i-=2; i >= 0; i-=2)
        kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
      kernel=(double **) RelinquishMagickMemory(kernel);
      edge_image=DestroyImage(edge_image);
      blur_image=DestroyImage(blur_image);
      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
    }
  /*
    Adaptively blur image.
  */
  for (y=0; y < (long) blur_image->rows; y++)
  {
    r=AcquireImagePixels(edge_image,0,y,edge_image->columns,1,exception);
    q=GetImagePixels(blur_image,0,y,blur_image->columns,1);
    if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
      break;
    indexes=GetIndexes(image);
    blur_indexes=GetIndexes(blur_image);
    for (x=0; x < (long) blur_image->columns; x++)
    {
      GetMagickPixelPacket(image,&pixel);
      gamma=0.0;
      i=(long) (width*QuantumScale*PixelIntensity(r)+0.5);
      if ((i & 0x01) != 0)
        i--;
      p=AcquireImagePixels(image,x-((long) (width-i)/2L),y-(long)
        ((width-i)/2L),width-i,width-i,exception);
      if (p == (const PixelPacket *) NULL)
        break;
      k=kernel[i];
      for (v=0; v < (long) (width-i); v++)
      {
        for (u=0; u < (long) (width-i); u++)
        {
          alpha=1.0;
          if (((channel & OpacityChannel) != 0) &&
              (image->matte != MagickFalse))
            alpha=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
          if ((channel & RedChannel) != 0)
            pixel.red+=(*k)*alpha*p->red;
          if ((channel & GreenChannel) != 0)
            pixel.green+=(*k)*alpha*p->green;
          if ((channel & BlueChannel) != 0)
            pixel.blue+=(*k)*alpha*p->blue;
          if ((channel & OpacityChannel) != 0)
            pixel.opacity+=(*k)*p->opacity;
          if (((channel & IndexChannel) != 0) &&
              (image->colorspace == CMYKColorspace))
            pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
          gamma+=(*k)*alpha;
          k++;
          p++;
        }
      }
      gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
      if ((channel & RedChannel) != 0)
        q->red=RoundToQuantum(gamma*pixel.red+image->bias);
      if ((channel & GreenChannel) != 0)
        q->green=RoundToQuantum(gamma*pixel.green+image->bias);
      if ((channel & BlueChannel) != 0)
        q->blue=RoundToQuantum(gamma*pixel.blue+image->bias);
      if ((channel & OpacityChannel) != 0)
        q->opacity=RoundToQuantum(pixel.opacity+image->bias);
      if (((channel & IndexChannel) != 0) &&
          (image->colorspace == CMYKColorspace))
        blur_indexes[x]=RoundToQuantum(gamma*pixel.index+image->bias);
      q++;
      r++;
    }
    if (SyncImagePixels(blur_image) == MagickFalse)
      break;
    if ((image->progress_monitor != (MagickProgressMonitor) NULL) &&
        (QuantumTick(y,image->rows) != MagickFalse))
      {
        status=image->progress_monitor(AdaptiveBlurImageTag,y,image->rows,
          image->client_data);
        if (status == MagickFalse)
          break;
      }
  }
  edge_image=DestroyImage(edge_image);
  for (i=0; i < (long) width;  i+=2)
    kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
  kernel=(double **) RelinquishMagickMemory(kernel);
  return(blur_image);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     A d a p t i v e S h a r p e n I m a g e                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AdaptiveSharpenImage() adaptively sharpens the image by sharpening more
%  intensely near image edges and less intensely far from edges. We sharpen the
%  image with a Gaussian operator of the given radius and standard deviation
%  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
%  radius of 0 and AdaptiveSharpenImage() selects a suitable radius for you.
%
%  The format of the AdaptiveSharpenImage method is:
%
%      Image *AdaptiveSharpenImage(const Image *image,const double radius,
%        const double sigma,ExceptionInfo *exception)
%      Image *AdaptiveSharpenImageChannel(const Image *image,
%        const ChannelType channel,double radius,const double sigma,
%        ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: The image.
%
%    o channel: The channel type.
%
%    o radius: The radius of the Gaussian, in pixels, not counting the center
%      pixel.
%
%    o sigma: The standard deviation of the Laplacian, in pixels.
%
%    o exception: Return any errors or warnings in this structure.
%
*/

MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
  const double sigma,ExceptionInfo *exception)
{
  Image
    *sharp_image;

  sharp_image=AdaptiveSharpenImageChannel(image,DefaultChannels,radius,sigma,
    exception);
  return(sharp_image);
}

MagickExport Image *AdaptiveSharpenImageChannel(const Image *image,
  const ChannelType channel,const double radius,const double sigma,
  ExceptionInfo *exception)
{
#define AdaptiveSharpenImageTag  "Convolve/Image"

  double
    **kernel;

  Image
    *blur_image,
    *edge_image,
    *sharp_image;

  IndexPacket
    *indexes,
    *sharp_indexes;

  long
    j,
    u,
    v,
    y;

  MagickBooleanType
    status;

  MagickPixelPacket
    pixel;

  MagickRealType
    alpha,
    gamma,
    normalize;

  register const double
    *k;

  register const PixelPacket
    *p,
    *r;

  register long
    i,
    x;

  register PixelPacket
    *q;

  unsigned long
    width;

  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);
  sharp_image=CloneImage(image,0,0,MagickTrue,exception);
  if (sharp_image == (Image *) NULL)
    return((Image *) NULL);
  if (SetImageStorageClass(sharp_image,DirectClass) == MagickFalse)
    {
      InheritException(exception,&sharp_image->exception);
      sharp_image=DestroyImage(sharp_image);
      return((Image *) NULL);
    }
  /*
    Edge detect the image brighness channel, level, blur, and level again.
  */
  edge_image=EdgeImage(image,radius,exception);
  if (edge_image == (Image *) NULL)
    {
      sharp_image=DestroyImage(sharp_image);
      return((Image *) NULL);
    }
  (void) LevelImage(edge_image,"20%,95%");
  blur_image=GaussianBlurImage(edge_image,radius,sigma,exception);
  if (blur_image != (Image *) NULL)
    {
      edge_image=DestroyImage(edge_image);
      edge_image=blur_image;
    }
  (void) LevelImage(edge_image,"10%,95%");
  /*
    Create a set of kernels from maximum (radius,sigma) to minimum.
  */
  width=GetOptimalKernelWidth2D(radius,sigma);
  kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
  if (kernel == (double **) NULL)
    {
      edge_image=DestroyImage(edge_image);
      sharp_image=DestroyImage(sharp_image);
      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
    }
  (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
  for (i=0; i < (long) width; i+=2)
  {
    kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
      sizeof(**kernel));
    if (kernel[i] == (double *) NULL)
      break;
    j=0;
    normalize=0.0;
    for (v=(-((long) (width-i)/2)); v <= (long) ((width-i)/2); v++)
    {
      for (u=(-((long) (width-i)/2)); u <= (long) ((width-i)/2); u++)
      {
        alpha=exp(-((double) u*u+v*v)/(2.0*sigma*sigma));
        kernel[i][j]=(double) (-alpha/(2.0*MagickPI*sigma*sigma));
        if (((width-i) < 3) || (u != 0) || (v != 0))
          normalize+=kernel[i][j];
        j++;
      }
    }
    kernel[i][j/2]=(double) ((-2.0)*normalize);
    normalize=0.0;
    for (j=0; j < (long) ((width-i)*(width-i)); j++)
      normalize+=kernel[i][j];
    if (fabs(normalize) <= MagickEpsilon)
      normalize=1.0;
    normalize=1.0/normalize;
    for (j=0; j < (long) ((width-i)*(width-i)); j++)
      kernel[i][j]=(double) (normalize*kernel[i][j]);
  }
  if (i < (long) width)
    {
      for (i-=2; i >= 0; i-=2)
        kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
      kernel=(double **) RelinquishMagickMemory(kernel);
      edge_image=DestroyImage(edge_image);
      sharp_image=DestroyImage(sharp_image);
      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
    }
  /*
    Adaptively sharpen image.
  */
  for (y=0; y < (long) sharp_image->rows; y++)
  {
    r=AcquireImagePixels(edge_image,0,y,edge_image->columns,1,exception);
    q=GetImagePixels(sharp_image,0,y,sharp_image->columns,1);
    if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
      break;
    indexes=GetIndexes(image);
    sharp_indexes=GetIndexes(sharp_image);
    for (x=0; x < (long) sharp_image->columns; x++)
    {
      GetMagickPixelPacket(image,&pixel);
      gamma=0.0;
      i=(long) (width*QuantumScale*(QuantumRange-PixelIntensity(r))+0.5);
      if ((i & 0x01) != 0)
        i--;
      p=AcquireImagePixels(image,x-((long) (width-i)/2L),y-(long)
        ((width-i)/2L),width-i,width-i,exception);
      if (p == (const PixelPacket *) NULL)
        break;
      k=kernel[i];
      for (v=0; v < (long) (width-i); v++)
      {
        for (u=0; u < (long) (width-i); u++)
        {
          alpha=1.0;
          if (((channel & OpacityChannel) != 0) &&
              (image->matte != MagickFalse))
            alpha=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
          if ((channel & RedChannel) != 0)
            pixel.red+=(*k)*alpha*p->red;
          if ((channel & GreenChannel) != 0)
            pixel.green+=(*k)*alpha*p->green;
          if ((channel & BlueChannel) != 0)
            pixel.blue+=(*k)*alpha*p->blue;
          if ((channel & OpacityChannel) != 0)
            pixel.opacity+=(*k)*p->opacity;
          if (((channel & IndexChannel) != 0) &&
              (image->colorspace == CMYKColorspace))
            pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
          gamma+=(*k)*alpha;
          k++;
          p++;
        }
      }
      gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
      if ((channel & RedChannel) != 0)
        q->red=RoundToQuantum(gamma*pixel.red+image->bias);
      if ((channel & GreenChannel) != 0)
        q->green=RoundToQuantum(gamma*pixel.green+image->bias);
      if ((channel & BlueChannel) != 0)
        q->blue=RoundToQuantum(gamma*pixel.blue+image->bias);
      if ((channel & OpacityChannel) != 0)
        q->opacity=RoundToQuantum(pixel.opacity+image->bias);
      if (((channel & IndexChannel) != 0) &&
          (image->colorspace == CMYKColorspace))
        sharp_indexes[x]=RoundToQuantum(gamma*pixel.index+image->bias);
      q++;
      r++;
    }
    if (SyncImagePixels(sharp_image) == MagickFalse)
      break;
    if ((image->progress_monitor != (MagickProgressMonitor) NULL) &&
        (QuantumTick(y,image->rows) != MagickFalse))
      {
        status=image->progress_monitor(AdaptiveSharpenImageTag,y,image->rows,
          image->client_data);
        if (status == MagickFalse)
          break;
      }
  }
  edge_image=DestroyImage(edge_image);
  for (i=0; i < (long) width;  i+=2)
    kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
  kernel=(double **) RelinquishMagickMemory(kernel);
  return(sharp_image);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     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=GetRandomValue();
  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=GetRandomValue();
      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=GetRandomValue();
      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=GetRandomValue();
        alpha=alpha*beta;
      }
      noise=(MagickRealType) ScaleCharToQuantum((unsigned char)
        (i/SigmaPoisson));
      break;
    }
    case RandomNoise:
    {
      noise=(MagickRealType) QuantumRange*GetRandomValue();
      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
    y;

  MagickBooleanType
    status;

  MagickRealType
    attenuate;

  register const IndexPacket
    *indexes;

  register const PixelPacket
    *pixels;

  register IndexPacket
    *noise_indexes;

  register long
    x;

  register PixelPacket
    *noise_pixels;

  /*
    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);
  for (y=0; y < (long) image->rows; y++)
  {
    pixels=AcquireImagePixels(image,0,y,image->columns,1,exception);
    noise_pixels=GetImagePixels(noise_image,0,y,noise_image->columns,1);
    if ((pixels == (PixelPacket *) NULL) ||
        (noise_pixels == (PixelPacket *) NULL))
      break;
    indexes=AcquireIndexes(image);
    noise_indexes=GetIndexes(noise_image);
    for (x=0; x < (long) image->columns; x++)
    {
      if ((channel & RedChannel) != 0)
        noise_pixels[x].red=GenerateNoise(pixels[x].red,noise_type,attenuate);
      if ((channel & GreenChannel) != 0)
        noise_pixels[x].green=GenerateNoise(pixels[x].green,noise_type,
          attenuate);
      if ((channel & BlueChannel) != 0)
        noise_pixels[x].blue=GenerateNoise(pixels[x].blue,noise_type,attenuate);
      if ((channel & OpacityChannel) != 0)
        noise_pixels[x].opacity=GenerateNoise(pixels[x].opacity,noise_type,
          attenuate);
      if (((channel & IndexChannel) != 0) &&
          (image->colorspace == CMYKColorspace))
        noise_indexes[x]=(IndexPacket) GenerateNoise(indexes[x],noise_type,
          attenuate);
    }
    if (SyncImagePixels(noise_image) == MagickFalse)
      break;
    if ((image->progress_monitor != (MagickProgressMonitor) NULL) &&
        (QuantumTick(y,image->rows) != MagickFalse))
      {
        status=image->progress_monitor(AddNoiseImageTag,y,image->rows,
          image->client_data);
        if (status == MagickFalse)
          break;
      }
  }
  return(noise_image);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     B l u r I m a g e                                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  BlurImage() blurs an image.  We convolve the image with a Gaussian operator
%  of the given radius and standard deviation (sigma).  For reasonable results,
%  the radius should be larger than sigma.  Use a radius of 0 and BlurImage()
%  selects a suitable radius for you.
%
%  BlurImage() differs from GaussianBlurImage() in that it uses a separable
%  kernel which is faster but mathematically equivalent to the non-separable
%  kernel.
%
%  The format of the BlurImage method is:
%
%      Image *BlurImage(const Image *image,const double radius,
%        const double sigma,ExceptionInfo *exception)
%      Image *BlurImageChannel(const Image *image,const ChannelType channel,
%        const double radius,const double sigma,ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: The image.
%
%    o channel: The channel type.
%
%    o radius: The radius of the Gaussian, in pixels, not counting the center
%      pixel.
%
%    o sigma: The standard deviation of the Gaussian, in pixels.
%
%    o exception: Return any errors or warnings in this structure.
%
*/

MagickExport Image *BlurImage(const Image *image,const double radius,
  const double sigma,ExceptionInfo *exception)
{
  Image
    *blur_image;

  blur_image=BlurImageChannel(image,DefaultChannels,radius,sigma,exception);
  return(blur_image);
}

static double *GetBlurKernel(unsigned long width,const MagickRealType sigma)
{
#define KernelRank 3

  double
    *kernel;

  long
    bias;

  MagickRealType
    alpha,
    normalize;

  register long
    i;

  /*
    Generate a 1-D convolution kernel.
  */
  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
  kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
  if (kernel == (double *) NULL)
    return(0);
  (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
  bias=KernelRank*(long) width/2;
  for (i=(-bias); i <= bias; i++)
  {
    alpha=exp((-((double) (i*i))/(double) (2.0*KernelRank*KernelRank*
      sigma*sigma)));
    kernel[(i+bias)/KernelRank]+=(double) (alpha/(MagickSQ2PI*sigma));
  }
  normalize=0.0;
  for (i=0; i < (long) width; i++)
    normalize+=kernel[i];
  for (i=0; i < (long) width; i++)
    kernel[i]/=normalize;
  return(kernel);
}

MagickExport Image *BlurImageChannel(const Image *image,
  const ChannelType channel,const double radius,const double sigma,
  ExceptionInfo *exception)
{
#define BlurImageTag  "Blur/Image"

  double
    *kernel;

  Image
    *blur_image;

  IndexPacket
    *blur_indexes;

  long
    y;

  MagickBooleanType
    status;

  MagickPixelPacket
    pixel;

  MagickOffsetType
    offset;

  MagickRealType
    alpha,
    gamma;

  register const IndexPacket
    *indexes;

  register const PixelPacket
    *pixels;

  register const double
    *k;

  register long
    i,
    x;

  register PixelPacket
    *blur_pixels;

  unsigned long
    width;

  ViewInfo
    *image_view;

  /*
    Initialize blur image attributes.
  */
  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickM