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

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