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