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-2010 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]=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
 241           (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]=(-exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
 558           (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]=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 **) RelinquishAlignedMemory(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 **) AcquireAlignedMemory(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,4) shared(status)
1865 #endif
1866   for (channel=0; channel <= 3; channel++)
1867   {
1868     int
1869       id;
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     id=GetOpenMPThreadId();
1886     pixel=pixels[id];
1887     (void) ResetMagickMemory(pixel,0,length*sizeof(*pixel));
1888     buffer=buffers[id];
1889     j=(ssize_t) image->columns+2;
1890     for (y=0; y < (ssize_t) image->rows; y++)
1891     {
1892       register const PixelPacket
1893         *restrict p;
1894
1895       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1896       if (p == (const PixelPacket *) NULL)
1897         break;
1898       j++;
1899       for (x=0; x < (ssize_t) image->columns; x++)
1900       {
1901         switch (channel)
1902         {
1903           case 0: pixel[j]=GetRedPixelComponent(p); break;
1904           case 1: pixel[j]=GetGreenPixelComponent(p); break;
1905           case 2: pixel[j]=GetBluePixelComponent(p); break;
1906           case 3: pixel[j]=GetOpacityPixelComponent(p); break;
1907           default: break;
1908         }
1909         p++;
1910         j++;
1911       }
1912       j++;
1913     }
1914     (void) ResetMagickMemory(buffer,0,length*sizeof(*buffer));
1915     for (i=0; i < 4; i++)
1916     {
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       Hull(X[i],Y[i],image->columns,image->rows,pixel,buffer,-1);
1921     }
1922     j=(ssize_t) image->columns+2;
1923     for (y=0; y < (ssize_t) image->rows; y++)
1924     {
1925       MagickBooleanType
1926         sync;
1927
1928       register PixelPacket
1929         *restrict q;
1930
1931       q=GetCacheViewAuthenticPixels(despeckle_view,0,y,despeckle_image->columns,
1932         1,exception);
1933       if (q == (PixelPacket *) NULL)
1934         break;
1935       j++;
1936       for (x=0; x < (ssize_t) image->columns; x++)
1937       {
1938         switch (channel)
1939         {
1940           case 0: q->red=pixel[j]; break;
1941           case 1: q->green=pixel[j]; break;
1942           case 2: q->blue=pixel[j]; break;
1943           case 3: q->opacity=pixel[j]; break;
1944           default: break;
1945         }
1946         q++;
1947         j++;
1948       }
1949       sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
1950       if (sync == MagickFalse)
1951         {
1952           status=MagickFalse;
1953           break;
1954         }
1955       j++;
1956     }
1957     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1958       {
1959         MagickBooleanType
1960           proceed;
1961
1962 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1963   #pragma omp critical (MagickCore_DespeckleImage)
1964 #endif
1965         proceed=SetImageProgress(image,DespeckleImageTag,(MagickOffsetType)
1966           channel,3);
1967         if (proceed == MagickFalse)
1968           status=MagickFalse;
1969       }
1970   }
1971   despeckle_view=DestroyCacheView(despeckle_view);
1972   image_view=DestroyCacheView(image_view);
1973   buffers=DestroyPixelThreadSet(buffers);
1974   pixels=DestroyPixelThreadSet(pixels);
1975   despeckle_image->type=image->type;
1976   if (status == MagickFalse)
1977     despeckle_image=DestroyImage(despeckle_image);
1978   return(despeckle_image);
1979 }
1980 \f
1981 /*
1982 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1983 %                                                                             %
1984 %                                                                             %
1985 %                                                                             %
1986 %     E d g e I m a g e                                                       %
1987 %                                                                             %
1988 %                                                                             %
1989 %                                                                             %
1990 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1991 %
1992 %  EdgeImage() finds edges in an image.  Radius defines the radius of the
1993 %  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
1994 %  radius for you.
1995 %
1996 %  The format of the EdgeImage method is:
1997 %
1998 %      Image *EdgeImage(const Image *image,const double radius,
1999 %        ExceptionInfo *exception)
2000 %
2001 %  A description of each parameter follows:
2002 %
2003 %    o image: the image.
2004 %
2005 %    o radius: the radius of the pixel neighborhood.
2006 %
2007 %    o exception: return any errors or warnings in this structure.
2008 %
2009 */
2010 MagickExport Image *EdgeImage(const Image *image,const double radius,
2011   ExceptionInfo *exception)
2012 {
2013   Image
2014     *edge_image;
2015
2016   double
2017     *kernel;
2018
2019   register ssize_t
2020     i;
2021
2022   size_t
2023     width;
2024
2025   assert(image != (const Image *) NULL);
2026   assert(image->signature == MagickSignature);
2027   if (image->debug != MagickFalse)
2028     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2029   assert(exception != (ExceptionInfo *) NULL);
2030   assert(exception->signature == MagickSignature);
2031   width=GetOptimalKernelWidth1D(radius,0.5);
2032   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2033   if (kernel == (double *) NULL)
2034     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2035   for (i=0; i < (ssize_t) (width*width); i++)
2036     kernel[i]=(-1.0);
2037   kernel[i/2]=(double) (width*width-1.0);
2038   edge_image=ConvolveImage(image,width,kernel,exception);
2039   kernel=(double *) RelinquishMagickMemory(kernel);
2040   return(edge_image);
2041 }
2042 \f
2043 /*
2044 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2045 %                                                                             %
2046 %                                                                             %
2047 %                                                                             %
2048 %     E m b o s s I m a g e                                                   %
2049 %                                                                             %
2050 %                                                                             %
2051 %                                                                             %
2052 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2053 %
2054 %  EmbossImage() returns a grayscale image with a three-dimensional effect.
2055 %  We convolve the image with a Gaussian operator of the given radius and
2056 %  standard deviation (sigma).  For reasonable results, radius should be
2057 %  larger than sigma.  Use a radius of 0 and Emboss() selects a suitable
2058 %  radius for you.
2059 %
2060 %  The format of the EmbossImage method is:
2061 %
2062 %      Image *EmbossImage(const Image *image,const double radius,
2063 %        const double sigma,ExceptionInfo *exception)
2064 %
2065 %  A description of each parameter follows:
2066 %
2067 %    o image: the image.
2068 %
2069 %    o radius: the radius of the pixel neighborhood.
2070 %
2071 %    o sigma: the standard deviation of the Gaussian, in pixels.
2072 %
2073 %    o exception: return any errors or warnings in this structure.
2074 %
2075 */
2076 MagickExport Image *EmbossImage(const Image *image,const double radius,
2077   const double sigma,ExceptionInfo *exception)
2078 {
2079   double
2080     *kernel;
2081
2082   Image
2083     *emboss_image;
2084
2085   register ssize_t
2086     i;
2087
2088   size_t
2089     width;
2090
2091   ssize_t
2092     j,
2093     k,
2094     u,
2095     v;
2096
2097   assert(image != (Image *) NULL);
2098   assert(image->signature == MagickSignature);
2099   if (image->debug != MagickFalse)
2100     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2101   assert(exception != (ExceptionInfo *) NULL);
2102   assert(exception->signature == MagickSignature);
2103   width=GetOptimalKernelWidth2D(radius,sigma);
2104   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2105   if (kernel == (double *) NULL)
2106     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2107   j=(ssize_t) width/2;
2108   k=j;
2109   i=0;
2110   for (v=(-j); v <= j; v++)
2111   {
2112     for (u=(-j); u <= j; u++)
2113     {
2114       kernel[i]=((u < 0) || (v < 0) ? -8.0 : 8.0)*
2115         exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
2116         (2.0*MagickPI*MagickSigma*MagickSigma);
2117       if (u != k)
2118         kernel[i]=0.0;
2119       i++;
2120     }
2121     k--;
2122   }
2123   emboss_image=ConvolveImage(image,width,kernel,exception);
2124   if (emboss_image != (Image *) NULL)
2125     (void) EqualizeImage(emboss_image);
2126   kernel=(double *) RelinquishMagickMemory(kernel);
2127   return(emboss_image);
2128 }
2129 \f
2130 /*
2131 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2132 %                                                                             %
2133 %                                                                             %
2134 %                                                                             %
2135 %     F i l t e r I m a g e                                                   %
2136 %                                                                             %
2137 %                                                                             %
2138 %                                                                             %
2139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2140 %
2141 %  FilterImage() applies a custom convolution kernel to the image.
2142 %
2143 %  The format of the FilterImage method is:
2144 %
2145 %      Image *FilterImage(const Image *image,const KernelInfo *kernel,
2146 %        ExceptionInfo *exception)
2147 %      Image *FilterImageChannel(const Image *image,const ChannelType channel,
2148 %        const KernelInfo *kernel,ExceptionInfo *exception)
2149 %
2150 %  A description of each parameter follows:
2151 %
2152 %    o image: the image.
2153 %
2154 %    o channel: the channel type.
2155 %
2156 %    o kernel: the filtering kernel.
2157 %
2158 %    o exception: return any errors or warnings in this structure.
2159 %
2160 */
2161
2162 MagickExport Image *FilterImage(const Image *image,const KernelInfo *kernel,
2163   ExceptionInfo *exception)
2164 {
2165   Image
2166     *filter_image;
2167
2168   filter_image=FilterImageChannel(image,DefaultChannels,kernel,exception);
2169   return(filter_image);
2170 }
2171
2172 MagickExport Image *FilterImageChannel(const Image *image,
2173   const ChannelType channel,const KernelInfo *kernel,ExceptionInfo *exception)
2174 {
2175 #define FilterImageTag  "Filter/Image"
2176
2177   CacheView
2178     *filter_view,
2179     *image_view;
2180
2181   Image
2182     *filter_image;
2183
2184   MagickBooleanType
2185     status;
2186
2187   MagickOffsetType
2188     progress;
2189
2190   MagickPixelPacket
2191     bias;
2192
2193   ssize_t
2194     y;
2195
2196   /*
2197     Initialize filter image attributes.
2198   */
2199   assert(image != (Image *) NULL);
2200   assert(image->signature == MagickSignature);
2201   if (image->debug != MagickFalse)
2202     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2203   assert(exception != (ExceptionInfo *) NULL);
2204   assert(exception->signature == MagickSignature);
2205   if ((kernel->width % 2) == 0)
2206     ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
2207   filter_image=CloneImage(image,0,0,MagickTrue,exception);
2208   if (filter_image == (Image *) NULL)
2209     return((Image *) NULL);
2210   if (SetImageStorageClass(filter_image,DirectClass) == MagickFalse)
2211     {
2212       InheritException(exception,&filter_image->exception);
2213       filter_image=DestroyImage(filter_image);
2214       return((Image *) NULL);
2215     }
2216   if (image->debug != MagickFalse)
2217     {
2218       char
2219         format[MaxTextExtent],
2220         *message;
2221
2222       register const double
2223         *k;
2224
2225       ssize_t
2226         u,
2227         v;
2228
2229       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
2230         "  FilterImage with %.20gx%.20g kernel:",(double) kernel->width,(double)
2231         kernel->height);
2232       message=AcquireString("");
2233       k=kernel->values;
2234       for (v=0; v < (ssize_t) kernel->height; v++)
2235       {
2236         *message='\0';
2237         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) v);
2238         (void) ConcatenateString(&message,format);
2239         for (u=0; u < (ssize_t) kernel->width; u++)
2240         {
2241           (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
2242           (void) ConcatenateString(&message,format);
2243         }
2244         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
2245       }
2246       message=DestroyString(message);
2247     }
2248   status=AccelerateConvolveImage(image,kernel,filter_image,exception);
2249   if (status == MagickTrue)
2250     return(filter_image);
2251   /*
2252     Filter image.
2253   */
2254   status=MagickTrue;
2255   progress=0;
2256   GetMagickPixelPacket(image,&bias);
2257   SetMagickPixelPacketBias(image,&bias);
2258   image_view=AcquireCacheView(image);
2259   filter_view=AcquireCacheView(filter_image);
2260 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2261   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2262 #endif
2263   for (y=0; y < (ssize_t) image->rows; y++)
2264   {
2265     MagickBooleanType
2266       sync;
2267
2268     register const IndexPacket
2269       *restrict indexes;
2270
2271     register const PixelPacket
2272       *restrict p;
2273
2274     register IndexPacket
2275       *restrict filter_indexes;
2276
2277     register PixelPacket
2278       *restrict q;
2279
2280     register ssize_t
2281       x;
2282
2283     if (status == MagickFalse)
2284       continue;
2285     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel->width/2L),
2286       y-(ssize_t) (kernel->height/2L),image->columns+kernel->width,
2287       kernel->height,exception);
2288     q=GetCacheViewAuthenticPixels(filter_view,0,y,filter_image->columns,1,
2289       exception);
2290     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2291       {
2292         status=MagickFalse;
2293         continue;
2294       }
2295     indexes=GetCacheViewVirtualIndexQueue(image_view);
2296     filter_indexes=GetCacheViewAuthenticIndexQueue(filter_view);
2297     for (x=0; x < (ssize_t) image->columns; x++)
2298     {
2299       MagickPixelPacket
2300         pixel;
2301
2302       register const double
2303         *restrict k;
2304
2305       register const PixelPacket
2306         *restrict kernel_pixels;
2307
2308       register ssize_t
2309         u;
2310
2311       ssize_t
2312         v;
2313
2314       pixel=bias;
2315       k=kernel->values;
2316       kernel_pixels=p;
2317       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
2318         {
2319           for (v=0; v < (ssize_t) kernel->width; v++)
2320           {
2321             for (u=0; u < (ssize_t) kernel->height; u++)
2322             {
2323               pixel.red+=(*k)*kernel_pixels[u].red;
2324               pixel.green+=(*k)*kernel_pixels[u].green;
2325               pixel.blue+=(*k)*kernel_pixels[u].blue;
2326               k++;
2327             }
2328             kernel_pixels+=image->columns+kernel->width;
2329           }
2330           if ((channel & RedChannel) != 0)
2331             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
2332           if ((channel & GreenChannel) != 0)
2333             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
2334           if ((channel & BlueChannel) != 0)
2335             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
2336           if ((channel & OpacityChannel) != 0)
2337             {
2338               k=kernel->values;
2339               kernel_pixels=p;
2340               for (v=0; v < (ssize_t) kernel->width; v++)
2341               {
2342                 for (u=0; u < (ssize_t) kernel->height; u++)
2343                 {
2344                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
2345                   k++;
2346                 }
2347                 kernel_pixels+=image->columns+kernel->width;
2348               }
2349               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2350             }
2351           if (((channel & IndexChannel) != 0) &&
2352               (image->colorspace == CMYKColorspace))
2353             {
2354               register const IndexPacket
2355                 *restrict kernel_indexes;
2356
2357               k=kernel->values;
2358               kernel_indexes=indexes;
2359               for (v=0; v < (ssize_t) kernel->width; v++)
2360               {
2361                 for (u=0; u < (ssize_t) kernel->height; u++)
2362                 {
2363                   pixel.index+=(*k)*kernel_indexes[u];
2364                   k++;
2365                 }
2366                 kernel_indexes+=image->columns+kernel->width;
2367               }
2368               filter_indexes[x]=ClampToQuantum(pixel.index);
2369             }
2370         }
2371       else
2372         {
2373           MagickRealType
2374             alpha,
2375             gamma;
2376
2377           gamma=0.0;
2378           for (v=0; v < (ssize_t) kernel->width; v++)
2379           {
2380             for (u=0; u < (ssize_t) kernel->height; u++)
2381             {
2382               alpha=(MagickRealType) (QuantumScale*(QuantumRange-
2383                 kernel_pixels[u].opacity));
2384               pixel.red+=(*k)*alpha*kernel_pixels[u].red;
2385               pixel.green+=(*k)*alpha*kernel_pixels[u].green;
2386               pixel.blue+=(*k)*alpha*kernel_pixels[u].blue;
2387               gamma+=(*k)*alpha;
2388               k++;
2389             }
2390             kernel_pixels+=image->columns+kernel->width;
2391           }
2392           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2393           if ((channel & RedChannel) != 0)
2394             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
2395           if ((channel & GreenChannel) != 0)
2396             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
2397           if ((channel & BlueChannel) != 0)
2398             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
2399           if ((channel & OpacityChannel) != 0)
2400             {
2401               k=kernel->values;
2402               kernel_pixels=p;
2403               for (v=0; v < (ssize_t) kernel->width; v++)
2404               {
2405                 for (u=0; u < (ssize_t) kernel->height; u++)
2406                 {
2407                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
2408                   k++;
2409                 }
2410                 kernel_pixels+=image->columns+kernel->width;
2411               }
2412               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2413             }
2414           if (((channel & IndexChannel) != 0) &&
2415               (image->colorspace == CMYKColorspace))
2416             {
2417               register const IndexPacket
2418                 *restrict kernel_indexes;
2419
2420               k=kernel->values;
2421               kernel_pixels=p;
2422               kernel_indexes=indexes;
2423               for (v=0; v < (ssize_t) kernel->width; v++)
2424               {
2425                 for (u=0; u < (ssize_t) kernel->height; u++)
2426                 {
2427                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
2428                     kernel_pixels[u].opacity));
2429                   pixel.index+=(*k)*alpha*kernel_indexes[u];
2430                   k++;
2431                 }
2432                 kernel_pixels+=image->columns+kernel->width;
2433                 kernel_indexes+=image->columns+kernel->width;
2434               }
2435               filter_indexes[x]=ClampToQuantum(gamma*
2436                 GetIndexPixelComponent(&pixel));
2437             }
2438         }
2439       p++;
2440       q++;
2441     }
2442     sync=SyncCacheViewAuthenticPixels(filter_view,exception);
2443     if (sync == MagickFalse)
2444       status=MagickFalse;
2445     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2446       {
2447         MagickBooleanType
2448           proceed;
2449
2450 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2451   #pragma omp critical (MagickCore_FilterImageChannel)
2452 #endif
2453         proceed=SetImageProgress(image,FilterImageTag,progress++,image->rows);
2454         if (proceed == MagickFalse)
2455           status=MagickFalse;
2456       }
2457   }
2458   filter_image->type=image->type;
2459   filter_view=DestroyCacheView(filter_view);
2460   image_view=DestroyCacheView(image_view);
2461   if (status == MagickFalse)
2462     filter_image=DestroyImage(filter_image);
2463   return(filter_image);
2464 }
2465 \f
2466 /*
2467 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2468 %                                                                             %
2469 %                                                                             %
2470 %                                                                             %
2471 %     G a u s s i a n B l u r I m a g e                                       %
2472 %                                                                             %
2473 %                                                                             %
2474 %                                                                             %
2475 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2476 %
2477 %  GaussianBlurImage() blurs an image.  We convolve the image with a
2478 %  Gaussian operator of the given radius and standard deviation (sigma).
2479 %  For reasonable results, the radius should be larger than sigma.  Use a
2480 %  radius of 0 and GaussianBlurImage() selects a suitable radius for you
2481 %
2482 %  The format of the GaussianBlurImage method is:
2483 %
2484 %      Image *GaussianBlurImage(const Image *image,onst double radius,
2485 %        const double sigma,ExceptionInfo *exception)
2486 %      Image *GaussianBlurImageChannel(const Image *image,
2487 %        const ChannelType channel,const double radius,const double sigma,
2488 %        ExceptionInfo *exception)
2489 %
2490 %  A description of each parameter follows:
2491 %
2492 %    o image: the image.
2493 %
2494 %    o channel: the channel type.
2495 %
2496 %    o radius: the radius of the Gaussian, in pixels, not counting the center
2497 %      pixel.
2498 %
2499 %    o sigma: the standard deviation of the Gaussian, in pixels.
2500 %
2501 %    o exception: return any errors or warnings in this structure.
2502 %
2503 */
2504
2505 MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
2506   const double sigma,ExceptionInfo *exception)
2507 {
2508   Image
2509     *blur_image;
2510
2511   blur_image=GaussianBlurImageChannel(image,DefaultChannels,radius,sigma,
2512     exception);
2513   return(blur_image);
2514 }
2515
2516 MagickExport Image *GaussianBlurImageChannel(const Image *image,
2517   const ChannelType channel,const double radius,const double sigma,
2518   ExceptionInfo *exception)
2519 {
2520   double
2521     *kernel;
2522
2523   Image
2524     *blur_image;
2525
2526   register ssize_t
2527     i;
2528
2529   size_t
2530     width;
2531
2532   ssize_t
2533     j,
2534     u,
2535     v;
2536
2537   assert(image != (const Image *) NULL);
2538   assert(image->signature == MagickSignature);
2539   if (image->debug != MagickFalse)
2540     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2541   assert(exception != (ExceptionInfo *) NULL);
2542   assert(exception->signature == MagickSignature);
2543   width=GetOptimalKernelWidth2D(radius,sigma);
2544   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2545   if (kernel == (double *) NULL)
2546     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2547   j=(ssize_t) width/2;
2548   i=0;
2549   for (v=(-j); v <= j; v++)
2550   {
2551     for (u=(-j); u <= j; u++)
2552       kernel[i++]=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
2553         (2.0*MagickPI*MagickSigma*MagickSigma);
2554   }
2555   blur_image=ConvolveImageChannel(image,channel,width,kernel,exception);
2556   kernel=(double *) RelinquishMagickMemory(kernel);
2557   return(blur_image);
2558 }
2559 \f
2560 /*
2561 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2562 %                                                                             %
2563 %                                                                             %
2564 %                                                                             %
2565 %     M e d i a n F i l t e r I m a g e                                       %
2566 %                                                                             %
2567 %                                                                             %
2568 %                                                                             %
2569 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2570 %
2571 %  MedianFilterImage() applies a digital filter that improves the quality
2572 %  of a noisy image.  Each pixel is replaced by the median in a set of
2573 %  neighboring pixels as defined by radius.
2574 %
2575 %  The algorithm was contributed by Mike Edmonds and implements an insertion
2576 %  sort for selecting median color-channel values.  For more on this algorithm
2577 %  see "Skip Lists: A probabilistic Alternative to Balanced Trees" by William
2578 %  Pugh in the June 1990 of Communications of the ACM.
2579 %
2580 %  The format of the MedianFilterImage method is:
2581 %
2582 %      Image *MedianFilterImage(const Image *image,const double radius,
2583 %        ExceptionInfo *exception)
2584 %
2585 %  A description of each parameter follows:
2586 %
2587 %    o image: the image.
2588 %
2589 %    o radius: the radius of the pixel neighborhood.
2590 %
2591 %    o exception: return any errors or warnings in this structure.
2592 %
2593 */
2594
2595 #define MedianListChannels  5
2596
2597 typedef struct _MedianListNode
2598 {
2599   size_t
2600     next[9],
2601     count,
2602     signature;
2603 } MedianListNode;
2604
2605 typedef struct _MedianSkipList
2606 {
2607   ssize_t
2608     level;
2609
2610   MedianListNode
2611     *nodes;
2612 } MedianSkipList;
2613
2614 typedef struct _MedianPixelList
2615 {
2616   size_t
2617     center,
2618     seed,
2619     signature;
2620
2621   MedianSkipList
2622     lists[MedianListChannels];
2623 } MedianPixelList;
2624
2625 static MedianPixelList *DestroyMedianPixelList(MedianPixelList *pixel_list)
2626 {
2627   register ssize_t
2628     i;
2629
2630   if (pixel_list == (MedianPixelList *) NULL)
2631     return((MedianPixelList *) NULL);
2632   for (i=0; i < MedianListChannels; i++)
2633     if (pixel_list->lists[i].nodes != (MedianListNode *) NULL)
2634       pixel_list->lists[i].nodes=(MedianListNode *) RelinquishMagickMemory(
2635         pixel_list->lists[i].nodes);
2636   pixel_list=(MedianPixelList *) RelinquishAlignedMemory(pixel_list);
2637   return(pixel_list);
2638 }
2639
2640 static MedianPixelList **DestroyMedianPixelListThreadSet(
2641   MedianPixelList **pixel_list)
2642 {
2643   register ssize_t
2644     i;
2645
2646   assert(pixel_list != (MedianPixelList **) NULL);
2647   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2648     if (pixel_list[i] != (MedianPixelList *) NULL)
2649       pixel_list[i]=DestroyMedianPixelList(pixel_list[i]);
2650   pixel_list=(MedianPixelList **) RelinquishAlignedMemory(pixel_list);
2651   return(pixel_list);
2652 }
2653
2654 static MedianPixelList *AcquireMedianPixelList(const size_t width)
2655 {
2656   MedianPixelList
2657     *pixel_list;
2658
2659   register ssize_t
2660     i;
2661
2662   pixel_list=(MedianPixelList *) AcquireAlignedMemory(1,sizeof(*pixel_list));
2663   if (pixel_list == (MedianPixelList *) NULL)
2664     return(pixel_list);
2665   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
2666   pixel_list->center=width*width/2;
2667   for (i=0; i < MedianListChannels; i++)
2668   {
2669     pixel_list->lists[i].nodes=(MedianListNode *) AcquireQuantumMemory(65537UL,
2670       sizeof(*pixel_list->lists[i].nodes));
2671     if (pixel_list->lists[i].nodes == (MedianListNode *) NULL)
2672       return(DestroyMedianPixelList(pixel_list));
2673     (void) ResetMagickMemory(pixel_list->lists[i].nodes,0,65537UL*
2674       sizeof(*pixel_list->lists[i].nodes));
2675   }
2676   pixel_list->signature=MagickSignature;
2677   return(pixel_list);
2678 }
2679
2680 static MedianPixelList **AcquireMedianPixelListThreadSet(const size_t width)
2681 {
2682   MedianPixelList
2683     **pixel_list;
2684
2685   register ssize_t
2686     i;
2687
2688   size_t
2689     number_threads;
2690
2691   number_threads=GetOpenMPMaximumThreads();
2692   pixel_list=(MedianPixelList **) AcquireAlignedMemory(number_threads,
2693     sizeof(*pixel_list));
2694   if (pixel_list == (MedianPixelList **) NULL)
2695     return((MedianPixelList **) NULL);
2696   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
2697   for (i=0; i < (ssize_t) number_threads; i++)
2698   {
2699     pixel_list[i]=AcquireMedianPixelList(width);
2700     if (pixel_list[i] == (MedianPixelList *) NULL)
2701       return(DestroyMedianPixelListThreadSet(pixel_list));
2702   }
2703   return(pixel_list);
2704 }
2705
2706 static void AddNodeMedianPixelList(MedianPixelList *pixel_list,
2707   const ssize_t channel,const size_t color)
2708 {
2709   register MedianSkipList
2710     *list;
2711
2712   register ssize_t
2713     level;
2714
2715   size_t
2716     search,
2717     update[9];
2718
2719   /*
2720     Initialize the node.
2721   */
2722   list=pixel_list->lists+channel;
2723   list->nodes[color].signature=pixel_list->signature;
2724   list->nodes[color].count=1;
2725   /*
2726     Determine where it belongs in the list.
2727   */
2728   search=65536UL;
2729   for (level=list->level; level >= 0; level--)
2730   {
2731     while (list->nodes[search].next[level] < color)
2732       search=list->nodes[search].next[level];
2733     update[level]=search;
2734   }
2735   /*
2736     Generate a pseudo-random level for this node.
2737   */
2738   for (level=0; ; level++)
2739   {
2740     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2741     if ((pixel_list->seed & 0x300) != 0x300)
2742       break;
2743   }
2744   if (level > 8)
2745     level=8;
2746   if (level > (list->level+2))
2747     level=list->level+2;
2748   /*
2749     If we're raising the list's level, link back to the root node.
2750   */
2751   while (level > list->level)
2752   {
2753     list->level++;
2754     update[list->level]=65536UL;
2755   }
2756   /*
2757     Link the node into the skip-list.
2758   */
2759   do
2760   {
2761     list->nodes[color].next[level]=list->nodes[update[level]].next[level];
2762     list->nodes[update[level]].next[level]=color;
2763   }
2764   while (level-- > 0);
2765 }
2766
2767 static MagickPixelPacket GetMedianPixelList(MedianPixelList *pixel_list)
2768 {
2769   MagickPixelPacket
2770     pixel;
2771
2772   register MedianSkipList
2773     *list;
2774
2775   register ssize_t
2776     channel;
2777
2778   size_t
2779     center,
2780     color,
2781     count;
2782
2783   unsigned short
2784     channels[MedianListChannels];
2785
2786   /*
2787     Find the median value for each of the color.
2788   */
2789   center=pixel_list->center;
2790   for (channel=0; channel < 5; channel++)
2791   {
2792     list=pixel_list->lists+channel;
2793     color=65536UL;
2794     count=0;
2795     do
2796     {
2797       color=list->nodes[color].next[0];
2798       count+=list->nodes[color].count;
2799     }
2800     while (count <= center);
2801     channels[channel]=(unsigned short) color;
2802   }
2803   GetMagickPixelPacket((const Image *) NULL,&pixel);
2804   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
2805   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
2806   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
2807   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
2808   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
2809   return(pixel);
2810 }
2811
2812 static inline void InsertMedianPixelList(const Image *image,
2813   const PixelPacket *pixel,const IndexPacket *indexes,
2814   MedianPixelList *pixel_list)
2815 {
2816   size_t
2817     signature;
2818
2819   unsigned short
2820     index;
2821
2822   index=ScaleQuantumToShort(pixel->red);
2823   signature=pixel_list->lists[0].nodes[index].signature;
2824   if (signature == pixel_list->signature)
2825     pixel_list->lists[0].nodes[index].count++;
2826   else
2827     AddNodeMedianPixelList(pixel_list,0,index);
2828   index=ScaleQuantumToShort(pixel->green);
2829   signature=pixel_list->lists[1].nodes[index].signature;
2830   if (signature == pixel_list->signature)
2831     pixel_list->lists[1].nodes[index].count++;
2832   else
2833     AddNodeMedianPixelList(pixel_list,1,index);
2834   index=ScaleQuantumToShort(pixel->blue);
2835   signature=pixel_list->lists[2].nodes[index].signature;
2836   if (signature == pixel_list->signature)
2837     pixel_list->lists[2].nodes[index].count++;
2838   else
2839     AddNodeMedianPixelList(pixel_list,2,index);
2840   index=ScaleQuantumToShort(pixel->opacity);
2841   signature=pixel_list->lists[3].nodes[index].signature;
2842   if (signature == pixel_list->signature)
2843     pixel_list->lists[3].nodes[index].count++;
2844   else
2845     AddNodeMedianPixelList(pixel_list,3,index);
2846   if (image->colorspace == CMYKColorspace)
2847     index=ScaleQuantumToShort(*indexes);
2848   signature=pixel_list->lists[4].nodes[index].signature;
2849   if (signature == pixel_list->signature)
2850     pixel_list->lists[4].nodes[index].count++;
2851   else
2852     AddNodeMedianPixelList(pixel_list,4,index);
2853 }
2854
2855 static void ResetMedianPixelList(MedianPixelList *pixel_list)
2856 {
2857   int
2858     level;
2859
2860   register MedianListNode
2861     *root;
2862
2863   register MedianSkipList
2864     *list;
2865
2866   register ssize_t
2867     channel;
2868
2869   /*
2870     Reset the skip-list.
2871   */
2872   for (channel=0; channel < 5; channel++)
2873   {
2874     list=pixel_list->lists+channel;
2875     root=list->nodes+65536UL;
2876     list->level=0;
2877     for (level=0; level < 9; level++)
2878       root->next[level]=65536UL;
2879   }
2880   pixel_list->seed=pixel_list->signature++;
2881 }
2882
2883 MagickExport Image *MedianFilterImage(const Image *image,const double radius,
2884   ExceptionInfo *exception)
2885 {
2886 #define MedianFilterImageTag  "MedianFilter/Image"
2887
2888   CacheView
2889     *image_view,
2890     *median_view;
2891
2892   Image
2893     *median_image;
2894
2895   MagickBooleanType
2896     status;
2897
2898   MagickOffsetType
2899     progress;
2900
2901   MedianPixelList
2902     **restrict pixel_list;
2903
2904   size_t
2905     width;
2906
2907   ssize_t
2908     y;
2909
2910   /*
2911     Initialize median image attributes.
2912   */
2913   assert(image != (Image *) NULL);
2914   assert(image->signature == MagickSignature);
2915   if (image->debug != MagickFalse)
2916     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2917   assert(exception != (ExceptionInfo *) NULL);
2918   assert(exception->signature == MagickSignature);
2919   width=GetOptimalKernelWidth2D(radius,0.5);
2920   if ((image->columns < width) || (image->rows < width))
2921     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
2922   median_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2923     exception);
2924   if (median_image == (Image *) NULL)
2925     return((Image *) NULL);
2926   if (SetImageStorageClass(median_image,DirectClass) == MagickFalse)
2927     {
2928       InheritException(exception,&median_image->exception);
2929       median_image=DestroyImage(median_image);
2930       return((Image *) NULL);
2931     }
2932   pixel_list=AcquireMedianPixelListThreadSet(width);
2933   if (pixel_list == (MedianPixelList **) NULL)
2934     {
2935       median_image=DestroyImage(median_image);
2936       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2937     }
2938   /*
2939     Median filter each image row.
2940   */
2941   status=MagickTrue;
2942   progress=0;
2943   image_view=AcquireCacheView(image);
2944   median_view=AcquireCacheView(median_image);
2945 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2946   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2947 #endif
2948   for (y=0; y < (ssize_t) median_image->rows; y++)
2949   {
2950     int
2951       id;
2952
2953     register const IndexPacket
2954       *restrict indexes;
2955
2956     register const PixelPacket
2957       *restrict p;
2958
2959     register IndexPacket
2960       *restrict median_indexes;
2961
2962     register PixelPacket
2963       *restrict q;
2964
2965     register ssize_t
2966       x;
2967
2968     if (status == MagickFalse)
2969       continue;
2970     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
2971       (width/2L),image->columns+width,width,exception);
2972     q=QueueCacheViewAuthenticPixels(median_view,0,y,median_image->columns,1,
2973       exception);
2974     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2975       {
2976         status=MagickFalse;
2977         continue;
2978       }
2979     indexes=GetCacheViewVirtualIndexQueue(image_view);
2980     median_indexes=GetCacheViewAuthenticIndexQueue(median_view);
2981     id=GetOpenMPThreadId();
2982     for (x=0; x < (ssize_t) median_image->columns; x++)
2983     {
2984       MagickPixelPacket
2985         pixel;
2986
2987       register const IndexPacket
2988         *restrict s;
2989
2990       register const PixelPacket
2991         *restrict r;
2992
2993       register ssize_t
2994         u,
2995         v;
2996
2997       r=p;
2998       s=indexes+x;
2999       ResetMedianPixelList(pixel_list[id]);
3000       for (v=0; v < (ssize_t) width; v++)
3001       {
3002         for (u=0; u < (ssize_t) width; u++)
3003           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
3004         r+=image->columns+width;
3005         s+=image->columns+width;
3006       }
3007       pixel=GetMedianPixelList(pixel_list[id]);
3008       SetPixelPacket(median_image,&pixel,q,median_indexes+x);
3009       p++;
3010       q++;
3011     }
3012     if (SyncCacheViewAuthenticPixels(median_view,exception) == MagickFalse)
3013       status=MagickFalse;
3014     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3015       {
3016         MagickBooleanType
3017           proceed;
3018
3019 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3020   #pragma omp critical (MagickCore_MedianFilterImage)
3021 #endif
3022         proceed=SetImageProgress(image,MedianFilterImageTag,progress++,
3023           image->rows);
3024         if (proceed == MagickFalse)
3025           status=MagickFalse;
3026       }
3027   }
3028   median_view=DestroyCacheView(median_view);
3029   image_view=DestroyCacheView(image_view);
3030   pixel_list=DestroyMedianPixelListThreadSet(pixel_list);
3031   return(median_image);
3032 }
3033 \f
3034 /*
3035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3036 %                                                                             %
3037 %                                                                             %
3038 %                                                                             %
3039 %     M o t i o n B l u r I m a g e                                           %
3040 %                                                                             %
3041 %                                                                             %
3042 %                                                                             %
3043 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3044 %
3045 %  MotionBlurImage() simulates motion blur.  We convolve the image with a
3046 %  Gaussian operator of the given radius and standard deviation (sigma).
3047 %  For reasonable results, radius should be larger than sigma.  Use a
3048 %  radius of 0 and MotionBlurImage() selects a suitable radius for you.
3049 %  Angle gives the angle of the blurring motion.
3050 %
3051 %  Andrew Protano contributed this effect.
3052 %
3053 %  The format of the MotionBlurImage method is:
3054 %
3055 %    Image *MotionBlurImage(const Image *image,const double radius,
3056 %      const double sigma,const double angle,ExceptionInfo *exception)
3057 %    Image *MotionBlurImageChannel(const Image *image,const ChannelType channel,
3058 %      const double radius,const double sigma,const double angle,
3059 %      ExceptionInfo *exception)
3060 %
3061 %  A description of each parameter follows:
3062 %
3063 %    o image: the image.
3064 %
3065 %    o channel: the channel type.
3066 %
3067 %    o radius: the radius of the Gaussian, in pixels, not counting the center
3068 %    o radius: the radius of the Gaussian, in pixels, not counting
3069 %      the center pixel.
3070 %
3071 %    o sigma: the standard deviation of the Gaussian, in pixels.
3072 %
3073 %    o angle: Apply the effect along this angle.
3074 %
3075 %    o exception: return any errors or warnings in this structure.
3076 %
3077 */
3078
3079 static double *GetMotionBlurKernel(const size_t width,const double sigma)
3080 {
3081   double
3082     *kernel,
3083     normalize;
3084
3085   register ssize_t
3086     i;
3087
3088   /*
3089    Generate a 1-D convolution kernel.
3090   */
3091   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
3092   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
3093   if (kernel == (double *) NULL)
3094     return(kernel);
3095   normalize=0.0;
3096   for (i=0; i < (ssize_t) width; i++)
3097   {
3098     kernel[i]=exp((-((double) i*i)/(double) (2.0*MagickSigma*MagickSigma)))/
3099       (MagickSQ2PI*MagickSigma);
3100     normalize+=kernel[i];
3101   }
3102   for (i=0; i < (ssize_t) width; i++)
3103     kernel[i]/=normalize;
3104   return(kernel);
3105 }
3106
3107 MagickExport Image *MotionBlurImage(const Image *image,const double radius,
3108   const double sigma,const double angle,ExceptionInfo *exception)
3109 {
3110   Image
3111     *motion_blur;
3112
3113   motion_blur=MotionBlurImageChannel(image,DefaultChannels,radius,sigma,angle,
3114     exception);
3115   return(motion_blur);
3116 }
3117
3118 MagickExport Image *MotionBlurImageChannel(const Image *image,
3119   const ChannelType channel,const double radius,const double sigma,
3120   const double angle,ExceptionInfo *exception)
3121 {
3122   CacheView
3123     *blur_view,
3124     *image_view;
3125
3126   double
3127     *kernel;
3128
3129   Image
3130     *blur_image;
3131
3132   MagickBooleanType
3133     status;
3134
3135   MagickOffsetType
3136     progress;
3137
3138   MagickPixelPacket
3139     bias;
3140
3141   OffsetInfo
3142     *offset;
3143
3144   PointInfo
3145     point;
3146
3147   register ssize_t
3148     i;
3149
3150   size_t
3151     width;
3152
3153   ssize_t
3154     y;
3155
3156   assert(image != (Image *) NULL);
3157   assert(image->signature == MagickSignature);
3158   if (image->debug != MagickFalse)
3159     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3160   assert(exception != (ExceptionInfo *) NULL);
3161   width=GetOptimalKernelWidth1D(radius,sigma);
3162   kernel=GetMotionBlurKernel(width,sigma);
3163   if (kernel == (double *) NULL)
3164     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3165   offset=(OffsetInfo *) AcquireQuantumMemory(width,sizeof(*offset));
3166   if (offset == (OffsetInfo *) NULL)
3167     {
3168       kernel=(double *) RelinquishMagickMemory(kernel);
3169       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3170     }
3171   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3172   if (blur_image == (Image *) NULL)
3173     {
3174       kernel=(double *) RelinquishMagickMemory(kernel);
3175       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
3176       return((Image *) NULL);
3177     }
3178   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
3179     {
3180       kernel=(double *) RelinquishMagickMemory(kernel);
3181       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
3182       InheritException(exception,&blur_image->exception);
3183       blur_image=DestroyImage(blur_image);
3184       return((Image *) NULL);
3185     }
3186   point.x=(double) width*sin(DegreesToRadians(angle));
3187   point.y=(double) width*cos(DegreesToRadians(angle));
3188   for (i=0; i < (ssize_t) width; i++)
3189   {
3190     offset[i].x=(ssize_t) ceil((double) (i*point.y)/hypot(point.x,point.y)-0.5);
3191     offset[i].y=(ssize_t) ceil((double) (i*point.x)/hypot(point.x,point.y)-0.5);
3192   }
3193   /*
3194     Motion blur image.
3195   */
3196   status=MagickTrue;
3197   progress=0;
3198   GetMagickPixelPacket(image,&bias);
3199   image_view=AcquireCacheView(image);
3200   blur_view=AcquireCacheView(blur_image);
3201 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3202   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
3203 #endif
3204   for (y=0; y < (ssize_t) image->rows; y++)
3205   {
3206     register IndexPacket
3207       *restrict blur_indexes;
3208
3209     register PixelPacket
3210       *restrict q;
3211
3212     register ssize_t
3213       x;
3214
3215     if (status == MagickFalse)
3216       continue;
3217     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
3218       exception);
3219     if (q == (PixelPacket *) NULL)
3220       {
3221         status=MagickFalse;
3222         continue;
3223       }
3224     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
3225     for (x=0; x < (ssize_t) image->columns; x++)
3226     {
3227       MagickPixelPacket
3228         qixel;
3229
3230       PixelPacket
3231         pixel;
3232
3233       register const IndexPacket
3234         *restrict indexes;
3235
3236       register double
3237         *restrict k;
3238
3239       register ssize_t
3240         i;
3241
3242       k=kernel;
3243       qixel=bias;
3244       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
3245         {
3246           for (i=0; i < (ssize_t) width; i++)
3247           {
3248             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
3249               offset[i].y,&pixel,exception);
3250             qixel.red+=(*k)*pixel.red;
3251             qixel.green+=(*k)*pixel.green;
3252             qixel.blue+=(*k)*pixel.blue;
3253             qixel.opacity+=(*k)*pixel.opacity;
3254             if (image->colorspace == CMYKColorspace)
3255               {
3256                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3257                 qixel.index+=(*k)*(*indexes);
3258               }
3259             k++;
3260           }
3261           if ((channel & RedChannel) != 0)
3262             q->red=ClampToQuantum(qixel.red);
3263           if ((channel & GreenChannel) != 0)
3264             q->green=ClampToQuantum(qixel.green);
3265           if ((channel & BlueChannel) != 0)
3266             q->blue=ClampToQuantum(qixel.blue);
3267           if ((channel & OpacityChannel) != 0)
3268             q->opacity=ClampToQuantum(qixel.opacity);
3269           if (((channel & IndexChannel) != 0) &&
3270               (image->colorspace == CMYKColorspace))
3271             blur_indexes[x]=(IndexPacket) ClampToQuantum(qixel.index);
3272         }
3273       else
3274         {
3275           MagickRealType
3276             alpha,
3277             gamma;
3278
3279           alpha=0.0;
3280           gamma=0.0;
3281           for (i=0; i < (ssize_t) width; i++)
3282           {
3283             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
3284               offset[i].y,&pixel,exception);
3285             alpha=(MagickRealType) (QuantumScale*
3286               GetAlphaPixelComponent(&pixel));
3287             qixel.red+=(*k)*alpha*pixel.red;
3288             qixel.green+=(*k)*alpha*pixel.green;
3289             qixel.blue+=(*k)*alpha*pixel.blue;
3290             qixel.opacity+=(*k)*pixel.opacity;
3291             if (image->colorspace == CMYKColorspace)
3292               {
3293                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3294                 qixel.index+=(*k)*alpha*(*indexes);
3295               }
3296             gamma+=(*k)*alpha;
3297             k++;
3298           }
3299           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
3300           if ((channel & RedChannel) != 0)
3301             q->red=ClampToQuantum(gamma*qixel.red);
3302           if ((channel & GreenChannel) != 0)
3303             q->green=ClampToQuantum(gamma*qixel.green);
3304           if ((channel & BlueChannel) != 0)
3305             q->blue=ClampToQuantum(gamma*qixel.blue);
3306           if ((channel & OpacityChannel) != 0)
3307             q->opacity=ClampToQuantum(qixel.opacity);
3308           if (((channel & IndexChannel) != 0) &&
3309               (image->colorspace == CMYKColorspace))
3310             blur_indexes[x]=(IndexPacket) ClampToQuantum(gamma*qixel.index);
3311         }
3312       q++;
3313     }
3314     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
3315       status=MagickFalse;
3316     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3317       {
3318         MagickBooleanType
3319           proceed;
3320
3321 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3322   #pragma omp critical (MagickCore_MotionBlurImageChannel)
3323 #endif
3324         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
3325         if (proceed == MagickFalse)
3326           status=MagickFalse;
3327       }
3328   }
3329   blur_view=DestroyCacheView(blur_view);
3330   image_view=DestroyCacheView(image_view);
3331   kernel=(double *) RelinquishMagickMemory(kernel);
3332   offset=(OffsetInfo *) RelinquishMagickMemory(offset);
3333   if (status == MagickFalse)
3334     blur_image=DestroyImage(blur_image);
3335   return(blur_image);
3336 }
3337 \f
3338 /*
3339 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3340 %                                                                             %
3341 %                                                                             %
3342 %                                                                             %
3343 %     P r e v i e w I m a g e                                                 %
3344 %                                                                             %
3345 %                                                                             %
3346 %                                                                             %
3347 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3348 %
3349 %  PreviewImage() tiles 9 thumbnails of the specified image with an image
3350 %  processing operation applied with varying parameters.  This may be helpful
3351 %  pin-pointing an appropriate parameter for a particular image processing
3352 %  operation.
3353 %
3354 %  The format of the PreviewImages method is:
3355 %
3356 %      Image *PreviewImages(const Image *image,const PreviewType preview,
3357 %        ExceptionInfo *exception)
3358 %
3359 %  A description of each parameter follows:
3360 %
3361 %    o image: the image.
3362 %
3363 %    o preview: the image processing operation.
3364 %
3365 %    o exception: return any errors or warnings in this structure.
3366 %
3367 */
3368 MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
3369   ExceptionInfo *exception)
3370 {
3371 #define NumberTiles  9
3372 #define PreviewImageTag  "Preview/Image"
3373 #define DefaultPreviewGeometry  "204x204+10+10"
3374
3375   char
3376     factor[MaxTextExtent],
3377     label[MaxTextExtent];
3378
3379   double
3380     degrees,
3381     gamma,
3382     percentage,
3383     radius,
3384     sigma,
3385     threshold;
3386
3387   Image
3388     *images,
3389     *montage_image,
3390     *preview_image,
3391     *thumbnail;
3392
3393   ImageInfo
3394     *preview_info;
3395
3396   MagickBooleanType
3397     proceed;
3398
3399   MontageInfo
3400     *montage_info;
3401
3402   QuantizeInfo
3403     quantize_info;
3404
3405   RectangleInfo
3406     geometry;
3407
3408   register ssize_t
3409     i,
3410     x;
3411
3412   size_t
3413     colors;
3414
3415   ssize_t
3416     y;
3417
3418   /*
3419     Open output image file.
3420   */
3421   assert(image != (Image *) NULL);
3422   assert(image->signature == MagickSignature);
3423   if (image->debug != MagickFalse)
3424     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3425   colors=2;
3426   degrees=0.0;
3427   gamma=(-0.2f);
3428   preview_info=AcquireImageInfo();
3429   SetGeometry(image,&geometry);
3430   (void) ParseMetaGeometry(DefaultPreviewGeometry,&geometry.x,&geometry.y,
3431     &geometry.width,&geometry.height);
3432   images=NewImageList();
3433   percentage=12.5;
3434   GetQuantizeInfo(&quantize_info);
3435   radius=0.0;
3436   sigma=1.0;
3437   threshold=0.0;
3438   x=0;
3439   y=0;
3440   for (i=0; i < NumberTiles; i++)
3441   {
3442     thumbnail=ThumbnailImage(image,geometry.width,geometry.height,exception);
3443     if (thumbnail == (Image *) NULL)
3444       break;
3445     (void) SetImageProgressMonitor(thumbnail,(MagickProgressMonitor) NULL,
3446       (void *) NULL);
3447     (void) SetImageProperty(thumbnail,"label",DefaultTileLabel);
3448     if (i == (NumberTiles/2))
3449       {
3450         (void) QueryColorDatabase("#dfdfdf",&thumbnail->matte_color,exception);
3451         AppendImageToList(&images,thumbnail);
3452         continue;
3453       }
3454     switch (preview)
3455     {
3456       case RotatePreview:
3457       {
3458         degrees+=45.0;
3459         preview_image=RotateImage(thumbnail,degrees,exception);
3460         (void) FormatMagickString(label,MaxTextExtent,"rotate %g",degrees);
3461         break;
3462       }
3463       case ShearPreview:
3464       {
3465         degrees+=5.0;
3466         preview_image=ShearImage(thumbnail,degrees,degrees,exception);
3467         (void) FormatMagickString(label,MaxTextExtent,"shear %gx%g",
3468           degrees,2.0*degrees);
3469         break;
3470       }
3471       case RollPreview:
3472       {
3473         x=(ssize_t) ((i+1)*thumbnail->columns)/NumberTiles;
3474         y=(ssize_t) ((i+1)*thumbnail->rows)/NumberTiles;
3475         preview_image=RollImage(thumbnail,x,y,exception);
3476         (void) FormatMagickString(label,MaxTextExtent,"roll %+.20gx%+.20g",
3477           (double) x,(double) y);
3478         break;
3479       }
3480       case HuePreview:
3481       {
3482         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3483         if (preview_image == (Image *) NULL)
3484           break;
3485         (void) FormatMagickString(factor,MaxTextExtent,"100,100,%g",
3486           2.0*percentage);
3487         (void) ModulateImage(preview_image,factor);
3488         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3489         break;
3490       }
3491       case SaturationPreview:
3492       {
3493         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3494         if (preview_image == (Image *) NULL)
3495           break;
3496         (void) FormatMagickString(factor,MaxTextExtent,"100,%g",
3497           2.0*percentage);
3498         (void) ModulateImage(preview_image,factor);
3499         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3500         break;
3501       }
3502       case BrightnessPreview:
3503       {
3504         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3505         if (preview_image == (Image *) NULL)
3506           break;
3507         (void) FormatMagickString(factor,MaxTextExtent,"%g",2.0*percentage);
3508         (void) ModulateImage(preview_image,factor);
3509         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3510         break;
3511       }
3512       case GammaPreview:
3513       default:
3514       {
3515         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3516         if (preview_image == (Image *) NULL)
3517           break;
3518         gamma+=0.4f;
3519         (void) GammaImageChannel(preview_image,DefaultChannels,gamma);
3520         (void) FormatMagickString(label,MaxTextExtent,"gamma %g",gamma);
3521         break;
3522       }
3523       case SpiffPreview:
3524       {
3525         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3526         if (preview_image != (Image *) NULL)
3527           for (x=0; x < i; x++)
3528             (void) ContrastImage(preview_image,MagickTrue);
3529         (void) FormatMagickString(label,MaxTextExtent,"contrast (%.20g)",
3530           (double) i+1);
3531         break;
3532       }
3533       case DullPreview:
3534       {
3535         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3536         if (preview_image == (Image *) NULL)
3537           break;
3538         for (x=0; x < i; x++)
3539           (void) ContrastImage(preview_image,MagickFalse);
3540         (void) FormatMagickString(label,MaxTextExtent,"+contrast (%.20g)",
3541           (double) i+1);
3542         break;
3543       }
3544       case GrayscalePreview:
3545       {
3546         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3547         if (preview_image == (Image *) NULL)
3548           break;
3549         colors<<=1;
3550         quantize_info.number_colors=colors;
3551         quantize_info.colorspace=GRAYColorspace;
3552         (void) QuantizeImage(&quantize_info,preview_image);
3553         (void) FormatMagickString(label,MaxTextExtent,
3554           "-colorspace gray -colors %.20g",(double) colors);
3555         break;
3556       }
3557       case QuantizePreview:
3558       {
3559         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3560         if (preview_image == (Image *) NULL)
3561           break;
3562         colors<<=1;
3563         quantize_info.number_colors=colors;
3564         (void) QuantizeImage(&quantize_info,preview_image);
3565         (void) FormatMagickString(label,MaxTextExtent,"colors %.20g",(double)
3566           colors);
3567         break;
3568       }
3569       case DespecklePreview:
3570       {
3571         for (x=0; x < (i-1); x++)
3572         {
3573           preview_image=DespeckleImage(thumbnail,exception);
3574           if (preview_image == (Image *) NULL)
3575             break;
3576           thumbnail=DestroyImage(thumbnail);
3577           thumbnail=preview_image;
3578         }
3579         preview_image=DespeckleImage(thumbnail,exception);
3580         if (preview_image == (Image *) NULL)
3581           break;
3582         (void) FormatMagickString(label,MaxTextExtent,"despeckle (%.20g)",
3583           (double) i+1);
3584         break;
3585       }
3586       case ReduceNoisePreview:
3587       {
3588         preview_image=ReduceNoiseImage(thumbnail,radius,exception);
3589         (void) FormatMagickString(label,MaxTextExtent,"noise %g",radius);
3590         break;
3591       }
3592       case AddNoisePreview:
3593       {
3594         switch ((int) i)
3595         {
3596           case 0:
3597           {
3598             (void) CopyMagickString(factor,"uniform",MaxTextExtent);
3599             break;
3600           }
3601           case 1:
3602           {
3603             (void) CopyMagickString(factor,"gaussian",MaxTextExtent);
3604             break;
3605           }
3606           case 2:
3607           {
3608             (void) CopyMagickString(factor,"multiplicative",MaxTextExtent);
3609             break;
3610           }
3611           case 3:
3612           {
3613             (void) CopyMagickString(factor,"impulse",MaxTextExtent);
3614             break;
3615           }
3616           case 4:
3617           {
3618             (void) CopyMagickString(factor,"laplacian",MaxTextExtent);
3619             break;
3620           }
3621           case 5:
3622           {
3623             (void) CopyMagickString(factor,"Poisson",MaxTextExtent);
3624             break;
3625           }
3626           default:
3627           {
3628             (void) CopyMagickString(thumbnail->magick,"NULL",MaxTextExtent);
3629             break;
3630           }
3631         }
3632         preview_image=ReduceNoiseImage(thumbnail,(double) i,exception);
3633         (void) FormatMagickString(label,MaxTextExtent,"+noise %s",factor);
3634         break;
3635       }
3636       case SharpenPreview:
3637       {
3638         preview_image=SharpenImage(thumbnail,radius,sigma,exception);
3639         (void) FormatMagickString(label,MaxTextExtent,"sharpen %gx%g",
3640           radius,sigma);
3641         break;
3642       }
3643       case BlurPreview:
3644       {
3645         preview_image=BlurImage(thumbnail,radius,sigma,exception);
3646         (void) FormatMagickString(label,MaxTextExtent,"blur %gx%g",radius,
3647           sigma);
3648         break;
3649       }
3650       case ThresholdPreview:
3651       {
3652         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3653         if (preview_image == (Image *) NULL)
3654           break;
3655         (void) BilevelImage(thumbnail,
3656           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
3657         (void) FormatMagickString(label,MaxTextExtent,"threshold %g",
3658           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
3659         break;
3660       }
3661       case EdgeDetectPreview:
3662       {
3663         preview_image=EdgeImage(thumbnail,radius,exception);
3664         (void) FormatMagickString(label,MaxTextExtent,"edge %g",radius);
3665         break;
3666       }
3667       case SpreadPreview:
3668       {
3669         preview_image=SpreadImage(thumbnail,radius,exception);
3670         (void) FormatMagickString(label,MaxTextExtent,"spread %g",
3671           radius+0.5);
3672         break;
3673       }
3674       case SolarizePreview:
3675       {
3676         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3677         if (preview_image == (Image *) NULL)
3678           break;
3679         (void) SolarizeImage(preview_image,(double) QuantumRange*
3680           percentage/100.0);
3681         (void) FormatMagickString(label,MaxTextExtent,"solarize %g",
3682           (QuantumRange*percentage)/100.0);
3683         break;
3684       }
3685       case ShadePreview:
3686       {
3687         degrees+=10.0;
3688         preview_image=ShadeImage(thumbnail,MagickTrue,degrees,degrees,
3689           exception);
3690         (void) FormatMagickString(label,MaxTextExtent,"shade %gx%g",
3691           degrees,degrees);
3692         break;
3693       }
3694       case RaisePreview:
3695       {
3696         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3697         if (preview_image == (Image *) NULL)
3698           break;
3699         geometry.width=(size_t) (2*i+2);
3700         geometry.height=(size_t) (2*i+2);
3701         geometry.x=i/2;
3702         geometry.y=i/2;
3703         (void) RaiseImage(preview_image,&geometry,MagickTrue);
3704         (void) FormatMagickString(label,MaxTextExtent,
3705           "raise %.20gx%.20g%+.20g%+.20g",(double) geometry.width,(double)
3706           geometry.height,(double) geometry.x,(double) geometry.y);
3707         break;
3708       }
3709       case SegmentPreview:
3710       {
3711         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3712         if (preview_image == (Image *) NULL)
3713           break;
3714         threshold+=0.4f;
3715         (void) SegmentImage(preview_image,RGBColorspace,MagickFalse,threshold,
3716           threshold);
3717         (void) FormatMagickString(label,MaxTextExtent,"segment %gx%g",
3718           threshold,threshold);
3719         break;
3720       }
3721       case SwirlPreview:
3722       {
3723         preview_image=SwirlImage(thumbnail,degrees,exception);
3724         (void) FormatMagickString(label,MaxTextExtent,"swirl %g",degrees);
3725         degrees+=45.0;
3726         break;
3727       }
3728       case ImplodePreview:
3729       {
3730         degrees+=0.1f;
3731         preview_image=ImplodeImage(thumbnail,degrees,exception);
3732         (void) FormatMagickString(label,MaxTextExtent,"implode %g",degrees);
3733         break;
3734       }
3735       case WavePreview:
3736       {
3737         degrees+=5.0f;
3738         preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,exception);
3739         (void) FormatMagickString(label,MaxTextExtent,"wave %gx%g",
3740           0.5*degrees,2.0*degrees);
3741         break;
3742       }
3743       case OilPaintPreview:
3744       {
3745         preview_image=OilPaintImage(thumbnail,(double) radius,exception);
3746         (void) FormatMagickString(label,MaxTextExtent,"paint %g",radius);
3747         break;
3748       }
3749       case CharcoalDrawingPreview:
3750       {
3751         preview_image=CharcoalImage(thumbnail,(double) radius,(double) sigma,
3752           exception);
3753         (void) FormatMagickString(label,MaxTextExtent,"charcoal %gx%g",
3754           radius,sigma);
3755         break;
3756       }
3757       case JPEGPreview:
3758       {
3759         char
3760           filename[MaxTextExtent];
3761
3762         int
3763           file;
3764
3765         MagickBooleanType
3766           status;
3767
3768         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3769         if (preview_image == (Image *) NULL)
3770           break;
3771         preview_info->quality=(size_t) percentage;
3772         (void) FormatMagickString(factor,MaxTextExtent,"%.20g",(double)
3773           preview_info->quality);
3774         file=AcquireUniqueFileResource(filename);
3775         if (file != -1)
3776           file=close(file)-1;
3777         (void) FormatMagickString(preview_image->filename,MaxTextExtent,
3778           "jpeg:%s",filename);
3779         status=WriteImage(preview_info,preview_image);
3780         if (status != MagickFalse)
3781           {
3782             Image
3783               *quality_image;
3784
3785             (void) CopyMagickString(preview_info->filename,
3786               preview_image->filename,MaxTextExtent);
3787             quality_image=ReadImage(preview_info,exception);
3788             if (quality_image != (Image *) NULL)
3789               {
3790                 preview_image=DestroyImage(preview_image);
3791                 preview_image=quality_image;
3792               }
3793           }
3794         (void) RelinquishUniqueFileResource(preview_image->filename);
3795         if ((GetBlobSize(preview_image)/1024) >= 1024)
3796           (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%gmb ",
3797             factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
3798             1024.0/1024.0);
3799         else
3800           if (GetBlobSize(preview_image) >= 1024)
3801             (void) FormatMagickString(label,MaxTextExtent,
3802               "quality %s\n%gkb ",factor,(double) ((MagickOffsetType)
3803               GetBlobSize(preview_image))/1024.0);
3804           else
3805             (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%.20gb ",
3806               factor,(double) GetBlobSize(thumbnail));
3807         break;
3808       }
3809     }
3810     thumbnail=DestroyImage(thumbnail);
3811     percentage+=12.5;
3812     radius+=0.5;
3813     sigma+=0.25;
3814     if (preview_image == (Image *) NULL)
3815       break;
3816     (void) DeleteImageProperty(preview_image,"label");
3817     (void) SetImageProperty(preview_image,"label",label);
3818     AppendImageToList(&images,preview_image);
3819     proceed=SetImageProgress(image,PreviewImageTag,(MagickOffsetType) i,
3820       NumberTiles);
3821     if (proceed == MagickFalse)
3822       break;
3823   }
3824   if (images == (Image *) NULL)
3825     {
3826       preview_info=DestroyImageInfo(preview_info);
3827       return((Image *) NULL);
3828     }
3829   /*
3830     Create the montage.
3831   */
3832   montage_info=CloneMontageInfo(preview_info,(MontageInfo *) NULL);
3833   (void) CopyMagickString(montage_info->filename,image->filename,MaxTextExtent);
3834   montage_info->shadow=MagickTrue;
3835   (void) CloneString(&montage_info->tile,"3x3");
3836   (void) CloneString(&montage_info->geometry,DefaultPreviewGeometry);
3837   (void) CloneString(&montage_info->frame,DefaultTileFrame);
3838   montage_image=MontageImages(images,montage_info,exception);
3839   montage_info=DestroyMontageInfo(montage_info);
3840   images=DestroyImageList(images);
3841   if (montage_image == (Image *) NULL)
3842     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3843   if (montage_image->montage != (char *) NULL)
3844     {
3845       /*
3846         Free image directory.
3847       */
3848       montage_image->montage=(char *) RelinquishMagickMemory(
3849         montage_image->montage);
3850       if (image->directory != (char *) NULL)
3851         montage_image->directory=(char *) RelinquishMagickMemory(
3852           montage_image->directory);
3853     }
3854   preview_info=DestroyImageInfo(preview_info);
3855   return(montage_image);
3856 }
3857 \f
3858 /*
3859 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3860 %                                                                             %
3861 %                                                                             %
3862 %                                                                             %
3863 %     R a d i a l B l u r I m a g e                                           %
3864 %                                                                             %
3865 %                                                                             %
3866 %                                                                             %
3867 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3868 %
3869 %  RadialBlurImage() applies a radial blur to the image.
3870 %
3871 %  Andrew Protano contributed this effect.
3872 %
3873 %  The format of the RadialBlurImage method is:
3874 %
3875 %    Image *RadialBlurImage(const Image *image,const double angle,
3876 %      ExceptionInfo *exception)
3877 %    Image *RadialBlurImageChannel(const Image *image,const ChannelType channel,
3878 %      const double angle,ExceptionInfo *exception)
3879 %
3880 %  A description of each parameter follows:
3881 %
3882 %    o image: the image.
3883 %
3884 %    o channel: the channel type.
3885 %
3886 %    o angle: the angle of the radial blur.
3887 %
3888 %    o exception: return any errors or warnings in this structure.
3889 %
3890 */
3891
3892 MagickExport Image *RadialBlurImage(const Image *image,const double angle,
3893   ExceptionInfo *exception)
3894 {
3895   Image
3896     *blur_image;
3897
3898   blur_image=RadialBlurImageChannel(image,DefaultChannels,angle,exception);
3899   return(blur_image);
3900 }
3901
3902 MagickExport Image *RadialBlurImageChannel(const Image *image,
3903   const ChannelType channel,const double angle,ExceptionInfo *exception)
3904 {
3905   CacheView
3906     *blur_view,
3907     *image_view;
3908
3909   Image
3910     *blur_image;
3911
3912   MagickBooleanType
3913     status;
3914
3915   MagickOffsetType
3916     progress;
3917
3918   MagickPixelPacket
3919     bias;
3920
3921   MagickRealType
3922     blur_radius,
3923     *cos_theta,
3924     offset,
3925     *sin_theta,
3926     theta;
3927
3928   PointInfo
3929     blur_center;
3930
3931   register ssize_t
3932     i;
3933
3934   size_t
3935     n;
3936
3937   ssize_t
3938     y;
3939
3940   /*
3941     Allocate blur image.
3942   */
3943   assert(image != (Image *) NULL);
3944   assert(image->signature == MagickSignature);
3945   if (image->debug != MagickFalse)
3946     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3947   assert(exception != (ExceptionInfo *) NULL);
3948   assert(exception->signature == MagickSignature);
3949   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3950   if (blur_image == (Image *) NULL)
3951     return((Image *) NULL);
3952   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
3953     {
3954       InheritException(exception,&blur_image->exception);
3955       blur_image=DestroyImage(blur_image);
3956       return((Image *) NULL);
3957     }
3958   blur_center.x=(double) image->columns/2.0;
3959   blur_center.y=(double) image->rows/2.0;
3960   blur_radius=hypot(blur_center.x,blur_center.y);
3961   n=(size_t) fabs(4.0*DegreesToRadians(angle)*sqrt((double) blur_radius)+2UL);
3962   theta=DegreesToRadians(angle)/(MagickRealType) (n-1);
3963   cos_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
3964     sizeof(*cos_theta));
3965   sin_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
3966     sizeof(*sin_theta));
3967   if ((cos_theta == (MagickRealType *) NULL) ||
3968       (sin_theta == (MagickRealType *) NULL))
3969     {
3970       blur_image=DestroyImage(blur_image);
3971       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3972     }
3973   offset=theta*(MagickRealType) (n-1)/2.0;
3974   for (i=0; i < (ssize_t) n; i++)
3975   {
3976     cos_theta[i]=cos((double) (theta*i-offset));
3977     sin_theta[i]=sin((double) (theta*i-offset));
3978   }
3979   /*
3980     Radial blur image.
3981   */
3982   status=MagickTrue;
3983   progress=0;
3984   GetMagickPixelPacket(image,&bias);
3985   image_view=AcquireCacheView(image);
3986   blur_view=AcquireCacheView(blur_image);
3987 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3988   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
3989 #endif
3990   for (y=0; y < (ssize_t) blur_image->rows; y++)
3991   {
3992     register const IndexPacket
3993       *restrict indexes;
3994
3995     register IndexPacket
3996       *restrict blur_indexes;
3997
3998     register PixelPacket
3999       *restrict q;
4000
4001     register ssize_t
4002       x;
4003
4004     if (status == MagickFalse)
4005       continue;
4006     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
4007       exception);
4008     if (q == (PixelPacket *) NULL)
4009       {
4010         status=MagickFalse;
4011         continue;
4012       }
4013     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
4014     for (x=0; x < (ssize_t) blur_image->columns; x++)
4015     {
4016       MagickPixelPacket
4017         qixel;
4018
4019       MagickRealType
4020         normalize,
4021         radius;
4022
4023       PixelPacket
4024         pixel;
4025
4026       PointInfo
4027         center;
4028
4029       register ssize_t
4030         i;
4031
4032       size_t
4033         step;
4034
4035       center.x=(double) x-blur_center.x;
4036       center.y=(double) y-blur_center.y;
4037       radius=hypot((double) center.x,center.y);
4038       if (radius == 0)
4039         step=1;
4040       else
4041         {
4042           step=(size_t) (blur_radius/radius);
4043           if (step == 0)
4044             step=1;
4045           else
4046             if (step >= n)
4047               step=n-1;
4048         }
4049       normalize=0.0;
4050       qixel=bias;
4051       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
4052         {
4053           for (i=0; i < (ssize_t) n; i+=(ssize_t) step)
4054           {
4055             (void) GetOneCacheViewVirtualPixel(image_view,(ssize_t)
4056               (blur_center.x+center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),
4057               (ssize_t) (blur_center.y+center.x*sin_theta[i]+center.y*
4058               cos_theta[i]+0.5),&pixel,exception);
4059             qixel.red+=pixel.red;
4060             qixel.green+=pixel.green;
4061             qixel.blue+=pixel.blue;
4062             qixel.opacity+=pixel.opacity;
4063             if (image->colorspace == CMYKColorspace)
4064               {
4065                 indexes=GetCacheViewVirtualIndexQueue(image_view);
4066                 qixel.index+=(*indexes);
4067               }
4068             normalize+=1.0;
4069           }
4070           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
4071             normalize);
4072           if ((channel & RedChannel) != 0)
4073             q->red=ClampToQuantum(normalize*qixel.red);
4074           if ((channel & GreenChannel) != 0)
4075             q->green=ClampToQuantum(normalize*qixel.green);
4076           if ((channel & BlueChannel) != 0)
4077             q->blue=ClampToQuantum(normalize*qixel.blue);
4078           if ((channel & OpacityChannel) != 0)
4079             q->opacity=ClampToQuantum(normalize*qixel.opacity);
4080           if (((channel & IndexChannel) != 0) &&
4081               (image->colorspace == CMYKColorspace))
4082             blur_indexes[x]=(IndexPacket) ClampToQuantum(normalize*qixel.index);
4083         }
4084       else
4085         {
4086           MagickRealType
4087             alpha,
4088             gamma;
4089
4090           alpha=1.0;
4091           gamma=0.0;
4092           for (i=0; i < (ssize_t) n; i+=(ssize_t) step)
4093           {
4094             (void) GetOneCacheViewVirtualPixel(image_view,(ssize_t)
4095               (blur_center.x+center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),
4096               (ssize_t) (blur_center.y+center.x*sin_theta[i]+center.y*
4097               cos_theta[i]+0.5),&pixel,exception);
4098             alpha=(MagickRealType) (QuantumScale*
4099               GetAlphaPixelComponent(&pixel));
4100             qixel.red+=alpha*pixel.red;
4101             qixel.green+=alpha*pixel.green;
4102             qixel.blue+=alpha*pixel.blue;
4103             qixel.opacity+=pixel.opacity;
4104             if (image->colorspace == CMYKColorspace)
4105               {
4106                 indexes=GetCacheViewVirtualIndexQueue(image_view);
4107                 qixel.index+=alpha*(*indexes);
4108               }
4109             gamma+=alpha;
4110             normalize+=1.0;
4111           }
4112           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4113           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
4114             normalize);
4115           if ((channel & RedChannel) != 0)
4116             q->red=ClampToQuantum(gamma*qixel.red);
4117           if ((channel & GreenChannel) != 0)
4118             q->green=ClampToQuantum(gamma*qixel.green);
4119           if ((channel & BlueChannel) != 0)
4120             q->blue=ClampToQuantum(gamma*qixel.blue);
4121           if ((channel & OpacityChannel) != 0)
4122             q->opacity=ClampToQuantum(normalize*qixel.opacity);
4123           if (((channel & IndexChannel) != 0) &&
4124               (image->colorspace == CMYKColorspace))
4125             blur_indexes[x]=(IndexPacket) ClampToQuantum(gamma*qixel.index);
4126         }
4127       q++;
4128     }
4129     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
4130       status=MagickFalse;
4131     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4132       {
4133         MagickBooleanType
4134           proceed;
4135
4136 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4137   #pragma omp critical (MagickCore_RadialBlurImageChannel)
4138 #endif
4139         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
4140         if (proceed == MagickFalse)
4141           status=MagickFalse;
4142       }
4143   }
4144   blur_view=DestroyCacheView(blur_view);
4145   image_view=DestroyCacheView(image_view);
4146   cos_theta=(MagickRealType *) RelinquishMagickMemory(cos_theta);
4147   sin_theta=(MagickRealType *) RelinquishMagickMemory(sin_theta);
4148   if (status == MagickFalse)
4149     blur_image=DestroyImage(blur_image);
4150   return(blur_image);
4151 }
4152 \f
4153 /*
4154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4155 %                                                                             %
4156 %                                                                             %
4157 %                                                                             %
4158 %     R e d u c e N o i s e I m a g e                                         %
4159 %                                                                             %
4160 %                                                                             %
4161 %                                                                             %
4162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4163 %
4164 %  ReduceNoiseImage() smooths the contours of an image while still preserving
4165 %  edge information.  The algorithm works by replacing each pixel with its
4166 %  neighbor closest in value.  A neighbor is defined by radius.  Use a radius
4167 %  of 0 and ReduceNoise() selects a suitable radius for you.
4168 %
4169 %  The format of the ReduceNoiseImage method is:
4170 %
4171 %      Image *ReduceNoiseImage(const Image *image,const double radius,
4172 %        ExceptionInfo *exception)
4173 %
4174 %  A description of each parameter follows:
4175 %
4176 %    o image: the image.
4177 %
4178 %    o radius: the radius of the pixel neighborhood.
4179 %
4180 %    o exception: return any errors or warnings in this structure.
4181 %
4182 */
4183
4184 static MagickPixelPacket GetNonpeakMedianPixelList(MedianPixelList *pixel_list)
4185 {
4186   MagickPixelPacket
4187     pixel;
4188
4189   register MedianSkipList
4190     *list;
4191
4192   register ssize_t
4193     channel;
4194
4195   size_t
4196     center,
4197     color,
4198     count,
4199     previous,
4200     next;
4201
4202   unsigned short
4203     channels[5];
4204
4205   /*
4206     Finds the median value for each of the color.
4207   */
4208   center=pixel_list->center;
4209   for (channel=0; channel < 5; channel++)
4210   {
4211     list=pixel_list->lists+channel;
4212     color=65536UL;
4213     next=list->nodes[color].next[0];
4214     count=0;
4215     do
4216     {
4217       previous=color;
4218       color=next;
4219       next=list->nodes[color].next[0];
4220       count+=list->nodes[color].count;
4221     }
4222     while (count <= center);
4223     if ((previous == 65536UL) && (next != 65536UL))
4224       color=next;
4225     else
4226       if ((previous != 65536UL) && (next == 65536UL))
4227         color=previous;
4228     channels[channel]=(unsigned short) color;
4229   }
4230   GetMagickPixelPacket((const Image *) NULL,&pixel);
4231   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
4232   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
4233   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
4234   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
4235   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
4236   return(pixel);
4237 }
4238
4239 MagickExport Image *ReduceNoiseImage(const Image *image,const double radius,
4240   ExceptionInfo *exception)
4241 {
4242 #define ReduceNoiseImageTag  "ReduceNoise/Image"
4243
4244   CacheView
4245     *image_view,
4246     *noise_view;
4247
4248   Image
4249     *noise_image;
4250
4251   MagickBooleanType
4252     status;
4253
4254   MagickOffsetType
4255     progress;
4256
4257   MedianPixelList
4258     **restrict pixel_list;
4259
4260   size_t
4261     width;
4262
4263   ssize_t
4264     y;
4265
4266   /*
4267     Initialize noise image attributes.
4268   */
4269   assert(image != (Image *) NULL);
4270   assert(image->signature == MagickSignature);
4271   if (image->debug != MagickFalse)
4272     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4273   assert(exception != (ExceptionInfo *) NULL);
4274   assert(exception->signature == MagickSignature);
4275   width=GetOptimalKernelWidth2D(radius,0.5);
4276   if ((image->columns < width) || (image->rows < width))
4277     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
4278   noise_image=CloneImage(image,image->columns,image->rows,MagickTrue,
4279     exception);
4280   if (noise_image == (Image *) NULL)
4281     return((Image *) NULL);
4282   if (SetImageStorageClass(noise_image,DirectClass) == MagickFalse)
4283     {
4284       InheritException(exception,&noise_image->exception);
4285       noise_image=DestroyImage(noise_image);
4286       return((Image *) NULL);
4287     }
4288   pixel_list=AcquireMedianPixelListThreadSet(width);
4289   if (pixel_list == (MedianPixelList **) NULL)
4290     {
4291       noise_image=DestroyImage(noise_image);
4292       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4293     }
4294   /*
4295     Reduce noise image.
4296   */
4297   status=MagickTrue;
4298   progress=0;
4299   image_view=AcquireCacheView(image);
4300   noise_view=AcquireCacheView(noise_image);
4301 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4302   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4303 #endif
4304   for (y=0; y < (ssize_t) noise_image->rows; y++)
4305   {
4306     int
4307       id;
4308
4309     register const IndexPacket
4310       *restrict indexes;
4311
4312     register const PixelPacket
4313       *restrict p;
4314
4315     register IndexPacket
4316       *restrict noise_indexes;
4317
4318     register PixelPacket
4319       *restrict q;
4320
4321     register ssize_t
4322       x;
4323
4324     if (status == MagickFalse)
4325       continue;
4326     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
4327       (width/2L),image->columns+width,width,exception);
4328     q=QueueCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
4329       exception);
4330     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4331       {
4332         status=MagickFalse;
4333         continue;
4334       }
4335     indexes=GetCacheViewVirtualIndexQueue(image_view);
4336     noise_indexes=GetCacheViewAuthenticIndexQueue(noise_view);
4337     id=GetOpenMPThreadId();
4338     for (x=0; x < (ssize_t) noise_image->columns; x++)
4339     {
4340       MagickPixelPacket
4341         pixel;
4342
4343       register const PixelPacket
4344         *restrict r;
4345
4346       register const IndexPacket
4347         *restrict s;
4348
4349       register ssize_t
4350         u,
4351         v;
4352
4353       r=p;
4354       s=indexes+x;
4355       ResetMedianPixelList(pixel_list[id]);
4356       for (v=0; v < (ssize_t) width; v++)
4357       {
4358         for (u=0; u < (ssize_t) width; u++)
4359           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
4360         r+=image->columns+width;
4361         s+=image->columns+width;
4362       }
4363       pixel=GetNonpeakMedianPixelList(pixel_list[id]);
4364       SetPixelPacket(noise_image,&pixel,q,noise_indexes+x);
4365       p++;
4366       q++;
4367     }
4368     if (SyncCacheViewAuthenticPixels(noise_view,exception) == MagickFalse)
4369       status=MagickFalse;
4370     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4371       {
4372         MagickBooleanType
4373           proceed;
4374
4375 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4376   #pragma omp critical (MagickCore_ReduceNoiseImage)
4377 #endif
4378         proceed=SetImageProgress(image,ReduceNoiseImageTag,progress++,
4379           image->rows);
4380         if (proceed == MagickFalse)
4381           status=MagickFalse;
4382       }
4383   }
4384   noise_view=DestroyCacheView(noise_view);
4385   image_view=DestroyCacheView(image_view);
4386   pixel_list=DestroyMedianPixelListThreadSet(pixel_list);
4387   return(noise_image);
4388 }
4389 \f
4390 /*
4391 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4392 %                                                                             %
4393 %                                                                             %
4394 %                                                                             %
4395 %     S e l e c t i v e B l u r I m a g e                                     %
4396 %                                                                             %
4397 %                                                                             %
4398 %                                                                             %
4399 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4400 %
4401 %  SelectiveBlurImage() selectively blur pixels within a contrast threshold.
4402 %  It is similar to the unsharpen mask that sharpens everything with contrast
4403 %  above a certain threshold.
4404 %
4405 %  The format of the SelectiveBlurImage method is:
4406 %
4407 %      Image *SelectiveBlurImage(const Image *image,const double radius,
4408 %        const double sigma,const double threshold,ExceptionInfo *exception)
4409 %      Image *SelectiveBlurImageChannel(const Image *image,
4410 %        const ChannelType channel,const double radius,const double sigma,
4411 %        const double threshold,ExceptionInfo *exception)
4412 %
4413 %  A description of each parameter follows:
4414 %
4415 %    o image: the image.
4416 %
4417 %    o channel: the channel type.
4418 %
4419 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4420 %      pixel.
4421 %
4422 %    o sigma: the standard deviation of the Gaussian, in pixels.
4423 %
4424 %    o threshold: only pixels within this contrast threshold are included
4425 %      in the blur operation.
4426 %
4427 %    o exception: return any errors or warnings in this structure.
4428 %
4429 */
4430
4431 static inline MagickBooleanType SelectiveContrast(const PixelPacket *p,
4432   const PixelPacket *q,const double threshold)
4433 {
4434   if (fabs(PixelIntensity(p)-PixelIntensity(q)) < threshold)
4435     return(MagickTrue);
4436   return(MagickFalse);
4437 }
4438
4439 MagickExport Image *SelectiveBlurImage(const Image *image,const double radius,
4440   const double sigma,const double threshold,ExceptionInfo *exception)
4441 {
4442   Image
4443     *blur_image;
4444
4445   blur_image=SelectiveBlurImageChannel(image,DefaultChannels,radius,sigma,
4446     threshold,exception);
4447   return(blur_image);
4448 }
4449
4450 MagickExport Image *SelectiveBlurImageChannel(const Image *image,
4451   const ChannelType channel,const double radius,const double sigma,
4452   const double threshold,ExceptionInfo *exception)
4453 {
4454 #define SelectiveBlurImageTag  "SelectiveBlur/Image"
4455
4456   CacheView
4457     *blur_view,
4458     *image_view;
4459
4460   double
4461     *kernel;
4462
4463   Image
4464     *blur_image;
4465
4466   MagickBooleanType
4467     status;
4468
4469   MagickOffsetType
4470     progress;
4471
4472   MagickPixelPacket
4473     bias;
4474
4475   register ssize_t
4476     i;
4477
4478   size_t
4479     width;
4480
4481   ssize_t
4482     j,
4483     u,
4484     v,
4485     y;
4486
4487   /*
4488     Initialize blur image attributes.
4489   */
4490   assert(image != (Image *) NULL);
4491   assert(image->signature == MagickSignature);
4492   if (image->debug != MagickFalse)
4493     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4494   assert(exception != (ExceptionInfo *) NULL);
4495   assert(exception->signature == MagickSignature);
4496   width=GetOptimalKernelWidth1D(radius,sigma);
4497   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
4498   if (kernel == (double *) NULL)
4499     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4500   j=(ssize_t) width/2;
4501   i=0;
4502   for (v=(-j); v <= j; v++)
4503   {
4504     for (u=(-j); u <= j; u++)
4505       kernel[i++]=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
4506         (2.0*MagickPI*MagickSigma*MagickSigma);
4507   }
4508   if (image->debug != MagickFalse)
4509     {
4510       char
4511         format[MaxTextExtent],
4512         *message;
4513
4514       register const double
4515         *k;
4516
4517       ssize_t
4518         u,
4519         v;
4520
4521       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
4522         "  SelectiveBlurImage with %.20gx%.20g kernel:",(double) width,(double)
4523         width);
4524       message=AcquireString("");
4525       k=kernel;
4526       for (v=0; v < (ssize_t) width; v++)
4527       {
4528         *message='\0';
4529         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) v);
4530         (void) ConcatenateString(&message,format);
4531         for (u=0; u < (ssize_t) width; u++)
4532         {
4533           (void) FormatMagickString(format,MaxTextExtent,"%+f ",*k++);
4534           (void) ConcatenateString(&message,format);
4535         }
4536         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
4537       }
4538       message=DestroyString(message);
4539     }
4540   blur_image=CloneImage(image,0,0,MagickTrue,exception);
4541   if (blur_image == (Image *) NULL)
4542     return((Image *) NULL);
4543   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
4544     {
4545       InheritException(exception,&blur_image->exception);
4546       blur_image=DestroyImage(blur_image);
4547       return((Image *) NULL);
4548     }
4549   /*
4550     Threshold blur image.
4551   */
4552   status=MagickTrue;
4553   progress=0;
4554   GetMagickPixelPacket(image,&bias);
4555   SetMagickPixelPacketBias(image,&bias);
4556   image_view=AcquireCacheView(image);
4557   blur_view=AcquireCacheView(blur_image);
4558 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4559   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4560 #endif
4561   for (y=0; y < (ssize_t) image->rows; y++)
4562   {
4563     MagickBooleanType
4564       sync;
4565
4566     MagickRealType
4567       gamma;
4568
4569     register const IndexPacket
4570       *restrict indexes;
4571
4572     register const PixelPacket
4573       *restrict p;
4574
4575     register IndexPacket
4576       *restrict blur_indexes;
4577
4578     register PixelPacket
4579       *restrict q;
4580
4581     register ssize_t
4582       x;
4583
4584     if (status == MagickFalse)
4585       continue;
4586     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
4587       (width/2L),image->columns+width,width,exception);
4588     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
4589       exception);
4590     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4591       {
4592         status=MagickFalse;
4593         continue;
4594       }
4595     indexes=GetCacheViewVirtualIndexQueue(image_view);
4596     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
4597     for (x=0; x < (ssize_t) image->columns; x++)
4598     {
4599       MagickPixelPacket
4600         pixel;
4601
4602       register const double
4603         *restrict k;
4604
4605       register ssize_t
4606         u;
4607
4608       ssize_t
4609         j,
4610         v;
4611
4612       pixel=bias;
4613       k=kernel;
4614       gamma=0.0;
4615       j=0;
4616       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
4617         {
4618           for (v=0; v < (ssize_t) width; v++)
4619           {
4620             for (u=0; u < (ssize_t) width; u++)
4621             {
4622               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4623                 {
4624                   pixel.red+=(*k)*(p+u+j)->red;
4625                   pixel.green+=(*k)*(p+u+j)->green;
4626                   pixel.blue+=(*k)*(p+u+j)->blue;
4627                   gamma+=(*k);
4628                   k++;
4629                 }
4630             }
4631             j+=(ssize_t) (image->columns+width);
4632           }
4633           if (gamma != 0.0)
4634             {
4635               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4636               if ((channel & RedChannel) != 0)
4637                 q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
4638               if ((channel & GreenChannel) != 0)
4639                 q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
4640               if ((channel & BlueChannel) != 0)
4641                 q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
4642             }
4643           if ((channel & OpacityChannel) != 0)
4644             {
4645               gamma=0.0;
4646               j=0;
4647               for (v=0; v < (ssize_t) width; v++)
4648               {
4649                 for (u=0; u < (ssize_t) width; u++)
4650                 {
4651                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4652                     {
4653                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4654                       gamma+=(*k);
4655                       k++;
4656                     }
4657                 }
4658                 j+=(ssize_t) (image->columns+width);
4659               }
4660               if (gamma != 0.0)
4661                 {
4662                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4663                     gamma);
4664                   SetOpacityPixelComponent(q,ClampToQuantum(gamma*
4665                     GetOpacityPixelComponent(&pixel)));
4666                 }
4667             }
4668           if (((channel & IndexChannel) != 0) &&
4669               (image->colorspace == CMYKColorspace))
4670             {
4671               gamma=0.0;
4672               j=0;
4673               for (v=0; v < (ssize_t) width; v++)
4674               {
4675                 for (u=0; u < (ssize_t) width; u++)
4676                 {
4677                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4678                     {
4679                       pixel.index+=(*k)*indexes[x+u+j];
4680                       gamma+=(*k);
4681                       k++;
4682                     }
4683                 }
4684                 j+=(ssize_t) (image->columns+width);
4685               }
4686               if (gamma != 0.0)
4687                 {
4688                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4689                     gamma);
4690                   blur_indexes[x]=ClampToQuantum(gamma*
4691                     GetIndexPixelComponent(&pixel));
4692                 }
4693             }
4694         }
4695       else
4696         {
4697           MagickRealType
4698             alpha;
4699
4700           for (v=0; v < (ssize_t) width; v++)
4701           {
4702             for (u=0; u < (ssize_t) width; u++)
4703             {
4704               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4705                 {
4706                   alpha=(MagickRealType) (QuantumScale*
4707                     GetAlphaPixelComponent(p+u+j));
4708                   pixel.red+=(*k)*alpha*(p+u+j)->red;
4709                   pixel.green+=(*k)*alpha*(p+u+j)->green;
4710                   pixel.blue+=(*k)*alpha*(p+u+j)->blue;
4711                   pixel.opacity+=(*k)*(p+u+j)->opacity;
4712                   gamma+=(*k)*alpha;
4713                   k++;
4714                 }
4715             }
4716             j+=(ssize_t) (image->columns+width);
4717           }
4718           if (gamma != 0.0)
4719             {
4720               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4721               if ((channel & RedChannel) != 0)
4722                 q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
4723               if ((channel & GreenChannel) != 0)
4724                 q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
4725               if ((channel & BlueChannel) != 0)
4726                 q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
4727             }
4728           if ((channel & OpacityChannel) != 0)
4729             {
4730               gamma=0.0;
4731               j=0;
4732               for (v=0; v < (ssize_t) width; v++)
4733               {
4734                 for (u=0; u < (ssize_t) width; u++)
4735                 {
4736                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4737                     {
4738                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4739                       gamma+=(*k);
4740                       k++;
4741                     }
4742                 }
4743                 j+=(ssize_t) (image->columns+width);
4744               }
4745               if (gamma != 0.0)
4746                 {
4747                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4748                     gamma);
4749                   SetOpacityPixelComponent(q,
4750                     ClampOpacityPixelComponent(&pixel));
4751                 }
4752             }
4753           if (((channel & IndexChannel) != 0) &&
4754               (image->colorspace == CMYKColorspace))
4755             {
4756               gamma=0.0;
4757               j=0;
4758               for (v=0; v < (ssize_t) width; v++)
4759               {
4760                 for (u=0; u < (ssize_t) width; u++)
4761                 {
4762                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4763                     {
4764                       alpha=(MagickRealType) (QuantumScale*
4765                         GetAlphaPixelComponent(p+u+j));
4766                       pixel.index+=(*k)*alpha*indexes[x+u+j];
4767                       gamma+=(*k);
4768                       k++;
4769                     }
4770                 }
4771                 j+=(ssize_t) (image->columns+width);
4772               }
4773               if (gamma != 0.0)
4774                 {
4775                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4776                     gamma);
4777                   blur_indexes[x]=ClampToQuantum(gamma*
4778                     GetIndexPixelComponent(&pixel));
4779                 }
4780             }
4781         }
4782       p++;
4783       q++;
4784     }
4785     sync=SyncCacheViewAuthenticPixels(blur_view,exception);
4786     if (sync == MagickFalse)
4787       status=MagickFalse;
4788     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4789       {
4790         MagickBooleanType
4791           proceed;
4792
4793 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4794   #pragma omp critical (MagickCore_SelectiveBlurImageChannel)
4795 #endif
4796         proceed=SetImageProgress(image,SelectiveBlurImageTag,progress++,
4797           image->rows);
4798         if (proceed == MagickFalse)
4799           status=MagickFalse;
4800       }
4801   }
4802   blur_image->type=image->type;
4803   blur_view=DestroyCacheView(blur_view);
4804   image_view=DestroyCacheView(image_view);
4805   kernel=(double *) RelinquishMagickMemory(kernel);
4806   if (status == MagickFalse)
4807     blur_image=DestroyImage(blur_image);
4808   return(blur_image);
4809 }
4810 \f
4811 /*
4812 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4813 %                                                                             %
4814 %                                                                             %
4815 %                                                                             %
4816 %     S h a d e I m a g e                                                     %
4817 %                                                                             %
4818 %                                                                             %
4819 %                                                                             %
4820 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4821 %
4822 %  ShadeImage() shines a distant light on an image to create a
4823 %  three-dimensional effect. You control the positioning of the light with
4824 %  azimuth and elevation; azimuth is measured in degrees off the x axis
4825 %  and elevation is measured in pixels above the Z axis.
4826 %
4827 %  The format of the ShadeImage method is:
4828 %
4829 %      Image *ShadeImage(const Image *image,const MagickBooleanType gray,
4830 %        const double azimuth,const double elevation,ExceptionInfo *exception)
4831 %
4832 %  A description of each parameter follows:
4833 %
4834 %    o image: the image.
4835 %
4836 %    o gray: A value other than zero shades the intensity of each pixel.
4837 %
4838 %    o azimuth, elevation:  Define the light source direction.
4839 %
4840 %    o exception: return any errors or warnings in this structure.
4841 %
4842 */
4843 MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
4844   const double azimuth,const double elevation,ExceptionInfo *exception)
4845 {
4846 #define ShadeImageTag  "Shade/Image"
4847
4848   CacheView
4849     *image_view,
4850     *shade_view;
4851
4852   Image
4853     *shade_image;
4854
4855   MagickBooleanType
4856     status;
4857
4858   MagickOffsetType
4859     progress;
4860
4861   PrimaryInfo
4862     light;
4863
4864   ssize_t
4865     y;
4866
4867   /*
4868     Initialize shaded image attributes.
4869   */
4870   assert(image != (const Image *) NULL);
4871   assert(image->signature == MagickSignature);
4872   if (image->debug != MagickFalse)
4873     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4874   assert(exception != (ExceptionInfo *) NULL);
4875   assert(exception->signature == MagickSignature);
4876   shade_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
4877   if (shade_image == (Image *) NULL)
4878     return((Image *) NULL);
4879   if (SetImageStorageClass(shade_image,DirectClass) == MagickFalse)
4880     {
4881       InheritException(exception,&shade_image->exception);
4882       shade_image=DestroyImage(shade_image);
4883       return((Image *) NULL);
4884     }
4885   /*
4886     Compute the light vector.
4887   */
4888   light.x=(double) QuantumRange*cos(DegreesToRadians(azimuth))*
4889     cos(DegreesToRadians(elevation));
4890   light.y=(double) QuantumRange*sin(DegreesToRadians(azimuth))*
4891     cos(DegreesToRadians(elevation));
4892   light.z=(double) QuantumRange*sin(DegreesToRadians(elevation));
4893   /*
4894     Shade image.
4895   */
4896   status=MagickTrue;
4897   progress=0;
4898   image_view=AcquireCacheView(image);
4899   shade_view=AcquireCacheView(shade_image);
4900 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4901   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4902 #endif
4903   for (y=0; y < (ssize_t) image->rows; y++)
4904   {
4905     MagickRealType
4906       distance,
4907       normal_distance,
4908       shade;
4909
4910     PrimaryInfo
4911       normal;
4912
4913     register const PixelPacket
4914       *restrict p,
4915       *restrict s0,
4916       *restrict s1,
4917       *restrict s2;
4918
4919     register PixelPacket
4920       *restrict q;
4921
4922     register ssize_t
4923       x;
4924
4925     if (status == MagickFalse)
4926       continue;
4927     p=GetCacheViewVirtualPixels(image_view,-1,y-1,image->columns+2,3,exception);
4928     q=QueueCacheViewAuthenticPixels(shade_view,0,y,shade_image->columns,1,
4929       exception);
4930     if ((p == (PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4931       {
4932         status=MagickFalse;
4933         continue;
4934       }
4935     /*
4936       Shade this row of pixels.
4937     */
4938     normal.z=2.0*(double) QuantumRange;  /* constant Z of surface normal */
4939     s0=p+1;
4940     s1=s0+image->columns+2;
4941     s2=s1+image->columns+2;
4942     for (x=0; x < (ssize_t) image->columns; x++)
4943     {
4944       /*
4945         Determine the surface normal and compute shading.
4946       */
4947       normal.x=(double) (PixelIntensity(s0-1)+PixelIntensity(s1-1)+
4948         PixelIntensity(s2-1)-PixelIntensity(s0+1)-PixelIntensity(s1+1)-
4949         PixelIntensity(s2+1));
4950       normal.y=(double) (PixelIntensity(s2-1)+PixelIntensity(s2)+
4951         PixelIntensity(s2+1)-PixelIntensity(s0-1)-PixelIntensity(s0)-
4952         PixelIntensity(s0+1));
4953       if ((normal.x == 0.0) && (normal.y == 0.0))
4954         shade=light.z;
4955       else
4956         {
4957           shade=0.0;
4958           distance=normal.x*light.x+normal.y*light.y+normal.z*light.z;
4959           if (distance > MagickEpsilon)
4960             {
4961               normal_distance=
4962                 normal.x*normal.x+normal.y*normal.y+normal.z*normal.z;
4963               if (normal_distance > (MagickEpsilon*MagickEpsilon))
4964                 shade=distance/sqrt((double) normal_distance);
4965             }
4966         }
4967       if (gray != MagickFalse)
4968         {
4969           q->red=(Quantum) shade;
4970           q->green=(Quantum) shade;
4971           q->blue=(Quantum) shade;
4972         }
4973       else
4974         {
4975           q->red=ClampToQuantum(QuantumScale*shade*s1->red);
4976           q->green=ClampToQuantum(QuantumScale*shade*s1->green);
4977           q->blue=ClampToQuantum(QuantumScale*shade*s1->blue);
4978         }
4979       q->opacity=s1->opacity;
4980       s0++;
4981       s1++;
4982       s2++;
4983       q++;
4984     }
4985     if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
4986       status=MagickFalse;
4987     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4988       {
4989         MagickBooleanType
4990           proceed;
4991
4992 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4993   #pragma omp critical (MagickCore_ShadeImage)
4994 #endif
4995         proceed=SetImageProgress(image,ShadeImageTag,progress++,image->rows);
4996         if (proceed == MagickFalse)
4997           status=MagickFalse;
4998       }
4999   }
5000   shade_view=DestroyCacheView(shade_view);
5001   image_view=DestroyCacheView(image_view);
5002   if (status == MagickFalse)
5003     shade_image=DestroyImage(shade_image);
5004   return(shade_image);
5005 }
5006 \f
5007 /*
5008 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5009 %                                                                             %
5010 %                                                                             %
5011 %                                                                             %
5012 %     S h a r p e n I m a g e                                                 %
5013 %                                                                             %
5014 %                                                                             %
5015 %                                                                             %
5016 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5017 %
5018 %  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
5019 %  operator of the given radius and standard deviation (sigma).  For
5020 %  reasonable results, radius should be larger than sigma.  Use a radius of 0
5021 %  and SharpenImage() selects a suitable radius for you.
5022 %
5023 %  Using a separable kernel would be faster, but the negative weights cancel
5024 %  out on the corners of the kernel producing often undesirable ringing in the
5025 %  filtered result; this can be avoided by using a 2D gaussian shaped image
5026 %  sharpening kernel instead.
5027 %
5028 %  The format of the SharpenImage method is:
5029 %
5030 %    Image *SharpenImage(const Image *image,const double radius,
5031 %      const double sigma,ExceptionInfo *exception)
5032 %    Image *SharpenImageChannel(const Image *image,const ChannelType channel,
5033 %      const double radius,const double sigma,ExceptionInfo *exception)
5034 %
5035 %  A description of each parameter follows:
5036 %
5037 %    o image: the image.
5038 %
5039 %    o channel: the channel type.
5040 %
5041 %    o radius: the radius of the Gaussian, in pixels, not counting the center
5042 %      pixel.
5043 %
5044 %    o sigma: the standard deviation of the Laplacian, in pixels.
5045 %
5046 %    o exception: return any errors or warnings in this structure.
5047 %
5048 */
5049
5050 MagickExport Image *SharpenImage(const Image *image,const double radius,
5051   const double sigma,ExceptionInfo *exception)
5052 {
5053   Image
5054     *sharp_image;
5055
5056   sharp_image=SharpenImageChannel(image,DefaultChannels,radius,sigma,exception);
5057   return(sharp_image);
5058 }
5059
5060 MagickExport Image *SharpenImageChannel(const Image *image,
5061   const ChannelType channel,const double radius,const double sigma,
5062   ExceptionInfo *exception)
5063 {
5064   double
5065     *kernel,
5066     normalize;
5067
5068   Image
5069     *sharp_image;
5070
5071   register ssize_t
5072     i;
5073
5074   size_t
5075     width;
5076
5077   ssize_t
5078     j,
5079     u,
5080     v;
5081
5082   assert(image != (const Image *) NULL);
5083   assert(image->signature == MagickSignature);
5084   if (image->debug != MagickFalse)
5085     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5086   assert(exception != (ExceptionInfo *) NULL);
5087   assert(exception->signature == MagickSignature);
5088   width=GetOptimalKernelWidth2D(radius,sigma);
5089   kernel=(double *) AcquireQuantumMemory((size_t) width*width,sizeof(*kernel));
5090   if (kernel == (double *) NULL)
5091     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
5092   normalize=0.0;
5093   j=(ssize_t) width/2;
5094   i=0;
5095   for (v=(-j); v <= j; v++)
5096   {
5097     for (u=(-j); u <= j; u++)
5098     {
5099       kernel[i]=(-exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
5100         (2.0*MagickPI*MagickSigma*MagickSigma));
5101       normalize+=kernel[i];
5102       i++;
5103     }
5104   }
5105   kernel[i/2]=(double) ((-2.0)*normalize);
5106   sharp_image=ConvolveImageChannel(image,channel,width,kernel,exception);
5107   kernel=(double *) RelinquishMagickMemory(kernel);
5108   return(sharp_image);
5109 }
5110 \f
5111 /*
5112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5113 %                                                                             %
5114 %                                                                             %
5115 %                                                                             %
5116 %     S p r e a d I m a g e                                                   %
5117 %                                                                             %
5118 %                                                                             %
5119 %                                                                             %
5120 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5121 %
5122 %  SpreadImage() is a special effects method that randomly displaces each
5123 %  pixel in a block defined by the radius parameter.
5124 %
5125 %  The format of the SpreadImage method is:
5126 %
5127 %      Image *SpreadImage(const Image *image,const double radius,
5128 %        ExceptionInfo *exception)
5129 %
5130 %  A description of each parameter follows:
5131 %
5132 %    o image: the image.
5133 %
5134 %    o radius:  Choose a random pixel in a neighborhood of this extent.
5135 %
5136 %    o exception: return any errors or warnings in this structure.
5137 %
5138 */
5139 MagickExport Image *SpreadImage(const Image *image,const double radius,
5140   ExceptionInfo *exception)
5141 {
5142 #define SpreadImageTag  "Spread/Image"
5143
5144   CacheView
5145     *image_view;
5146
5147   Image
5148     *spread_image;
5149
5150   MagickBooleanType
5151     status;
5152
5153   MagickOffsetType
5154     progress;
5155
5156   MagickPixelPacket
5157     bias;
5158
5159   RandomInfo
5160     **restrict random_info;
5161
5162   ResampleFilter
5163     **restrict resample_filter;
5164
5165   size_t
5166     width;
5167
5168   ssize_t
5169     y;
5170
5171   /*
5172     Initialize spread image attributes.
5173   */
5174   assert(image != (Image *) NULL);
5175   assert(image->signature == MagickSignature);
5176   if (image->debug != MagickFalse)
5177     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5178   assert(exception != (ExceptionInfo *) NULL);
5179   assert(exception->signature == MagickSignature);
5180   spread_image=CloneImage(image,image->columns,image->rows,MagickTrue,
5181     exception);
5182   if (spread_image == (Image *) NULL)
5183     return((Image *) NULL);
5184   if (SetImageStorageClass(spread_image,DirectClass) == MagickFalse)
5185     {
5186       InheritException(exception,&spread_image->exception);
5187       spread_image=DestroyImage(spread_image);
5188       return((Image *) NULL);
5189     }
5190   /*
5191     Spread image.
5192   */
5193   status=MagickTrue;
5194   progress=0;
5195   GetMagickPixelPacket(spread_image,&bias);
5196   width=GetOptimalKernelWidth1D(radius,0.5);
5197   resample_filter=AcquireResampleFilterThreadSet(image,
5198     UndefinedVirtualPixelMethod,MagickTrue,exception);
5199   random_info=AcquireRandomInfoThreadSet();
5200   image_view=AcquireCacheView(spread_image);
5201 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5202   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
5203 #endif
5204   for (y=0; y < (ssize_t) spread_image->rows; y++)
5205   {
5206     int
5207       id;
5208
5209     MagickPixelPacket
5210       pixel;
5211
5212     register IndexPacket
5213       *restrict indexes;
5214
5215     register PixelPacket
5216       *restrict q;
5217
5218     register ssize_t
5219       x;
5220
5221     if (status == MagickFalse)
5222       continue;
5223     q=QueueCacheViewAuthenticPixels(image_view,0,y,spread_image->columns,1,
5224       exception);
5225     if (q == (PixelPacket *) NULL)
5226       {
5227         status=MagickFalse;
5228         continue;
5229       }
5230     indexes=GetCacheViewAuthenticIndexQueue(image_view);
5231     pixel=bias;
5232     id=GetOpenMPThreadId();
5233     for (x=0; x < (ssize_t) spread_image->columns; x++)
5234     {
5235       (void) ResamplePixelColor(resample_filter[id],(double) x+width*
5236         (GetPseudoRandomValue(random_info[id])-0.5),(double) y+width*
5237         (GetPseudoRandomValue(random_info[id])-0.5),&pixel);
5238       SetPixelPacket(spread_image,&pixel,q,indexes+x);
5239       q++;
5240     }
5241     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
5242       status=MagickFalse;
5243     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5244       {
5245         MagickBooleanType
5246           proceed;
5247
5248 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5249   #pragma omp critical (MagickCore_SpreadImage)
5250 #endif
5251         proceed=SetImageProgress(image,SpreadImageTag,progress++,image->rows);
5252         if (proceed == MagickFalse)
5253           status=MagickFalse;
5254       }
5255   }
5256   image_view=DestroyCacheView(image_view);
5257   random_info=DestroyRandomInfoThreadSet(random_info);
5258   resample_filter=DestroyResampleFilterThreadSet(resample_filter);
5259   return(spread_image);
5260 }
5261 \f
5262 /*
5263 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5264 %                                                                             %
5265 %                                                                             %
5266 %                                                                             %
5267 %     U n s h a r p M a s k I m a g e                                         %
5268 %                                                                             %
5269 %                                                                             %
5270 %                                                                             %
5271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5272 %
5273 %  UnsharpMaskImage() sharpens one or more image channels.  We convolve the
5274 %  image with a Gaussian operator of the given radius and standard deviation
5275 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
5276 %  radius of 0 and UnsharpMaskImage() selects a suitable radius for you.
5277 %
5278 %  The format of the UnsharpMaskImage method is:
5279 %
5280 %    Image *UnsharpMaskImage(const Image *image,const double radius,
5281 %      const double sigma,const double amount,const double threshold,
5282 %      ExceptionInfo *exception)
5283 %    Image *UnsharpMaskImageChannel(const Image *image,
5284 %      const ChannelType channel,const double radius,const double sigma,
5285 %      const double amount,const double threshold,ExceptionInfo *exception)
5286 %
5287 %  A description of each parameter follows:
5288 %
5289 %    o image: the image.
5290 %
5291 %    o channel: the channel type.
5292 %
5293 %    o radius: the radius of the Gaussian, in pixels, not counting the center
5294 %      pixel.
5295 %
5296 %    o sigma: the standard deviation of the Gaussian, in pixels.
5297 %
5298 %    o amount: the percentage of the difference between the original and the
5299 %      blur image that is added back into the original.
5300 %
5301 %    o threshold: the threshold in pixels needed to apply the diffence amount.
5302 %
5303 %    o exception: return any errors or warnings in this structure.
5304 %
5305 */
5306
5307 MagickExport Image *UnsharpMaskImage(const Image *image,const double radius,
5308   const double sigma,const double amount,const double threshold,
5309   ExceptionInfo *exception)
5310 {
5311   Image
5312     *sharp_image;
5313
5314   sharp_image=UnsharpMaskImageChannel(image,DefaultChannels,radius,sigma,amount,
5315     threshold,exception);
5316   return(sharp_image);
5317 }
5318
5319 MagickExport Image *UnsharpMaskImageChannel(const Image *image,
5320   const ChannelType channel,const double radius,const double sigma,
5321   const double amount,const double threshold,ExceptionInfo *exception)
5322 {
5323 #define SharpenImageTag  "Sharpen/Image"
5324
5325   CacheView
5326     *image_view,
5327     *unsharp_view;
5328
5329   Image
5330     *unsharp_image;
5331
5332   MagickBooleanType
5333     status;
5334
5335   MagickOffsetType
5336     progress;
5337
5338   MagickPixelPacket
5339     bias;
5340
5341   MagickRealType
5342     quantum_threshold;
5343
5344   ssize_t
5345     y;
5346
5347   assert(image != (const Image *) NULL);
5348   assert(image->signature == MagickSignature);
5349   if (image->debug != MagickFalse)
5350     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5351   assert(exception != (ExceptionInfo *) NULL);
5352   unsharp_image=BlurImageChannel(image,channel,radius,sigma,exception);
5353   if (unsharp_image == (Image *) NULL)
5354     return((Image *) NULL);
5355   quantum_threshold=(MagickRealType) QuantumRange*threshold;
5356   /*
5357     Unsharp-mask image.
5358   */
5359   status=MagickTrue;
5360   progress=0;
5361   GetMagickPixelPacket(image,&bias);
5362   image_view=AcquireCacheView(image);
5363   unsharp_view=AcquireCacheView(unsharp_image);
5364 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5365   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
5366 #endif
5367   for (y=0; y < (ssize_t) image->rows; y++)
5368   {
5369     MagickPixelPacket
5370       pixel;
5371
5372     register const IndexPacket
5373       *restrict indexes;
5374
5375     register const PixelPacket
5376       *restrict p;
5377
5378     register IndexPacket
5379       *restrict unsharp_indexes;
5380
5381     register PixelPacket
5382       *restrict q;
5383
5384     register ssize_t
5385       x;
5386
5387     if (status == MagickFalse)
5388       continue;
5389     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
5390     q=GetCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
5391       exception);
5392     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
5393       {
5394         status=MagickFalse;
5395         continue;
5396       }
5397     indexes=GetCacheViewVirtualIndexQueue(image_view);
5398     unsharp_indexes=GetCacheViewAuthenticIndexQueue(unsharp_view);
5399     pixel=bias;
5400     for (x=0; x < (ssize_t) image->columns; x++)
5401     {
5402       if ((channel & RedChannel) != 0)
5403         {
5404           pixel.red=p->red-(MagickRealType) q->red;
5405           if (fabs(2.0*pixel.red) < quantum_threshold)
5406             pixel.red=(MagickRealType) GetRedPixelComponent(p);
5407           else
5408             pixel.red=(MagickRealType) p->red+(pixel.red*amount);
5409           SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
5410         }
5411       if ((channel & GreenChannel) != 0)
5412         {
5413           pixel.green=p->green-(MagickRealType) q->green;
5414           if (fabs(2.0*pixel.green) < quantum_threshold)
5415             pixel.green=(MagickRealType) GetGreenPixelComponent(p);
5416           else
5417             pixel.green=(MagickRealType) p->green+(pixel.green*amount);
5418           SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
5419         }
5420       if ((channel & BlueChannel) != 0)
5421         {
5422           pixel.blue=p->blue-(MagickRealType) q->blue;
5423           if (fabs(2.0*pixel.blue) < quantum_threshold)
5424             pixel.blue=(MagickRealType) GetBluePixelComponent(p);
5425           else
5426             pixel.blue=(MagickRealType) p->blue+(pixel.blue*amount);
5427           SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
5428         }
5429       if ((channel & OpacityChannel) != 0)
5430         {
5431           pixel.opacity=p->opacity-(MagickRealType) q->opacity;
5432           if (fabs(2.0*pixel.opacity) < quantum_threshold)
5433             pixel.opacity=(MagickRealType) GetOpacityPixelComponent(p);
5434           else
5435             pixel.opacity=p->opacity+(pixel.opacity*amount);
5436           SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
5437         }
5438       if (((channel & IndexChannel) != 0) &&
5439           (image->colorspace == CMYKColorspace))
5440         {
5441           pixel.index=unsharp_indexes[x]-(MagickRealType) indexes[x];
5442           if (fabs(2.0*pixel.index) < quantum_threshold)
5443             pixel.index=(MagickRealType) unsharp_indexes[x];
5444           else
5445             pixel.index=(MagickRealType) unsharp_indexes[x]+(pixel.index*
5446               amount);
5447           unsharp_indexes[x]=ClampToQuantum(pixel.index);
5448         }
5449       p++;
5450       q++;
5451     }
5452     if (SyncCacheViewAuthenticPixels(unsharp_view,exception) == MagickFalse)
5453       status=MagickFalse;
5454     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5455       {
5456         MagickBooleanType
5457           proceed;
5458
5459 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5460   #pragma omp critical (MagickCore_UnsharpMaskImageChannel)
5461 #endif
5462         proceed=SetImageProgress(image,SharpenImageTag,progress++,image->rows);
5463         if (proceed == MagickFalse)
5464           status=MagickFalse;
5465       }
5466   }
5467   unsharp_image->type=image->type;
5468   unsharp_view=DestroyCacheView(unsharp_view);
5469   image_view=DestroyCacheView(image_view);
5470   if (status == MagickFalse)
5471     unsharp_image=DestroyImage(unsharp_image);
5472   return(unsharp_image);
5473 }