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