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

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