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