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