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