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