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