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

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