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