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