root / ImageMagick / trunk / 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-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/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/monitor-private.h"
#include "magick/montage.h"
#include "magick/paint.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-private.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"
#define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)

  double
    **kernel;

  Image
    *blur_image,
    *edge_image,
    *gaussian_image;

  long
    j,
    progress,
    y;

  MagickPixelPacket
    zero;

  MagickRealType
    alpha,
    normalize;

  register long
    i,
    u,
    v;

  unsigned long
    width;

  ViewInfo
    **blur_view,
    **edge_view,
    **image_view;

  volatile MagickBooleanType
    status;

  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,image->columns,image->rows,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*MagickSigma*MagickSigma));
        kernel[i][j]=(double) (alpha/(2.0*MagickPI*MagickSigma*MagickSigma));
        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.
  */
  status=MagickTrue;
  progress=0;
  GetMagickPixelPacket(image,&zero);
  image_view=AcquireCacheViewThreadSet(image);
  edge_view=AcquireCacheViewThreadSet(edge_image);
  blur_view=AcquireCacheViewThreadSet(blur_image);
  #pragma omp parallel for schedule(static,64)
  for (y=0; y < (long) blur_image->rows; y++)
  {
    IndexPacket
      *blur_indexes,
      *indexes;

    register const PixelPacket
      *p,
      *r;

    register long
      id,
      x;

    register PixelPacket
      *q;

    if (status == MagickFalse)
      continue;
    id=GetCacheViewThreadId();
    r=AcquireCacheViewPixels(edge_view[id],0,y,edge_image->columns,1,
      exception);
    q=SetCacheViewPixels(blur_view[id],0,y,blur_image->columns,1);
    if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    blur_indexes=GetCacheViewIndexes(blur_view[id]);
    for (x=0; x < (long) blur_image->columns; x++)
    {
      MagickPixelPacket
        pixel;

      MagickRealType
        alpha,
        gamma;

      register const double
        *k;

      register long
        i,
        u,
        v;

      gamma=0.0;
      i=(long) (width*QuantumScale*PixelIntensity(r)+0.5);
      if (i < 0)
        i=0;
      else
        if (i > (long) width)
          i=(long) width;
      if ((i & 0x01) != 0)
        i--;
      p=AcquireCacheViewPixels(image_view[id],x-((long) (width-i)/2L),y-(long)
        ((width-i)/2L),width-i,width-i,exception);
      if (p == (const PixelPacket *) NULL)
        break;
      indexes=GetCacheViewIndexes(image_view[id]);
      pixel=zero;
      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=(MagickRealType) (QuantumScale*(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 (SyncCacheViewPixels(blur_view[id]) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

        #pragma omp critical
        proceed=SetImageProgress(image,AdaptiveBlurImageTag,progress++,
          image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  blur_image->type=image->type;
  blur_view=DestroyCacheViewThreadSet(blur_view);
  edge_view=DestroyCacheViewThreadSet(edge_view);
  image_view=DestroyCacheViewThreadSet(image_view);
  edge_image=DestroyImage(edge_image);
  for (i=0; i < (long) width;  i+=2)
    kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
  kernel=(double **) RelinquishMagickMemory(kernel);
  if (status == MagickFalse)
    blur_image=DestroyImage(blur_image);
  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"
#define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)

  double
    **kernel;

  Image
    *sharp_image,
    *edge_image,
    *gaussian_image;

  long
    j,
    progress,
    y;

  MagickPixelPacket
    zero;

  MagickRealType
    alpha,
    normalize;

  register long
    i,
    u,
    v;

  unsigned long
    width;

  ViewInfo
    **sharp_view,
    **edge_view,
    **image_view;

  volatile MagickBooleanType
    status;

  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,image->columns,image->rows,MagickTrue,exception);
  if (sharp_image == (Image *) NULL)
    return((Image *) NULL);
  if (fabs(sigma) <= MagickEpsilon)
    return(sharp_image);
  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, sharp, 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%");
  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);
      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*MagickSigma*MagickSigma));
        kernel[i][j]=(double) (-alpha/(2.0*MagickPI*MagickSigma*MagickSigma));
        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 sharp image.
  */
  status=MagickTrue;
  progress=0;
  GetMagickPixelPacket(image,&zero);
  image_view=AcquireCacheViewThreadSet(image);
  edge_view=AcquireCacheViewThreadSet(edge_image);
  sharp_view=AcquireCacheViewThreadSet(sharp_image);
  #pragma omp parallel for schedule(static,64)
  for (y=0; y < (long) sharp_image->rows; y++)
  {
    IndexPacket
      *sharp_indexes,
      *indexes;

    register const PixelPacket
      *p,
      *r;

    register long
      id,
      x;

    register PixelPacket
      *q;

    if (status == MagickFalse)
      continue;
    id=GetCacheViewThreadId();
    r=AcquireCacheViewPixels(edge_view[id],0,y,edge_image->columns,1,
      exception);
    q=SetCacheViewPixels(sharp_view[id],0,y,sharp_image->columns,1);
    if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    sharp_indexes=GetCacheViewIndexes(sharp_view[id]);
    for (x=0; x < (long) sharp_image->columns; x++)
    {
      MagickPixelPacket
        pixel;

      MagickRealType
        alpha,
        gamma;

      register const double
        *k;

      register long
        i,
        u,
        v;

      gamma=0.0;
      i=(long) (width*(QuantumRange-QuantumScale*PixelIntensity(r))+0.5);
      if (i < 0)
        i=0;
      else
        if (i > (long) width)
          i=(long) width;
      if ((i & 0x01) != 0)
        i--;
      p=AcquireCacheViewPixels(image_view[id],x-((long) (width-i)/2L),y-(long)
        ((width-i)/2L),width-i,width-i,exception);
      if (p == (const PixelPacket *) NULL)
        break;
      indexes=GetCacheViewIndexes(image_view[id]);
      k=kernel[i];
      pixel=zero;
      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=(MagickRealType) (QuantumScale*(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 (SyncCacheViewPixels(sharp_view[id]) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

        #pragma omp critical
        proceed=SetImageProgress(image,AdaptiveSharpenImageTag,progress++,
          image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  sharp_image->type=image->type;
  sharp_view=DestroyCacheViewThreadSet(sharp_view);
  edge_view=DestroyCacheViewThreadSet(edge_view);
  image_view=DestroyCacheViewThreadSet(image_view);
  edge_image=DestroyImage(edge_image);
  for (i=0; i < (long) width;  i+=2)
    kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
  kernel=(double **) RelinquishMagickMemory(kernel);
  if (status == MagickFalse)
    sharp_image=DestroyImage(sharp_image);
  return(sharp_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*
      MagickSigma*MagickSigma)));
    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;

  long
    progress,
    x,
    y;

  MagickPixelPacket
    zero;

  MagickRealType
    bias;

  register long
    i;

  unsigned long
    width;

  ViewInfo
    **blur_view,
    **image_view;

  volatile MagickBooleanType
    status;

  /*
    Initialize blur 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);
  blur_image=CloneImage(image,image->columns,image->rows,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);
    }
  width=GetOptimalKernelWidth1D(radius,sigma);
  kernel=GetBlurKernel(width,sigma);
  if (kernel == (double *) NULL)
    {
      blur_image=DestroyImage(blur_image);
      return((Image *) NULL);
      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
    }
  if (image->debug != MagickFalse)
    {
      char
        format[MaxTextExtent],
        *message;

      register const double
        *k;

      (void) LogMagickEvent(TransformEvent,GetMagickModule(),
        "  BlurImage with %ld kernel:",width);
      message=AcquireString("");
      k=kernel;
      for (i=0; i < (long) width; i++)
      {
        *message='\0';
        (void) FormatMagickString(format,MaxTextExtent,"%ld: ",i);
        (void) ConcatenateString(&message,format);
        (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
        (void) ConcatenateString(&message,format);
        (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
      }
      message=DestroyString(message);
    }
  /*
    Blur rows.
  */
  status=MagickTrue;
  progress=0;
  GetMagickPixelPacket(image,&zero);
  bias=image->bias;
  image_view=AcquireCacheViewThreadSet(image);
  blur_view=AcquireCacheViewThreadSet(blur_image);
  #pragma omp parallel for schedule(static,64)
  for (y=0; y < (long) blur_image->rows; y++)
  {
    const IndexPacket
      *indexes;

    IndexPacket
      *blur_indexes;

    register const PixelPacket
      *p;

    register long
      id,
      x;

    register PixelPacket
      *q;

    if (status == MagickFalse)
      continue;
    id=GetCacheViewThreadId();
    p=AcquireCacheViewPixels(image_view[id],-((long) width/2L),y,image->columns+
      width,1,exception);
    q=SetCacheViewPixels(blur_view[id],0,y,blur_image->columns,1);
    if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    indexes=AcquireCacheViewIndexes(image_view[id]);
    blur_indexes=GetCacheViewIndexes(blur_view[id]);
    for (x=0; x < (long) blur_image->columns; x++)
    {
      MagickPixelPacket
        pixel;

      register const double
        *k;

      register long
        i;

      pixel=zero;
      k=kernel;
      if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
        {
          for (i=0; i < (long) width; i++)
          {
            pixel.red+=(*k)*(p+i)->red;
            pixel.green+=(*k)*(p+i)->green;
            pixel.blue+=(*k)*(p+i)->blue;
            pixel.opacity+=(*k)*(p+i)->opacity;
            k++;
          }
          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=kernel;
              for (i=0; i < (long) width; i++)
              {
                pixel.index+=(*k)*indexes[x+i];
                k++;
              }
              blur_indexes[x]=RoundToQuantum(pixel.index+bias);
            }
        }
      else
        {
          MagickRealType
            alpha,
            gamma;

          gamma=0.0;
          for (i=0; i < (long) width; i++)
          {
            alpha=(MagickRealType) (QuantumScale*(QuantumRange-(p+i)->opacity));
            pixel.red+=(*k)*alpha*(p+i)->red;
            pixel.green+=(*k)*alpha*(p+i)->green;
            pixel.blue+=(*k)*alpha*(p+i)->blue;
            pixel.opacity+=(*k)*(p+i)->opacity;
            gamma+=(*k)*alpha;
            k++;
          }
          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=kernel;
              for (i=0; i < (long) width; i++)
              {
                alpha=(MagickRealType) (QuantumScale*(QuantumRange-
                  (p+i)->opacity));
                pixel.index+=(*k)*alpha*indexes[x+i];
                k++;
              }
              blur_indexes[x]=RoundToQuantum(gamma*pixel.index+bias);
            }
        }
      p++;
      q++;
    }
    if (SyncCacheViewPixels(blur_view[id]) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

        #pragma omp critical
        proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
          blur_image->columns);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  /*
    Blur columns.
  */
  image_view=DestroyCacheViewThreadSet(image_view);
  image_view=AcquireCacheViewThreadSet(blur_image);