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

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