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