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