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