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