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