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

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