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