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