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