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