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