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