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

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