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

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