granicus.if.org Git - imagemagick/blob - magick/effect.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %                   EEEEE  FFFFF  FFFFF  EEEEE  CCCC  TTTTT                   %
   7 %                   E      F      F      E     C        T                     %
   8 %                   EEE    FFF    FFF    EEE   C        T                     %
   9 %                   E      F      F      E     C        T                     %
  10 %                   EEEEE  F      F      EEEEE  CCCC    T                     %
  11 %                                                                             %
  12 %                                                                             %
  13 %                       MagickCore Image Effects Methods                      %
  14 %                                                                             %
  15 %                               Software Design                               %
  16 %                                 John Cristy                                 %
  17 %                                 October 1996                                %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "magick/studio.h"
  44 #include "magick/accelerate.h"
  45 #include "magick/blob.h"
  46 #include "magick/cache-view.h"
  47 #include "magick/color.h"
  48 #include "magick/color-private.h"
  49 #include "magick/colorspace.h"
  50 #include "magick/constitute.h"
  51 #include "magick/decorate.h"
  52 #include "magick/draw.h"
  53 #include "magick/enhance.h"
  54 #include "magick/exception.h"
  55 #include "magick/exception-private.h"
  56 #include "magick/effect.h"
  57 #include "magick/fx.h"
  58 #include "magick/gem.h"
  59 #include "magick/geometry.h"
  60 #include "magick/image-private.h"
  61 #include "magick/list.h"
  62 #include "magick/log.h"
  63 #include "magick/memory_.h"
  64 #include "magick/monitor.h"
  65 #include "magick/monitor-private.h"
  66 #include "magick/montage.h"
  67 #include "magick/morphology.h"
  68 #include "magick/paint.h"
  69 #include "magick/pixel-private.h"
  70 #include "magick/property.h"
  71 #include "magick/quantize.h"
  72 #include "magick/quantum.h"
  73 #include "magick/random_.h"
  74 #include "magick/random-private.h"
  75 #include "magick/resample.h"
  76 #include "magick/resample-private.h"
  77 #include "magick/resize.h"
  78 #include "magick/resource_.h"
  79 #include "magick/segment.h"
  80 #include "magick/shear.h"
  81 #include "magick/signature-private.h"
  82 #include "magick/string_.h"
  83 #include "magick/thread-private.h"
  84 #include "magick/transform.h"
  85 #include "magick/threshold.h"
  86 \f
  87 /*
  88 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  89 %                                                                             %
  90 %                                                                             %
  91 %                                                                             %
  92 %     A d a p t i v e B l u r I m a g e                                       %
  93 %                                                                             %
  94 %                                                                             %
  95 %                                                                             %
  96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  97 %
  98 %  AdaptiveBlurImage() adaptively blurs the image by blurring less
  99 %  intensely near image edges and more intensely far from edges.  We blur the
 100 %  image with a Gaussian operator of the given radius and standard deviation
 101 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
 102 %  radius of 0 and AdaptiveBlurImage() selects a suitable radius for you.
 103 %
 104 %  The format of the AdaptiveBlurImage method is:
 105 %
 106 %      Image *AdaptiveBlurImage(const Image *image,const double radius,
 107 %        const double sigma,ExceptionInfo *exception)
 108 %      Image *AdaptiveBlurImageChannel(const Image *image,
 109 %        const ChannelType channel,double radius,const double sigma,
 110 %        ExceptionInfo *exception)
 111 %
 112 %  A description of each parameter follows:
 113 %
 114 %    o image: the image.
 115 %
 116 %    o channel: the channel type.
 117 %
 118 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 119 %      pixel.
 120 %
 121 %    o sigma: the standard deviation of the Laplacian, in pixels.
 122 %
 123 %    o exception: return any errors or warnings in this structure.
 124 %
 125 */
 126
 127 MagickExport Image *AdaptiveBlurImage(const Image *image,const double radius,
 128   const double sigma,ExceptionInfo *exception)
 129 {
 130   Image
 131     *blur_image;
 132
 133   blur_image=AdaptiveBlurImageChannel(image,DefaultChannels,radius,sigma,
 134     exception);
 135   return(blur_image);
 136 }
 137
 138 MagickExport Image *AdaptiveBlurImageChannel(const Image *image,
 139   const ChannelType channel,const double radius,const double sigma,
 140   ExceptionInfo *exception)
 141 {
 142 #define AdaptiveBlurImageTag  "Convolve/Image"
 143 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
 144
 145   CacheView
 146     *blur_view,
 147     *edge_view,
 148     *image_view;
 149
 150   double
 151     **kernel,
 152     normalize;
 153
 154   Image
 155     *blur_image,
 156     *edge_image,
 157     *gaussian_image;
 158
 159   MagickBooleanType
 160     status;
 161
 162   MagickOffsetType
 163     progress;
 164
 165   MagickPixelPacket
 166     bias;
 167
 168   register ssize_t
 169     i;
 170
 171   size_t
 172     width;
 173
 174   ssize_t
 175     j,
 176     k,
 177     u,
 178     v,
 179     y;
 180
 181   assert(image != (const Image *) NULL);
 182   assert(image->signature == MagickSignature);
 183   if (image->debug != MagickFalse)
 184     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 185   assert(exception != (ExceptionInfo *) NULL);
 186   assert(exception->signature == MagickSignature);
 187   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
 188   if (blur_image == (Image *) NULL)
 189     return((Image *) NULL);
 190   if (fabs(sigma) <= MagickEpsilon)
 191     return(blur_image);
 192   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
 193     {
 194       InheritException(exception,&blur_image->exception);
 195       blur_image=DestroyImage(blur_image);
 196       return((Image *) NULL);
 197     }
 198   /*
 199     Edge detect the image brighness channel, level, blur, and level again.
 200   */
 201   edge_image=EdgeImage(image,radius,exception);
 202   if (edge_image == (Image *) NULL)
 203     {
 204       blur_image=DestroyImage(blur_image);
 205       return((Image *) NULL);
 206     }
 207   (void) LevelImage(edge_image,"20%,95%");
 208   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
 209   if (gaussian_image != (Image *) NULL)
 210     {
 211       edge_image=DestroyImage(edge_image);
 212       edge_image=gaussian_image;
 213     }
 214   (void) LevelImage(edge_image,"10%,95%");
 215   /*
 216     Create a set of kernels from maximum (radius,sigma) to minimum.
 217   */
 218   width=GetOptimalKernelWidth2D(radius,sigma);
 219   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 220   if (kernel == (double **) NULL)
 221     {
 222       edge_image=DestroyImage(edge_image);
 223       blur_image=DestroyImage(blur_image);
 224       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 225     }
 226   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
 227   for (i=0; i < (ssize_t) width; i+=2)
 228   {
 229     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
 230       sizeof(**kernel));
 231     if (kernel[i] == (double *) NULL)
 232       break;
 233     normalize=0.0;
 234     j=(ssize_t) (width-i)/2;
 235     k=0;
 236     for (v=(-j); v <= j; v++)
 237     {
 238       for (u=(-j); u <= j; u++)
 239       {
 240         kernel[i][k]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
 241           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
 242         normalize+=kernel[i][k];
 243         k++;
 244       }
 245     }
 246     if (fabs(normalize) <= MagickEpsilon)
 247       normalize=1.0;
 248     normalize=1.0/normalize;
 249     for (k=0; k < (j*j); k++)
 250       kernel[i][k]=normalize*kernel[i][k];
 251   }
 252   if (i < (ssize_t) width)
 253     {
 254       for (i-=2; i >= 0; i-=2)
 255         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 256       kernel=(double **) RelinquishMagickMemory(kernel);
 257       edge_image=DestroyImage(edge_image);
 258       blur_image=DestroyImage(blur_image);
 259       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 260     }
 261   /*
 262     Adaptively blur image.
 263   */
 264   status=MagickTrue;
 265   progress=0;
 266   GetMagickPixelPacket(image,&bias);
 267   SetMagickPixelPacketBias(image,&bias);
 268   image_view=AcquireCacheView(image);
 269   edge_view=AcquireCacheView(edge_image);
 270   blur_view=AcquireCacheView(blur_image);
 271 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 272   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 273 #endif
 274   for (y=0; y < (ssize_t) blur_image->rows; y++)
 275   {
 276     register const IndexPacket
 277       *restrict indexes;
 278
 279     register const PixelPacket
 280       *restrict p,
 281       *restrict r;
 282
 283     register IndexPacket
 284       *restrict blur_indexes;
 285
 286     register PixelPacket
 287       *restrict q;
 288
 289     register ssize_t
 290       x;
 291
 292     if (status == MagickFalse)
 293       continue;
 294     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
 295     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
 296       exception);
 297     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 298       {
 299         status=MagickFalse;
 300         continue;
 301       }
 302     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
 303     for (x=0; x < (ssize_t) blur_image->columns; x++)
 304     {
 305       MagickPixelPacket
 306         pixel;
 307
 308       MagickRealType
 309         alpha,
 310         gamma;
 311
 312       register const double
 313         *restrict k;
 314
 315       register ssize_t
 316         i,
 317         u,
 318         v;
 319
 320       gamma=0.0;
 321       i=(ssize_t) ceil((double) width*QuantumScale*PixelIntensity(r)-0.5);
 322       if (i < 0)
 323         i=0;
 324       else
 325         if (i > (ssize_t) width)
 326           i=(ssize_t) width;
 327       if ((i & 0x01) != 0)
 328         i--;
 329       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-i)/2L),y-
 330         (ssize_t) ((width-i)/2L),width-i,width-i,exception);
 331       if (p == (const PixelPacket *) NULL)
 332         break;
 333       indexes=GetCacheViewVirtualIndexQueue(image_view);
 334       pixel=bias;
 335       k=kernel[i];
 336       for (v=0; v < (ssize_t) (width-i); v++)
 337       {
 338         for (u=0; u < (ssize_t) (width-i); u++)
 339         {
 340           alpha=1.0;
 341           if (((channel & OpacityChannel) != 0) &&
 342               (image->matte != MagickFalse))
 343             alpha=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(p));
 344           if ((channel & RedChannel) != 0)
 345             pixel.red+=(*k)*alpha*GetRedPixelComponent(p);
 346           if ((channel & GreenChannel) != 0)
 347             pixel.green+=(*k)*alpha*GetGreenPixelComponent(p);
 348           if ((channel & BlueChannel) != 0)
 349             pixel.blue+=(*k)*alpha*GetBluePixelComponent(p);
 350           if ((channel & OpacityChannel) != 0)
 351             pixel.opacity+=(*k)*GetOpacityPixelComponent(p);
 352           if (((channel & IndexChannel) != 0) &&
 353               (image->colorspace == CMYKColorspace))
 354             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
 355           gamma+=(*k)*alpha;
 356           k++;
 357           p++;
 358         }
 359       }
 360       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 361       if ((channel & RedChannel) != 0)
 362         q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
 363       if ((channel & GreenChannel) != 0)
 364         q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
 365       if ((channel & BlueChannel) != 0)
 366         q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
 367       if ((channel & OpacityChannel) != 0)
 368         SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
 369       if (((channel & IndexChannel) != 0) &&
 370           (image->colorspace == CMYKColorspace))
 371         blur_indexes[x]=ClampToQuantum(gamma*GetIndexPixelComponent(&pixel));
 372       q++;
 373       r++;
 374     }
 375     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
 376       status=MagickFalse;
 377     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 378       {
 379         MagickBooleanType
 380           proceed;
 381
 382 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 383   #pragma omp critical (MagickCore_AdaptiveBlurImageChannel)
 384 #endif
 385         proceed=SetImageProgress(image,AdaptiveBlurImageTag,progress++,
 386           image->rows);
 387         if (proceed == MagickFalse)
 388           status=MagickFalse;
 389       }
 390   }
 391   blur_image->type=image->type;
 392   blur_view=DestroyCacheView(blur_view);
 393   edge_view=DestroyCacheView(edge_view);
 394   image_view=DestroyCacheView(image_view);
 395   edge_image=DestroyImage(edge_image);
 396   for (i=0; i < (ssize_t) width;  i+=2)
 397     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 398   kernel=(double **) RelinquishMagickMemory(kernel);
 399   if (status == MagickFalse)
 400     blur_image=DestroyImage(blur_image);
 401   return(blur_image);
 402 }
 403 \f
 404 /*
 405 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 406 %                                                                             %
 407 %                                                                             %
 408 %                                                                             %
 409 %     A d a p t i v e S h a r p e n I m a g e                                 %
 410 %                                                                             %
 411 %                                                                             %
 412 %                                                                             %
 413 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 414 %
 415 %  AdaptiveSharpenImage() adaptively sharpens the image by sharpening more
 416 %  intensely near image edges and less intensely far from edges. We sharpen the
 417 %  image with a Gaussian operator of the given radius and standard deviation
 418 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
 419 %  radius of 0 and AdaptiveSharpenImage() selects a suitable radius for you.
 420 %
 421 %  The format of the AdaptiveSharpenImage method is:
 422 %
 423 %      Image *AdaptiveSharpenImage(const Image *image,const double radius,
 424 %        const double sigma,ExceptionInfo *exception)
 425 %      Image *AdaptiveSharpenImageChannel(const Image *image,
 426 %        const ChannelType channel,double radius,const double sigma,
 427 %        ExceptionInfo *exception)
 428 %
 429 %  A description of each parameter follows:
 430 %
 431 %    o image: the image.
 432 %
 433 %    o channel: the channel type.
 434 %
 435 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 436 %      pixel.
 437 %
 438 %    o sigma: the standard deviation of the Laplacian, in pixels.
 439 %
 440 %    o exception: return any errors or warnings in this structure.
 441 %
 442 */
 443
 444 MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
 445   const double sigma,ExceptionInfo *exception)
 446 {
 447   Image
 448     *sharp_image;
 449
 450   sharp_image=AdaptiveSharpenImageChannel(image,DefaultChannels,radius,sigma,
 451     exception);
 452   return(sharp_image);
 453 }
 454
 455 MagickExport Image *AdaptiveSharpenImageChannel(const Image *image,
 456   const ChannelType channel,const double radius,const double sigma,
 457   ExceptionInfo *exception)
 458 {
 459 #define AdaptiveSharpenImageTag  "Convolve/Image"
 460 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
 461
 462   CacheView
 463     *sharp_view,
 464     *edge_view,
 465     *image_view;
 466
 467   double
 468     **kernel,
 469     normalize;
 470
 471   Image
 472     *sharp_image,
 473     *edge_image,
 474     *gaussian_image;
 475
 476   MagickBooleanType
 477     status;
 478
 479   MagickOffsetType
 480     progress;
 481
 482   MagickPixelPacket
 483     bias;
 484
 485   register ssize_t
 486     i;
 487
 488   size_t
 489     width;
 490
 491   ssize_t
 492     j,
 493     k,
 494     u,
 495     v,
 496     y;
 497
 498   assert(image != (const Image *) NULL);
 499   assert(image->signature == MagickSignature);
 500   if (image->debug != MagickFalse)
 501     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 502   assert(exception != (ExceptionInfo *) NULL);
 503   assert(exception->signature == MagickSignature);
 504   sharp_image=CloneImage(image,0,0,MagickTrue,exception);
 505   if (sharp_image == (Image *) NULL)
 506     return((Image *) NULL);
 507   if (fabs(sigma) <= MagickEpsilon)
 508     return(sharp_image);
 509   if (SetImageStorageClass(sharp_image,DirectClass) == MagickFalse)
 510     {
 511       InheritException(exception,&sharp_image->exception);
 512       sharp_image=DestroyImage(sharp_image);
 513       return((Image *) NULL);
 514     }
 515   /*
 516     Edge detect the image brighness channel, level, sharp, and level again.
 517   */
 518   edge_image=EdgeImage(image,radius,exception);
 519   if (edge_image == (Image *) NULL)
 520     {
 521       sharp_image=DestroyImage(sharp_image);
 522       return((Image *) NULL);
 523     }
 524   (void) LevelImage(edge_image,"20%,95%");
 525   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
 526   if (gaussian_image != (Image *) NULL)
 527     {
 528       edge_image=DestroyImage(edge_image);
 529       edge_image=gaussian_image;
 530     }
 531   (void) LevelImage(edge_image,"10%,95%");
 532   /*
 533     Create a set of kernels from maximum (radius,sigma) to minimum.
 534   */
 535   width=GetOptimalKernelWidth2D(radius,sigma);
 536   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 537   if (kernel == (double **) NULL)
 538     {
 539       edge_image=DestroyImage(edge_image);
 540       sharp_image=DestroyImage(sharp_image);
 541       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 542     }
 543   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
 544   for (i=0; i < (ssize_t) width; i+=2)
 545   {
 546     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
 547       sizeof(**kernel));
 548     if (kernel[i] == (double *) NULL)
 549       break;
 550     normalize=0.0;
 551     j=(ssize_t) (width-i)/2;
 552     k=0;
 553     for (v=(-j); v <= j; v++)
 554     {
 555       for (u=(-j); u <= j; u++)
 556       {
 557         kernel[i][k]=(double) (-exp(-((double) u*u+v*v)/(2.0*MagickSigma*
 558           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
 559         normalize+=kernel[i][k];
 560         k++;
 561       }
 562     }
 563     if (fabs(normalize) <= MagickEpsilon)
 564       normalize=1.0;
 565     normalize=1.0/normalize;
 566     for (k=0; k < (j*j); k++)
 567       kernel[i][k]=normalize*kernel[i][k];
 568   }
 569   if (i < (ssize_t) width)
 570     {
 571       for (i-=2; i >= 0; i-=2)
 572         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 573       kernel=(double **) RelinquishMagickMemory(kernel);
 574       edge_image=DestroyImage(edge_image);
 575       sharp_image=DestroyImage(sharp_image);
 576       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 577     }
 578   /*
 579     Adaptively sharpen image.
 580   */
 581   status=MagickTrue;
 582   progress=0;
 583   GetMagickPixelPacket(image,&bias);
 584   SetMagickPixelPacketBias(image,&bias);
 585   image_view=AcquireCacheView(image);
 586   edge_view=AcquireCacheView(edge_image);
 587   sharp_view=AcquireCacheView(sharp_image);
 588 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 589   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 590 #endif
 591   for (y=0; y < (ssize_t) sharp_image->rows; y++)
 592   {
 593     register const IndexPacket
 594       *restrict indexes;
 595
 596     register const PixelPacket
 597       *restrict p,
 598       *restrict r;
 599
 600     register IndexPacket
 601       *restrict sharp_indexes;
 602
 603     register PixelPacket
 604       *restrict q;
 605
 606     register ssize_t
 607       x;
 608
 609     if (status == MagickFalse)
 610       continue;
 611     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
 612     q=QueueCacheViewAuthenticPixels(sharp_view,0,y,sharp_image->columns,1,
 613       exception);
 614     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 615       {
 616         status=MagickFalse;
 617         continue;
 618       }
 619     sharp_indexes=GetCacheViewAuthenticIndexQueue(sharp_view);
 620     for (x=0; x < (ssize_t) sharp_image->columns; x++)
 621     {
 622       MagickPixelPacket
 623         pixel;
 624
 625       MagickRealType
 626         alpha,
 627         gamma;
 628
 629       register const double
 630         *restrict k;
 631
 632       register ssize_t
 633         i,
 634         u,
 635         v;
 636
 637       gamma=0.0;
 638       i=(ssize_t) ceil((double) width*(QuantumRange-QuantumScale*
 639         PixelIntensity(r))-0.5);
 640       if (i < 0)
 641         i=0;
 642       else
 643         if (i > (ssize_t) width)
 644           i=(ssize_t) width;
 645       if ((i & 0x01) != 0)
 646         i--;
 647       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-i)/2L),y-
 648         (ssize_t) ((width-i)/2L),width-i,width-i,exception);
 649       if (p == (const PixelPacket *) NULL)
 650         break;
 651       indexes=GetCacheViewVirtualIndexQueue(image_view);
 652       k=kernel[i];
 653       pixel=bias;
 654       for (v=0; v < (ssize_t) (width-i); v++)
 655       {
 656         for (u=0; u < (ssize_t) (width-i); u++)
 657         {
 658           alpha=1.0;
 659           if (((channel & OpacityChannel) != 0) &&
 660               (image->matte != MagickFalse))
 661             alpha=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(p));
 662           if ((channel & RedChannel) != 0)
 663             pixel.red+=(*k)*alpha*GetRedPixelComponent(p);
 664           if ((channel & GreenChannel) != 0)
 665             pixel.green+=(*k)*alpha*GetGreenPixelComponent(p);
 666           if ((channel & BlueChannel) != 0)
 667             pixel.blue+=(*k)*alpha*GetBluePixelComponent(p);
 668           if ((channel & OpacityChannel) != 0)
 669             pixel.opacity+=(*k)*GetOpacityPixelComponent(p);
 670           if (((channel & IndexChannel) != 0) &&
 671               (image->colorspace == CMYKColorspace))
 672             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
 673           gamma+=(*k)*alpha;
 674           k++;
 675           p++;
 676         }
 677       }
 678       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 679       if ((channel & RedChannel) != 0)
 680         q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
 681       if ((channel & GreenChannel) != 0)
 682         q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
 683       if ((channel & BlueChannel) != 0)
 684         q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
 685       if ((channel & OpacityChannel) != 0)
 686         SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
 687       if (((channel & IndexChannel) != 0) &&
 688           (image->colorspace == CMYKColorspace))
 689         sharp_indexes[x]=ClampToQuantum(gamma*GetIndexPixelComponent(&pixel));
 690       q++;
 691       r++;
 692     }
 693     if (SyncCacheViewAuthenticPixels(sharp_view,exception) == MagickFalse)
 694       status=MagickFalse;
 695     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 696       {
 697         MagickBooleanType
 698           proceed;
 699
 700 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 701   #pragma omp critical (MagickCore_AdaptiveSharpenImageChannel)
 702 #endif
 703         proceed=SetImageProgress(image,AdaptiveSharpenImageTag,progress++,
 704           image->rows);
 705         if (proceed == MagickFalse)
 706           status=MagickFalse;
 707       }
 708   }
 709   sharp_image->type=image->type;
 710   sharp_view=DestroyCacheView(sharp_view);
 711   edge_view=DestroyCacheView(edge_view);
 712   image_view=DestroyCacheView(image_view);
 713   edge_image=DestroyImage(edge_image);
 714   for (i=0; i < (ssize_t) width;  i+=2)
 715     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 716   kernel=(double **) RelinquishMagickMemory(kernel);
 717   if (status == MagickFalse)
 718     sharp_image=DestroyImage(sharp_image);
 719   return(sharp_image);
 720 }
 721 \f
 722 /*
 723 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 724 %                                                                             %
 725 %                                                                             %
 726 %                                                                             %
 727 %     B l u r I m a g e                                                       %
 728 %                                                                             %
 729 %                                                                             %
 730 %                                                                             %
 731 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 732 %
 733 %  BlurImage() blurs an image.  We convolve the image with a Gaussian operator
 734 %  of the given radius and standard deviation (sigma).  For reasonable results,
 735 %  the radius should be larger than sigma.  Use a radius of 0 and BlurImage()
 736 %  selects a suitable radius for you.
 737 %
 738 %  BlurImage() differs from GaussianBlurImage() in that it uses a separable
 739 %  kernel which is faster but mathematically equivalent to the non-separable
 740 %  kernel.
 741 %
 742 %  The format of the BlurImage method is:
 743 %
 744 %      Image *BlurImage(const Image *image,const double radius,
 745 %        const double sigma,ExceptionInfo *exception)
 746 %      Image *BlurImageChannel(const Image *image,const ChannelType channel,
 747 %        const double radius,const double sigma,ExceptionInfo *exception)
 748 %
 749 %  A description of each parameter follows:
 750 %
 751 %    o image: the image.
 752 %
 753 %    o channel: the channel type.
 754 %
 755 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 756 %      pixel.
 757 %
 758 %    o sigma: the standard deviation of the Gaussian, in pixels.
 759 %
 760 %    o exception: return any errors or warnings in this structure.
 761 %
 762 */
 763
 764 MagickExport Image *BlurImage(const Image *image,const double radius,
 765   const double sigma,ExceptionInfo *exception)
 766 {
 767   Image
 768     *blur_image;
 769
 770   blur_image=BlurImageChannel(image,DefaultChannels,radius,sigma,exception);
 771   return(blur_image);
 772 }
 773
 774 static double *GetBlurKernel(const size_t width,const double sigma)
 775 {
 776   double
 777     *kernel,
 778     normalize;
 779
 780   register ssize_t
 781     i;
 782
 783   ssize_t
 784     j,
 785     k;
 786
 787   /*
 788     Generate a 1-D convolution kernel.
 789   */
 790   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
 791   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 792   if (kernel == (double *) NULL)
 793     return(0);
 794   normalize=0.0;
 795   j=(ssize_t) width/2;
 796   i=0;
 797   for (k=(-j); k <= j; k++)
 798   {
 799     kernel[i]=(double) (exp(-((double) k*k)/(2.0*MagickSigma*MagickSigma))/
 800       (MagickSQ2PI*MagickSigma));
 801     normalize+=kernel[i];
 802     i++;
 803   }
 804   for (i=0; i < (ssize_t) width; i++)
 805     kernel[i]/=normalize;
 806   return(kernel);
 807 }
 808
 809 MagickExport Image *BlurImageChannel(const Image *image,
 810   const ChannelType channel,const double radius,const double sigma,
 811   ExceptionInfo *exception)
 812 {
 813 #define BlurImageTag  "Blur/Image"
 814
 815   CacheView
 816     *blur_view,
 817     *image_view;
 818
 819   double
 820     *kernel;
 821
 822   Image
 823     *blur_image;
 824
 825   MagickBooleanType
 826     status;
 827
 828   MagickOffsetType
 829     progress;
 830
 831   MagickPixelPacket
 832     bias;
 833
 834   register ssize_t
 835     i;
 836
 837   size_t
 838     width;
 839
 840   ssize_t
 841     x,
 842     y;
 843
 844   /*
 845     Initialize blur image attributes.
 846   */
 847   assert(image != (Image *) NULL);
 848   assert(image->signature == MagickSignature);
 849   if (image->debug != MagickFalse)
 850     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 851   assert(exception != (ExceptionInfo *) NULL);
 852   assert(exception->signature == MagickSignature);
 853   blur_image=CloneImage(image,0,0,MagickTrue,exception);
 854   if (blur_image == (Image *) NULL)
 855     return((Image *) NULL);
 856   if (fabs(sigma) <= MagickEpsilon)
 857     return(blur_image);
 858   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
 859     {
 860       InheritException(exception,&blur_image->exception);
 861       blur_image=DestroyImage(blur_image);
 862       return((Image *) NULL);
 863     }
 864   width=GetOptimalKernelWidth1D(radius,sigma);
 865   kernel=GetBlurKernel(width,sigma);
 866   if (kernel == (double *) NULL)
 867     {
 868       blur_image=DestroyImage(blur_image);
 869       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 870     }
 871   if (image->debug != MagickFalse)
 872     {
 873       char
 874         format[MaxTextExtent],
 875         *message;
 876
 877       register const double
 878         *k;
 879
 880       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
 881         "  BlurImage with %.20g kernel:",(double) width);
 882       message=AcquireString("");
 883       k=kernel;
 884       for (i=0; i < (ssize_t) width; i++)
 885       {
 886         *message='\0';
 887         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) i);
 888         (void) ConcatenateString(&message,format);
 889         (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
 890         (void) ConcatenateString(&message,format);
 891         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
 892       }
 893       message=DestroyString(message);
 894     }
 895   /*
 896     Blur rows.
 897   */
 898   status=MagickTrue;
 899   progress=0;
 900   GetMagickPixelPacket(image,&bias);
 901   SetMagickPixelPacketBias(image,&bias);
 902   image_view=AcquireCacheView(image);
 903   blur_view=AcquireCacheView(blur_image);
 904 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 905   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 906 #endif
 907   for (y=0; y < (ssize_t) blur_image->rows; y++)
 908   {
 909     register const IndexPacket
 910       *restrict indexes;
 911
 912     register const PixelPacket
 913       *restrict p;
 914
 915     register IndexPacket
 916       *restrict blur_indexes;
 917
 918     register PixelPacket
 919       *restrict q;
 920
 921     register ssize_t
 922       x;
 923
 924     if (status == MagickFalse)
 925       continue;
 926     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y,
 927       image->columns+width,1,exception);
 928     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
 929       exception);
 930     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 931       {
 932         status=MagickFalse;
 933         continue;
 934       }
 935     indexes=GetCacheViewVirtualIndexQueue(image_view);
 936     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
 937     for (x=0; x < (ssize_t) blur_image->columns; x++)
 938     {
 939       MagickPixelPacket
 940         pixel;
 941
 942       register const double
 943         *restrict k;
 944
 945       register const PixelPacket
 946         *restrict kernel_pixels;
 947
 948       register ssize_t
 949         i;
 950
 951       pixel=bias;
 952       k=kernel;
 953       kernel_pixels=p;
 954       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
 955         {
 956           for (i=0; i < (ssize_t) width; i++)
 957           {
 958             pixel.red+=(*k)*kernel_pixels->red;
 959             pixel.green+=(*k)*kernel_pixels->green;
 960             pixel.blue+=(*k)*kernel_pixels->blue;
 961             k++;
 962             kernel_pixels++;
 963           }
 964           if ((channel & RedChannel) != 0)
 965             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
 966           if ((channel & GreenChannel) != 0)
 967             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
 968           if ((channel & BlueChannel) != 0)
 969             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
 970           if ((channel & OpacityChannel) != 0)
 971             {
 972               k=kernel;
 973               kernel_pixels=p;
 974               for (i=0; i < (ssize_t) width; i++)
 975               {
 976                 pixel.opacity+=(*k)*kernel_pixels->opacity;
 977                 k++;
 978                 kernel_pixels++;
 979               }
 980               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
 981             }
 982           if (((channel & IndexChannel) != 0) &&
 983               (image->colorspace == CMYKColorspace))
 984             {
 985               register const IndexPacket
 986                 *restrict kernel_indexes;
 987
 988               k=kernel;
 989               kernel_indexes=indexes;
 990               for (i=0; i < (ssize_t) width; i++)
 991               {
 992                 pixel.index+=(*k)*(*kernel_indexes);
 993                 k++;
 994                 kernel_indexes++;
 995               }
 996               blur_indexes[x]=ClampToQuantum(pixel.index);
 997             }
 998         }
 999       else
1000         {
1001           MagickRealType
1002             alpha,
1003             gamma;
1004
1005           gamma=0.0;
1006           for (i=0; i < (ssize_t) width; i++)
1007           {
1008             alpha=(MagickRealType) (QuantumScale*
1009               GetAlphaPixelComponent(kernel_pixels));
1010             pixel.red+=(*k)*alpha*kernel_pixels->red;
1011             pixel.green+=(*k)*alpha*kernel_pixels->green;
1012             pixel.blue+=(*k)*alpha*kernel_pixels->blue;
1013             gamma+=(*k)*alpha;
1014             k++;
1015             kernel_pixels++;
1016           }
1017           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1018           if ((channel & RedChannel) != 0)
1019             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
1020           if ((channel & GreenChannel) != 0)
1021             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
1022           if ((channel & BlueChannel) != 0)
1023             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
1024           if ((channel & OpacityChannel) != 0)
1025             {
1026               k=kernel;
1027               kernel_pixels=p;
1028               for (i=0; i < (ssize_t) width; i++)
1029               {
1030                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1031                 k++;
1032                 kernel_pixels++;
1033               }
1034               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1035             }
1036           if (((channel & IndexChannel) != 0) &&
1037               (image->colorspace == CMYKColorspace))
1038             {
1039               register const IndexPacket
1040                 *restrict kernel_indexes;
1041
1042               k=kernel;
1043               kernel_pixels=p;
1044               kernel_indexes=indexes;
1045               for (i=0; i < (ssize_t) width; i++)
1046               {
1047                 alpha=(MagickRealType) (QuantumScale*
1048                   GetAlphaPixelComponent(kernel_pixels));
1049                 pixel.index+=(*k)*alpha*(*kernel_indexes);
1050                 k++;
1051                 kernel_pixels++;
1052                 kernel_indexes++;
1053               }
1054               blur_indexes[x]=ClampToQuantum(gamma*
1055                 GetIndexPixelComponent(&pixel));
1056             }
1057         }
1058       p++;
1059       q++;
1060     }
1061     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
1062       status=MagickFalse;
1063     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1064       {
1065         MagickBooleanType
1066           proceed;
1067
1068 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1069   #pragma omp critical (MagickCore_BlurImageChannel)
1070 #endif
1071         proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
1072           blur_image->columns);
1073         if (proceed == MagickFalse)
1074           status=MagickFalse;
1075       }
1076   }
1077   blur_view=DestroyCacheView(blur_view);
1078   image_view=DestroyCacheView(image_view);
1079   /*
1080     Blur columns.
1081   */
1082   image_view=AcquireCacheView(blur_image);
1083   blur_view=AcquireCacheView(blur_image);
1084 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1085   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1086 #endif
1087   for (x=0; x < (ssize_t) blur_image->columns; x++)
1088   {
1089     register const IndexPacket
1090       *restrict indexes;
1091
1092     register const PixelPacket
1093       *restrict p;
1094
1095     register IndexPacket
1096       *restrict blur_indexes;
1097
1098     register PixelPacket
1099       *restrict q;
1100
1101     register ssize_t
1102       y;
1103
1104     if (status == MagickFalse)
1105       continue;
1106     p=GetCacheViewVirtualPixels(image_view,x,-((ssize_t) width/2L),1,
1107       image->rows+width,exception);
1108     q=GetCacheViewAuthenticPixels(blur_view,x,0,1,blur_image->rows,exception);
1109     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1110       {
1111         status=MagickFalse;
1112         continue;
1113       }
1114     indexes=GetCacheViewVirtualIndexQueue(image_view);
1115     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
1116     for (y=0; y < (ssize_t) blur_image->rows; y++)
1117     {
1118       MagickPixelPacket
1119         pixel;
1120
1121       register const double
1122         *restrict k;
1123
1124       register const PixelPacket
1125         *restrict kernel_pixels;
1126
1127       register ssize_t
1128         i;
1129
1130       pixel=bias;
1131       k=kernel;
1132       kernel_pixels=p;
1133       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
1134         {
1135           for (i=0; i < (ssize_t) width; i++)
1136           {
1137             pixel.red+=(*k)*kernel_pixels->red;
1138             pixel.green+=(*k)*kernel_pixels->green;
1139             pixel.blue+=(*k)*kernel_pixels->blue;
1140             k++;
1141             kernel_pixels++;
1142           }
1143           if ((channel & RedChannel) != 0)
1144             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
1145           if ((channel & GreenChannel) != 0)
1146             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
1147           if ((channel & BlueChannel) != 0)
1148             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
1149           if ((channel & OpacityChannel) != 0)
1150             {
1151               k=kernel;
1152               kernel_pixels=p;
1153               for (i=0; i < (ssize_t) width; i++)
1154               {
1155                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1156                 k++;
1157                 kernel_pixels++;
1158               }
1159               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1160             }
1161           if (((channel & IndexChannel) != 0) &&
1162               (image->colorspace == CMYKColorspace))
1163             {
1164               register const IndexPacket
1165                 *restrict kernel_indexes;
1166
1167               k=kernel;
1168               kernel_indexes=indexes;
1169               for (i=0; i < (ssize_t) width; i++)
1170               {
1171                 pixel.index+=(*k)*(*kernel_indexes);
1172                 k++;
1173                 kernel_indexes++;
1174               }
1175               blur_indexes[y]=ClampToQuantum(pixel.index);
1176             }
1177         }
1178       else
1179         {
1180           MagickRealType
1181             alpha,
1182             gamma;
1183
1184           gamma=0.0;
1185           for (i=0; i < (ssize_t) width; i++)
1186           {
1187             alpha=(MagickRealType) (QuantumScale*
1188               GetAlphaPixelComponent(kernel_pixels));
1189             pixel.red+=(*k)*alpha*kernel_pixels->red;
1190             pixel.green+=(*k)*alpha*kernel_pixels->green;
1191             pixel.blue+=(*k)*alpha*kernel_pixels->blue;
1192             gamma+=(*k)*alpha;
1193             k++;
1194             kernel_pixels++;
1195           }
1196           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1197           if ((channel & RedChannel) != 0)
1198             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
1199           if ((channel & GreenChannel) != 0)
1200             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
1201           if ((channel & BlueChannel) != 0)
1202             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
1203           if ((channel & OpacityChannel) != 0)
1204             {
1205               k=kernel;
1206               kernel_pixels=p;
1207               for (i=0; i < (ssize_t) width; i++)
1208               {
1209                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1210                 k++;
1211                 kernel_pixels++;
1212               }
1213               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1214             }
1215           if (((channel & IndexChannel) != 0) &&
1216               (image->colorspace == CMYKColorspace))
1217             {
1218               register const IndexPacket
1219                 *restrict kernel_indexes;
1220
1221               k=kernel;
1222               kernel_pixels=p;
1223               kernel_indexes=indexes;
1224               for (i=0; i < (ssize_t) width; i++)
1225               {
1226                 alpha=(MagickRealType) (QuantumScale*
1227                   GetAlphaPixelComponent(kernel_pixels));
1228                 pixel.index+=(*k)*alpha*(*kernel_indexes);
1229                 k++;
1230                 kernel_pixels++;
1231                 kernel_indexes++;
1232               }
1233               blur_indexes[y]=ClampToQuantum(gamma*
1234                 GetIndexPixelComponent(&pixel));
1235             }
1236         }
1237       p++;
1238       q++;
1239     }
1240     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
1241       status=MagickFalse;
1242     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1243       {
1244         MagickBooleanType
1245           proceed;
1246
1247 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1248   #pragma omp critical (MagickCore_BlurImageChannel)
1249 #endif
1250         proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
1251           blur_image->columns);
1252         if (proceed == MagickFalse)
1253           status=MagickFalse;
1254       }
1255   }
1256   blur_view=DestroyCacheView(blur_view);
1257   image_view=DestroyCacheView(image_view);
1258   kernel=(double *) RelinquishMagickMemory(kernel);
1259   if (status == MagickFalse)
1260     blur_image=DestroyImage(blur_image);
1261   blur_image->type=image->type;
1262   return(blur_image);
1263 }
1264 \f
1265 /*
1266 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1267 %                                                                             %
1268 %                                                                             %
1269 %                                                                             %
1270 %     C o n v o l v e I m a g e                                               %
1271 %                                                                             %
1272 %                                                                             %
1273 %                                                                             %
1274 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1275 %
1276 %  ConvolveImage() applies a custom convolution kernel to the image.
1277 %
1278 %  The format of the ConvolveImage method is:
1279 %
1280 %      Image *ConvolveImage(const Image *image,const size_t order,
1281 %        const double *kernel,ExceptionInfo *exception)
1282 %      Image *ConvolveImageChannel(const Image *image,const ChannelType channel,
1283 %        const size_t order,const double *kernel,ExceptionInfo *exception)
1284 %
1285 %  A description of each parameter follows:
1286 %
1287 %    o image: the image.
1288 %
1289 %    o channel: the channel type.
1290 %
1291 %    o order: the number of columns and rows in the filter kernel.
1292 %
1293 %    o kernel: An array of double representing the convolution kernel.
1294 %
1295 %    o exception: return any errors or warnings in this structure.
1296 %
1297 */
1298
1299 MagickExport Image *ConvolveImage(const Image *image,const size_t order,
1300   const double *kernel,ExceptionInfo *exception)
1301 {
1302   Image
1303     *convolve_image;
1304
1305   convolve_image=ConvolveImageChannel(image,DefaultChannels,order,kernel,
1306     exception);
1307   return(convolve_image);
1308 }
1309
1310 MagickExport Image *ConvolveImageChannel(const Image *image,
1311   const ChannelType channel,const size_t order,const double *kernel,
1312   ExceptionInfo *exception)
1313 {
1314 #define ConvolveImageTag  "Convolve/Image"
1315
1316   CacheView
1317     *convolve_view,
1318     *image_view;
1319
1320   double
1321     *normal_kernel;
1322
1323   Image
1324     *convolve_image;
1325
1326   MagickBooleanType
1327     status;
1328
1329   MagickOffsetType
1330     progress;
1331
1332   MagickPixelPacket
1333     bias;
1334
1335   MagickRealType
1336     gamma;
1337
1338   register ssize_t
1339     i;
1340
1341   size_t
1342     width;
1343
1344   ssize_t
1345     y;
1346
1347   /*
1348     Initialize convolve image attributes.
1349   */
1350   assert(image != (Image *) NULL);
1351   assert(image->signature == MagickSignature);
1352   if (image->debug != MagickFalse)
1353     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1354   assert(exception != (ExceptionInfo *) NULL);
1355   assert(exception->signature == MagickSignature);
1356   width=order;
1357   if ((width % 2) == 0)
1358     ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
1359   convolve_image=CloneImage(image,0,0,MagickTrue,exception);
1360   if (convolve_image == (Image *) NULL)
1361     return((Image *) NULL);
1362   if (SetImageStorageClass(convolve_image,DirectClass) == MagickFalse)
1363     {
1364       InheritException(exception,&convolve_image->exception);
1365       convolve_image=DestroyImage(convolve_image);
1366       return((Image *) NULL);
1367     }
1368   if (image->debug != MagickFalse)
1369     {
1370       char
1371         format[MaxTextExtent],
1372         *message;
1373
1374       register const double
1375         *k;
1376
1377       ssize_t
1378         u,
1379         v;
1380
1381       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1382         "  ConvolveImage with %.20gx%.20g kernel:",(double) width,(double)
1383         width);
1384       message=AcquireString("");
1385       k=kernel;
1386       for (v=0; v < (ssize_t) width; v++)
1387       {
1388         *message='\0';
1389         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) v);
1390         (void) ConcatenateString(&message,format);
1391         for (u=0; u < (ssize_t) width; u++)
1392         {
1393           (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
1394           (void) ConcatenateString(&message,format);
1395         }
1396         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
1397       }
1398       message=DestroyString(message);
1399     }
1400   /*
1401     Normalize kernel.
1402   */
1403   normal_kernel=(double *) AcquireQuantumMemory(width*width,
1404     sizeof(*normal_kernel));
1405   if (normal_kernel == (double *) NULL)
1406     {
1407       convolve_image=DestroyImage(convolve_image);
1408       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1409     }
1410   gamma=0.0;
1411   for (i=0; i < (ssize_t) (width*width); i++)
1412     gamma+=kernel[i];
1413   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1414   for (i=0; i < (ssize_t) (width*width); i++)
1415     normal_kernel[i]=gamma*kernel[i];
1416   /*
1417     Convolve image.
1418   */
1419   status=MagickTrue;
1420   progress=0;
1421   GetMagickPixelPacket(image,&bias);
1422   SetMagickPixelPacketBias(image,&bias);
1423   image_view=AcquireCacheView(image);
1424   convolve_view=AcquireCacheView(convolve_image);
1425 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1426   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1427 #endif
1428   for (y=0; y < (ssize_t) image->rows; y++)
1429   {
1430     MagickBooleanType
1431       sync;
1432
1433     register const IndexPacket
1434       *restrict indexes;
1435
1436     register const PixelPacket
1437       *restrict p;
1438
1439     register IndexPacket
1440       *restrict convolve_indexes;
1441
1442     register PixelPacket
1443       *restrict q;
1444
1445     register ssize_t
1446       x;
1447
1448     if (status == MagickFalse)
1449       continue;
1450     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
1451       (width/2L),image->columns+width,width,exception);
1452     q=GetCacheViewAuthenticPixels(convolve_view,0,y,convolve_image->columns,1,
1453       exception);
1454     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1455       {
1456         status=MagickFalse;
1457         continue;
1458       }
1459     indexes=GetCacheViewVirtualIndexQueue(image_view);
1460     convolve_indexes=GetCacheViewAuthenticIndexQueue(convolve_view);
1461     for (x=0; x < (ssize_t) image->columns; x++)
1462     {
1463       MagickPixelPacket
1464         pixel;
1465
1466       register const double
1467         *restrict k;
1468
1469       register const PixelPacket
1470         *restrict kernel_pixels;
1471
1472       register ssize_t
1473         u;
1474
1475       ssize_t
1476         v;
1477
1478       pixel=bias;
1479       k=normal_kernel;
1480       kernel_pixels=p;
1481       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
1482         {
1483           for (v=0; v < (ssize_t) width; v++)
1484           {
1485             for (u=0; u < (ssize_t) width; u++)
1486             {
1487               pixel.red+=(*k)*kernel_pixels[u].red;
1488               pixel.green+=(*k)*kernel_pixels[u].green;
1489               pixel.blue+=(*k)*kernel_pixels[u].blue;
1490               k++;
1491             }
1492             kernel_pixels+=image->columns+width;
1493           }
1494           if ((channel & RedChannel) != 0)
1495             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
1496           if ((channel & GreenChannel) != 0)
1497             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
1498           if ((channel & BlueChannel) != 0)
1499             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
1500           if ((channel & OpacityChannel) != 0)
1501             {
1502               k=normal_kernel;
1503               kernel_pixels=p;
1504               for (v=0; v < (ssize_t) width; v++)
1505               {
1506                 for (u=0; u < (ssize_t) width; u++)
1507                 {
1508                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
1509                   k++;
1510                 }
1511                 kernel_pixels+=image->columns+width;
1512               }
1513               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1514             }
1515           if (((channel & IndexChannel) != 0) &&
1516               (image->colorspace == CMYKColorspace))
1517             {
1518               register const IndexPacket
1519                 *restrict kernel_indexes;
1520
1521               k=normal_kernel;
1522               kernel_indexes=indexes;
1523               for (v=0; v < (ssize_t) width; v++)
1524               {
1525                 for (u=0; u < (ssize_t) width; u++)
1526                 {
1527                   pixel.index+=(*k)*kernel_indexes[u];
1528                   k++;
1529                 }
1530                 kernel_indexes+=image->columns+width;
1531               }
1532               convolve_indexes[x]=ClampToQuantum(pixel.index);
1533             }
1534         }
1535       else
1536         {
1537           MagickRealType
1538             alpha,
1539             gamma;
1540
1541           gamma=0.0;
1542           for (v=0; v < (ssize_t) width; v++)
1543           {
1544             for (u=0; u < (ssize_t) width; u++)
1545             {
1546               alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1547                 kernel_pixels[u].opacity));
1548               pixel.red+=(*k)*alpha*kernel_pixels[u].red;
1549               pixel.green+=(*k)*alpha*kernel_pixels[u].green;
1550               pixel.blue+=(*k)*alpha*kernel_pixels[u].blue;
1551               gamma+=(*k)*alpha;
1552               k++;
1553             }
1554             kernel_pixels+=image->columns+width;
1555           }
1556           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1557           if ((channel & RedChannel) != 0)
1558             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
1559           if ((channel & GreenChannel) != 0)
1560             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
1561           if ((channel & BlueChannel) != 0)
1562             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
1563           if ((channel & OpacityChannel) != 0)
1564             {
1565               k=normal_kernel;
1566               kernel_pixels=p;
1567               for (v=0; v < (ssize_t) width; v++)
1568               {
1569                 for (u=0; u < (ssize_t) width; u++)
1570                 {
1571                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
1572                   k++;
1573                 }
1574                 kernel_pixels+=image->columns+width;
1575               }
1576               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1577             }
1578           if (((channel & IndexChannel) != 0) &&
1579               (image->colorspace == CMYKColorspace))
1580             {
1581               register const IndexPacket
1582                 *restrict kernel_indexes;
1583
1584               k=normal_kernel;
1585               kernel_pixels=p;
1586               kernel_indexes=indexes;
1587               for (v=0; v < (ssize_t) width; v++)
1588               {
1589                 for (u=0; u < (ssize_t) width; u++)
1590                 {
1591                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1592                     kernel_pixels[u].opacity));
1593                   pixel.index+=(*k)*alpha*kernel_indexes[u];
1594                   k++;
1595                 }
1596                 kernel_pixels+=image->columns+width;
1597                 kernel_indexes+=image->columns+width;
1598               }
1599               convolve_indexes[x]=ClampToQuantum(gamma*
1600                 GetIndexPixelComponent(&pixel));
1601             }
1602         }
1603       p++;
1604       q++;
1605     }
1606     sync=SyncCacheViewAuthenticPixels(convolve_view,exception);
1607     if (sync == MagickFalse)
1608       status=MagickFalse;
1609     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1610       {
1611         MagickBooleanType
1612           proceed;
1613
1614 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1615   #pragma omp critical (MagickCore_ConvolveImageChannel)
1616 #endif
1617         proceed=SetImageProgress(image,ConvolveImageTag,progress++,image->rows);
1618         if (proceed == MagickFalse)
1619           status=MagickFalse;
1620       }
1621   }
1622   convolve_image->type=image->type;
1623   convolve_view=DestroyCacheView(convolve_view);
1624   image_view=DestroyCacheView(image_view);
1625   normal_kernel=(double *) RelinquishMagickMemory(normal_kernel);
1626   if (status == MagickFalse)
1627     convolve_image=DestroyImage(convolve_image);
1628   return(convolve_image);
1629 }
1630 \f
1631 /*
1632 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1633 %                                                                             %
1634 %                                                                             %
1635 %                                                                             %
1636 %     D e s p e c k l e I m a g e                                             %
1637 %                                                                             %
1638 %                                                                             %
1639 %                                                                             %
1640 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1641 %
1642 %  DespeckleImage() reduces the speckle noise in an image while perserving the
1643 %  edges of the original image.
1644 %
1645 %  The format of the DespeckleImage method is:
1646 %
1647 %      Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
1648 %
1649 %  A description of each parameter follows:
1650 %
1651 %    o image: the image.
1652 %
1653 %    o exception: return any errors or warnings in this structure.
1654 %
1655 */
1656
1657 static Quantum **DestroyPixelThreadSet(Quantum **pixels)
1658 {
1659   register ssize_t
1660     i;
1661
1662   assert(pixels != (Quantum **) NULL);
1663   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
1664     if (pixels[i] != (Quantum *) NULL)
1665       pixels[i]=(Quantum *) RelinquishMagickMemory(pixels[i]);
1666   pixels=(Quantum **) RelinquishMagickMemory(pixels);
1667   return(pixels);
1668 }
1669
1670 static Quantum **AcquirePixelThreadSet(const size_t count)
1671 {
1672   Quantum
1673     **pixels;
1674
1675   register ssize_t
1676     i;
1677
1678   size_t
1679     number_threads;
1680
1681   number_threads=GetOpenMPMaximumThreads();
1682   pixels=(Quantum **) AcquireQuantumMemory(number_threads,sizeof(*pixels));
1683   if (pixels == (Quantum **) NULL)
1684     return((Quantum **) NULL);
1685   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
1686   for (i=0; i < (ssize_t) number_threads; i++)
1687   {
1688     pixels[i]=(Quantum *) AcquireQuantumMemory(count,sizeof(**pixels));
1689     if (pixels[i] == (Quantum *) NULL)
1690       return(DestroyPixelThreadSet(pixels));
1691   }
1692   return(pixels);
1693 }
1694
1695 static void Hull(const ssize_t x_offset,const ssize_t y_offset,
1696   const size_t columns,const size_t rows,Quantum *f,Quantum *g,
1697   const int polarity)
1698 {
1699   MagickRealType
1700     v;
1701
1702   register Quantum
1703     *p,
1704     *q,
1705     *r,
1706     *s;
1707
1708   register ssize_t
1709     x;
1710
1711   ssize_t
1712     y;
1713
1714   assert(f != (Quantum *) NULL);
1715   assert(g != (Quantum *) NULL);
1716   p=f+(columns+2);
1717   q=g+(columns+2);
1718   r=p+(y_offset*((ssize_t) columns+2)+x_offset);
1719   for (y=0; y < (ssize_t) rows; y++)
1720   {
1721     p++;
1722     q++;
1723     r++;
1724     if (polarity > 0)
1725       for (x=(ssize_t) columns; x != 0; x--)
1726       {
1727         v=(MagickRealType) (*p);
1728         if ((MagickRealType) *r >= (v+(MagickRealType) ScaleCharToQuantum(2)))
1729           v+=ScaleCharToQuantum(1);
1730         *q=(Quantum) v;
1731         p++;
1732         q++;
1733         r++;
1734       }
1735     else
1736       for (x=(ssize_t) columns; x != 0; x--)
1737       {
1738         v=(MagickRealType) (*p);
1739         if ((MagickRealType) *r <= (v-(MagickRealType) ScaleCharToQuantum(2)))
1740           v-=(ssize_t) ScaleCharToQuantum(1);
1741         *q=(Quantum) v;
1742         p++;
1743         q++;
1744         r++;
1745       }
1746     p++;
1747     q++;
1748     r++;
1749   }
1750   p=f+(columns+2);
1751   q=g+(columns+2);
1752   r=q+(y_offset*((ssize_t) columns+2)+x_offset);
1753   s=q-(y_offset*((ssize_t) columns+2)+x_offset);
1754   for (y=0; y < (ssize_t) rows; y++)
1755   {
1756     p++;
1757     q++;
1758     r++;
1759     s++;
1760     if (polarity > 0)
1761       for (x=(ssize_t) columns; x != 0; x--)
1762       {
1763         v=(MagickRealType) (*q);
1764         if (((MagickRealType) *s >=
1765              (v+(MagickRealType) ScaleCharToQuantum(2))) &&
1766             ((MagickRealType) *r > v))
1767           v+=ScaleCharToQuantum(1);
1768         *p=(Quantum) v;
1769         p++;
1770         q++;
1771         r++;
1772         s++;
1773       }
1774     else
1775       for (x=(ssize_t) columns; x != 0; x--)
1776       {
1777         v=(MagickRealType) (*q);
1778         if (((MagickRealType) *s <=
1779              (v-(MagickRealType) ScaleCharToQuantum(2))) &&
1780             ((MagickRealType) *r < v))
1781           v-=(MagickRealType) ScaleCharToQuantum(1);
1782         *p=(Quantum) v;
1783         p++;
1784         q++;
1785         r++;
1786         s++;
1787       }
1788     p++;
1789     q++;
1790     r++;
1791     s++;
1792   }
1793 }
1794
1795 MagickExport Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
1796 {
1797 #define DespeckleImageTag  "Despeckle/Image"
1798
1799   CacheView
1800     *despeckle_view,
1801     *image_view;
1802
1803   Image
1804     *despeckle_image;
1805
1806   MagickBooleanType
1807     status;
1808
1809   Quantum
1810     **restrict buffers,
1811     **restrict pixels;
1812
1813   size_t
1814     length;
1815
1816   ssize_t
1817     channel;
1818
1819   static const ssize_t
1820     X[4] = {0, 1, 1,-1},
1821     Y[4] = {1, 0, 1, 1};
1822
1823   /*
1824     Allocate despeckled image.
1825   */
1826   assert(image != (const Image *) NULL);
1827   assert(image->signature == MagickSignature);
1828   if (image->debug != MagickFalse)
1829     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1830   assert(exception != (ExceptionInfo *) NULL);
1831   assert(exception->signature == MagickSignature);
1832   despeckle_image=CloneImage(image,image->columns,image->rows,MagickTrue,
1833     exception);
1834   if (despeckle_image == (Image *) NULL)
1835     return((Image *) NULL);
1836   if (SetImageStorageClass(despeckle_image,DirectClass) == MagickFalse)
1837     {
1838       InheritException(exception,&despeckle_image->exception);
1839       despeckle_image=DestroyImage(despeckle_image);
1840       return((Image *) NULL);
1841     }
1842   /*
1843     Allocate image buffers.
1844   */
1845   length=(size_t) ((image->columns+2)*(image->rows+2));
1846   pixels=AcquirePixelThreadSet(length);
1847   buffers=AcquirePixelThreadSet(length);
1848   if ((pixels == (Quantum **) NULL) || (buffers == (Quantum **) NULL))
1849     {
1850       if (buffers != (Quantum **) NULL)
1851         buffers=DestroyPixelThreadSet(buffers);
1852       if (pixels != (Quantum **) NULL)
1853         pixels=DestroyPixelThreadSet(pixels);
1854       despeckle_image=DestroyImage(despeckle_image);
1855       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1856     }
1857   /*
1858     Reduce speckle in the image.
1859   */
1860   status=MagickTrue;
1861   image_view=AcquireCacheView(image);
1862   despeckle_view=AcquireCacheView(despeckle_image);
1863 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1864   #pragma omp parallel for schedule(dynamic,1) shared(status)
1865 #endif
1866   for (channel=0; channel <= 3; channel++)
1867   {
1868     const int
1869       id = GetOpenMPThreadId();
1870
1871     register ssize_t
1872       i,
1873       x;
1874
1875     register Quantum
1876       *buffer,
1877       *pixel;
1878
1879     ssize_t
1880       j,
1881       y;
1882
1883     if (status == MagickFalse)
1884       continue;
1885     pixel=pixels[id];
1886     (void) ResetMagickMemory(pixel,0,length*sizeof(*pixel));
1887     buffer=buffers[id];
1888     j=(ssize_t) image->columns+2;
1889     for (y=0; y < (ssize_t) image->rows; y++)
1890     {
1891       register const PixelPacket
1892         *restrict p;
1893
1894       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1895       if (p == (const PixelPacket *) NULL)
1896         break;
1897       j++;
1898       for (x=0; x < (ssize_t) image->columns; x++)
1899       {
1900         switch (channel)
1901         {
1902           case 0: pixel[j]=GetRedPixelComponent(p); break;
1903           case 1: pixel[j]=GetGreenPixelComponent(p); break;
1904           case 2: pixel[j]=GetBluePixelComponent(p); break;
1905           case 3: pixel[j]=GetOpacityPixelComponent(p); break;
1906           default: break;
1907         }
1908         p++;
1909         j++;
1910       }
1911       j++;
1912     }
1913     (void) ResetMagickMemory(buffer,0,length*sizeof(*buffer));
1914     for (i=0; i < 4; i++)
1915     {
1916       Hull(X[i],Y[i],image->columns,image->rows,pixel,buffer,1);
1917       Hull(-X[i],-Y[i],image->columns,image->rows,pixel,buffer,1);
1918       Hull(-X[i],-Y[i],image->columns,image->rows,pixel,buffer,-1);
1919       Hull(X[i],Y[i],image->columns,image->rows,pixel,buffer,-1);
1920     }
1921     j=(ssize_t) image->columns+2;
1922     for (y=0; y < (ssize_t) image->rows; y++)
1923     {
1924       MagickBooleanType
1925         sync;
1926
1927       register PixelPacket
1928         *restrict q;
1929
1930       q=GetCacheViewAuthenticPixels(despeckle_view,0,y,despeckle_image->columns,
1931         1,exception);
1932       if (q == (PixelPacket *) NULL)
1933         break;
1934       j++;
1935       for (x=0; x < (ssize_t) image->columns; x++)
1936       {
1937         switch (channel)
1938         {
1939           case 0: q->red=pixel[j]; break;
1940           case 1: q->green=pixel[j]; break;
1941           case 2: q->blue=pixel[j]; break;
1942           case 3: q->opacity=pixel[j]; break;
1943           default: break;
1944         }
1945         q++;
1946         j++;
1947       }
1948       sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
1949       if (sync == MagickFalse)
1950         {
1951           status=MagickFalse;
1952           break;
1953         }
1954       j++;
1955     }
1956     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1957       {
1958         MagickBooleanType
1959           proceed;
1960
1961 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1962   #pragma omp critical (MagickCore_DespeckleImage)
1963 #endif
1964         proceed=SetImageProgress(image,DespeckleImageTag,(MagickOffsetType)
1965           channel,3);
1966         if (proceed == MagickFalse)
1967           status=MagickFalse;
1968       }
1969   }
1970   despeckle_view=DestroyCacheView(despeckle_view);
1971   image_view=DestroyCacheView(image_view);
1972   buffers=DestroyPixelThreadSet(buffers);
1973   pixels=DestroyPixelThreadSet(pixels);
1974   despeckle_image->type=image->type;
1975   if (status == MagickFalse)
1976     despeckle_image=DestroyImage(despeckle_image);
1977   return(despeckle_image);
1978 }
1979 \f
1980 /*
1981 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1982 %                                                                             %
1983 %                                                                             %
1984 %                                                                             %
1985 %     E d g e I m a g e                                                       %
1986 %                                                                             %
1987 %                                                                             %
1988 %                                                                             %
1989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1990 %
1991 %  EdgeImage() finds edges in an image.  Radius defines the radius of the
1992 %  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
1993 %  radius for you.
1994 %
1995 %  The format of the EdgeImage method is:
1996 %
1997 %      Image *EdgeImage(const Image *image,const double radius,
1998 %        ExceptionInfo *exception)
1999 %
2000 %  A description of each parameter follows:
2001 %
2002 %    o image: the image.
2003 %
2004 %    o radius: the radius of the pixel neighborhood.
2005 %
2006 %    o exception: return any errors or warnings in this structure.
2007 %
2008 */
2009 MagickExport Image *EdgeImage(const Image *image,const double radius,
2010   ExceptionInfo *exception)
2011 {
2012   Image
2013     *edge_image;
2014
2015   double
2016     *kernel;
2017
2018   register ssize_t
2019     i;
2020
2021   size_t
2022     width;
2023
2024   assert(image != (const Image *) NULL);
2025   assert(image->signature == MagickSignature);
2026   if (image->debug != MagickFalse)
2027     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2028   assert(exception != (ExceptionInfo *) NULL);
2029   assert(exception->signature == MagickSignature);
2030   width=GetOptimalKernelWidth1D(radius,0.5);
2031   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2032   if (kernel == (double *) NULL)
2033     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2034   for (i=0; i < (ssize_t) (width*width); i++)
2035     kernel[i]=(-1.0);
2036   kernel[i/2]=(double) (width*width-1.0);
2037   edge_image=ConvolveImage(image,width,kernel,exception);
2038   kernel=(double *) RelinquishMagickMemory(kernel);
2039   return(edge_image);
2040 }
2041 \f
2042 /*
2043 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2044 %                                                                             %
2045 %                                                                             %
2046 %                                                                             %
2047 %     E m b o s s I m a g e                                                   %
2048 %                                                                             %
2049 %                                                                             %
2050 %                                                                             %
2051 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2052 %
2053 %  EmbossImage() returns a grayscale image with a three-dimensional effect.
2054 %  We convolve the image with a Gaussian operator of the given radius and
2055 %  standard deviation (sigma).  For reasonable results, radius should be
2056 %  larger than sigma.  Use a radius of 0 and Emboss() selects a suitable
2057 %  radius for you.
2058 %
2059 %  The format of the EmbossImage method is:
2060 %
2061 %      Image *EmbossImage(const Image *image,const double radius,
2062 %        const double sigma,ExceptionInfo *exception)
2063 %
2064 %  A description of each parameter follows:
2065 %
2066 %    o image: the image.
2067 %
2068 %    o radius: the radius of the pixel neighborhood.
2069 %
2070 %    o sigma: the standard deviation of the Gaussian, in pixels.
2071 %
2072 %    o exception: return any errors or warnings in this structure.
2073 %
2074 */
2075 MagickExport Image *EmbossImage(const Image *image,const double radius,
2076   const double sigma,ExceptionInfo *exception)
2077 {
2078   double
2079     *kernel;
2080
2081   Image
2082     *emboss_image;
2083
2084   register ssize_t
2085     i;
2086
2087   size_t
2088     width;
2089
2090   ssize_t
2091     j,
2092     k,
2093     u,
2094     v;
2095
2096   assert(image != (Image *) NULL);
2097   assert(image->signature == MagickSignature);
2098   if (image->debug != MagickFalse)
2099     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2100   assert(exception != (ExceptionInfo *) NULL);
2101   assert(exception->signature == MagickSignature);
2102   width=GetOptimalKernelWidth2D(radius,sigma);
2103   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2104   if (kernel == (double *) NULL)
2105     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2106   j=(ssize_t) width/2;
2107   k=j;
2108   i=0;
2109   for (v=(-j); v <= j; v++)
2110   {
2111     for (u=(-j); u <= j; u++)
2112     {
2113       kernel[i]=(double) (((u < 0) || (v < 0) ? -8.0 : 8.0)*
2114         exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
2115         (2.0*MagickPI*MagickSigma*MagickSigma));
2116       if (u != k)
2117         kernel[i]=0.0;
2118       i++;
2119     }
2120     k--;
2121   }
2122   emboss_image=ConvolveImage(image,width,kernel,exception);
2123   if (emboss_image != (Image *) NULL)
2124     (void) EqualizeImage(emboss_image);
2125   kernel=(double *) RelinquishMagickMemory(kernel);
2126   return(emboss_image);
2127 }
2128 \f
2129 /*
2130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2131 %                                                                             %
2132 %                                                                             %
2133 %                                                                             %
2134 %     F i l t e r I m a g e                                                   %
2135 %                                                                             %
2136 %                                                                             %
2137 %                                                                             %
2138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2139 %
2140 %  FilterImage() applies a custom convolution kernel to the image.
2141 %
2142 %  The format of the FilterImage method is:
2143 %
2144 %      Image *FilterImage(const Image *image,const KernelInfo *kernel,
2145 %        ExceptionInfo *exception)
2146 %      Image *FilterImageChannel(const Image *image,const ChannelType channel,
2147 %        const KernelInfo *kernel,ExceptionInfo *exception)
2148 %
2149 %  A description of each parameter follows:
2150 %
2151 %    o image: the image.
2152 %
2153 %    o channel: the channel type.
2154 %
2155 %    o kernel: the filtering kernel.
2156 %
2157 %    o exception: return any errors or warnings in this structure.
2158 %
2159 */
2160
2161 MagickExport Image *FilterImage(const Image *image,const KernelInfo *kernel,
2162   ExceptionInfo *exception)
2163 {
2164   Image
2165     *filter_image;
2166
2167   filter_image=FilterImageChannel(image,DefaultChannels,kernel,exception);
2168   return(filter_image);
2169 }
2170
2171 MagickExport Image *FilterImageChannel(const Image *image,
2172   const ChannelType channel,const KernelInfo *kernel,ExceptionInfo *exception)
2173 {
2174 #define FilterImageTag  "Filter/Image"
2175
2176   CacheView
2177     *filter_view,
2178     *image_view;
2179
2180   Image
2181     *filter_image;
2182
2183   MagickBooleanType
2184     status;
2185
2186   MagickOffsetType
2187     progress;
2188
2189   MagickPixelPacket
2190     bias;
2191
2192   ssize_t
2193     y;
2194
2195   /*
2196     Initialize filter image attributes.
2197   */
2198   assert(image != (Image *) NULL);
2199   assert(image->signature == MagickSignature);
2200   if (image->debug != MagickFalse)
2201     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2202   assert(exception != (ExceptionInfo *) NULL);
2203   assert(exception->signature == MagickSignature);
2204   if ((kernel->width % 2) == 0)
2205     ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
2206   filter_image=CloneImage(image,0,0,MagickTrue,exception);
2207   if (filter_image == (Image *) NULL)
2208     return((Image *) NULL);
2209   if (SetImageStorageClass(filter_image,DirectClass) == MagickFalse)
2210     {
2211       InheritException(exception,&filter_image->exception);
2212       filter_image=DestroyImage(filter_image);
2213       return((Image *) NULL);
2214     }
2215   if (image->debug != MagickFalse)
2216     {
2217       char
2218         format[MaxTextExtent],
2219         *message;
2220
2221       register const double
2222         *k;
2223
2224       ssize_t
2225         u,
2226         v;
2227
2228       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
2229         "  FilterImage with %.20gx%.20g kernel:",(double) kernel->width,(double)
2230         kernel->height);
2231       message=AcquireString("");
2232       k=kernel->values;
2233       for (v=0; v < (ssize_t) kernel->height; v++)
2234       {
2235         *message='\0';
2236         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) v);
2237         (void) ConcatenateString(&message,format);
2238         for (u=0; u < (ssize_t) kernel->width; u++)
2239         {
2240           (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
2241           (void) ConcatenateString(&message,format);
2242         }
2243         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
2244       }
2245       message=DestroyString(message);
2246     }
2247   status=AccelerateConvolveImage(image,kernel,filter_image,exception);
2248   if (status == MagickTrue)
2249     return(filter_image);
2250   /*
2251     Filter image.
2252   */
2253   status=MagickTrue;
2254   progress=0;
2255   GetMagickPixelPacket(image,&bias);
2256   SetMagickPixelPacketBias(image,&bias);
2257   image_view=AcquireCacheView(image);
2258   filter_view=AcquireCacheView(filter_image);
2259 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2260   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2261 #endif
2262   for (y=0; y < (ssize_t) image->rows; y++)
2263   {
2264     MagickBooleanType
2265       sync;
2266
2267     register const IndexPacket
2268       *restrict indexes;
2269
2270     register const PixelPacket
2271       *restrict p;
2272
2273     register IndexPacket
2274       *restrict filter_indexes;
2275
2276     register PixelPacket
2277       *restrict q;
2278
2279     register ssize_t
2280       x;
2281
2282     if (status == MagickFalse)
2283       continue;
2284     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel->width/2L),
2285       y-(ssize_t) (kernel->height/2L),image->columns+kernel->width,
2286       kernel->height,exception);
2287     q=GetCacheViewAuthenticPixels(filter_view,0,y,filter_image->columns,1,
2288       exception);
2289     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2290       {
2291         status=MagickFalse;
2292         continue;
2293       }
2294     indexes=GetCacheViewVirtualIndexQueue(image_view);
2295     filter_indexes=GetCacheViewAuthenticIndexQueue(filter_view);
2296     for (x=0; x < (ssize_t) image->columns; x++)
2297     {
2298       MagickPixelPacket
2299         pixel;
2300
2301       register const double
2302         *restrict k;
2303
2304       register const PixelPacket
2305         *restrict kernel_pixels;
2306
2307       register ssize_t
2308         u;
2309
2310       ssize_t
2311         v;
2312
2313       pixel=bias;
2314       k=kernel->values;
2315       kernel_pixels=p;
2316       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
2317         {
2318           for (v=0; v < (ssize_t) kernel->width; v++)
2319           {
2320             for (u=0; u < (ssize_t) kernel->height; u++)
2321             {
2322               pixel.red+=(*k)*kernel_pixels[u].red;
2323               pixel.green+=(*k)*kernel_pixels[u].green;
2324               pixel.blue+=(*k)*kernel_pixels[u].blue;
2325               k++;
2326             }
2327             kernel_pixels+=image->columns+kernel->width;
2328           }
2329           if ((channel & RedChannel) != 0)
2330             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
2331           if ((channel & GreenChannel) != 0)
2332             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
2333           if ((channel & BlueChannel) != 0)
2334             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
2335           if ((channel & OpacityChannel) != 0)
2336             {
2337               k=kernel->values;
2338               kernel_pixels=p;
2339               for (v=0; v < (ssize_t) kernel->width; v++)
2340               {
2341                 for (u=0; u < (ssize_t) kernel->height; u++)
2342                 {
2343                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
2344                   k++;
2345                 }
2346                 kernel_pixels+=image->columns+kernel->width;
2347               }
2348               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2349             }
2350           if (((channel & IndexChannel) != 0) &&
2351               (image->colorspace == CMYKColorspace))
2352             {
2353               register const IndexPacket
2354                 *restrict kernel_indexes;
2355
2356               k=kernel->values;
2357               kernel_indexes=indexes;
2358               for (v=0; v < (ssize_t) kernel->width; v++)
2359               {
2360                 for (u=0; u < (ssize_t) kernel->height; u++)
2361                 {
2362                   pixel.index+=(*k)*kernel_indexes[u];
2363                   k++;
2364                 }
2365                 kernel_indexes+=image->columns+kernel->width;
2366               }
2367               filter_indexes[x]=ClampToQuantum(pixel.index);
2368             }
2369         }
2370       else
2371         {
2372           MagickRealType
2373             alpha,
2374             gamma;
2375
2376           gamma=0.0;
2377           for (v=0; v < (ssize_t) kernel->width; v++)
2378           {
2379             for (u=0; u < (ssize_t) kernel->height; u++)
2380             {
2381               alpha=(MagickRealType) (QuantumScale*(QuantumRange-
2382                 kernel_pixels[u].opacity));
2383               pixel.red+=(*k)*alpha*kernel_pixels[u].red;
2384               pixel.green+=(*k)*alpha*kernel_pixels[u].green;
2385               pixel.blue+=(*k)*alpha*kernel_pixels[u].blue;
2386               gamma+=(*k)*alpha;
2387               k++;
2388             }
2389             kernel_pixels+=image->columns+kernel->width;
2390           }
2391           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2392           if ((channel & RedChannel) != 0)
2393             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
2394           if ((channel & GreenChannel) != 0)
2395             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
2396           if ((channel & BlueChannel) != 0)
2397             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
2398           if ((channel & OpacityChannel) != 0)
2399             {
2400               k=kernel->values;
2401               kernel_pixels=p;
2402               for (v=0; v < (ssize_t) kernel->width; v++)
2403               {
2404                 for (u=0; u < (ssize_t) kernel->height; u++)
2405                 {
2406                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
2407                   k++;
2408                 }
2409                 kernel_pixels+=image->columns+kernel->width;
2410               }
2411               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2412             }
2413           if (((channel & IndexChannel) != 0) &&
2414               (image->colorspace == CMYKColorspace))
2415             {
2416               register const IndexPacket
2417                 *restrict kernel_indexes;
2418
2419               k=kernel->values;
2420               kernel_pixels=p;
2421               kernel_indexes=indexes;
2422               for (v=0; v < (ssize_t) kernel->width; v++)
2423               {
2424                 for (u=0; u < (ssize_t) kernel->height; u++)
2425                 {
2426                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
2427                     kernel_pixels[u].opacity));
2428                   pixel.index+=(*k)*alpha*kernel_indexes[u];
2429                   k++;
2430                 }
2431                 kernel_pixels+=image->columns+kernel->width;
2432                 kernel_indexes+=image->columns+kernel->width;
2433               }
2434               filter_indexes[x]=ClampToQuantum(gamma*
2435                 GetIndexPixelComponent(&pixel));
2436             }
2437         }
2438       p++;
2439       q++;
2440     }
2441     sync=SyncCacheViewAuthenticPixels(filter_view,exception);
2442     if (sync == MagickFalse)
2443       status=MagickFalse;
2444     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2445       {
2446         MagickBooleanType
2447           proceed;
2448
2449 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2450   #pragma omp critical (MagickCore_FilterImageChannel)
2451 #endif
2452         proceed=SetImageProgress(image,FilterImageTag,progress++,image->rows);
2453         if (proceed == MagickFalse)
2454           status=MagickFalse;
2455       }
2456   }
2457   filter_image->type=image->type;
2458   filter_view=DestroyCacheView(filter_view);
2459   image_view=DestroyCacheView(image_view);
2460   if (status == MagickFalse)
2461     filter_image=DestroyImage(filter_image);
2462   return(filter_image);
2463 }
2464 \f
2465 /*
2466 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2467 %                                                                             %
2468 %                                                                             %
2469 %                                                                             %
2470 %     G a u s s i a n B l u r I m a g e                                       %
2471 %                                                                             %
2472 %                                                                             %
2473 %                                                                             %
2474 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2475 %
2476 %  GaussianBlurImage() blurs an image.  We convolve the image with a
2477 %  Gaussian operator of the given radius and standard deviation (sigma).
2478 %  For reasonable results, the radius should be larger than sigma.  Use a
2479 %  radius of 0 and GaussianBlurImage() selects a suitable radius for you
2480 %
2481 %  The format of the GaussianBlurImage method is:
2482 %
2483 %      Image *GaussianBlurImage(const Image *image,onst double radius,
2484 %        const double sigma,ExceptionInfo *exception)
2485 %      Image *GaussianBlurImageChannel(const Image *image,
2486 %        const ChannelType channel,const double radius,const double sigma,
2487 %        ExceptionInfo *exception)
2488 %
2489 %  A description of each parameter follows:
2490 %
2491 %    o image: the image.
2492 %
2493 %    o channel: the channel type.
2494 %
2495 %    o radius: the radius of the Gaussian, in pixels, not counting the center
2496 %      pixel.
2497 %
2498 %    o sigma: the standard deviation of the Gaussian, in pixels.
2499 %
2500 %    o exception: return any errors or warnings in this structure.
2501 %
2502 */
2503
2504 MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
2505   const double sigma,ExceptionInfo *exception)
2506 {
2507   Image
2508     *blur_image;
2509
2510   blur_image=GaussianBlurImageChannel(image,DefaultChannels,radius,sigma,
2511     exception);
2512   return(blur_image);
2513 }
2514
2515 MagickExport Image *GaussianBlurImageChannel(const Image *image,
2516   const ChannelType channel,const double radius,const double sigma,
2517   ExceptionInfo *exception)
2518 {
2519   double
2520     *kernel;
2521
2522   Image
2523     *blur_image;
2524
2525   register ssize_t
2526     i;
2527
2528   size_t
2529     width;
2530
2531   ssize_t
2532     j,
2533     u,
2534     v;
2535
2536   assert(image != (const Image *) NULL);
2537   assert(image->signature == MagickSignature);
2538   if (image->debug != MagickFalse)
2539     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2540   assert(exception != (ExceptionInfo *) NULL);
2541   assert(exception->signature == MagickSignature);
2542   width=GetOptimalKernelWidth2D(radius,sigma);
2543   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2544   if (kernel == (double *) NULL)
2545     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2546   j=(ssize_t) width/2;
2547   i=0;
2548   for (v=(-j); v <= j; v++)
2549   {
2550     for (u=(-j); u <= j; u++)
2551       kernel[i++]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
2552         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
2553   }
2554   blur_image=ConvolveImageChannel(image,channel,width,kernel,exception);
2555   kernel=(double *) RelinquishMagickMemory(kernel);
2556   return(blur_image);
2557 }
2558 \f
2559 /*
2560 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2561 %                                                                             %
2562 %                                                                             %
2563 %                                                                             %
2564 %     M e d i a n F i l t e r I m a g e                                       %
2565 %                                                                             %
2566 %                                                                             %
2567 %                                                                             %
2568 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2569 %
2570 %  MedianFilterImage() applies a digital filter that improves the quality
2571 %  of a noisy image.  Each pixel is replaced by the median in a set of
2572 %  neighboring pixels as defined by radius.
2573 %
2574 %  The algorithm was contributed by Mike Edmonds and implements an insertion
2575 %  sort for selecting median color-channel values.  For more on this algorithm
2576 %  see "Skip Lists: A probabilistic Alternative to Balanced Trees" by William
2577 %  Pugh in the June 1990 of Communications of the ACM.
2578 %
2579 %  The format of the MedianFilterImage method is:
2580 %
2581 %      Image *MedianFilterImage(const Image *image,const double radius,
2582 %        ExceptionInfo *exception)
2583 %
2584 %  A description of each parameter follows:
2585 %
2586 %    o image: the image.
2587 %
2588 %    o radius: the radius of the pixel neighborhood.
2589 %
2590 %    o exception: return any errors or warnings in this structure.
2591 %
2592 */
2593
2594 #define MedianListChannels  5
2595
2596 typedef struct _MedianListNode
2597 {
2598   size_t
2599     next[9],
2600     count,
2601     signature;
2602 } MedianListNode;
2603
2604 typedef struct _MedianSkipList
2605 {
2606   ssize_t
2607     level;
2608
2609   MedianListNode
2610     *nodes;
2611 } MedianSkipList;
2612
2613 typedef struct _MedianPixelList
2614 {
2615   size_t
2616     center,
2617     seed,
2618     signature;
2619
2620   MedianSkipList
2621     lists[MedianListChannels];
2622 } MedianPixelList;
2623
2624 static MedianPixelList *DestroyMedianPixelList(MedianPixelList *pixel_list)
2625 {
2626   register ssize_t
2627     i;
2628
2629   if (pixel_list == (MedianPixelList *) NULL)
2630     return((MedianPixelList *) NULL);
2631   for (i=0; i < MedianListChannels; i++)
2632     if (pixel_list->lists[i].nodes != (MedianListNode *) NULL)
2633       pixel_list->lists[i].nodes=(MedianListNode *) RelinquishMagickMemory(
2634         pixel_list->lists[i].nodes);
2635   pixel_list=(MedianPixelList *) RelinquishMagickMemory(pixel_list);
2636   return(pixel_list);
2637 }
2638
2639 static MedianPixelList **DestroyMedianPixelListThreadSet(
2640   MedianPixelList **pixel_list)
2641 {
2642   register ssize_t
2643     i;
2644
2645   assert(pixel_list != (MedianPixelList **) NULL);
2646   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2647     if (pixel_list[i] != (MedianPixelList *) NULL)
2648       pixel_list[i]=DestroyMedianPixelList(pixel_list[i]);
2649   pixel_list=(MedianPixelList **) RelinquishMagickMemory(pixel_list);
2650   return(pixel_list);
2651 }
2652
2653 static MedianPixelList *AcquireMedianPixelList(const size_t width)
2654 {
2655   MedianPixelList
2656     *pixel_list;
2657
2658   register ssize_t
2659     i;
2660
2661   pixel_list=(MedianPixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2662   if (pixel_list == (MedianPixelList *) NULL)
2663     return(pixel_list);
2664   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
2665   pixel_list->center=width*width/2;
2666   for (i=0; i < MedianListChannels; i++)
2667   {
2668     pixel_list->lists[i].nodes=(MedianListNode *) AcquireQuantumMemory(65537UL,
2669       sizeof(*pixel_list->lists[i].nodes));
2670     if (pixel_list->lists[i].nodes == (MedianListNode *) NULL)
2671       return(DestroyMedianPixelList(pixel_list));
2672     (void) ResetMagickMemory(pixel_list->lists[i].nodes,0,65537UL*
2673       sizeof(*pixel_list->lists[i].nodes));
2674   }
2675   pixel_list->signature=MagickSignature;
2676   return(pixel_list);
2677 }
2678
2679 static MedianPixelList **AcquireMedianPixelListThreadSet(const size_t width)
2680 {
2681   MedianPixelList
2682     **pixel_list;
2683
2684   register ssize_t
2685     i;
2686
2687   size_t
2688     number_threads;
2689
2690   number_threads=GetOpenMPMaximumThreads();
2691   pixel_list=(MedianPixelList **) AcquireQuantumMemory(number_threads,
2692     sizeof(*pixel_list));
2693   if (pixel_list == (MedianPixelList **) NULL)
2694     return((MedianPixelList **) NULL);
2695   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
2696   for (i=0; i < (ssize_t) number_threads; i++)
2697   {
2698     pixel_list[i]=AcquireMedianPixelList(width);
2699     if (pixel_list[i] == (MedianPixelList *) NULL)
2700       return(DestroyMedianPixelListThreadSet(pixel_list));
2701   }
2702   return(pixel_list);
2703 }
2704
2705 static void AddNodeMedianPixelList(MedianPixelList *pixel_list,
2706   const ssize_t channel,const size_t color)
2707 {
2708   register MedianSkipList
2709     *list;
2710
2711   register ssize_t
2712     level;
2713
2714   size_t
2715     search,
2716     update[9];
2717
2718   /*
2719     Initialize the node.
2720   */
2721   list=pixel_list->lists+channel;
2722   list->nodes[color].signature=pixel_list->signature;
2723   list->nodes[color].count=1;
2724   /*
2725     Determine where it belongs in the list.
2726   */
2727   search=65536UL;
2728   for (level=list->level; level >= 0; level--)
2729   {
2730     while (list->nodes[search].next[level] < color)
2731       search=list->nodes[search].next[level];
2732     update[level]=search;
2733   }
2734   /*
2735     Generate a pseudo-random level for this node.
2736   */
2737   for (level=0; ; level++)
2738   {
2739     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2740     if ((pixel_list->seed & 0x300) != 0x300)
2741       break;
2742   }
2743   if (level > 8)
2744     level=8;
2745   if (level > (list->level+2))
2746     level=list->level+2;
2747   /*
2748     If we're raising the list's level, link back to the root node.
2749   */
2750   while (level > list->level)
2751   {
2752     list->level++;
2753     update[list->level]=65536UL;
2754   }
2755   /*
2756     Link the node into the skip-list.
2757   */
2758   do
2759   {
2760     list->nodes[color].next[level]=list->nodes[update[level]].next[level];
2761     list->nodes[update[level]].next[level]=color;
2762   }
2763   while (level-- > 0);
2764 }
2765
2766 static MagickPixelPacket GetMedianPixelList(MedianPixelList *pixel_list)
2767 {
2768   MagickPixelPacket
2769     pixel;
2770
2771   register MedianSkipList
2772     *list;
2773
2774   register ssize_t
2775     channel;
2776
2777   size_t
2778     center,
2779     color,
2780     count;
2781
2782   unsigned short
2783     channels[MedianListChannels];
2784
2785   /*
2786     Find the median value for each of the color.
2787   */
2788   center=pixel_list->center;
2789   for (channel=0; channel < 5; channel++)
2790   {
2791     list=pixel_list->lists+channel;
2792     color=65536UL;
2793     count=0;
2794     do
2795     {
2796       color=list->nodes[color].next[0];
2797       count+=list->nodes[color].count;
2798     }
2799     while (count <= center);
2800     channels[channel]=(unsigned short) color;
2801   }
2802   GetMagickPixelPacket((const Image *) NULL,&pixel);
2803   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
2804   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
2805   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
2806   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
2807   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
2808   return(pixel);
2809 }
2810
2811 static inline void InsertMedianPixelList(const Image *image,
2812   const PixelPacket *pixel,const IndexPacket *indexes,
2813   MedianPixelList *pixel_list)
2814 {
2815   size_t
2816     signature;
2817
2818   unsigned short
2819     index;
2820
2821   index=ScaleQuantumToShort(pixel->red);
2822   signature=pixel_list->lists[0].nodes[index].signature;
2823   if (signature == pixel_list->signature)
2824     pixel_list->lists[0].nodes[index].count++;
2825   else
2826     AddNodeMedianPixelList(pixel_list,0,index);
2827   index=ScaleQuantumToShort(pixel->green);
2828   signature=pixel_list->lists[1].nodes[index].signature;
2829   if (signature == pixel_list->signature)
2830     pixel_list->lists[1].nodes[index].count++;
2831   else
2832     AddNodeMedianPixelList(pixel_list,1,index);
2833   index=ScaleQuantumToShort(pixel->blue);
2834   signature=pixel_list->lists[2].nodes[index].signature;
2835   if (signature == pixel_list->signature)
2836     pixel_list->lists[2].nodes[index].count++;
2837   else
2838     AddNodeMedianPixelList(pixel_list,2,index);
2839   index=ScaleQuantumToShort(pixel->opacity);
2840   signature=pixel_list->lists[3].nodes[index].signature;
2841   if (signature == pixel_list->signature)
2842     pixel_list->lists[3].nodes[index].count++;
2843   else
2844     AddNodeMedianPixelList(pixel_list,3,index);
2845   if (image->colorspace == CMYKColorspace)
2846     index=ScaleQuantumToShort(*indexes);
2847   signature=pixel_list->lists[4].nodes[index].signature;
2848   if (signature == pixel_list->signature)
2849     pixel_list->lists[4].nodes[index].count++;
2850   else
2851     AddNodeMedianPixelList(pixel_list,4,index);
2852 }
2853
2854 static void ResetMedianPixelList(MedianPixelList *pixel_list)
2855 {
2856   int
2857     level;
2858
2859   register MedianListNode
2860     *root;
2861
2862   register MedianSkipList
2863     *list;
2864
2865   register ssize_t
2866     channel;
2867
2868   /*
2869     Reset the skip-list.
2870   */
2871   for (channel=0; channel < 5; channel++)
2872   {
2873     list=pixel_list->lists+channel;
2874     root=list->nodes+65536UL;
2875     list->level=0;
2876     for (level=0; level < 9; level++)
2877       root->next[level]=65536UL;
2878   }
2879   pixel_list->seed=pixel_list->signature++;
2880 }
2881
2882 MagickExport Image *MedianFilterImage(const Image *image,const double radius,
2883   ExceptionInfo *exception)
2884 {
2885 #define MedianFilterImageTag  "MedianFilter/Image"
2886
2887   CacheView
2888     *image_view,
2889     *median_view;
2890
2891   Image
2892     *median_image;
2893
2894   MagickBooleanType
2895     status;
2896
2897   MagickOffsetType
2898     progress;
2899
2900   MedianPixelList
2901     **restrict pixel_list;
2902
2903   size_t
2904     width;
2905
2906   ssize_t
2907     y;
2908
2909   /*
2910     Initialize median image attributes.
2911   */
2912   assert(image != (Image *) NULL);
2913   assert(image->signature == MagickSignature);
2914   if (image->debug != MagickFalse)
2915     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2916   assert(exception != (ExceptionInfo *) NULL);
2917   assert(exception->signature == MagickSignature);
2918   width=GetOptimalKernelWidth2D(radius,0.5);
2919   if ((image->columns < width) || (image->rows < width))
2920     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
2921   median_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2922     exception);
2923   if (median_image == (Image *) NULL)
2924     return((Image *) NULL);
2925   if (SetImageStorageClass(median_image,DirectClass) == MagickFalse)
2926     {
2927       InheritException(exception,&median_image->exception);
2928       median_image=DestroyImage(median_image);
2929       return((Image *) NULL);
2930     }
2931   pixel_list=AcquireMedianPixelListThreadSet(width);
2932   if (pixel_list == (MedianPixelList **) NULL)
2933     {
2934       median_image=DestroyImage(median_image);
2935       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2936     }
2937   /*
2938     Median filter each image row.
2939   */
2940   status=MagickTrue;
2941   progress=0;
2942   image_view=AcquireCacheView(image);
2943   median_view=AcquireCacheView(median_image);
2944 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2945   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2946 #endif
2947   for (y=0; y < (ssize_t) median_image->rows; y++)
2948   {
2949     const int
2950       id = GetOpenMPThreadId();
2951
2952     register const IndexPacket
2953       *restrict indexes;
2954
2955     register const PixelPacket
2956       *restrict p;
2957
2958     register IndexPacket
2959       *restrict median_indexes;
2960
2961     register PixelPacket
2962       *restrict q;
2963
2964     register ssize_t
2965       x;
2966
2967     if (status == MagickFalse)
2968       continue;
2969     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
2970       (width/2L),image->columns+width,width,exception);
2971     q=QueueCacheViewAuthenticPixels(median_view,0,y,median_image->columns,1,
2972       exception);
2973     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2974       {
2975         status=MagickFalse;
2976         continue;
2977       }
2978     indexes=GetCacheViewVirtualIndexQueue(image_view);
2979     median_indexes=GetCacheViewAuthenticIndexQueue(median_view);
2980     for (x=0; x < (ssize_t) median_image->columns; x++)
2981     {
2982       MagickPixelPacket
2983         pixel;
2984
2985       register const IndexPacket
2986         *restrict s;
2987
2988       register const PixelPacket
2989         *restrict r;
2990
2991       register ssize_t
2992         u,
2993         v;
2994
2995       r=p;
2996       s=indexes+x;
2997       ResetMedianPixelList(pixel_list[id]);
2998       for (v=0; v < (ssize_t) width; v++)
2999       {
3000         for (u=0; u < (ssize_t) width; u++)
3001           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
3002         r+=image->columns+width;
3003         s+=image->columns+width;
3004       }
3005       pixel=GetMedianPixelList(pixel_list[id]);
3006       SetPixelPacket(median_image,&pixel,q,median_indexes+x);
3007       p++;
3008       q++;
3009     }
3010     if (SyncCacheViewAuthenticPixels(median_view,exception) == MagickFalse)
3011       status=MagickFalse;
3012     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3013       {
3014         MagickBooleanType
3015           proceed;
3016
3017 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3018   #pragma omp critical (MagickCore_MedianFilterImage)
3019 #endif
3020         proceed=SetImageProgress(image,MedianFilterImageTag,progress++,
3021           image->rows);
3022         if (proceed == MagickFalse)
3023           status=MagickFalse;
3024       }
3025   }
3026   median_view=DestroyCacheView(median_view);
3027   image_view=DestroyCacheView(image_view);
3028   pixel_list=DestroyMedianPixelListThreadSet(pixel_list);
3029   return(median_image);
3030 }
3031 \f
3032 /*
3033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3034 %                                                                             %
3035 %                                                                             %
3036 %                                                                             %
3037 %     M o t i o n B l u r I m a g e                                           %
3038 %                                                                             %
3039 %                                                                             %
3040 %                                                                             %
3041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3042 %
3043 %  MotionBlurImage() simulates motion blur.  We convolve the image with a
3044 %  Gaussian operator of the given radius and standard deviation (sigma).
3045 %  For reasonable results, radius should be larger than sigma.  Use a
3046 %  radius of 0 and MotionBlurImage() selects a suitable radius for you.
3047 %  Angle gives the angle of the blurring motion.
3048 %
3049 %  Andrew Protano contributed this effect.
3050 %
3051 %  The format of the MotionBlurImage method is:
3052 %
3053 %    Image *MotionBlurImage(const Image *image,const double radius,
3054 %      const double sigma,const double angle,ExceptionInfo *exception)
3055 %    Image *MotionBlurImageChannel(const Image *image,const ChannelType channel,
3056 %      const double radius,const double sigma,const double angle,
3057 %      ExceptionInfo *exception)
3058 %
3059 %  A description of each parameter follows:
3060 %
3061 %    o image: the image.
3062 %
3063 %    o channel: the channel type.
3064 %
3065 %    o radius: the radius of the Gaussian, in pixels, not counting the center
3066 %    o radius: the radius of the Gaussian, in pixels, not counting
3067 %      the center pixel.
3068 %
3069 %    o sigma: the standard deviation of the Gaussian, in pixels.
3070 %
3071 %    o angle: Apply the effect along this angle.
3072 %
3073 %    o exception: return any errors or warnings in this structure.
3074 %
3075 */
3076
3077 static double *GetMotionBlurKernel(const size_t width,const double sigma)
3078 {
3079   double
3080     *kernel,
3081     normalize;
3082
3083   register ssize_t
3084     i;
3085
3086   /*
3087    Generate a 1-D convolution kernel.
3088   */
3089   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
3090   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
3091   if (kernel == (double *) NULL)
3092     return(kernel);
3093   normalize=0.0;
3094   for (i=0; i < (ssize_t) width; i++)
3095   {
3096     kernel[i]=(double) (exp((-((double) i*i)/(double) (2.0*MagickSigma*
3097       MagickSigma)))/(MagickSQ2PI*MagickSigma));
3098     normalize+=kernel[i];
3099   }
3100   for (i=0; i < (ssize_t) width; i++)
3101     kernel[i]/=normalize;
3102   return(kernel);
3103 }
3104
3105 MagickExport Image *MotionBlurImage(const Image *image,const double radius,
3106   const double sigma,const double angle,ExceptionInfo *exception)
3107 {
3108   Image
3109     *motion_blur;
3110
3111   motion_blur=MotionBlurImageChannel(image,DefaultChannels,radius,sigma,angle,
3112     exception);
3113   return(motion_blur);
3114 }
3115
3116 MagickExport Image *MotionBlurImageChannel(const Image *image,
3117   const ChannelType channel,const double radius,const double sigma,
3118   const double angle,ExceptionInfo *exception)
3119 {
3120   CacheView
3121     *blur_view,
3122     *image_view;
3123
3124   double
3125     *kernel;
3126
3127   Image
3128     *blur_image;
3129
3130   MagickBooleanType
3131     status;
3132
3133   MagickOffsetType
3134     progress;
3135
3136   MagickPixelPacket
3137     bias;
3138
3139   OffsetInfo
3140     *offset;
3141
3142   PointInfo
3143     point;
3144
3145   register ssize_t
3146     i;
3147
3148   size_t
3149     width;
3150
3151   ssize_t
3152     y;
3153
3154   assert(image != (Image *) NULL);
3155   assert(image->signature == MagickSignature);
3156   if (image->debug != MagickFalse)
3157     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3158   assert(exception != (ExceptionInfo *) NULL);
3159   width=GetOptimalKernelWidth1D(radius,sigma);
3160   kernel=GetMotionBlurKernel(width,sigma);
3161   if (kernel == (double *) NULL)
3162     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3163   offset=(OffsetInfo *) AcquireQuantumMemory(width,sizeof(*offset));
3164   if (offset == (OffsetInfo *) NULL)
3165     {
3166       kernel=(double *) RelinquishMagickMemory(kernel);
3167       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3168     }
3169   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3170   if (blur_image == (Image *) NULL)
3171     {
3172       kernel=(double *) RelinquishMagickMemory(kernel);
3173       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
3174       return((Image *) NULL);
3175     }
3176   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
3177     {
3178       kernel=(double *) RelinquishMagickMemory(kernel);
3179       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
3180       InheritException(exception,&blur_image->exception);
3181       blur_image=DestroyImage(blur_image);
3182       return((Image *) NULL);
3183     }
3184   point.x=(double) width*sin(DegreesToRadians(angle));
3185   point.y=(double) width*cos(DegreesToRadians(angle));
3186   for (i=0; i < (ssize_t) width; i++)
3187   {
3188     offset[i].x=(ssize_t) ceil((double) (i*point.y)/hypot(point.x,point.y)-0.5);
3189     offset[i].y=(ssize_t) ceil((double) (i*point.x)/hypot(point.x,point.y)-0.5);
3190   }
3191   /*
3192     Motion blur image.
3193   */
3194   status=MagickTrue;
3195   progress=0;
3196   GetMagickPixelPacket(image,&bias);
3197   image_view=AcquireCacheView(image);
3198   blur_view=AcquireCacheView(blur_image);
3199 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3200   #pragma omp parallel for schedule(dynamic,4) shared(progress,status) omp_throttle(1)
3201 #endif
3202   for (y=0; y < (ssize_t) image->rows; y++)
3203   {
3204     register IndexPacket
3205       *restrict blur_indexes;
3206
3207     register PixelPacket
3208       *restrict q;
3209
3210     register ssize_t
3211       x;
3212
3213     if (status == MagickFalse)
3214       continue;
3215     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
3216       exception);
3217     if (q == (PixelPacket *) NULL)
3218       {
3219         status=MagickFalse;
3220         continue;
3221       }
3222     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
3223     for (x=0; x < (ssize_t) image->columns; x++)
3224     {
3225       MagickPixelPacket
3226         qixel;
3227
3228       PixelPacket
3229         pixel;
3230
3231       register const IndexPacket
3232         *restrict indexes;
3233
3234       register double
3235         *restrict k;
3236
3237       register ssize_t
3238         i;
3239
3240       k=kernel;
3241       qixel=bias;
3242       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
3243         {
3244           for (i=0; i < (ssize_t) width; i++)
3245           {
3246             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
3247               offset[i].y,&pixel,exception);
3248             qixel.red+=(*k)*pixel.red;
3249             qixel.green+=(*k)*pixel.green;
3250             qixel.blue+=(*k)*pixel.blue;
3251             qixel.opacity+=(*k)*pixel.opacity;
3252             if (image->colorspace == CMYKColorspace)
3253               {
3254                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3255                 qixel.index+=(*k)*(*indexes);
3256               }
3257             k++;
3258           }
3259           if ((channel & RedChannel) != 0)
3260             q->red=ClampToQuantum(qixel.red);
3261           if ((channel & GreenChannel) != 0)
3262             q->green=ClampToQuantum(qixel.green);
3263           if ((channel & BlueChannel) != 0)
3264             q->blue=ClampToQuantum(qixel.blue);
3265           if ((channel & OpacityChannel) != 0)
3266             q->opacity=ClampToQuantum(qixel.opacity);
3267           if (((channel & IndexChannel) != 0) &&
3268               (image->colorspace == CMYKColorspace))
3269             blur_indexes[x]=(IndexPacket) ClampToQuantum(qixel.index);
3270         }
3271       else
3272         {
3273           MagickRealType
3274             alpha,
3275             gamma;
3276
3277           alpha=0.0;
3278           gamma=0.0;
3279           for (i=0; i < (ssize_t) width; i++)
3280           {
3281             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
3282               offset[i].y,&pixel,exception);
3283             alpha=(MagickRealType) (QuantumScale*
3284               GetAlphaPixelComponent(&pixel));
3285             qixel.red+=(*k)*alpha*pixel.red;
3286             qixel.green+=(*k)*alpha*pixel.green;
3287             qixel.blue+=(*k)*alpha*pixel.blue;
3288             qixel.opacity+=(*k)*pixel.opacity;
3289             if (image->colorspace == CMYKColorspace)
3290               {
3291                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3292                 qixel.index+=(*k)*alpha*(*indexes);
3293               }
3294             gamma+=(*k)*alpha;
3295             k++;
3296           }
3297           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
3298           if ((channel & RedChannel) != 0)
3299             q->red=ClampToQuantum(gamma*qixel.red);
3300           if ((channel & GreenChannel) != 0)
3301             q->green=ClampToQuantum(gamma*qixel.green);
3302           if ((channel & BlueChannel) != 0)
3303             q->blue=ClampToQuantum(gamma*qixel.blue);
3304           if ((channel & OpacityChannel) != 0)
3305             q->opacity=ClampToQuantum(qixel.opacity);
3306           if (((channel & IndexChannel) != 0) &&
3307               (image->colorspace == CMYKColorspace))
3308             blur_indexes[x]=(IndexPacket) ClampToQuantum(gamma*qixel.index);
3309         }
3310       q++;
3311     }
3312     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
3313       status=MagickFalse;
3314     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3315       {
3316         MagickBooleanType
3317           proceed;
3318
3319 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3320   #pragma omp critical (MagickCore_MotionBlurImageChannel)
3321 #endif
3322         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
3323         if (proceed == MagickFalse)
3324           status=MagickFalse;
3325       }
3326   }
3327   blur_view=DestroyCacheView(blur_view);
3328   image_view=DestroyCacheView(image_view);
3329   kernel=(double *) RelinquishMagickMemory(kernel);
3330   offset=(OffsetInfo *) RelinquishMagickMemory(offset);
3331   if (status == MagickFalse)
3332     blur_image=DestroyImage(blur_image);
3333   return(blur_image);
3334 }
3335 \f
3336 /*
3337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3338 %                                                                             %
3339 %                                                                             %
3340 %                                                                             %
3341 %     P r e v i e w I m a g e                                                 %
3342 %                                                                             %
3343 %                                                                             %
3344 %                                                                             %
3345 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3346 %
3347 %  PreviewImage() tiles 9 thumbnails of the specified image with an image
3348 %  processing operation applied with varying parameters.  This may be helpful
3349 %  pin-pointing an appropriate parameter for a particular image processing
3350 %  operation.
3351 %
3352 %  The format of the PreviewImages method is:
3353 %
3354 %      Image *PreviewImages(const Image *image,const PreviewType preview,
3355 %        ExceptionInfo *exception)
3356 %
3357 %  A description of each parameter follows:
3358 %
3359 %    o image: the image.
3360 %
3361 %    o preview: the image processing operation.
3362 %
3363 %    o exception: return any errors or warnings in this structure.
3364 %
3365 */
3366 MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
3367   ExceptionInfo *exception)
3368 {
3369 #define NumberTiles  9
3370 #define PreviewImageTag  "Preview/Image"
3371 #define DefaultPreviewGeometry  "204x204+10+10"
3372
3373   char
3374     factor[MaxTextExtent],
3375     label[MaxTextExtent];
3376
3377   double
3378     degrees,
3379     gamma,
3380     percentage,
3381     radius,
3382     sigma,
3383     threshold;
3384
3385   Image
3386     *images,
3387     *montage_image,
3388     *preview_image,
3389     *thumbnail;
3390
3391   ImageInfo
3392     *preview_info;
3393
3394   MagickBooleanType
3395     proceed;
3396
3397   MontageInfo
3398     *montage_info;
3399
3400   QuantizeInfo
3401     quantize_info;
3402
3403   RectangleInfo
3404     geometry;
3405
3406   register ssize_t
3407     i,
3408     x;
3409
3410   size_t
3411     colors;
3412
3413   ssize_t
3414     y;
3415
3416   /*
3417     Open output image file.
3418   */
3419   assert(image != (Image *) NULL);
3420   assert(image->signature == MagickSignature);
3421   if (image->debug != MagickFalse)
3422     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3423   colors=2;
3424   degrees=0.0;
3425   gamma=(-0.2f);
3426   preview_info=AcquireImageInfo();
3427   SetGeometry(image,&geometry);
3428   (void) ParseMetaGeometry(DefaultPreviewGeometry,&geometry.x,&geometry.y,
3429     &geometry.width,&geometry.height);
3430   images=NewImageList();
3431   percentage=12.5;
3432   GetQuantizeInfo(&quantize_info);
3433   radius=0.0;
3434   sigma=1.0;
3435   threshold=0.0;
3436   x=0;
3437   y=0;
3438   for (i=0; i < NumberTiles; i++)
3439   {
3440     thumbnail=ThumbnailImage(image,geometry.width,geometry.height,exception);
3441     if (thumbnail == (Image *) NULL)
3442       break;
3443     (void) SetImageProgressMonitor(thumbnail,(MagickProgressMonitor) NULL,
3444       (void *) NULL);
3445     (void) SetImageProperty(thumbnail,"label",DefaultTileLabel);
3446     if (i == (NumberTiles/2))
3447       {
3448         (void) QueryColorDatabase("#dfdfdf",&thumbnail->matte_color,exception);
3449         AppendImageToList(&images,thumbnail);
3450         continue;
3451       }
3452     switch (preview)
3453     {
3454       case RotatePreview:
3455       {
3456         degrees+=45.0;
3457         preview_image=RotateImage(thumbnail,degrees,exception);
3458         (void) FormatMagickString(label,MaxTextExtent,"rotate %g",degrees);
3459         break;
3460       }
3461       case ShearPreview:
3462       {
3463         degrees+=5.0;
3464         preview_image=ShearImage(thumbnail,degrees,degrees,exception);
3465         (void) FormatMagickString(label,MaxTextExtent,"shear %gx%g",
3466           degrees,2.0*degrees);
3467         break;
3468       }
3469       case RollPreview:
3470       {
3471         x=(ssize_t) ((i+1)*thumbnail->columns)/NumberTiles;
3472         y=(ssize_t) ((i+1)*thumbnail->rows)/NumberTiles;
3473         preview_image=RollImage(thumbnail,x,y,exception);
3474         (void) FormatMagickString(label,MaxTextExtent,"roll %+.20gx%+.20g",
3475           (double) x,(double) y);
3476         break;
3477       }
3478       case HuePreview:
3479       {
3480         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3481         if (preview_image == (Image *) NULL)
3482           break;
3483         (void) FormatMagickString(factor,MaxTextExtent,"100,100,%g",
3484           2.0*percentage);
3485         (void) ModulateImage(preview_image,factor);
3486         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3487         break;
3488       }
3489       case SaturationPreview:
3490       {
3491         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3492         if (preview_image == (Image *) NULL)
3493           break;
3494         (void) FormatMagickString(factor,MaxTextExtent,"100,%g",
3495           2.0*percentage);
3496         (void) ModulateImage(preview_image,factor);
3497         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3498         break;
3499       }
3500       case BrightnessPreview:
3501       {
3502         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3503         if (preview_image == (Image *) NULL)
3504           break;
3505         (void) FormatMagickString(factor,MaxTextExtent,"%g",2.0*percentage);
3506         (void) ModulateImage(preview_image,factor);
3507         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3508         break;
3509       }
3510       case GammaPreview:
3511       default:
3512       {
3513         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3514         if (preview_image == (Image *) NULL)
3515           break;
3516         gamma+=0.4f;
3517         (void) GammaImageChannel(preview_image,DefaultChannels,gamma);
3518         (void) FormatMagickString(label,MaxTextExtent,"gamma %g",gamma);
3519         break;
3520       }
3521       case SpiffPreview:
3522       {
3523         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3524         if (preview_image != (Image *) NULL)
3525           for (x=0; x < i; x++)
3526             (void) ContrastImage(preview_image,MagickTrue);
3527         (void) FormatMagickString(label,MaxTextExtent,"contrast (%.20g)",
3528           (double) i+1);
3529         break;
3530       }
3531       case DullPreview:
3532       {
3533         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3534         if (preview_image == (Image *) NULL)
3535           break;
3536         for (x=0; x < i; x++)
3537           (void) ContrastImage(preview_image,MagickFalse);
3538         (void) FormatMagickString(label,MaxTextExtent,"+contrast (%.20g)",
3539           (double) i+1);
3540         break;
3541       }
3542       case GrayscalePreview:
3543       {
3544         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3545         if (preview_image == (Image *) NULL)
3546           break;
3547         colors<<=1;
3548         quantize_info.number_colors=colors;
3549         quantize_info.colorspace=GRAYColorspace;
3550         (void) QuantizeImage(&quantize_info,preview_image);
3551         (void) FormatMagickString(label,MaxTextExtent,
3552           "-colorspace gray -colors %.20g",(double) colors);
3553         break;
3554       }
3555       case QuantizePreview:
3556       {
3557         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3558         if (preview_image == (Image *) NULL)
3559           break;
3560         colors<<=1;
3561         quantize_info.number_colors=colors;
3562         (void) QuantizeImage(&quantize_info,preview_image);
3563         (void) FormatMagickString(label,MaxTextExtent,"colors %.20g",(double)
3564           colors);
3565         break;
3566       }
3567       case DespecklePreview:
3568       {
3569         for (x=0; x < (i-1); x++)
3570         {
3571           preview_image=DespeckleImage(thumbnail,exception);
3572           if (preview_image == (Image *) NULL)
3573             break;
3574           thumbnail=DestroyImage(thumbnail);
3575           thumbnail=preview_image;
3576         }
3577         preview_image=DespeckleImage(thumbnail,exception);
3578         if (preview_image == (Image *) NULL)
3579           break;
3580         (void) FormatMagickString(label,MaxTextExtent,"despeckle (%.20g)",
3581           (double) i+1);
3582         break;
3583       }
3584       case ReduceNoisePreview:
3585       {
3586         preview_image=ReduceNoiseImage(thumbnail,radius,exception);
3587         (void) FormatMagickString(label,MaxTextExtent,"noise %g",radius);
3588         break;
3589       }
3590       case AddNoisePreview:
3591       {
3592         switch ((int) i)
3593         {
3594           case 0:
3595           {
3596             (void) CopyMagickString(factor,"uniform",MaxTextExtent);
3597             break;
3598           }
3599           case 1:
3600           {
3601             (void) CopyMagickString(factor,"gaussian",MaxTextExtent);
3602             break;
3603           }
3604           case 2:
3605           {
3606             (void) CopyMagickString(factor,"multiplicative",MaxTextExtent);
3607             break;
3608           }
3609           case 3:
3610           {
3611             (void) CopyMagickString(factor,"impulse",MaxTextExtent);
3612             break;
3613           }
3614           case 4:
3615           {
3616             (void) CopyMagickString(factor,"laplacian",MaxTextExtent);
3617             break;
3618           }
3619           case 5:
3620           {
3621             (void) CopyMagickString(factor,"Poisson",MaxTextExtent);
3622             break;
3623           }
3624           default:
3625           {
3626             (void) CopyMagickString(thumbnail->magick,"NULL",MaxTextExtent);
3627             break;
3628           }
3629         }
3630         preview_image=ReduceNoiseImage(thumbnail,(double) i,exception);
3631         (void) FormatMagickString(label,MaxTextExtent,"+noise %s",factor);
3632         break;
3633       }
3634       case SharpenPreview:
3635       {
3636         preview_image=SharpenImage(thumbnail,radius,sigma,exception);
3637         (void) FormatMagickString(label,MaxTextExtent,"sharpen %gx%g",
3638           radius,sigma);
3639         break;
3640       }
3641       case BlurPreview:
3642       {
3643         preview_image=BlurImage(thumbnail,radius,sigma,exception);
3644         (void) FormatMagickString(label,MaxTextExtent,"blur %gx%g",radius,
3645           sigma);
3646         break;
3647       }
3648       case ThresholdPreview:
3649       {
3650         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3651         if (preview_image == (Image *) NULL)
3652           break;
3653         (void) BilevelImage(thumbnail,
3654           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
3655         (void) FormatMagickString(label,MaxTextExtent,"threshold %g",
3656           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
3657         break;
3658       }
3659       case EdgeDetectPreview:
3660       {
3661         preview_image=EdgeImage(thumbnail,radius,exception);
3662         (void) FormatMagickString(label,MaxTextExtent,"edge %g",radius);
3663         break;
3664       }
3665       case SpreadPreview:
3666       {
3667         preview_image=SpreadImage(thumbnail,radius,exception);
3668         (void) FormatMagickString(label,MaxTextExtent,"spread %g",
3669           radius+0.5);
3670         break;
3671       }
3672       case SolarizePreview:
3673       {
3674         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3675         if (preview_image == (Image *) NULL)
3676           break;
3677         (void) SolarizeImage(preview_image,(double) QuantumRange*
3678           percentage/100.0);
3679         (void) FormatMagickString(label,MaxTextExtent,"solarize %g",
3680           (QuantumRange*percentage)/100.0);
3681         break;
3682       }
3683       case ShadePreview:
3684       {
3685         degrees+=10.0;
3686         preview_image=ShadeImage(thumbnail,MagickTrue,degrees,degrees,
3687           exception);
3688         (void) FormatMagickString(label,MaxTextExtent,"shade %gx%g",
3689           degrees,degrees);
3690         break;
3691       }
3692       case RaisePreview:
3693       {
3694         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3695         if (preview_image == (Image *) NULL)
3696           break;
3697         geometry.width=(size_t) (2*i+2);
3698         geometry.height=(size_t) (2*i+2);
3699         geometry.x=i/2;
3700         geometry.y=i/2;
3701         (void) RaiseImage(preview_image,&geometry,MagickTrue);
3702         (void) FormatMagickString(label,MaxTextExtent,
3703           "raise %.20gx%.20g%+.20g%+.20g",(double) geometry.width,(double)
3704           geometry.height,(double) geometry.x,(double) geometry.y);
3705         break;
3706       }
3707       case SegmentPreview:
3708       {
3709         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3710         if (preview_image == (Image *) NULL)
3711           break;
3712         threshold+=0.4f;
3713         (void) SegmentImage(preview_image,RGBColorspace,MagickFalse,threshold,
3714           threshold);
3715         (void) FormatMagickString(label,MaxTextExtent,"segment %gx%g",
3716           threshold,threshold);
3717         break;
3718       }
3719       case SwirlPreview:
3720       {
3721         preview_image=SwirlImage(thumbnail,degrees,exception);
3722         (void) FormatMagickString(label,MaxTextExtent,"swirl %g",degrees);
3723         degrees+=45.0;
3724         break;
3725       }
3726       case ImplodePreview:
3727       {
3728         degrees+=0.1f;
3729         preview_image=ImplodeImage(thumbnail,degrees,exception);
3730         (void) FormatMagickString(label,MaxTextExtent,"implode %g",degrees);
3731         break;
3732       }
3733       case WavePreview:
3734       {
3735         degrees+=5.0f;
3736         preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,exception);
3737         (void) FormatMagickString(label,MaxTextExtent,"wave %gx%g",
3738           0.5*degrees,2.0*degrees);
3739         break;
3740       }
3741       case OilPaintPreview:
3742       {
3743         preview_image=OilPaintImage(thumbnail,(double) radius,exception);
3744         (void) FormatMagickString(label,MaxTextExtent,"paint %g",radius);
3745         break;
3746       }
3747       case CharcoalDrawingPreview:
3748       {
3749         preview_image=CharcoalImage(thumbnail,(double) radius,(double) sigma,
3750           exception);
3751         (void) FormatMagickString(label,MaxTextExtent,"charcoal %gx%g",
3752           radius,sigma);
3753         break;
3754       }
3755       case JPEGPreview:
3756       {
3757         char
3758           filename[MaxTextExtent];
3759
3760         int
3761           file;
3762
3763         MagickBooleanType
3764           status;
3765
3766         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3767         if (preview_image == (Image *) NULL)
3768           break;
3769         preview_info->quality=(size_t) percentage;
3770         (void) FormatMagickString(factor,MaxTextExtent,"%.20g",(double)
3771           preview_info->quality);
3772         file=AcquireUniqueFileResource(filename);
3773         if (file != -1)
3774           file=close(file)-1;
3775         (void) FormatMagickString(preview_image->filename,MaxTextExtent,
3776           "jpeg:%s",filename);
3777         status=WriteImage(preview_info,preview_image);
3778         if (status != MagickFalse)
3779           {
3780             Image
3781               *quality_image;
3782
3783             (void) CopyMagickString(preview_info->filename,
3784               preview_image->filename,MaxTextExtent);
3785             quality_image=ReadImage(preview_info,exception);
3786             if (quality_image != (Image *) NULL)
3787               {
3788                 preview_image=DestroyImage(preview_image);
3789                 preview_image=quality_image;
3790               }
3791           }
3792         (void) RelinquishUniqueFileResource(preview_image->filename);
3793         if ((GetBlobSize(preview_image)/1024) >= 1024)
3794           (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%gmb ",
3795             factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
3796             1024.0/1024.0);
3797         else
3798           if (GetBlobSize(preview_image) >= 1024)
3799             (void) FormatMagickString(label,MaxTextExtent,
3800               "quality %s\n%gkb ",factor,(double) ((MagickOffsetType)
3801               GetBlobSize(preview_image))/1024.0);
3802           else
3803             (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%.20gb ",
3804               factor,(double) GetBlobSize(thumbnail));
3805         break;
3806       }
3807     }
3808     thumbnail=DestroyImage(thumbnail);
3809     percentage+=12.5;
3810     radius+=0.5;
3811     sigma+=0.25;
3812     if (preview_image == (Image *) NULL)
3813       break;
3814     (void) DeleteImageProperty(preview_image,"label");
3815     (void) SetImageProperty(preview_image,"label",label);
3816     AppendImageToList(&images,preview_image);
3817     proceed=SetImageProgress(image,PreviewImageTag,(MagickOffsetType) i,
3818       NumberTiles);
3819     if (proceed == MagickFalse)
3820       break;
3821   }
3822   if (images == (Image *) NULL)
3823     {
3824       preview_info=DestroyImageInfo(preview_info);
3825       return((Image *) NULL);
3826     }
3827   /*
3828     Create the montage.
3829   */
3830   montage_info=CloneMontageInfo(preview_info,(MontageInfo *) NULL);
3831   (void) CopyMagickString(montage_info->filename,image->filename,MaxTextExtent);
3832   montage_info->shadow=MagickTrue;
3833   (void) CloneString(&montage_info->tile,"3x3");
3834   (void) CloneString(&montage_info->geometry,DefaultPreviewGeometry);
3835   (void) CloneString(&montage_info->frame,DefaultTileFrame);
3836   montage_image=MontageImages(images,montage_info,exception);
3837   montage_info=DestroyMontageInfo(montage_info);
3838   images=DestroyImageList(images);
3839   if (montage_image == (Image *) NULL)
3840     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3841   if (montage_image->montage != (char *) NULL)
3842     {
3843       /*
3844         Free image directory.
3845       */
3846       montage_image->montage=(char *) RelinquishMagickMemory(
3847         montage_image->montage);
3848       if (image->directory != (char *) NULL)
3849         montage_image->directory=(char *) RelinquishMagickMemory(
3850           montage_image->directory);
3851     }
3852   preview_info=DestroyImageInfo(preview_info);
3853   return(montage_image);
3854 }
3855 \f
3856 /*
3857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3858 %                                                                             %
3859 %                                                                             %
3860 %                                                                             %
3861 %     R a d i a l B l u r I m a g e                                           %
3862 %                                                                             %
3863 %                                                                             %
3864 %                                                                             %
3865 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3866 %
3867 %  RadialBlurImage() applies a radial blur to the image.
3868 %
3869 %  Andrew Protano contributed this effect.
3870 %
3871 %  The format of the RadialBlurImage method is:
3872 %
3873 %    Image *RadialBlurImage(const Image *image,const double angle,
3874 %      ExceptionInfo *exception)
3875 %    Image *RadialBlurImageChannel(const Image *image,const ChannelType channel,
3876 %      const double angle,ExceptionInfo *exception)
3877 %
3878 %  A description of each parameter follows:
3879 %
3880 %    o image: the image.
3881 %
3882 %    o channel: the channel type.
3883 %
3884 %    o angle: the angle of the radial blur.
3885 %
3886 %    o exception: return any errors or warnings in this structure.
3887 %
3888 */
3889
3890 MagickExport Image *RadialBlurImage(const Image *image,const double angle,
3891   ExceptionInfo *exception)
3892 {
3893   Image
3894     *blur_image;
3895
3896   blur_image=RadialBlurImageChannel(image,DefaultChannels,angle,exception);
3897   return(blur_image);
3898 }
3899
3900 MagickExport Image *RadialBlurImageChannel(const Image *image,
3901   const ChannelType channel,const double angle,ExceptionInfo *exception)
3902 {
3903   CacheView
3904     *blur_view,
3905     *image_view;
3906
3907   Image
3908     *blur_image;
3909
3910   MagickBooleanType
3911     status;
3912
3913   MagickOffsetType
3914     progress;
3915
3916   MagickPixelPacket
3917     bias;
3918
3919   MagickRealType
3920     blur_radius,
3921     *cos_theta,
3922     offset,
3923     *sin_theta,
3924     theta;
3925
3926   PointInfo
3927     blur_center;
3928
3929   register ssize_t
3930     i;
3931
3932   size_t
3933     n;
3934
3935   ssize_t
3936     y;
3937
3938   /*
3939     Allocate blur image.
3940   */
3941   assert(image != (Image *) NULL);
3942   assert(image->signature == MagickSignature);
3943   if (image->debug != MagickFalse)
3944     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3945   assert(exception != (ExceptionInfo *) NULL);
3946   assert(exception->signature == MagickSignature);
3947   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3948   if (blur_image == (Image *) NULL)
3949     return((Image *) NULL);
3950   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
3951     {
3952       InheritException(exception,&blur_image->exception);
3953       blur_image=DestroyImage(blur_image);
3954       return((Image *) NULL);
3955     }
3956   blur_center.x=(double) image->columns/2.0;
3957   blur_center.y=(double) image->rows/2.0;
3958   blur_radius=hypot(blur_center.x,blur_center.y);
3959   n=(size_t) fabs(4.0*DegreesToRadians(angle)*sqrt((double) blur_radius)+2UL);
3960   theta=DegreesToRadians(angle)/(MagickRealType) (n-1);
3961   cos_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
3962     sizeof(*cos_theta));
3963   sin_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
3964     sizeof(*sin_theta));
3965   if ((cos_theta == (MagickRealType *) NULL) ||
3966       (sin_theta == (MagickRealType *) NULL))
3967     {
3968       blur_image=DestroyImage(blur_image);
3969       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3970     }
3971   offset=theta*(MagickRealType) (n-1)/2.0;
3972   for (i=0; i < (ssize_t) n; i++)
3973   {
3974     cos_theta[i]=cos((double) (theta*i-offset));
3975     sin_theta[i]=sin((double) (theta*i-offset));
3976   }
3977   /*
3978     Radial blur image.
3979   */
3980   status=MagickTrue;
3981   progress=0;
3982   GetMagickPixelPacket(image,&bias);
3983   image_view=AcquireCacheView(image);
3984   blur_view=AcquireCacheView(blur_image);
3985 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3986   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
3987 #endif
3988   for (y=0; y < (ssize_t) blur_image->rows; y++)
3989   {
3990     register const IndexPacket
3991       *restrict indexes;
3992
3993     register IndexPacket
3994       *restrict blur_indexes;
3995
3996     register PixelPacket
3997       *restrict q;
3998
3999     register ssize_t
4000       x;
4001
4002     if (status == MagickFalse)
4003       continue;
4004     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
4005       exception);
4006     if (q == (PixelPacket *) NULL)
4007       {
4008         status=MagickFalse;
4009         continue;
4010       }
4011     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
4012     for (x=0; x < (ssize_t) blur_image->columns; x++)
4013     {
4014       MagickPixelPacket
4015         qixel;
4016
4017       MagickRealType
4018         normalize,
4019         radius;
4020
4021       PixelPacket
4022         pixel;
4023
4024       PointInfo
4025         center;
4026
4027       register ssize_t
4028         i;
4029
4030       size_t
4031         step;
4032
4033       center.x=(double) x-blur_center.x;
4034       center.y=(double) y-blur_center.y;
4035       radius=hypot((double) center.x,center.y);
4036       if (radius == 0)
4037         step=1;
4038       else
4039         {
4040           step=(size_t) (blur_radius/radius);
4041           if (step == 0)
4042             step=1;
4043           else
4044             if (step >= n)
4045               step=n-1;
4046         }
4047       normalize=0.0;
4048       qixel=bias;
4049       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
4050         {
4051           for (i=0; i < (ssize_t) n; i+=(ssize_t) step)
4052           {
4053             (void) GetOneCacheViewVirtualPixel(image_view,(ssize_t)
4054               (blur_center.x+center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),
4055               (ssize_t) (blur_center.y+center.x*sin_theta[i]+center.y*
4056               cos_theta[i]+0.5),&pixel,exception);
4057             qixel.red+=pixel.red;
4058             qixel.green+=pixel.green;
4059             qixel.blue+=pixel.blue;
4060             qixel.opacity+=pixel.opacity;
4061             if (image->colorspace == CMYKColorspace)
4062               {
4063                 indexes=GetCacheViewVirtualIndexQueue(image_view);
4064                 qixel.index+=(*indexes);
4065               }
4066             normalize+=1.0;
4067           }
4068           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
4069             normalize);
4070           if ((channel & RedChannel) != 0)
4071             q->red=ClampToQuantum(normalize*qixel.red);
4072           if ((channel & GreenChannel) != 0)
4073             q->green=ClampToQuantum(normalize*qixel.green);
4074           if ((channel & BlueChannel) != 0)
4075             q->blue=ClampToQuantum(normalize*qixel.blue);
4076           if ((channel & OpacityChannel) != 0)
4077             q->opacity=ClampToQuantum(normalize*qixel.opacity);
4078           if (((channel & IndexChannel) != 0) &&
4079               (image->colorspace == CMYKColorspace))
4080             blur_indexes[x]=(IndexPacket) ClampToQuantum(normalize*qixel.index);
4081         }
4082       else
4083         {
4084           MagickRealType
4085             alpha,
4086             gamma;
4087
4088           alpha=1.0;
4089           gamma=0.0;
4090           for (i=0; i < (ssize_t) n; i+=(ssize_t) step)
4091           {
4092             (void) GetOneCacheViewVirtualPixel(image_view,(ssize_t)
4093               (blur_center.x+center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),
4094               (ssize_t) (blur_center.y+center.x*sin_theta[i]+center.y*
4095               cos_theta[i]+0.5),&pixel,exception);
4096             alpha=(MagickRealType) (QuantumScale*
4097               GetAlphaPixelComponent(&pixel));
4098             qixel.red+=alpha*pixel.red;
4099             qixel.green+=alpha*pixel.green;
4100             qixel.blue+=alpha*pixel.blue;
4101             qixel.opacity+=pixel.opacity;
4102             if (image->colorspace == CMYKColorspace)
4103               {
4104                 indexes=GetCacheViewVirtualIndexQueue(image_view);
4105                 qixel.index+=alpha*(*indexes);
4106               }
4107             gamma+=alpha;
4108             normalize+=1.0;
4109           }
4110           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4111           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
4112             normalize);
4113           if ((channel & RedChannel) != 0)
4114             q->red=ClampToQuantum(gamma*qixel.red);
4115           if ((channel & GreenChannel) != 0)
4116             q->green=ClampToQuantum(gamma*qixel.green);
4117           if ((channel & BlueChannel) != 0)
4118             q->blue=ClampToQuantum(gamma*qixel.blue);
4119           if ((channel & OpacityChannel) != 0)
4120             q->opacity=ClampToQuantum(normalize*qixel.opacity);
4121           if (((channel & IndexChannel) != 0) &&
4122               (image->colorspace == CMYKColorspace))
4123             blur_indexes[x]=(IndexPacket) ClampToQuantum(gamma*qixel.index);
4124         }
4125       q++;
4126     }
4127     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
4128       status=MagickFalse;
4129     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4130       {
4131         MagickBooleanType
4132           proceed;
4133
4134 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4135   #pragma omp critical (MagickCore_RadialBlurImageChannel)
4136 #endif
4137         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
4138         if (proceed == MagickFalse)
4139           status=MagickFalse;
4140       }
4141   }
4142   blur_view=DestroyCacheView(blur_view);
4143   image_view=DestroyCacheView(image_view);
4144   cos_theta=(MagickRealType *) RelinquishMagickMemory(cos_theta);
4145   sin_theta=(MagickRealType *) RelinquishMagickMemory(sin_theta);
4146   if (status == MagickFalse)
4147     blur_image=DestroyImage(blur_image);
4148   return(blur_image);
4149 }
4150 \f
4151 /*
4152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4153 %                                                                             %
4154 %                                                                             %
4155 %                                                                             %
4156 %     R e d u c e N o i s e I m a g e                                         %
4157 %                                                                             %
4158 %                                                                             %
4159 %                                                                             %
4160 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4161 %
4162 %  ReduceNoiseImage() smooths the contours of an image while still preserving
4163 %  edge information.  The algorithm works by replacing each pixel with its
4164 %  neighbor closest in value.  A neighbor is defined by radius.  Use a radius
4165 %  of 0 and ReduceNoise() selects a suitable radius for you.
4166 %
4167 %  The format of the ReduceNoiseImage method is:
4168 %
4169 %      Image *ReduceNoiseImage(const Image *image,const double radius,
4170 %        ExceptionInfo *exception)
4171 %
4172 %  A description of each parameter follows:
4173 %
4174 %    o image: the image.
4175 %
4176 %    o radius: the radius of the pixel neighborhood.
4177 %
4178 %    o exception: return any errors or warnings in this structure.
4179 %
4180 */
4181
4182 static MagickPixelPacket GetNonpeakMedianPixelList(MedianPixelList *pixel_list)
4183 {
4184   MagickPixelPacket
4185     pixel;
4186
4187   register MedianSkipList
4188     *list;
4189
4190   register ssize_t
4191     channel;
4192
4193   size_t
4194     center,
4195     color,
4196     count,
4197     previous,
4198     next;
4199
4200   unsigned short
4201     channels[5];
4202
4203   /*
4204     Finds the median value for each of the color.
4205   */
4206   center=pixel_list->center;
4207   for (channel=0; channel < 5; channel++)
4208   {
4209     list=pixel_list->lists+channel;
4210     color=65536UL;
4211     next=list->nodes[color].next[0];
4212     count=0;
4213     do
4214     {
4215       previous=color;
4216       color=next;
4217       next=list->nodes[color].next[0];
4218       count+=list->nodes[color].count;
4219     }
4220     while (count <= center);
4221     if ((previous == 65536UL) && (next != 65536UL))
4222       color=next;
4223     else
4224       if ((previous != 65536UL) && (next == 65536UL))
4225         color=previous;
4226     channels[channel]=(unsigned short) color;
4227   }
4228   GetMagickPixelPacket((const Image *) NULL,&pixel);
4229   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
4230   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
4231   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
4232   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
4233   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
4234   return(pixel);
4235 }
4236
4237 MagickExport Image *ReduceNoiseImage(const Image *image,const double radius,
4238   ExceptionInfo *exception)
4239 {
4240 #define ReduceNoiseImageTag  "ReduceNoise/Image"
4241
4242   CacheView
4243     *image_view,
4244     *noise_view;
4245
4246   Image
4247     *noise_image;
4248
4249   MagickBooleanType
4250     status;
4251
4252   MagickOffsetType
4253     progress;
4254
4255   MedianPixelList
4256     **restrict pixel_list;
4257
4258   size_t
4259     width;
4260
4261   ssize_t
4262     y;
4263
4264   /*
4265     Initialize noise image attributes.
4266   */
4267   assert(image != (Image *) NULL);
4268   assert(image->signature == MagickSignature);
4269   if (image->debug != MagickFalse)
4270     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4271   assert(exception != (ExceptionInfo *) NULL);
4272   assert(exception->signature == MagickSignature);
4273   width=GetOptimalKernelWidth2D(radius,0.5);
4274   if ((image->columns < width) || (image->rows < width))
4275     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
4276   noise_image=CloneImage(image,image->columns,image->rows,MagickTrue,
4277     exception);
4278   if (noise_image == (Image *) NULL)
4279     return((Image *) NULL);
4280   if (SetImageStorageClass(noise_image,DirectClass) == MagickFalse)
4281     {
4282       InheritException(exception,&noise_image->exception);
4283       noise_image=DestroyImage(noise_image);
4284       return((Image *) NULL);
4285     }
4286   pixel_list=AcquireMedianPixelListThreadSet(width);
4287   if (pixel_list == (MedianPixelList **) NULL)
4288     {
4289       noise_image=DestroyImage(noise_image);
4290       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4291     }
4292   /*
4293     Reduce noise image.
4294   */
4295   status=MagickTrue;
4296   progress=0;
4297   image_view=AcquireCacheView(image);
4298   noise_view=AcquireCacheView(noise_image);
4299 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4300   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4301 #endif
4302   for (y=0; y < (ssize_t) noise_image->rows; y++)
4303   {
4304     const int
4305       id = GetOpenMPThreadId();
4306
4307     register const IndexPacket
4308       *restrict indexes;
4309
4310     register const PixelPacket
4311       *restrict p;
4312
4313     register IndexPacket
4314       *restrict noise_indexes;
4315
4316     register PixelPacket
4317       *restrict q;
4318
4319     register ssize_t
4320       x;
4321
4322     if (status == MagickFalse)
4323       continue;
4324     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
4325       (width/2L),image->columns+width,width,exception);
4326     q=QueueCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
4327       exception);
4328     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4329       {
4330         status=MagickFalse;
4331         continue;
4332       }
4333     indexes=GetCacheViewVirtualIndexQueue(image_view);
4334     noise_indexes=GetCacheViewAuthenticIndexQueue(noise_view);
4335     for (x=0; x < (ssize_t) noise_image->columns; x++)
4336     {
4337       MagickPixelPacket
4338         pixel;
4339
4340       register const PixelPacket
4341         *restrict r;
4342
4343       register const IndexPacket
4344         *restrict s;
4345
4346       register ssize_t
4347         u,
4348         v;
4349
4350       r=p;
4351       s=indexes+x;
4352       ResetMedianPixelList(pixel_list[id]);
4353       for (v=0; v < (ssize_t) width; v++)
4354       {
4355         for (u=0; u < (ssize_t) width; u++)
4356           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
4357         r+=image->columns+width;
4358         s+=image->columns+width;
4359       }
4360       pixel=GetNonpeakMedianPixelList(pixel_list[id]);
4361       SetPixelPacket(noise_image,&pixel,q,noise_indexes+x);
4362       p++;
4363       q++;
4364     }
4365     if (SyncCacheViewAuthenticPixels(noise_view,exception) == MagickFalse)
4366       status=MagickFalse;
4367     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4368       {
4369         MagickBooleanType
4370           proceed;
4371
4372 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4373   #pragma omp critical (MagickCore_ReduceNoiseImage)
4374 #endif
4375         proceed=SetImageProgress(image,ReduceNoiseImageTag,progress++,
4376           image->rows);
4377         if (proceed == MagickFalse)
4378           status=MagickFalse;
4379       }
4380   }
4381   noise_view=DestroyCacheView(noise_view);
4382   image_view=DestroyCacheView(image_view);
4383   pixel_list=DestroyMedianPixelListThreadSet(pixel_list);
4384   return(noise_image);
4385 }
4386 \f
4387 /*
4388 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4389 %                                                                             %
4390 %                                                                             %
4391 %                                                                             %
4392 %     S e l e c t i v e B l u r I m a g e                                     %
4393 %                                                                             %
4394 %                                                                             %
4395 %                                                                             %
4396 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4397 %
4398 %  SelectiveBlurImage() selectively blur pixels within a contrast threshold.
4399 %  It is similar to the unsharpen mask that sharpens everything with contrast
4400 %  above a certain threshold.
4401 %
4402 %  The format of the SelectiveBlurImage method is:
4403 %
4404 %      Image *SelectiveBlurImage(const Image *image,const double radius,
4405 %        const double sigma,const double threshold,ExceptionInfo *exception)
4406 %      Image *SelectiveBlurImageChannel(const Image *image,
4407 %        const ChannelType channel,const double radius,const double sigma,
4408 %        const double threshold,ExceptionInfo *exception)
4409 %
4410 %  A description of each parameter follows:
4411 %
4412 %    o image: the image.
4413 %
4414 %    o channel: the channel type.
4415 %
4416 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4417 %      pixel.
4418 %
4419 %    o sigma: the standard deviation of the Gaussian, in pixels.
4420 %
4421 %    o threshold: only pixels within this contrast threshold are included
4422 %      in the blur operation.
4423 %
4424 %    o exception: return any errors or warnings in this structure.
4425 %
4426 */
4427
4428 static inline MagickBooleanType SelectiveContrast(const PixelPacket *p,
4429   const PixelPacket *q,const double threshold)
4430 {
4431   if (fabs(PixelIntensity(p)-PixelIntensity(q)) < threshold)
4432     return(MagickTrue);
4433   return(MagickFalse);
4434 }
4435
4436 MagickExport Image *SelectiveBlurImage(const Image *image,const double radius,
4437   const double sigma,const double threshold,ExceptionInfo *exception)
4438 {
4439   Image
4440     *blur_image;
4441
4442   blur_image=SelectiveBlurImageChannel(image,DefaultChannels,radius,sigma,
4443     threshold,exception);
4444   return(blur_image);
4445 }
4446
4447 MagickExport Image *SelectiveBlurImageChannel(const Image *image,
4448   const ChannelType channel,const double radius,const double sigma,
4449   const double threshold,ExceptionInfo *exception)
4450 {
4451 #define SelectiveBlurImageTag  "SelectiveBlur/Image"
4452
4453   CacheView
4454     *blur_view,
4455     *image_view;
4456
4457   double
4458     *kernel;
4459
4460   Image
4461     *blur_image;
4462
4463   MagickBooleanType
4464     status;
4465
4466   MagickOffsetType
4467     progress;
4468
4469   MagickPixelPacket
4470     bias;
4471
4472   register ssize_t
4473     i;
4474
4475   size_t
4476     width;
4477
4478   ssize_t
4479     j,
4480     u,
4481     v,
4482     y;
4483
4484   /*
4485     Initialize blur image attributes.
4486   */
4487   assert(image != (Image *) NULL);
4488   assert(image->signature == MagickSignature);
4489   if (image->debug != MagickFalse)
4490     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4491   assert(exception != (ExceptionInfo *) NULL);
4492   assert(exception->signature == MagickSignature);
4493   width=GetOptimalKernelWidth1D(radius,sigma);
4494   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
4495   if (kernel == (double *) NULL)
4496     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4497   j=(ssize_t) width/2;
4498   i=0;
4499   for (v=(-j); v <= j; v++)
4500   {
4501     for (u=(-j); u <= j; u++)
4502       kernel[i++]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
4503         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
4504   }
4505   if (image->debug != MagickFalse)
4506     {
4507       char
4508         format[MaxTextExtent],
4509         *message;
4510
4511       register const double
4512         *k;
4513
4514       ssize_t
4515         u,
4516         v;
4517
4518       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
4519         "  SelectiveBlurImage with %.20gx%.20g kernel:",(double) width,(double)
4520         width);
4521       message=AcquireString("");
4522       k=kernel;
4523       for (v=0; v < (ssize_t) width; v++)
4524       {
4525         *message='\0';
4526         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) v);
4527         (void) ConcatenateString(&message,format);
4528         for (u=0; u < (ssize_t) width; u++)
4529         {
4530           (void) FormatMagickString(format,MaxTextExtent,"%+f ",*k++);
4531           (void) ConcatenateString(&message,format);
4532         }
4533         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
4534       }
4535       message=DestroyString(message);
4536     }
4537   blur_image=CloneImage(image,0,0,MagickTrue,exception);
4538   if (blur_image == (Image *) NULL)
4539     return((Image *) NULL);
4540   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
4541     {
4542       InheritException(exception,&blur_image->exception);
4543       blur_image=DestroyImage(blur_image);
4544       return((Image *) NULL);
4545     }
4546   /*
4547     Threshold blur image.
4548   */
4549   status=MagickTrue;
4550   progress=0;
4551   GetMagickPixelPacket(image,&bias);
4552   SetMagickPixelPacketBias(image,&bias);
4553   image_view=AcquireCacheView(image);
4554   blur_view=AcquireCacheView(blur_image);
4555 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4556   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4557 #endif
4558   for (y=0; y < (ssize_t) image->rows; y++)
4559   {
4560     MagickBooleanType
4561       sync;
4562
4563     MagickRealType
4564       gamma;
4565
4566     register const IndexPacket
4567       *restrict indexes;
4568
4569     register const PixelPacket
4570       *restrict p;
4571
4572     register IndexPacket
4573       *restrict blur_indexes;
4574
4575     register PixelPacket
4576       *restrict q;
4577
4578     register ssize_t
4579       x;
4580
4581     if (status == MagickFalse)
4582       continue;
4583     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
4584       (width/2L),image->columns+width,width,exception);
4585     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
4586       exception);
4587     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4588       {
4589         status=MagickFalse;
4590         continue;
4591       }
4592     indexes=GetCacheViewVirtualIndexQueue(image_view);
4593     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
4594     for (x=0; x < (ssize_t) image->columns; x++)
4595     {
4596       MagickPixelPacket
4597         pixel;
4598
4599       register const double
4600         *restrict k;
4601
4602       register ssize_t
4603         u;
4604
4605       ssize_t
4606         j,
4607         v;
4608
4609       pixel=bias;
4610       k=kernel;
4611       gamma=0.0;
4612       j=0;
4613       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
4614         {
4615           for (v=0; v < (ssize_t) width; v++)
4616           {
4617             for (u=0; u < (ssize_t) width; u++)
4618             {
4619               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4620                 {
4621                   pixel.red+=(*k)*(p+u+j)->red;
4622                   pixel.green+=(*k)*(p+u+j)->green;
4623                   pixel.blue+=(*k)*(p+u+j)->blue;
4624                   gamma+=(*k);
4625                   k++;
4626                 }
4627             }
4628             j+=(ssize_t) (image->columns+width);
4629           }
4630           if (gamma != 0.0)
4631             {
4632               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4633               if ((channel & RedChannel) != 0)
4634                 q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
4635               if ((channel & GreenChannel) != 0)
4636                 q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
4637               if ((channel & BlueChannel) != 0)
4638                 q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
4639             }
4640           if ((channel & OpacityChannel) != 0)
4641             {
4642               gamma=0.0;
4643               j=0;
4644               for (v=0; v < (ssize_t) width; v++)
4645               {
4646                 for (u=0; u < (ssize_t) width; u++)
4647                 {
4648                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4649                     {
4650                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4651                       gamma+=(*k);
4652                       k++;
4653                     }
4654                 }
4655                 j+=(ssize_t) (image->columns+width);
4656               }
4657               if (gamma != 0.0)
4658                 {
4659                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4660                     gamma);
4661                   SetOpacityPixelComponent(q,ClampToQuantum(gamma*
4662                     GetOpacityPixelComponent(&pixel)));
4663                 }
4664             }
4665           if (((channel & IndexChannel) != 0) &&
4666               (image->colorspace == CMYKColorspace))
4667             {
4668               gamma=0.0;
4669               j=0;
4670               for (v=0; v < (ssize_t) width; v++)
4671               {
4672                 for (u=0; u < (ssize_t) width; u++)
4673                 {
4674                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4675                     {
4676                       pixel.index+=(*k)*indexes[x+u+j];
4677                       gamma+=(*k);
4678                       k++;
4679                     }
4680                 }
4681                 j+=(ssize_t) (image->columns+width);
4682               }
4683               if (gamma != 0.0)
4684                 {
4685                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4686                     gamma);
4687                   blur_indexes[x]=ClampToQuantum(gamma*
4688                     GetIndexPixelComponent(&pixel));
4689                 }
4690             }
4691         }
4692       else
4693         {
4694           MagickRealType
4695             alpha;
4696
4697           for (v=0; v < (ssize_t) width; v++)
4698           {
4699             for (u=0; u < (ssize_t) width; u++)
4700             {
4701               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4702                 {
4703                   alpha=(MagickRealType) (QuantumScale*
4704                     GetAlphaPixelComponent(p+u+j));
4705                   pixel.red+=(*k)*alpha*(p+u+j)->red;
4706                   pixel.green+=(*k)*alpha*(p+u+j)->green;
4707                   pixel.blue+=(*k)*alpha*(p+u+j)->blue;
4708                   pixel.opacity+=(*k)*(p+u+j)->opacity;
4709                   gamma+=(*k)*alpha;
4710                   k++;
4711                 }
4712             }
4713             j+=(ssize_t) (image->columns+width);
4714           }
4715           if (gamma != 0.0)
4716             {
4717               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4718               if ((channel & RedChannel) != 0)
4719                 q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
4720               if ((channel & GreenChannel) != 0)
4721                 q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
4722               if ((channel & BlueChannel) != 0)
4723                 q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
4724             }
4725           if ((channel & OpacityChannel) != 0)
4726             {
4727               gamma=0.0;
4728               j=0;
4729               for (v=0; v < (ssize_t) width; v++)
4730               {
4731                 for (u=0; u < (ssize_t) width; u++)
4732                 {
4733                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4734                     {
4735                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4736                       gamma+=(*k);
4737                       k++;
4738                     }
4739                 }
4740                 j+=(ssize_t) (image->columns+width);
4741               }
4742               if (gamma != 0.0)
4743                 {
4744                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4745                     gamma);
4746                   SetOpacityPixelComponent(q,
4747                     ClampOpacityPixelComponent(&pixel));
4748                 }
4749             }
4750           if (((channel & IndexChannel) != 0) &&
4751               (image->colorspace == CMYKColorspace))
4752             {
4753               gamma=0.0;
4754               j=0;
4755               for (v=0; v < (ssize_t) width; v++)
4756               {
4757                 for (u=0; u < (ssize_t) width; u++)
4758                 {
4759                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4760                     {
4761                       alpha=(MagickRealType) (QuantumScale*
4762                         GetAlphaPixelComponent(p+u+j));
4763                       pixel.index+=(*k)*alpha*indexes[x+u+j];
4764                       gamma+=(*k);
4765                       k++;
4766                     }
4767                 }
4768                 j+=(ssize_t) (image->columns+width);
4769               }
4770               if (gamma != 0.0)
4771                 {
4772                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4773                     gamma);
4774                   blur_indexes[x]=ClampToQuantum(gamma*
4775                     GetIndexPixelComponent(&pixel));
4776                 }
4777             }
4778         }
4779       p++;
4780       q++;
4781     }
4782     sync=SyncCacheViewAuthenticPixels(blur_view,exception);
4783     if (sync == MagickFalse)
4784       status=MagickFalse;
4785     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4786       {
4787         MagickBooleanType
4788           proceed;
4789
4790 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4791   #pragma omp critical (MagickCore_SelectiveBlurImageChannel)
4792 #endif
4793         proceed=SetImageProgress(image,SelectiveBlurImageTag,progress++,
4794           image->rows);
4795         if (proceed == MagickFalse)
4796           status=MagickFalse;
4797       }
4798   }
4799   blur_image->type=image->type;
4800   blur_view=DestroyCacheView(blur_view);
4801   image_view=DestroyCacheView(image_view);
4802   kernel=(double *) RelinquishMagickMemory(kernel);
4803   if (status == MagickFalse)
4804     blur_image=DestroyImage(blur_image);
4805   return(blur_image);
4806 }
4807 \f
4808 /*
4809 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4810 %                                                                             %
4811 %                                                                             %
4812 %                                                                             %
4813 %     S h a d e I m a g e                                                     %
4814 %                                                                             %
4815 %                                                                             %
4816 %                                                                             %
4817 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4818 %
4819 %  ShadeImage() shines a distant light on an image to create a
4820 %  three-dimensional effect. You control the positioning of the light with
4821 %  azimuth and elevation; azimuth is measured in degrees off the x axis
4822 %  and elevation is measured in pixels above the Z axis.
4823 %
4824 %  The format of the ShadeImage method is:
4825 %
4826 %      Image *ShadeImage(const Image *image,const MagickBooleanType gray,
4827 %        const double azimuth,const double elevation,ExceptionInfo *exception)
4828 %
4829 %  A description of each parameter follows:
4830 %
4831 %    o image: the image.
4832 %
4833 %    o gray: A value other than zero shades the intensity of each pixel.
4834 %
4835 %    o azimuth, elevation:  Define the light source direction.
4836 %
4837 %    o exception: return any errors or warnings in this structure.
4838 %
4839 */
4840 MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
4841   const double azimuth,const double elevation,ExceptionInfo *exception)
4842 {
4843 #define ShadeImageTag  "Shade/Image"
4844
4845   CacheView
4846     *image_view,
4847     *shade_view;
4848
4849   Image
4850     *shade_image;
4851
4852   MagickBooleanType
4853     status;
4854
4855   MagickOffsetType
4856     progress;
4857
4858   PrimaryInfo
4859     light;
4860
4861   ssize_t
4862     y;
4863
4864   /*
4865     Initialize shaded image attributes.
4866   */
4867   assert(image != (const Image *) NULL);
4868   assert(image->signature == MagickSignature);
4869   if (image->debug != MagickFalse)
4870     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4871   assert(exception != (ExceptionInfo *) NULL);
4872   assert(exception->signature == MagickSignature);
4873   shade_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
4874   if (shade_image == (Image *) NULL)
4875     return((Image *) NULL);
4876   if (SetImageStorageClass(shade_image,DirectClass) == MagickFalse)
4877     {
4878       InheritException(exception,&shade_image->exception);
4879       shade_image=DestroyImage(shade_image);
4880       return((Image *) NULL);
4881     }
4882   /*
4883     Compute the light vector.
4884   */
4885   light.x=(double) QuantumRange*cos(DegreesToRadians(azimuth))*
4886     cos(DegreesToRadians(elevation));
4887   light.y=(double) QuantumRange*sin(DegreesToRadians(azimuth))*
4888     cos(DegreesToRadians(elevation));
4889   light.z=(double) QuantumRange*sin(DegreesToRadians(elevation));
4890   /*
4891     Shade image.
4892   */
4893   status=MagickTrue;
4894   progress=0;
4895   image_view=AcquireCacheView(image);
4896   shade_view=AcquireCacheView(shade_image);
4897 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4898   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4899 #endif
4900   for (y=0; y < (ssize_t) image->rows; y++)
4901   {
4902     MagickRealType
4903       distance,
4904       normal_distance,
4905       shade;
4906
4907     PrimaryInfo
4908       normal;
4909
4910     register const PixelPacket
4911       *restrict p,
4912       *restrict s0,
4913       *restrict s1,
4914       *restrict s2;
4915
4916     register PixelPacket
4917       *restrict q;
4918
4919     register ssize_t
4920       x;
4921
4922     if (status == MagickFalse)
4923       continue;
4924     p=GetCacheViewVirtualPixels(image_view,-1,y-1,image->columns+2,3,exception);
4925     q=QueueCacheViewAuthenticPixels(shade_view,0,y,shade_image->columns,1,
4926       exception);
4927     if ((p == (PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4928       {
4929         status=MagickFalse;
4930         continue;
4931       }
4932     /*
4933       Shade this row of pixels.
4934     */
4935     normal.z=2.0*(double) QuantumRange;  /* constant Z of surface normal */
4936     s0=p+1;
4937     s1=s0+image->columns+2;
4938     s2=s1+image->columns+2;
4939     for (x=0; x < (ssize_t) image->columns; x++)
4940     {
4941       /*
4942         Determine the surface normal and compute shading.
4943       */
4944       normal.x=(double) (PixelIntensity(s0-1)+PixelIntensity(s1-1)+
4945         PixelIntensity(s2-1)-PixelIntensity(s0+1)-PixelIntensity(s1+1)-
4946         PixelIntensity(s2+1));
4947       normal.y=(double) (PixelIntensity(s2-1)+PixelIntensity(s2)+
4948         PixelIntensity(s2+1)-PixelIntensity(s0-1)-PixelIntensity(s0)-
4949         PixelIntensity(s0+1));
4950       if ((normal.x == 0.0) && (normal.y == 0.0))
4951         shade=light.z;
4952       else
4953         {
4954           shade=0.0;
4955           distance=normal.x*light.x+normal.y*light.y+normal.z*light.z;
4956           if (distance > MagickEpsilon)
4957             {
4958               normal_distance=
4959                 normal.x*normal.x+normal.y*normal.y+normal.z*normal.z;
4960               if (normal_distance > (MagickEpsilon*MagickEpsilon))
4961                 shade=distance/sqrt((double) normal_distance);
4962             }
4963         }
4964       if (gray != MagickFalse)
4965         {
4966           q->red=(Quantum) shade;
4967           q->green=(Quantum) shade;
4968           q->blue=(Quantum) shade;
4969         }
4970       else
4971         {
4972           q->red=ClampToQuantum(QuantumScale*shade*s1->red);
4973           q->green=ClampToQuantum(QuantumScale*shade*s1->green);
4974           q->blue=ClampToQuantum(QuantumScale*shade*s1->blue);
4975         }
4976       q->opacity=s1->opacity;
4977       s0++;
4978       s1++;
4979       s2++;
4980       q++;
4981     }
4982     if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
4983       status=MagickFalse;
4984     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4985       {
4986         MagickBooleanType
4987           proceed;
4988
4989 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4990   #pragma omp critical (MagickCore_ShadeImage)
4991 #endif
4992         proceed=SetImageProgress(image,ShadeImageTag,progress++,image->rows);
4993         if (proceed == MagickFalse)
4994           status=MagickFalse;
4995       }
4996   }
4997   shade_view=DestroyCacheView(shade_view);
4998   image_view=DestroyCacheView(image_view);
4999   if (status == MagickFalse)
5000     shade_image=DestroyImage(shade_image);
5001   return(shade_image);
5002 }
5003 \f
5004 /*
5005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5006 %                                                                             %
5007 %                                                                             %
5008 %                                                                             %
5009 %     S h a r p e n I m a g e                                                 %
5010 %                                                                             %
5011 %                                                                             %
5012 %                                                                             %
5013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5014 %
5015 %  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
5016 %  operator of the given radius and standard deviation (sigma).  For
5017 %  reasonable results, radius should be larger than sigma.  Use a radius of 0
5018 %  and SharpenImage() selects a suitable radius for you.
5019 %
5020 %  Using a separable kernel would be faster, but the negative weights cancel
5021 %  out on the corners of the kernel producing often undesirable ringing in the
5022 %  filtered result; this can be avoided by using a 2D gaussian shaped image
5023 %  sharpening kernel instead.
5024 %
5025 %  The format of the SharpenImage method is:
5026 %
5027 %    Image *SharpenImage(const Image *image,const double radius,
5028 %      const double sigma,ExceptionInfo *exception)
5029 %    Image *SharpenImageChannel(const Image *image,const ChannelType channel,
5030 %      const double radius,const double sigma,ExceptionInfo *exception)
5031 %
5032 %  A description of each parameter follows:
5033 %
5034 %    o image: the image.
5035 %
5036 %    o channel: the channel type.
5037 %
5038 %    o radius: the radius of the Gaussian, in pixels, not counting the center
5039 %      pixel.
5040 %
5041 %    o sigma: the standard deviation of the Laplacian, in pixels.
5042 %
5043 %    o exception: return any errors or warnings in this structure.
5044 %
5045 */
5046
5047 MagickExport Image *SharpenImage(const Image *image,const double radius,
5048   const double sigma,ExceptionInfo *exception)
5049 {
5050   Image
5051     *sharp_image;
5052
5053   sharp_image=SharpenImageChannel(image,DefaultChannels,radius,sigma,exception);
5054   return(sharp_image);
5055 }
5056
5057 MagickExport Image *SharpenImageChannel(const Image *image,
5058   const ChannelType channel,const double radius,const double sigma,
5059   ExceptionInfo *exception)
5060 {
5061   double
5062     *kernel,
5063     normalize;
5064
5065   Image
5066     *sharp_image;
5067
5068   register ssize_t
5069     i;
5070
5071   size_t
5072     width;
5073
5074   ssize_t
5075     j,
5076     u,
5077     v;
5078
5079   assert(image != (const Image *) NULL);
5080   assert(image->signature == MagickSignature);
5081   if (image->debug != MagickFalse)
5082     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5083   assert(exception != (ExceptionInfo *) NULL);
5084   assert(exception->signature == MagickSignature);
5085   width=GetOptimalKernelWidth2D(radius,sigma);
5086   kernel=(double *) AcquireQuantumMemory((size_t) width*width,sizeof(*kernel));
5087   if (kernel == (double *) NULL)
5088     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
5089   normalize=0.0;
5090   j=(ssize_t) width/2;
5091   i=0;
5092   for (v=(-j); v <= j; v++)
5093   {
5094     for (u=(-j); u <= j; u++)
5095     {
5096       kernel[i]=(double) (-exp(-((double) u*u+v*v)/(2.0*MagickSigma*
5097         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
5098       normalize+=kernel[i];
5099       i++;
5100     }
5101   }
5102   kernel[i/2]=(double) ((-2.0)*normalize);
5103   sharp_image=ConvolveImageChannel(image,channel,width,kernel,exception);
5104   kernel=(double *) RelinquishMagickMemory(kernel);
5105   return(sharp_image);
5106 }
5107 \f
5108 /*
5109 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5110 %                                                                             %
5111 %                                                                             %
5112 %                                                                             %
5113 %     S p r e a d I m a g e                                                   %
5114 %                                                                             %
5115 %                                                                             %
5116 %                                                                             %
5117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5118 %
5119 %  SpreadImage() is a special effects method that randomly displaces each
5120 %  pixel in a block defined by the radius parameter.
5121 %
5122 %  The format of the SpreadImage method is:
5123 %
5124 %      Image *SpreadImage(const Image *image,const double radius,
5125 %        ExceptionInfo *exception)
5126 %
5127 %  A description of each parameter follows:
5128 %
5129 %    o image: the image.
5130 %
5131 %    o radius:  Choose a random pixel in a neighborhood of this extent.
5132 %
5133 %    o exception: return any errors or warnings in this structure.
5134 %
5135 */
5136 MagickExport Image *SpreadImage(const Image *image,const double radius,
5137   ExceptionInfo *exception)
5138 {
5139 #define SpreadImageTag  "Spread/Image"
5140
5141   CacheView
5142     *image_view;
5143
5144   Image
5145     *spread_image;
5146
5147   MagickBooleanType
5148     status;
5149
5150   MagickOffsetType
5151     progress;
5152
5153   MagickPixelPacket
5154     bias;
5155
5156   RandomInfo
5157     **restrict random_info;
5158
5159   ResampleFilter
5160     **restrict resample_filter;
5161
5162   size_t
5163     width;
5164
5165   ssize_t
5166     y;
5167
5168   /*
5169     Initialize spread image attributes.
5170   */
5171   assert(image != (Image *) NULL);
5172   assert(image->signature == MagickSignature);
5173   if (image->debug != MagickFalse)
5174     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5175   assert(exception != (ExceptionInfo *) NULL);
5176   assert(exception->signature == MagickSignature);
5177   spread_image=CloneImage(image,image->columns,image->rows,MagickTrue,
5178     exception);
5179   if (spread_image == (Image *) NULL)
5180     return((Image *) NULL);
5181   if (SetImageStorageClass(spread_image,DirectClass) == MagickFalse)
5182     {
5183       InheritException(exception,&spread_image->exception);
5184       spread_image=DestroyImage(spread_image);
5185       return((Image *) NULL);
5186     }
5187   /*
5188     Spread image.
5189   */
5190   status=MagickTrue;
5191   progress=0;
5192   GetMagickPixelPacket(spread_image,&bias);
5193   width=GetOptimalKernelWidth1D(radius,0.5);
5194   resample_filter=AcquireResampleFilterThreadSet(image,
5195     UndefinedVirtualPixelMethod,MagickTrue,exception);
5196   random_info=AcquireRandomInfoThreadSet();
5197   image_view=AcquireCacheView(spread_image);
5198 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5199   #pragma omp parallel for schedule(dynamic,4) shared(progress,status) omp_throttle(1)
5200 #endif
5201   for (y=0; y < (ssize_t) spread_image->rows; y++)
5202   {
5203     const int
5204       id = GetOpenMPThreadId();
5205
5206     MagickPixelPacket
5207       pixel;
5208
5209     register IndexPacket
5210       *restrict indexes;
5211
5212     register PixelPacket
5213       *restrict q;
5214
5215     register ssize_t
5216       x;
5217
5218     if (status == MagickFalse)
5219       continue;
5220     q=QueueCacheViewAuthenticPixels(image_view,0,y,spread_image->columns,1,
5221       exception);
5222     if (q == (PixelPacket *) NULL)
5223       {
5224         status=MagickFalse;
5225         continue;
5226       }
5227     indexes=GetCacheViewAuthenticIndexQueue(image_view);
5228     pixel=bias;
5229     for (x=0; x < (ssize_t) spread_image->columns; x++)
5230     {
5231       (void) ResamplePixelColor(resample_filter[id],(double) x+width*
5232         (GetPseudoRandomValue(random_info[id])-0.5),(double) y+width*
5233         (GetPseudoRandomValue(random_info[id])-0.5),&pixel);
5234       SetPixelPacket(spread_image,&pixel,q,indexes+x);
5235       q++;
5236     }
5237     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
5238       status=MagickFalse;
5239     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5240       {
5241         MagickBooleanType
5242           proceed;
5243
5244 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5245   #pragma omp critical (MagickCore_SpreadImage)
5246 #endif
5247         proceed=SetImageProgress(image,SpreadImageTag,progress++,image->rows);
5248         if (proceed == MagickFalse)
5249           status=MagickFalse;
5250       }
5251   }
5252   image_view=DestroyCacheView(image_view);
5253   random_info=DestroyRandomInfoThreadSet(random_info);
5254   resample_filter=DestroyResampleFilterThreadSet(resample_filter);
5255   return(spread_image);
5256 }
5257 \f
5258 /*
5259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5260 %                                                                             %
5261 %                                                                             %
5262 %                                                                             %
5263 %     U n s h a r p M a s k I m a g e                                         %
5264 %                                                                             %
5265 %                                                                             %
5266 %                                                                             %
5267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5268 %
5269 %  UnsharpMaskImage() sharpens one or more image channels.  We convolve the
5270 %  image with a Gaussian operator of the given radius and standard deviation
5271 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
5272 %  radius of 0 and UnsharpMaskImage() selects a suitable radius for you.
5273 %
5274 %  The format of the UnsharpMaskImage method is:
5275 %
5276 %    Image *UnsharpMaskImage(const Image *image,const double radius,
5277 %      const double sigma,const double amount,const double threshold,
5278 %      ExceptionInfo *exception)
5279 %    Image *UnsharpMaskImageChannel(const Image *image,
5280 %      const ChannelType channel,const double radius,const double sigma,
5281 %      const double amount,const double threshold,ExceptionInfo *exception)
5282 %
5283 %  A description of each parameter follows:
5284 %
5285 %    o image: the image.
5286 %
5287 %    o channel: the channel type.
5288 %
5289 %    o radius: the radius of the Gaussian, in pixels, not counting the center
5290 %      pixel.
5291 %
5292 %    o sigma: the standard deviation of the Gaussian, in pixels.
5293 %
5294 %    o amount: the percentage of the difference between the original and the
5295 %      blur image that is added back into the original.
5296 %
5297 %    o threshold: the threshold in pixels needed to apply the diffence amount.
5298 %
5299 %    o exception: return any errors or warnings in this structure.
5300 %
5301 */
5302
5303 MagickExport Image *UnsharpMaskImage(const Image *image,const double radius,
5304   const double sigma,const double amount,const double threshold,
5305   ExceptionInfo *exception)
5306 {
5307   Image
5308     *sharp_image;
5309
5310   sharp_image=UnsharpMaskImageChannel(image,DefaultChannels,radius,sigma,amount,
5311     threshold,exception);
5312   return(sharp_image);
5313 }
5314
5315 MagickExport Image *UnsharpMaskImageChannel(const Image *image,
5316   const ChannelType channel,const double radius,const double sigma,
5317   const double amount,const double threshold,ExceptionInfo *exception)
5318 {
5319 #define SharpenImageTag  "Sharpen/Image"
5320
5321   CacheView
5322     *image_view,
5323     *unsharp_view;
5324
5325   Image
5326     *unsharp_image;
5327
5328   MagickBooleanType
5329     status;
5330
5331   MagickOffsetType
5332     progress;
5333
5334   MagickPixelPacket
5335     bias;
5336
5337   MagickRealType
5338     quantum_threshold;
5339
5340   ssize_t
5341     y;
5342
5343   assert(image != (const Image *) NULL);
5344   assert(image->signature == MagickSignature);
5345   if (image->debug != MagickFalse)
5346     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5347   assert(exception != (ExceptionInfo *) NULL);
5348   unsharp_image=BlurImageChannel(image,channel,radius,sigma,exception);
5349   if (unsharp_image == (Image *) NULL)
5350     return((Image *) NULL);
5351   quantum_threshold=(MagickRealType) QuantumRange*threshold;
5352   /*
5353     Unsharp-mask image.
5354   */
5355   status=MagickTrue;
5356   progress=0;
5357   GetMagickPixelPacket(image,&bias);
5358   image_view=AcquireCacheView(image);
5359   unsharp_view=AcquireCacheView(unsharp_image);
5360 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5361   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
5362 #endif
5363   for (y=0; y < (ssize_t) image->rows; y++)
5364   {
5365     MagickPixelPacket
5366       pixel;
5367
5368     register const IndexPacket
5369       *restrict indexes;
5370
5371     register const PixelPacket
5372       *restrict p;
5373
5374     register IndexPacket
5375       *restrict unsharp_indexes;
5376
5377     register PixelPacket
5378       *restrict q;
5379
5380     register ssize_t
5381       x;
5382
5383     if (status == MagickFalse)
5384       continue;
5385     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
5386     q=GetCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
5387       exception);
5388     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
5389       {
5390         status=MagickFalse;
5391         continue;
5392       }
5393     indexes=GetCacheViewVirtualIndexQueue(image_view);
5394     unsharp_indexes=GetCacheViewAuthenticIndexQueue(unsharp_view);
5395     pixel=bias;
5396     for (x=0; x < (ssize_t) image->columns; x++)
5397     {
5398       if ((channel & RedChannel) != 0)
5399         {
5400           pixel.red=p->red-(MagickRealType) q->red;
5401           if (fabs(2.0*pixel.red) < quantum_threshold)
5402             pixel.red=(MagickRealType) GetRedPixelComponent(p);
5403           else
5404             pixel.red=(MagickRealType) p->red+(pixel.red*amount);
5405           SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
5406         }
5407       if ((channel & GreenChannel) != 0)
5408         {
5409           pixel.green=p->green-(MagickRealType) q->green;
5410           if (fabs(2.0*pixel.green) < quantum_threshold)
5411             pixel.green=(MagickRealType) GetGreenPixelComponent(p);
5412           else
5413             pixel.green=(MagickRealType) p->green+(pixel.green*amount);
5414           SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
5415         }
5416       if ((channel & BlueChannel) != 0)
5417         {
5418           pixel.blue=p->blue-(MagickRealType) q->blue;
5419           if (fabs(2.0*pixel.blue) < quantum_threshold)
5420             pixel.blue=(MagickRealType) GetBluePixelComponent(p);
5421           else
5422             pixel.blue=(MagickRealType) p->blue+(pixel.blue*amount);
5423           SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
5424         }
5425       if ((channel & OpacityChannel) != 0)
5426         {
5427           pixel.opacity=p->opacity-(MagickRealType) q->opacity;
5428           if (fabs(2.0*pixel.opacity) < quantum_threshold)
5429             pixel.opacity=(MagickRealType) GetOpacityPixelComponent(p);
5430           else
5431             pixel.opacity=p->opacity+(pixel.opacity*amount);
5432           SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
5433         }
5434       if (((channel & IndexChannel) != 0) &&
5435           (image->colorspace == CMYKColorspace))
5436         {
5437           pixel.index=unsharp_indexes[x]-(MagickRealType) indexes[x];
5438           if (fabs(2.0*pixel.index) < quantum_threshold)
5439             pixel.index=(MagickRealType) unsharp_indexes[x];
5440           else
5441             pixel.index=(MagickRealType) unsharp_indexes[x]+(pixel.index*
5442               amount);
5443           unsharp_indexes[x]=ClampToQuantum(pixel.index);
5444         }
5445       p++;
5446       q++;
5447     }
5448     if (SyncCacheViewAuthenticPixels(unsharp_view,exception) == MagickFalse)
5449       status=MagickFalse;
5450     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5451       {
5452         MagickBooleanType
5453           proceed;
5454
5455 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5456   #pragma omp critical (MagickCore_UnsharpMaskImageChannel)
5457 #endif
5458         proceed=SetImageProgress(image,SharpenImageTag,progress++,image->rows);
5459         if (proceed == MagickFalse)
5460           status=MagickFalse;
5461       }
5462   }
5463   unsharp_image->type=image->type;
5464   unsharp_view=DestroyCacheView(unsharp_view);
5465   image_view=DestroyCacheView(image_view);
5466   if (status == MagickFalse)
5467     unsharp_image=DestroyImage(unsharp_image);
5468   return(unsharp_image);
5469 }