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

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