granicus.if.org Git - imagemagick/blob - MagickCore/effect.c

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