granicus.if.org Git - imagemagick/blob - magick/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-2010 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 "magick/studio.h"
  44 #include "magick/property.h"
  45 #include "magick/blob.h"
  46 #include "magick/cache-view.h"
  47 #include "magick/color.h"
  48 #include "magick/color-private.h"
  49 #include "magick/colorspace.h"
  50 #include "magick/constitute.h"
  51 #include "magick/decorate.h"
  52 #include "magick/draw.h"
  53 #include "magick/enhance.h"
  54 #include "magick/exception.h"
  55 #include "magick/exception-private.h"
  56 #include "magick/effect.h"
  57 #include "magick/fx.h"
  58 #include "magick/gem.h"
  59 #include "magick/geometry.h"
  60 #include "magick/image-private.h"
  61 #include "magick/list.h"
  62 #include "magick/log.h"
  63 #include "magick/memory_.h"
  64 #include "magick/monitor.h"
  65 #include "magick/monitor-private.h"
  66 #include "magick/montage.h"
  67 #include "magick/paint.h"
  68 #include "magick/pixel-private.h"
  69 #include "magick/property.h"
  70 #include "magick/quantize.h"
  71 #include "magick/quantum.h"
  72 #include "magick/random_.h"
  73 #include "magick/random-private.h"
  74 #include "magick/resample.h"
  75 #include "magick/resample-private.h"
  76 #include "magick/resize.h"
  77 #include "magick/resource_.h"
  78 #include "magick/segment.h"
  79 #include "magick/shear.h"
  80 #include "magick/signature-private.h"
  81 #include "magick/string_.h"
  82 #include "magick/thread-private.h"
  83 #include "magick/transform.h"
  84 #include "magick/threshold.h"
  85 \f
  86 /*
  87 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  88 %                                                                             %
  89 %                                                                             %
  90 %                                                                             %
  91 %     A d a p t i v e B l u r I m a g e                                       %
  92 %                                                                             %
  93 %                                                                             %
  94 %                                                                             %
  95 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  96 %
  97 %  AdaptiveBlurImage() adaptively blurs the image by blurring less
  98 %  intensely near image edges and more intensely far from edges.  We blur the
  99 %  image with a Gaussian operator of the given radius and standard deviation
 100 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
 101 %  radius of 0 and AdaptiveBlurImage() selects a suitable radius for you.
 102 %
 103 %  The format of the AdaptiveBlurImage method is:
 104 %
 105 %      Image *AdaptiveBlurImage(const Image *image,const double radius,
 106 %        const double sigma,ExceptionInfo *exception)
 107 %      Image *AdaptiveBlurImageChannel(const Image *image,
 108 %        const ChannelType channel,double radius,const double sigma,
 109 %        ExceptionInfo *exception)
 110 %
 111 %  A description of each parameter follows:
 112 %
 113 %    o image: the image.
 114 %
 115 %    o channel: the channel type.
 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 exception: return any errors or warnings in this structure.
 123 %
 124 */
 125
 126 MagickExport Image *AdaptiveBlurImage(const Image *image,const double radius,
 127   const double sigma,ExceptionInfo *exception)
 128 {
 129   Image
 130     *blur_image;
 131
 132   blur_image=AdaptiveBlurImageChannel(image,DefaultChannels,radius,sigma,
 133     exception);
 134   return(blur_image);
 135 }
 136
 137 MagickExport Image *AdaptiveBlurImageChannel(const Image *image,
 138   const ChannelType channel,const double radius,const double sigma,
 139   ExceptionInfo *exception)
 140 {
 141 #define AdaptiveBlurImageTag  "Convolve/Image"
 142 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
 143
 144   double
 145     **kernel,
 146     normalize;
 147
 148   Image
 149     *blur_image,
 150     *edge_image,
 151     *gaussian_image;
 152
 153   long
 154     j,
 155     k,
 156     progress,
 157     u,
 158     v,
 159     y;
 160
 161   MagickBooleanType
 162     status;
 163
 164   MagickPixelPacket
 165     bias;
 166
 167   register long
 168     i;
 169
 170   unsigned long
 171     width;
 172
 173   CacheView
 174     *blur_view,
 175     *edge_view,
 176     *image_view;
 177
 178   assert(image != (const Image *) NULL);
 179   assert(image->signature == MagickSignature);
 180   if (image->debug != MagickFalse)
 181     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 182   assert(exception != (ExceptionInfo *) NULL);
 183   assert(exception->signature == MagickSignature);
 184   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
 185   if (blur_image == (Image *) NULL)
 186     return((Image *) NULL);
 187   if (fabs(sigma) <= MagickEpsilon)
 188     return(blur_image);
 189   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
 190     {
 191       InheritException(exception,&blur_image->exception);
 192       blur_image=DestroyImage(blur_image);
 193       return((Image *) NULL);
 194     }
 195   /*
 196     Edge detect the image brighness channel, level, blur, and level again.
 197   */
 198   edge_image=EdgeImage(image,radius,exception);
 199   if (edge_image == (Image *) NULL)
 200     {
 201       blur_image=DestroyImage(blur_image);
 202       return((Image *) NULL);
 203     }
 204   (void) LevelImage(edge_image,"20%,95%");
 205   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
 206   if (gaussian_image != (Image *) NULL)
 207     {
 208       edge_image=DestroyImage(edge_image);
 209       edge_image=gaussian_image;
 210     }
 211   (void) LevelImage(edge_image,"10%,95%");
 212   /*
 213     Create a set of kernels from maximum (radius,sigma) to minimum.
 214   */
 215   width=GetOptimalKernelWidth2D(radius,sigma);
 216   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 217   if (kernel == (double **) NULL)
 218     {
 219       edge_image=DestroyImage(edge_image);
 220       blur_image=DestroyImage(blur_image);
 221       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 222     }
 223   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
 224   for (i=0; i < (long) width; i+=2)
 225   {
 226     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
 227       sizeof(**kernel));
 228     if (kernel[i] == (double *) NULL)
 229       break;
 230     normalize=0.0;
 231     j=(long) (width-i)/2;
 232     k=0;
 233     for (v=(-j); v <= j; v++)
 234     {
 235       for (u=(-j); u <= j; u++)
 236       {
 237         kernel[i][k]=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
 238           (2.0*MagickPI*MagickSigma*MagickSigma);
 239         normalize+=kernel[i][k];
 240         k++;
 241       }
 242     }
 243     if (fabs(normalize) <= MagickEpsilon)
 244       normalize=1.0;
 245     normalize=1.0/normalize;
 246     for (k=0; k < (j*j); k++)
 247       kernel[i][k]=normalize*kernel[i][k];
 248   }
 249   if (i < (long) width)
 250     {
 251       for (i-=2; i >= 0; i-=2)
 252         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 253       kernel=(double **) RelinquishMagickMemory(kernel);
 254       edge_image=DestroyImage(edge_image);
 255       blur_image=DestroyImage(blur_image);
 256       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 257     }
 258   /*
 259     Adaptively blur image.
 260   */
 261   status=MagickTrue;
 262   progress=0;
 263   GetMagickPixelPacket(image,&bias);
 264   SetMagickPixelPacketBias(image,&bias);
 265   image_view=AcquireCacheView(image);
 266   edge_view=AcquireCacheView(edge_image);
 267   blur_view=AcquireCacheView(blur_image);
 268 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 269   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 270 #endif
 271   for (y=0; y < (long) blur_image->rows; y++)
 272   {
 273     register const IndexPacket
 274       *restrict indexes;
 275
 276     register const PixelPacket
 277       *restrict p,
 278       *restrict r;
 279
 280     register IndexPacket
 281       *restrict blur_indexes;
 282
 283     register long
 284       x;
 285
 286     register PixelPacket
 287       *restrict q;
 288
 289     if (status == MagickFalse)
 290       continue;
 291     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
 292     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
 293       exception);
 294     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 295       {
 296         status=MagickFalse;
 297         continue;
 298       }
 299     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
 300     for (x=0; x < (long) blur_image->columns; x++)
 301     {
 302       MagickPixelPacket
 303         pixel;
 304
 305       MagickRealType
 306         alpha,
 307         gamma;
 308
 309       register const double
 310         *restrict k;
 311
 312       register long
 313         i,
 314         u,
 315         v;
 316
 317       gamma=0.0;
 318       i=(long) (width*QuantumScale*PixelIntensity(r)+0.5);
 319       if (i < 0)
 320         i=0;
 321       else
 322         if (i > (long) width)
 323           i=(long) width;
 324       if ((i & 0x01) != 0)
 325         i--;
 326       p=GetCacheViewVirtualPixels(image_view,x-((long) (width-i)/2L),y-(long)
 327         ((width-i)/2L),width-i,width-i,exception);
 328       if (p == (const PixelPacket *) NULL)
 329         break;
 330       indexes=GetCacheViewVirtualIndexQueue(image_view);
 331       pixel=bias;
 332       k=kernel[i];
 333       for (v=0; v < (long) (width-i); v++)
 334       {
 335         for (u=0; u < (long) (width-i); u++)
 336         {
 337           alpha=1.0;
 338           if (((channel & OpacityChannel) != 0) &&
 339               (image->matte != MagickFalse))
 340             alpha=(MagickRealType) (QuantumScale*(QuantumRange-p->opacity));
 341           if ((channel & RedChannel) != 0)
 342             pixel.red+=(*k)*alpha*p->red;
 343           if ((channel & GreenChannel) != 0)
 344             pixel.green+=(*k)*alpha*p->green;
 345           if ((channel & BlueChannel) != 0)
 346             pixel.blue+=(*k)*alpha*p->blue;
 347           if ((channel & OpacityChannel) != 0)
 348             pixel.opacity+=(*k)*p->opacity;
 349           if (((channel & IndexChannel) != 0) &&
 350               (image->colorspace == CMYKColorspace))
 351             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
 352           gamma+=(*k)*alpha;
 353           k++;
 354           p++;
 355         }
 356       }
 357       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 358       if ((channel & RedChannel) != 0)
 359         q->red=RoundToQuantum(gamma*pixel.red);
 360       if ((channel & GreenChannel) != 0)
 361         q->green=RoundToQuantum(gamma*pixel.green);
 362       if ((channel & BlueChannel) != 0)
 363         q->blue=RoundToQuantum(gamma*pixel.blue);
 364       if ((channel & OpacityChannel) != 0)
 365         q->opacity=RoundToQuantum(pixel.opacity);
 366       if (((channel & IndexChannel) != 0) &&
 367           (image->colorspace == CMYKColorspace))
 368         blur_indexes[x]=RoundToQuantum(gamma*pixel.index);
 369       q++;
 370       r++;
 371     }
 372     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
 373       status=MagickFalse;
 374     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 375       {
 376         MagickBooleanType
 377           proceed;
 378
 379 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 380   #pragma omp critical (MagickCore_AdaptiveBlurImageChannel)
 381 #endif
 382         proceed=SetImageProgress(image,AdaptiveBlurImageTag,progress++,
 383           image->rows);
 384         if (proceed == MagickFalse)
 385           status=MagickFalse;
 386       }
 387   }
 388   blur_image->type=image->type;
 389   blur_view=DestroyCacheView(blur_view);
 390   edge_view=DestroyCacheView(edge_view);
 391   image_view=DestroyCacheView(image_view);
 392   edge_image=DestroyImage(edge_image);
 393   for (i=0; i < (long) width;  i+=2)
 394     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 395   kernel=(double **) RelinquishMagickMemory(kernel);
 396   if (status == MagickFalse)
 397     blur_image=DestroyImage(blur_image);
 398   return(blur_image);
 399 }
 400 \f
 401 /*
 402 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 403 %                                                                             %
 404 %                                                                             %
 405 %                                                                             %
 406 %     A d a p t i v e S h a r p e n I m a g e                                 %
 407 %                                                                             %
 408 %                                                                             %
 409 %                                                                             %
 410 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 411 %
 412 %  AdaptiveSharpenImage() adaptively sharpens the image by sharpening more
 413 %  intensely near image edges and less intensely far from edges. We sharpen the
 414 %  image with a Gaussian operator of the given radius and standard deviation
 415 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
 416 %  radius of 0 and AdaptiveSharpenImage() selects a suitable radius for you.
 417 %
 418 %  The format of the AdaptiveSharpenImage method is:
 419 %
 420 %      Image *AdaptiveSharpenImage(const Image *image,const double radius,
 421 %        const double sigma,ExceptionInfo *exception)
 422 %      Image *AdaptiveSharpenImageChannel(const Image *image,
 423 %        const ChannelType channel,double radius,const double sigma,
 424 %        ExceptionInfo *exception)
 425 %
 426 %  A description of each parameter follows:
 427 %
 428 %    o image: the image.
 429 %
 430 %    o channel: the channel type.
 431 %
 432 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 433 %      pixel.
 434 %
 435 %    o sigma: the standard deviation of the Laplacian, in pixels.
 436 %
 437 %    o exception: return any errors or warnings in this structure.
 438 %
 439 */
 440
 441 MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
 442   const double sigma,ExceptionInfo *exception)
 443 {
 444   Image
 445     *sharp_image;
 446
 447   sharp_image=AdaptiveSharpenImageChannel(image,DefaultChannels,radius,sigma,
 448     exception);
 449   return(sharp_image);
 450 }
 451
 452 MagickExport Image *AdaptiveSharpenImageChannel(const Image *image,
 453   const ChannelType channel,const double radius,const double sigma,
 454   ExceptionInfo *exception)
 455 {
 456 #define AdaptiveSharpenImageTag  "Convolve/Image"
 457 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
 458
 459   double
 460     **kernel,
 461     normalize;
 462
 463   Image
 464     *sharp_image,
 465     *edge_image,
 466     *gaussian_image;
 467
 468   long
 469     j,
 470     k,
 471     progress,
 472     u,
 473     v,
 474     y;
 475
 476   MagickBooleanType
 477     status;
 478
 479   MagickPixelPacket
 480     bias;
 481
 482   register long
 483     i;
 484
 485   unsigned long
 486     width;
 487
 488   CacheView
 489     *sharp_view,
 490     *edge_view,
 491     *image_view;
 492
 493   assert(image != (const Image *) NULL);
 494   assert(image->signature == MagickSignature);
 495   if (image->debug != MagickFalse)
 496     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 497   assert(exception != (ExceptionInfo *) NULL);
 498   assert(exception->signature == MagickSignature);
 499   sharp_image=CloneImage(image,0,0,MagickTrue,exception);
 500   if (sharp_image == (Image *) NULL)
 501     return((Image *) NULL);
 502   if (fabs(sigma) <= MagickEpsilon)
 503     return(sharp_image);
 504   if (SetImageStorageClass(sharp_image,DirectClass) == MagickFalse)
 505     {
 506       InheritException(exception,&sharp_image->exception);
 507       sharp_image=DestroyImage(sharp_image);
 508       return((Image *) NULL);
 509     }
 510   /*
 511     Edge detect the image brighness channel, level, sharp, and level again.
 512   */
 513   edge_image=EdgeImage(image,radius,exception);
 514   if (edge_image == (Image *) NULL)
 515     {
 516       sharp_image=DestroyImage(sharp_image);
 517       return((Image *) NULL);
 518     }
 519   (void) LevelImage(edge_image,"20%,95%");
 520   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
 521   if (gaussian_image != (Image *) NULL)
 522     {
 523       edge_image=DestroyImage(edge_image);
 524       edge_image=gaussian_image;
 525     }
 526   (void) LevelImage(edge_image,"10%,95%");
 527   /*
 528     Create a set of kernels from maximum (radius,sigma) to minimum.
 529   */
 530   width=GetOptimalKernelWidth2D(radius,sigma);
 531   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 532   if (kernel == (double **) NULL)
 533     {
 534       edge_image=DestroyImage(edge_image);
 535       sharp_image=DestroyImage(sharp_image);
 536       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 537     }
 538   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
 539   for (i=0; i < (long) width; i+=2)
 540   {
 541     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
 542       sizeof(**kernel));
 543     if (kernel[i] == (double *) NULL)
 544       break;
 545     normalize=0.0;
 546     j=(long) (width-i)/2;
 547     k=0;
 548     for (v=(-j); v <= j; v++)
 549     {
 550       for (u=(-j); u <= j; u++)
 551       {
 552         kernel[i][k]=(-exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
 553           (2.0*MagickPI*MagickSigma*MagickSigma));
 554         normalize+=kernel[i][k];
 555         k++;
 556       }
 557     }
 558     if (fabs(normalize) <= MagickEpsilon)
 559       normalize=1.0;
 560     normalize=1.0/normalize;
 561     for (k=0; k < (j*j); k++)
 562       kernel[i][k]=normalize*kernel[i][k];
 563   }
 564   if (i < (long) width)
 565     {
 566       for (i-=2; i >= 0; i-=2)
 567         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 568       kernel=(double **) RelinquishMagickMemory(kernel);
 569       edge_image=DestroyImage(edge_image);
 570       sharp_image=DestroyImage(sharp_image);
 571       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 572     }
 573   /*
 574     Adaptively sharpen image.
 575   */
 576   status=MagickTrue;
 577   progress=0;
 578   GetMagickPixelPacket(image,&bias);
 579   SetMagickPixelPacketBias(image,&bias);
 580   image_view=AcquireCacheView(image);
 581   edge_view=AcquireCacheView(edge_image);
 582   sharp_view=AcquireCacheView(sharp_image);
 583 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 584   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 585 #endif
 586   for (y=0; y < (long) sharp_image->rows; y++)
 587   {
 588     register const IndexPacket
 589       *restrict indexes;
 590
 591     register const PixelPacket
 592       *restrict p,
 593       *restrict r;
 594
 595     register IndexPacket
 596       *restrict sharp_indexes;
 597
 598     register long
 599       x;
 600
 601     register PixelPacket
 602       *restrict q;
 603
 604     if (status == MagickFalse)
 605       continue;
 606     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
 607     q=QueueCacheViewAuthenticPixels(sharp_view,0,y,sharp_image->columns,1,
 608       exception);
 609     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 610       {
 611         status=MagickFalse;
 612         continue;
 613       }
 614     sharp_indexes=GetCacheViewAuthenticIndexQueue(sharp_view);
 615     for (x=0; x < (long) sharp_image->columns; x++)
 616     {
 617       MagickPixelPacket
 618         pixel;
 619
 620       MagickRealType
 621         alpha,
 622         gamma;
 623
 624       register const double
 625         *restrict k;
 626
 627       register long
 628         i,
 629         u,
 630         v;
 631
 632       gamma=0.0;
 633       i=(long) (width*(QuantumRange-QuantumScale*PixelIntensity(r))+0.5);
 634       if (i < 0)
 635         i=0;
 636       else
 637         if (i > (long) width)
 638           i=(long) width;
 639       if ((i & 0x01) != 0)
 640         i--;
 641       p=GetCacheViewVirtualPixels(image_view,x-((long) (width-i)/2L),y-(long)
 642         ((width-i)/2L),width-i,width-i,exception);
 643       if (p == (const PixelPacket *) NULL)
 644         break;
 645       indexes=GetCacheViewVirtualIndexQueue(image_view);
 646       k=kernel[i];
 647       pixel=bias;
 648       for (v=0; v < (long) (width-i); v++)
 649       {
 650         for (u=0; u < (long) (width-i); u++)
 651         {
 652           alpha=1.0;
 653           if (((channel & OpacityChannel) != 0) &&
 654               (image->matte != MagickFalse))
 655             alpha=(MagickRealType) (QuantumScale*(QuantumRange-p->opacity));
 656           if ((channel & RedChannel) != 0)
 657             pixel.red+=(*k)*alpha*p->red;
 658           if ((channel & GreenChannel) != 0)
 659             pixel.green+=(*k)*alpha*p->green;
 660           if ((channel & BlueChannel) != 0)
 661             pixel.blue+=(*k)*alpha*p->blue;
 662           if ((channel & OpacityChannel) != 0)
 663             pixel.opacity+=(*k)*p->opacity;
 664           if (((channel & IndexChannel) != 0) &&
 665               (image->colorspace == CMYKColorspace))
 666             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
 667           gamma+=(*k)*alpha;
 668           k++;
 669           p++;
 670         }
 671       }
 672       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 673       if ((channel & RedChannel) != 0)
 674         q->red=RoundToQuantum(gamma*pixel.red);
 675       if ((channel & GreenChannel) != 0)
 676         q->green=RoundToQuantum(gamma*pixel.green);
 677       if ((channel & BlueChannel) != 0)
 678         q->blue=RoundToQuantum(gamma*pixel.blue);
 679       if ((channel & OpacityChannel) != 0)
 680         q->opacity=RoundToQuantum(pixel.opacity);
 681       if (((channel & IndexChannel) != 0) &&
 682           (image->colorspace == CMYKColorspace))
 683         sharp_indexes[x]=RoundToQuantum(gamma*pixel.index);
 684       q++;
 685       r++;
 686     }
 687     if (SyncCacheViewAuthenticPixels(sharp_view,exception) == MagickFalse)
 688       status=MagickFalse;
 689     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 690       {
 691         MagickBooleanType
 692           proceed;
 693
 694 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 695   #pragma omp critical (MagickCore_AdaptiveSharpenImageChannel)
 696 #endif
 697         proceed=SetImageProgress(image,AdaptiveSharpenImageTag,progress++,
 698           image->rows);
 699         if (proceed == MagickFalse)
 700           status=MagickFalse;
 701       }
 702   }
 703   sharp_image->type=image->type;
 704   sharp_view=DestroyCacheView(sharp_view);
 705   edge_view=DestroyCacheView(edge_view);
 706   image_view=DestroyCacheView(image_view);
 707   edge_image=DestroyImage(edge_image);
 708   for (i=0; i < (long) width;  i+=2)
 709     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 710   kernel=(double **) RelinquishMagickMemory(kernel);
 711   if (status == MagickFalse)
 712     sharp_image=DestroyImage(sharp_image);
 713   return(sharp_image);
 714 }
 715 \f
 716 /*
 717 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 718 %                                                                             %
 719 %                                                                             %
 720 %                                                                             %
 721 %     B l u r I m a g e                                                       %
 722 %                                                                             %
 723 %                                                                             %
 724 %                                                                             %
 725 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 726 %
 727 %  BlurImage() blurs an image.  We convolve the image with a Gaussian operator
 728 %  of the given radius and standard deviation (sigma).  For reasonable results,
 729 %  the radius should be larger than sigma.  Use a radius of 0 and BlurImage()
 730 %  selects a suitable radius for you.
 731 %
 732 %  BlurImage() differs from GaussianBlurImage() in that it uses a separable
 733 %  kernel which is faster but mathematically equivalent to the non-separable
 734 %  kernel.
 735 %
 736 %  The format of the BlurImage method is:
 737 %
 738 %      Image *BlurImage(const Image *image,const double radius,
 739 %        const double sigma,ExceptionInfo *exception)
 740 %      Image *BlurImageChannel(const Image *image,const ChannelType channel,
 741 %        const double radius,const double sigma,ExceptionInfo *exception)
 742 %
 743 %  A description of each parameter follows:
 744 %
 745 %    o image: the image.
 746 %
 747 %    o channel: the channel type.
 748 %
 749 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 750 %      pixel.
 751 %
 752 %    o sigma: the standard deviation of the Gaussian, in pixels.
 753 %
 754 %    o exception: return any errors or warnings in this structure.
 755 %
 756 */
 757
 758 MagickExport Image *BlurImage(const Image *image,const double radius,
 759   const double sigma,ExceptionInfo *exception)
 760 {
 761   Image
 762     *blur_image;
 763
 764   blur_image=BlurImageChannel(image,DefaultChannels,radius,sigma,exception);
 765   return(blur_image);
 766 }
 767
 768 static double *GetBlurKernel(const unsigned long width,const double sigma)
 769 {
 770   double
 771     *kernel,
 772     normalize;
 773
 774   long
 775     j,
 776     k;
 777
 778   register long
 779     i;
 780
 781   /*
 782     Generate a 1-D convolution kernel.
 783   */
 784   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
 785   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 786   if (kernel == (double *) NULL)
 787     return(0);
 788   normalize=0.0;
 789   j=(long) width/2;
 790   i=0;
 791   for (k=(-j); k <= j; k++)
 792   {
 793     kernel[i]=exp(-((double) k*k)/(2.0*MagickSigma*MagickSigma))/
 794       (MagickSQ2PI*MagickSigma);
 795     normalize+=kernel[i];
 796     i++;
 797   }
 798   for (i=0; i < (long) width; i++)
 799     kernel[i]/=normalize;
 800   return(kernel);
 801 }
 802
 803 MagickExport Image *BlurImageChannel(const Image *image,
 804   const ChannelType channel,const double radius,const double sigma,
 805   ExceptionInfo *exception)
 806 {
 807 #define BlurImageTag  "Blur/Image"
 808
 809   double
 810     *kernel;
 811
 812   Image
 813     *blur_image;
 814
 815   long
 816     progress,
 817     x,
 818     y;
 819
 820   MagickBooleanType
 821     status;
 822
 823   MagickPixelPacket
 824     bias;
 825
 826   register long
 827     i;
 828
 829   unsigned long
 830     width;
 831
 832   CacheView
 833     *blur_view,
 834     *image_view;
 835
 836   /*
 837     Initialize blur image attributes.
 838   */
 839   assert(image != (Image *) NULL);
 840   assert(image->signature == MagickSignature);
 841   if (image->debug != MagickFalse)
 842     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 843   assert(exception != (ExceptionInfo *) NULL);
 844   assert(exception->signature == MagickSignature);
 845   blur_image=CloneImage(image,0,0,MagickTrue,exception);
 846   if (blur_image == (Image *) NULL)
 847     return((Image *) NULL);
 848   if (fabs(sigma) <= MagickEpsilon)
 849     return(blur_image);
 850   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
 851     {
 852       InheritException(exception,&blur_image->exception);
 853       blur_image=DestroyImage(blur_image);
 854       return((Image *) NULL);
 855     }
 856   width=GetOptimalKernelWidth1D(radius,sigma);
 857   kernel=GetBlurKernel(width,sigma);
 858   if (kernel == (double *) NULL)
 859     {
 860       blur_image=DestroyImage(blur_image);
 861       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 862     }
 863   if (image->debug != MagickFalse)
 864     {
 865       char
 866         format[MaxTextExtent],
 867         *message;
 868
 869       register const double
 870         *k;
 871
 872       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
 873         "  BlurImage with %ld kernel:",width);
 874       message=AcquireString("");
 875       k=kernel;
 876       for (i=0; i < (long) width; i++)
 877       {
 878         *message='\0';
 879         (void) FormatMagickString(format,MaxTextExtent,"%ld: ",i);
 880         (void) ConcatenateString(&message,format);
 881         (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
 882         (void) ConcatenateString(&message,format);
 883         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
 884       }
 885       message=DestroyString(message);
 886     }
 887   /*
 888     Blur rows.
 889   */
 890   status=MagickTrue;
 891   progress=0;
 892   GetMagickPixelPacket(image,&bias);
 893   SetMagickPixelPacketBias(image,&bias);
 894   image_view=AcquireCacheView(image);
 895   blur_view=AcquireCacheView(blur_image);
 896 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 897   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 898 #endif
 899   for (y=0; y < (long) blur_image->rows; y++)
 900   {
 901     register const IndexPacket
 902       *restrict indexes;
 903
 904     register const PixelPacket
 905       *restrict p;
 906
 907     register IndexPacket
 908       *restrict blur_indexes;
 909
 910     register long
 911       x;
 912
 913     register PixelPacket
 914       *restrict q;
 915
 916     if (status == MagickFalse)
 917       continue;
 918     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y,image->columns+
 919       width,1,exception);
 920     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
 921       exception);
 922     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 923       {
 924         status=MagickFalse;
 925         continue;
 926       }
 927     indexes=GetCacheViewVirtualIndexQueue(image_view);
 928     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
 929     for (x=0; x < (long) blur_image->columns; x++)
 930     {
 931       MagickPixelPacket
 932         pixel;
 933
 934       register const double
 935         *restrict k;
 936
 937       register const PixelPacket
 938         *restrict kernel_pixels;
 939
 940       register long
 941         i;
 942
 943       pixel=bias;
 944       k=kernel;
 945       kernel_pixels=p;
 946       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
 947         {
 948           for (i=0; i < (long) width; i++)
 949           {
 950             pixel.red+=(*k)*kernel_pixels->red;
 951             pixel.green+=(*k)*kernel_pixels->green;
 952             pixel.blue+=(*k)*kernel_pixels->blue;
 953             k++;
 954             kernel_pixels++;
 955           }
 956           if ((channel & RedChannel) != 0)
 957             q->red=RoundToQuantum(pixel.red);
 958           if ((channel & GreenChannel) != 0)
 959             q->green=RoundToQuantum(pixel.green);
 960           if ((channel & BlueChannel) != 0)
 961             q->blue=RoundToQuantum(pixel.blue);
 962           if ((channel & OpacityChannel) != 0)
 963             {
 964               k=kernel;
 965               kernel_pixels=p;
 966               for (i=0; i < (long) width; i++)
 967               {
 968                 pixel.opacity+=(*k)*kernel_pixels->opacity;
 969                 k++;
 970                 kernel_pixels++;
 971               }
 972               q->opacity=RoundToQuantum(pixel.opacity);
 973             }
 974           if (((channel & IndexChannel) != 0) &&
 975               (image->colorspace == CMYKColorspace))
 976             {
 977               register const IndexPacket
 978                 *restrict kernel_indexes;
 979
 980               k=kernel;
 981               kernel_indexes=indexes;
 982               for (i=0; i < (long) width; i++)
 983               {
 984                 pixel.index+=(*k)*(*kernel_indexes);
 985                 k++;
 986                 kernel_indexes++;
 987               }
 988               blur_indexes[x]=RoundToQuantum(pixel.index);
 989             }
 990         }
 991       else
 992         {
 993           MagickRealType
 994             alpha,
 995             gamma;
 996
 997           gamma=0.0;
 998           for (i=0; i < (long) width; i++)
 999           {
1000             alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1001               kernel_pixels->opacity));
1002             pixel.red+=(*k)*alpha*kernel_pixels->red;
1003             pixel.green+=(*k)*alpha*kernel_pixels->green;
1004             pixel.blue+=(*k)*alpha*kernel_pixels->blue;
1005             gamma+=(*k)*alpha;
1006             k++;
1007             kernel_pixels++;
1008           }
1009           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1010           if ((channel & RedChannel) != 0)
1011             q->red=RoundToQuantum(gamma*pixel.red);
1012           if ((channel & GreenChannel) != 0)
1013             q->green=RoundToQuantum(gamma*pixel.green);
1014           if ((channel & BlueChannel) != 0)
1015             q->blue=RoundToQuantum(gamma*pixel.blue);
1016           if ((channel & OpacityChannel) != 0)
1017             {
1018               k=kernel;
1019               kernel_pixels=p;
1020               for (i=0; i < (long) width; i++)
1021               {
1022                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1023                 k++;
1024                 kernel_pixels++;
1025               }
1026               q->opacity=RoundToQuantum(pixel.opacity);
1027             }
1028           if (((channel & IndexChannel) != 0) &&
1029               (image->colorspace == CMYKColorspace))
1030             {
1031               register const IndexPacket
1032                 *restrict kernel_indexes;
1033
1034               k=kernel;
1035               kernel_pixels=p;
1036               kernel_indexes=indexes;
1037               for (i=0; i < (long) width; i++)
1038               {
1039                 alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1040                   kernel_pixels->opacity));
1041                 pixel.index+=(*k)*alpha*(*kernel_indexes);
1042                 k++;
1043                 kernel_pixels++;
1044                 kernel_indexes++;
1045               }
1046               blur_indexes[x]=RoundToQuantum(gamma*pixel.index);
1047             }
1048         }
1049       p++;
1050       q++;
1051     }
1052     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
1053       status=MagickFalse;
1054     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1055       {
1056         MagickBooleanType
1057           proceed;
1058
1059 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1060   #pragma omp critical (MagickCore_BlurImageChannel)
1061 #endif
1062         proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
1063           blur_image->columns);
1064         if (proceed == MagickFalse)
1065           status=MagickFalse;
1066       }
1067   }
1068   blur_view=DestroyCacheView(blur_view);
1069   image_view=DestroyCacheView(image_view);
1070   /*
1071     Blur columns.
1072   */
1073   image_view=AcquireCacheView(blur_image);
1074   blur_view=AcquireCacheView(blur_image);
1075 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1076   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1077 #endif
1078   for (x=0; x < (long) blur_image->columns; x++)
1079   {
1080     register const IndexPacket
1081       *restrict indexes;
1082
1083     register const PixelPacket
1084       *restrict p;
1085
1086     register IndexPacket
1087       *restrict blur_indexes;
1088
1089     register long
1090       y;
1091
1092     register PixelPacket
1093       *restrict q;
1094
1095     if (status == MagickFalse)
1096       continue;
1097     p=GetCacheViewVirtualPixels(image_view,x,-((long) width/2L),1,image->rows+
1098       width,exception);
1099     q=GetCacheViewAuthenticPixels(blur_view,x,0,1,blur_image->rows,exception);
1100     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1101       {
1102         status=MagickFalse;
1103         continue;
1104       }
1105     indexes=GetCacheViewVirtualIndexQueue(image_view);
1106     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
1107     for (y=0; y < (long) blur_image->rows; y++)
1108     {
1109       MagickPixelPacket
1110         pixel;
1111
1112       register const double
1113         *restrict k;
1114
1115       register const PixelPacket
1116         *restrict kernel_pixels;
1117
1118       register long
1119         i;
1120
1121       pixel=bias;
1122       k=kernel;
1123       kernel_pixels=p;
1124       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
1125         {
1126           for (i=0; i < (long) width; i++)
1127           {
1128             pixel.red+=(*k)*kernel_pixels->red;
1129             pixel.green+=(*k)*kernel_pixels->green;
1130             pixel.blue+=(*k)*kernel_pixels->blue;
1131             k++;
1132             kernel_pixels++;
1133           }
1134           if ((channel & RedChannel) != 0)
1135             q->red=RoundToQuantum(pixel.red);
1136           if ((channel & GreenChannel) != 0)
1137             q->green=RoundToQuantum(pixel.green);
1138           if ((channel & BlueChannel) != 0)
1139             q->blue=RoundToQuantum(pixel.blue);
1140           if ((channel & OpacityChannel) != 0)
1141             {
1142               k=kernel;
1143               kernel_pixels=p;
1144               for (i=0; i < (long) width; i++)
1145               {
1146                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1147                 k++;
1148                 kernel_pixels++;
1149               }
1150               q->opacity=RoundToQuantum(pixel.opacity);
1151             }
1152           if (((channel & IndexChannel) != 0) &&
1153               (image->colorspace == CMYKColorspace))
1154             {
1155               register const IndexPacket
1156                 *restrict kernel_indexes;
1157
1158               k=kernel;
1159               kernel_indexes=indexes;
1160               for (i=0; i < (long) width; i++)
1161               {
1162                 pixel.index+=(*k)*(*kernel_indexes);
1163                 k++;
1164                 kernel_indexes++;
1165               }
1166               blur_indexes[y]=RoundToQuantum(pixel.index);
1167             }
1168         }
1169       else
1170         {
1171           MagickRealType
1172             alpha,
1173             gamma;
1174
1175           gamma=0.0;
1176           for (i=0; i < (long) width; i++)
1177           {
1178             alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1179               kernel_pixels->opacity));
1180             pixel.red+=(*k)*alpha*kernel_pixels->red;
1181             pixel.green+=(*k)*alpha*kernel_pixels->green;
1182             pixel.blue+=(*k)*alpha*kernel_pixels->blue;
1183             gamma+=(*k)*alpha;
1184             k++;
1185             kernel_pixels++;
1186           }
1187           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1188           if ((channel & RedChannel) != 0)
1189             q->red=RoundToQuantum(gamma*pixel.red);
1190           if ((channel & GreenChannel) != 0)
1191             q->green=RoundToQuantum(gamma*pixel.green);
1192           if ((channel & BlueChannel) != 0)
1193             q->blue=RoundToQuantum(gamma*pixel.blue);
1194           if ((channel & OpacityChannel) != 0)
1195             {
1196               k=kernel;
1197               kernel_pixels=p;
1198               for (i=0; i < (long) width; i++)
1199               {
1200                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1201                 k++;
1202                 kernel_pixels++;
1203               }
1204               q->opacity=RoundToQuantum(pixel.opacity);
1205             }
1206           if (((channel & IndexChannel) != 0) &&
1207               (image->colorspace == CMYKColorspace))
1208             {
1209               register const IndexPacket
1210                 *restrict kernel_indexes;
1211
1212               k=kernel;
1213               kernel_pixels=p;
1214               kernel_indexes=indexes;
1215               for (i=0; i < (long) width; i++)
1216               {
1217                 alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1218                   kernel_pixels->opacity));
1219                 pixel.index+=(*k)*alpha*(*kernel_indexes);
1220                 k++;
1221                 kernel_pixels++;
1222                 kernel_indexes++;
1223               }
1224               blur_indexes[y]=RoundToQuantum(gamma*pixel.index);
1225             }
1226         }
1227       p++;
1228       q++;
1229     }
1230     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
1231       status=MagickFalse;
1232     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1233       {
1234         MagickBooleanType
1235           proceed;
1236
1237 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1238   #pragma omp critical (MagickCore_BlurImageChannel)
1239 #endif
1240         proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
1241           blur_image->columns);
1242         if (proceed == MagickFalse)
1243           status=MagickFalse;
1244       }
1245   }
1246   blur_view=DestroyCacheView(blur_view);
1247   image_view=DestroyCacheView(image_view);
1248   kernel=(double *) RelinquishMagickMemory(kernel);
1249   if (status == MagickFalse)
1250     blur_image=DestroyImage(blur_image);
1251   blur_image->type=image->type;
1252   return(blur_image);
1253 }
1254 \f
1255 /*
1256 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1257 %                                                                             %
1258 %                                                                             %
1259 %                                                                             %
1260 %     C o n v o l v e I m a g e                                               %
1261 %                                                                             %
1262 %                                                                             %
1263 %                                                                             %
1264 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1265 %
1266 %  ConvolveImage() applies a custom convolution kernel to the image.
1267 %
1268 %  The format of the ConvolveImage method is:
1269 %
1270 %      Image *ConvolveImage(const Image *image,const unsigned long order,
1271 %        const double *kernel,ExceptionInfo *exception)
1272 %      Image *ConvolveImageChannel(const Image *image,const ChannelType channel,
1273 %        const unsigned long order,const double *kernel,
1274 %        ExceptionInfo *exception)
1275 %
1276 %  A description of each parameter follows:
1277 %
1278 %    o image: the image.
1279 %
1280 %    o channel: the channel type.
1281 %
1282 %    o order: the number of columns and rows in the filter kernel.
1283 %
1284 %    o kernel: An array of double representing the convolution kernel.
1285 %
1286 %    o exception: return any errors or warnings in this structure.
1287 %
1288 */
1289
1290 MagickExport Image *ConvolveImage(const Image *image,const unsigned long order,
1291   const double *kernel,ExceptionInfo *exception)
1292 {
1293   Image
1294     *convolve_image;
1295
1296   convolve_image=ConvolveImageChannel(image,DefaultChannels,order,kernel,
1297     exception);
1298   return(convolve_image);
1299 }
1300
1301 MagickExport Image *ConvolveImageChannel(const Image *image,
1302   const ChannelType channel,const unsigned long order,const double *kernel,
1303   ExceptionInfo *exception)
1304 {
1305 #define ConvolveImageTag  "Convolve/Image"
1306
1307   double
1308     *normal_kernel;
1309
1310   Image
1311     *convolve_image;
1312
1313   long
1314     progress,
1315     y;
1316
1317   MagickBooleanType
1318     status;
1319
1320   MagickPixelPacket
1321     bias;
1322
1323   MagickRealType
1324     gamma;
1325
1326   register long
1327     i;
1328
1329   unsigned long
1330     width;
1331
1332   CacheView
1333     *convolve_view,
1334     *image_view;
1335
1336   /*
1337     Initialize convolve image attributes.
1338   */
1339   assert(image != (Image *) NULL);
1340   assert(image->signature == MagickSignature);
1341   if (image->debug != MagickFalse)
1342     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1343   assert(exception != (ExceptionInfo *) NULL);
1344   assert(exception->signature == MagickSignature);
1345   width=order;
1346   if ((width % 2) == 0)
1347     ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
1348   convolve_image=CloneImage(image,0,0,MagickTrue,exception);
1349   if (convolve_image == (Image *) NULL)
1350     return((Image *) NULL);
1351   if (SetImageStorageClass(convolve_image,DirectClass) == MagickFalse)
1352     {
1353       InheritException(exception,&convolve_image->exception);
1354       convolve_image=DestroyImage(convolve_image);
1355       return((Image *) NULL);
1356     }
1357   if (image->debug != MagickFalse)
1358     {
1359       char
1360         format[MaxTextExtent],
1361         *message;
1362
1363       long
1364         u,
1365         v;
1366
1367       register const double
1368         *k;
1369
1370       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1371         "  ConvolveImage with %ldx%ld kernel:",width,width);
1372       message=AcquireString("");
1373       k=kernel;
1374       for (v=0; v < (long) width; v++)
1375       {
1376         *message='\0';
1377         (void) FormatMagickString(format,MaxTextExtent,"%ld: ",v);
1378         (void) ConcatenateString(&message,format);
1379         for (u=0; u < (long) width; u++)
1380         {
1381           (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
1382           (void) ConcatenateString(&message,format);
1383         }
1384         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
1385       }
1386       message=DestroyString(message);
1387     }
1388   /*
1389     Normalize kernel.
1390   */
1391   normal_kernel=(double *) AcquireQuantumMemory(width*width,
1392     sizeof(*normal_kernel));
1393   if (normal_kernel == (double *) NULL)
1394     {
1395       convolve_image=DestroyImage(convolve_image);
1396       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1397     }
1398   gamma=0.0;
1399   for (i=0; i < (long) (width*width); i++)
1400     gamma+=kernel[i];
1401   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1402   for (i=0; i < (long) (width*width); i++)
1403     normal_kernel[i]=gamma*kernel[i];
1404   /*
1405     Convolve image.
1406   */
1407   status=MagickTrue;
1408   progress=0;
1409   GetMagickPixelPacket(image,&bias);
1410   SetMagickPixelPacketBias(image,&bias);
1411   image_view=AcquireCacheView(image);
1412   convolve_view=AcquireCacheView(convolve_image);
1413 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1414   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1415 #endif
1416   for (y=0; y < (long) image->rows; y++)
1417   {
1418     MagickBooleanType
1419       sync;
1420
1421     register const IndexPacket
1422       *restrict indexes;
1423
1424     register const PixelPacket
1425       *restrict p;
1426
1427     register IndexPacket
1428       *restrict convolve_indexes;
1429
1430     register long
1431       x;
1432
1433     register PixelPacket
1434       *restrict q;
1435
1436     if (status == MagickFalse)
1437       continue;
1438     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
1439       2L),image->columns+width,width,exception);
1440     q=GetCacheViewAuthenticPixels(convolve_view,0,y,convolve_image->columns,1,
1441       exception);
1442     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1443       {
1444         status=MagickFalse;
1445         continue;
1446       }
1447     indexes=GetCacheViewVirtualIndexQueue(image_view);
1448     convolve_indexes=GetCacheViewAuthenticIndexQueue(convolve_view);
1449     for (x=0; x < (long) image->columns; x++)
1450     {
1451       long
1452         v;
1453
1454       MagickPixelPacket
1455         pixel;
1456
1457       register const double
1458         *restrict k;
1459
1460       register const PixelPacket
1461         *restrict kernel_pixels;
1462
1463       register long
1464         u;
1465
1466       pixel=bias;
1467       k=normal_kernel;
1468       kernel_pixels=p;
1469       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
1470         {
1471           for (v=0; v < (long) width; v++)
1472           {
1473             for (u=0; u < (long) width; u++)
1474             {
1475               pixel.red+=(*k)*kernel_pixels[u].red;
1476               pixel.green+=(*k)*kernel_pixels[u].green;
1477               pixel.blue+=(*k)*kernel_pixels[u].blue;
1478               k++;
1479             }
1480             kernel_pixels+=image->columns+width;
1481           }
1482           if ((channel & RedChannel) != 0)
1483             q->red=RoundToQuantum(pixel.red);
1484           if ((channel & GreenChannel) != 0)
1485             q->green=RoundToQuantum(pixel.green);
1486           if ((channel & BlueChannel) != 0)
1487             q->blue=RoundToQuantum(pixel.blue);
1488           if ((channel & OpacityChannel) != 0)
1489             {
1490               k=normal_kernel;
1491               kernel_pixels=p;
1492               for (v=0; v < (long) width; v++)
1493               {
1494                 for (u=0; u < (long) width; u++)
1495                 {
1496                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
1497                   k++;
1498                 }
1499                 kernel_pixels+=image->columns+width;
1500               }
1501               q->opacity=RoundToQuantum(pixel.opacity);
1502             }
1503           if (((channel & IndexChannel) != 0) &&
1504               (image->colorspace == CMYKColorspace))
1505             {
1506               register const IndexPacket
1507                 *restrict kernel_indexes;
1508
1509               k=normal_kernel;
1510               kernel_indexes=indexes;
1511               for (v=0; v < (long) width; v++)
1512               {
1513                 for (u=0; u < (long) width; u++)
1514                 {
1515                   pixel.index+=(*k)*kernel_indexes[u];
1516                   k++;
1517                 }
1518                 kernel_indexes+=image->columns+width;
1519               }
1520               convolve_indexes[x]=RoundToQuantum(pixel.index);
1521             }
1522         }
1523       else
1524         {
1525           MagickRealType
1526             alpha,
1527             gamma;
1528
1529           gamma=0.0;
1530           for (v=0; v < (long) width; v++)
1531           {
1532             for (u=0; u < (long) width; u++)
1533             {
1534               alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1535                 kernel_pixels[u].opacity));
1536               pixel.red+=(*k)*alpha*kernel_pixels[u].red;
1537               pixel.green+=(*k)*alpha*kernel_pixels[u].green;
1538               pixel.blue+=(*k)*alpha*kernel_pixels[u].blue;
1539               pixel.opacity+=(*k)*kernel_pixels[u].opacity;
1540               gamma+=(*k)*alpha;
1541               k++;
1542             }
1543             kernel_pixels+=image->columns+width;
1544           }
1545           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1546           if ((channel & RedChannel) != 0)
1547             q->red=RoundToQuantum(gamma*pixel.red);
1548           if ((channel & GreenChannel) != 0)
1549             q->green=RoundToQuantum(gamma*pixel.green);
1550           if ((channel & BlueChannel) != 0)
1551             q->blue=RoundToQuantum(gamma*pixel.blue);
1552           if ((channel & OpacityChannel) != 0)
1553             {
1554               k=normal_kernel;
1555               kernel_pixels=p;
1556               for (v=0; v < (long) width; v++)
1557               {
1558                 for (u=0; u < (long) width; u++)
1559                 {
1560                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
1561                   k++;
1562                 }
1563                 kernel_pixels+=image->columns+width;
1564               }
1565               q->opacity=RoundToQuantum(pixel.opacity);
1566             }
1567           if (((channel & IndexChannel) != 0) &&
1568               (image->colorspace == CMYKColorspace))
1569             {
1570               register const IndexPacket
1571                 *restrict kernel_indexes;
1572
1573               k=normal_kernel;
1574               kernel_pixels=p;
1575               kernel_indexes=indexes;
1576               for (v=0; v < (long) width; v++)
1577               {
1578                 for (u=0; u < (long) width; u++)
1579                 {
1580                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1581                     kernel_pixels[u].opacity));
1582                   pixel.index+=(*k)*alpha*kernel_indexes[u];
1583                   k++;
1584                 }
1585                 kernel_pixels+=image->columns+width;
1586                 kernel_indexes+=image->columns+width;
1587               }
1588               convolve_indexes[x]=RoundToQuantum(gamma*pixel.index);
1589             }
1590         }
1591       p++;
1592       q++;
1593     }
1594     sync=SyncCacheViewAuthenticPixels(convolve_view,exception);
1595     if (sync == MagickFalse)
1596       status=MagickFalse;
1597     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1598       {
1599         MagickBooleanType
1600           proceed;
1601
1602 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1603   #pragma omp critical (MagickCore_ConvolveImageChannel)
1604 #endif
1605         proceed=SetImageProgress(image,ConvolveImageTag,progress++,image->rows);
1606         if (proceed == MagickFalse)
1607           status=MagickFalse;
1608       }
1609   }
1610   convolve_image->type=image->type;
1611   convolve_view=DestroyCacheView(convolve_view);
1612   image_view=DestroyCacheView(image_view);
1613   normal_kernel=(double *) RelinquishMagickMemory(normal_kernel);
1614   if (status == MagickFalse)
1615     convolve_image=DestroyImage(convolve_image);
1616   return(convolve_image);
1617 }
1618 \f
1619 /*
1620 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1621 %                                                                             %
1622 %                                                                             %
1623 %                                                                             %
1624 %     D e s p e c k l e I m a g e                                             %
1625 %                                                                             %
1626 %                                                                             %
1627 %                                                                             %
1628 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1629 %
1630 %  DespeckleImage() reduces the speckle noise in an image while perserving the
1631 %  edges of the original image.
1632 %
1633 %  The format of the DespeckleImage method is:
1634 %
1635 %      Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
1636 %
1637 %  A description of each parameter follows:
1638 %
1639 %    o image: the image.
1640 %
1641 %    o exception: return any errors or warnings in this structure.
1642 %
1643 */
1644
1645 static Quantum **DestroyPixelThreadSet(Quantum **pixels)
1646 {
1647   register long
1648     i;
1649
1650   assert(pixels != (Quantum **) NULL);
1651   for (i=0; i < (long) GetOpenMPMaximumThreads(); i++)
1652     if (pixels[i] != (Quantum *) NULL)
1653       pixels[i]=(Quantum *) RelinquishMagickMemory(pixels[i]);
1654   pixels=(Quantum **) RelinquishAlignedMemory(pixels);
1655   return(pixels);
1656 }
1657
1658 static Quantum **AcquirePixelThreadSet(const size_t count)
1659 {
1660   register long
1661     i;
1662
1663   Quantum
1664     **pixels;
1665
1666   unsigned long
1667     number_threads;
1668
1669   number_threads=GetOpenMPMaximumThreads();
1670   pixels=(Quantum **) AcquireAlignedMemory(number_threads,sizeof(*pixels));
1671   if (pixels == (Quantum **) NULL)
1672     return((Quantum **) NULL);
1673   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
1674   for (i=0; i < (long) number_threads; i++)
1675   {
1676     pixels[i]=(Quantum *) AcquireQuantumMemory(count,sizeof(**pixels));
1677     if (pixels[i] == (Quantum *) NULL)
1678       return(DestroyPixelThreadSet(pixels));
1679   }
1680   return(pixels);
1681 }
1682
1683 static void Hull(const long x_offset,const long y_offset,
1684   const unsigned long columns,const unsigned long rows,Quantum *f,Quantum *g,
1685   const int polarity)
1686 {
1687   long
1688     y;
1689
1690   MagickRealType
1691     v;
1692
1693   register long
1694     x;
1695
1696   register Quantum
1697     *p,
1698     *q,
1699     *r,
1700     *s;
1701
1702   assert(f != (Quantum *) NULL);
1703   assert(g != (Quantum *) NULL);
1704   p=f+(columns+2);
1705   q=g+(columns+2);
1706   r=p+(y_offset*((long) columns+2)+x_offset);
1707   for (y=0; y < (long) rows; y++)
1708   {
1709     p++;
1710     q++;
1711     r++;
1712     if (polarity > 0)
1713       for (x=(long) columns; x != 0; x--)
1714       {
1715         v=(MagickRealType) (*p);
1716         if ((MagickRealType) *r >= (v+(MagickRealType) ScaleCharToQuantum(2)))
1717           v+=ScaleCharToQuantum(1);
1718         *q=(Quantum) v;
1719         p++;
1720         q++;
1721         r++;
1722       }
1723     else
1724       for (x=(long) columns; x != 0; x--)
1725       {
1726         v=(MagickRealType) (*p);
1727         if ((MagickRealType) *r <= (v-(MagickRealType) ScaleCharToQuantum(2)))
1728           v-=(long) ScaleCharToQuantum(1);
1729         *q=(Quantum) v;
1730         p++;
1731         q++;
1732         r++;
1733       }
1734     p++;
1735     q++;
1736     r++;
1737   }
1738   p=f+(columns+2);
1739   q=g+(columns+2);
1740   r=q+(y_offset*((long) columns+2)+x_offset);
1741   s=q-(y_offset*((long) columns+2)+x_offset);
1742   for (y=0; y < (long) rows; y++)
1743   {
1744     p++;
1745     q++;
1746     r++;
1747     s++;
1748     if (polarity > 0)
1749       for (x=(long) columns; x != 0; x--)
1750       {
1751         v=(MagickRealType) (*q);
1752         if (((MagickRealType) *s >=
1753              (v+(MagickRealType) ScaleCharToQuantum(2))) &&
1754             ((MagickRealType) *r > v))
1755           v+=ScaleCharToQuantum(1);
1756         *p=(Quantum) v;
1757         p++;
1758         q++;
1759         r++;
1760         s++;
1761       }
1762     else
1763       for (x=(long) columns; x != 0; x--)
1764       {
1765         v=(MagickRealType) (*q);
1766         if (((MagickRealType) *s <=
1767              (v-(MagickRealType) ScaleCharToQuantum(2))) &&
1768             ((MagickRealType) *r < v))
1769           v-=(MagickRealType) ScaleCharToQuantum(1);
1770         *p=(Quantum) v;
1771         p++;
1772         q++;
1773         r++;
1774         s++;
1775       }
1776     p++;
1777     q++;
1778     r++;
1779     s++;
1780   }
1781 }
1782
1783 MagickExport Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
1784 {
1785 #define DespeckleImageTag  "Despeckle/Image"
1786
1787   CacheView
1788     *despeckle_view,
1789     *image_view;
1790
1791   Image
1792     *despeckle_image;
1793
1794   long
1795     channel;
1796
1797   MagickBooleanType
1798     status;
1799
1800   Quantum
1801     **buffers,
1802     **pixels;
1803
1804   size_t
1805     length;
1806
1807   static const int
1808     X[4] = {0, 1, 1,-1},
1809     Y[4] = {1, 0, 1, 1};
1810
1811   /*
1812     Allocate despeckled image.
1813   */
1814   assert(image != (const Image *) NULL);
1815   assert(image->signature == MagickSignature);
1816   if (image->debug != MagickFalse)
1817     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1818   assert(exception != (ExceptionInfo *) NULL);
1819   assert(exception->signature == MagickSignature);
1820   despeckle_image=CloneImage(image,image->columns,image->rows,MagickTrue,
1821     exception);
1822   if (despeckle_image == (Image *) NULL)
1823     return((Image *) NULL);
1824   if (SetImageStorageClass(despeckle_image,DirectClass) == MagickFalse)
1825     {
1826       InheritException(exception,&despeckle_image->exception);
1827       despeckle_image=DestroyImage(despeckle_image);
1828       return((Image *) NULL);
1829     }
1830   /*
1831     Allocate image buffers.
1832   */
1833   length=(size_t) ((image->columns+2)*(image->rows+2));
1834   pixels=AcquirePixelThreadSet(length);
1835   buffers=AcquirePixelThreadSet(length);
1836   if ((pixels == (Quantum **) NULL) || (buffers == (Quantum **) NULL))
1837     {
1838       if (buffers != (Quantum **) NULL)
1839         buffers=DestroyPixelThreadSet(buffers);
1840       if (pixels != (Quantum **) NULL)
1841         pixels=DestroyPixelThreadSet(pixels);
1842       despeckle_image=DestroyImage(despeckle_image);
1843       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1844     }
1845   /*
1846     Reduce speckle in the image.
1847   */
1848   status=MagickTrue;
1849   image_view=AcquireCacheView(image);
1850   despeckle_view=AcquireCacheView(despeckle_image);
1851 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1852   #pragma omp parallel for schedule(dynamic,4) shared(status)
1853 #endif
1854   for (channel=0; channel <= 3; channel++)
1855   {
1856     long
1857       j,
1858       y;
1859
1860     register long
1861       i,
1862       id,
1863       x;
1864
1865     register Quantum
1866       *buffer,
1867       *pixel;
1868
1869     if (status == MagickFalse)
1870       continue;
1871     id=GetOpenMPThreadId();
1872     pixel=pixels[id];
1873     (void) ResetMagickMemory(pixel,0,length*sizeof(*pixel));
1874     buffer=buffers[id];
1875     j=(long) image->columns+2;
1876     for (y=0; y < (long) image->rows; y++)
1877     {
1878       register const PixelPacket
1879         *restrict p;
1880
1881       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1882       if (p == (const PixelPacket *) NULL)
1883         break;
1884       j++;
1885       for (x=0; x < (long) image->columns; x++)
1886       {
1887         switch (channel)
1888         {
1889           case 0: pixel[j]=p->red; break;
1890           case 1: pixel[j]=p->green; break;
1891           case 2: pixel[j]=p->blue; break;
1892           case 3: pixel[j]=p->opacity; break;
1893           default: break;
1894         }
1895         p++;
1896         j++;
1897       }
1898       j++;
1899     }
1900     (void) ResetMagickMemory(buffer,0,length*sizeof(*buffer));
1901     for (i=0; i < 4; i++)
1902     {
1903       Hull(X[i],Y[i],image->columns,image->rows,pixel,buffer,1);
1904       Hull(-X[i],-Y[i],image->columns,image->rows,pixel,buffer,1);
1905       Hull(-X[i],-Y[i],image->columns,image->rows,pixel,buffer,-1);
1906       Hull(X[i],Y[i],image->columns,image->rows,pixel,buffer,-1);
1907     }
1908     j=(long) image->columns+2;
1909     for (y=0; y < (long) image->rows; y++)
1910     {
1911       MagickBooleanType
1912         sync;
1913
1914       register PixelPacket
1915         *restrict q;
1916
1917       q=GetCacheViewAuthenticPixels(despeckle_view,0,y,despeckle_image->columns,
1918         1,exception);
1919       if (q == (PixelPacket *) NULL)
1920         break;
1921       j++;
1922       for (x=0; x < (long) image->columns; x++)
1923       {
1924         switch (channel)
1925         {
1926           case 0: q->red=pixel[j]; break;
1927           case 1: q->green=pixel[j]; break;
1928           case 2: q->blue=pixel[j]; break;
1929           case 3: q->opacity=pixel[j]; break;
1930           default: break;
1931         }
1932         q++;
1933         j++;
1934       }
1935       sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
1936       if (sync == MagickFalse)
1937         {
1938           status=MagickFalse;
1939           break;
1940         }
1941       j++;
1942     }
1943     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1944       {
1945         MagickBooleanType
1946           proceed;
1947
1948 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1949   #pragma omp critical (MagickCore_DespeckleImage)
1950 #endif
1951         proceed=SetImageProgress(image,DespeckleImageTag,channel,3);
1952         if (proceed == MagickFalse)
1953           status=MagickFalse;
1954       }
1955   }
1956   despeckle_view=DestroyCacheView(despeckle_view);
1957   image_view=DestroyCacheView(image_view);
1958   buffers=DestroyPixelThreadSet(buffers);
1959   pixels=DestroyPixelThreadSet(pixels);
1960   despeckle_image->type=image->type;
1961   if (status == MagickFalse)
1962     despeckle_image=DestroyImage(despeckle_image);
1963   return(despeckle_image);
1964 }
1965 \f
1966 /*
1967 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1968 %                                                                             %
1969 %                                                                             %
1970 %                                                                             %
1971 %     E d g e I m a g e                                                       %
1972 %                                                                             %
1973 %                                                                             %
1974 %                                                                             %
1975 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1976 %
1977 %  EdgeImage() finds edges in an image.  Radius defines the radius of the
1978 %  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
1979 %  radius for you.
1980 %
1981 %  The format of the EdgeImage method is:
1982 %
1983 %      Image *EdgeImage(const Image *image,const double radius,
1984 %        ExceptionInfo *exception)
1985 %
1986 %  A description of each parameter follows:
1987 %
1988 %    o image: the image.
1989 %
1990 %    o radius: the radius of the pixel neighborhood.
1991 %
1992 %    o exception: return any errors or warnings in this structure.
1993 %
1994 */
1995 MagickExport Image *EdgeImage(const Image *image,const double radius,
1996   ExceptionInfo *exception)
1997 {
1998   Image
1999     *edge_image;
2000
2001   double
2002     *kernel;
2003
2004   register long
2005     i;
2006
2007   unsigned long
2008     width;
2009
2010   assert(image != (const Image *) NULL);
2011   assert(image->signature == MagickSignature);
2012   if (image->debug != MagickFalse)
2013     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2014   assert(exception != (ExceptionInfo *) NULL);
2015   assert(exception->signature == MagickSignature);
2016   width=GetOptimalKernelWidth1D(radius,0.5);
2017   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2018   if (kernel == (double *) NULL)
2019     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2020   for (i=0; i < (long) (width*width); i++)
2021     kernel[i]=(-1.0);
2022   kernel[i/2]=(double) (width*width-1.0);
2023   edge_image=ConvolveImage(image,width,kernel,exception);
2024   kernel=(double *) RelinquishMagickMemory(kernel);
2025   return(edge_image);
2026 }
2027 \f
2028 /*
2029 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2030 %                                                                             %
2031 %                                                                             %
2032 %                                                                             %
2033 %     E m b o s s I m a g e                                                   %
2034 %                                                                             %
2035 %                                                                             %
2036 %                                                                             %
2037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2038 %
2039 %  EmbossImage() returns a grayscale image with a three-dimensional effect.
2040 %  We convolve the image with a Gaussian operator of the given radius and
2041 %  standard deviation (sigma).  For reasonable results, radius should be
2042 %  larger than sigma.  Use a radius of 0 and Emboss() selects a suitable
2043 %  radius for you.
2044 %
2045 %  The format of the EmbossImage method is:
2046 %
2047 %      Image *EmbossImage(const Image *image,const double radius,
2048 %        const double sigma,ExceptionInfo *exception)
2049 %
2050 %  A description of each parameter follows:
2051 %
2052 %    o image: the image.
2053 %
2054 %    o radius: the radius of the pixel neighborhood.
2055 %
2056 %    o sigma: the standard deviation of the Gaussian, in pixels.
2057 %
2058 %    o exception: return any errors or warnings in this structure.
2059 %
2060 */
2061 MagickExport Image *EmbossImage(const Image *image,const double radius,
2062   const double sigma,ExceptionInfo *exception)
2063 {
2064   double
2065     *kernel;
2066
2067   Image
2068     *emboss_image;
2069
2070   long
2071     j,
2072     k,
2073     u,
2074     v;
2075
2076   register long
2077     i;
2078
2079   unsigned long
2080     width;
2081
2082   assert(image != (Image *) NULL);
2083   assert(image->signature == MagickSignature);
2084   if (image->debug != MagickFalse)
2085     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2086   assert(exception != (ExceptionInfo *) NULL);
2087   assert(exception->signature == MagickSignature);
2088   width=GetOptimalKernelWidth2D(radius,sigma);
2089   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2090   if (kernel == (double *) NULL)
2091     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2092   j=(long) width/2;
2093   k=j;
2094   i=0;
2095   for (v=(-j); v <= j; v++)
2096   {
2097     for (u=(-j); u <= j; u++)
2098     {
2099       kernel[i]=((u < 0) || (v < 0) ? -8.0 : 8.0)*
2100         exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
2101         (2.0*MagickPI*MagickSigma*MagickSigma);
2102       if (u != k)
2103         kernel[i]=0.0;
2104       i++;
2105     }
2106     k--;
2107   }
2108   emboss_image=ConvolveImage(image,width,kernel,exception);
2109   if (emboss_image != (Image *) NULL)
2110     (void) EqualizeImage(emboss_image);
2111   kernel=(double *) RelinquishMagickMemory(kernel);
2112   return(emboss_image);
2113 }
2114 \f
2115 /*
2116 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2117 %                                                                             %
2118 %                                                                             %
2119 %                                                                             %
2120 %     G a u s s i a n B l u r I m a g e                                       %
2121 %                                                                             %
2122 %                                                                             %
2123 %                                                                             %
2124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2125 %
2126 %  GaussianBlurImage() blurs an image.  We convolve the image with a
2127 %  Gaussian operator of the given radius and standard deviation (sigma).
2128 %  For reasonable results, the radius should be larger than sigma.  Use a
2129 %  radius of 0 and GaussianBlurImage() selects a suitable radius for you
2130 %
2131 %  The format of the GaussianBlurImage method is:
2132 %
2133 %      Image *GaussianBlurImage(const Image *image,onst double radius,
2134 %        const double sigma,ExceptionInfo *exception)
2135 %      Image *GaussianBlurImageChannel(const Image *image,
2136 %        const ChannelType channel,const double radius,const double sigma,
2137 %        ExceptionInfo *exception)
2138 %
2139 %  A description of each parameter follows:
2140 %
2141 %    o image: the image.
2142 %
2143 %    o channel: the channel type.
2144 %
2145 %    o radius: the radius of the Gaussian, in pixels, not counting the center
2146 %      pixel.
2147 %
2148 %    o sigma: the standard deviation of the Gaussian, in pixels.
2149 %
2150 %    o exception: return any errors or warnings in this structure.
2151 %
2152 */
2153
2154 MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
2155   const double sigma,ExceptionInfo *exception)
2156 {
2157   Image
2158     *blur_image;
2159
2160   blur_image=GaussianBlurImageChannel(image,DefaultChannels,radius,sigma,
2161     exception);
2162   return(blur_image);
2163 }
2164
2165 MagickExport Image *GaussianBlurImageChannel(const Image *image,
2166   const ChannelType channel,const double radius,const double sigma,
2167   ExceptionInfo *exception)
2168 {
2169   double
2170     *kernel;
2171
2172   Image
2173     *blur_image;
2174
2175   long
2176     j,
2177     u,
2178     v;
2179
2180   register long
2181     i;
2182
2183   unsigned long
2184     width;
2185
2186   assert(image != (const Image *) NULL);
2187   assert(image->signature == MagickSignature);
2188   if (image->debug != MagickFalse)
2189     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2190   assert(exception != (ExceptionInfo *) NULL);
2191   assert(exception->signature == MagickSignature);
2192   width=GetOptimalKernelWidth2D(radius,sigma);
2193   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2194   if (kernel == (double *) NULL)
2195     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2196   j=(long) width/2;
2197   i=0;
2198   for (v=(-j); v <= j; v++)
2199   {
2200     for (u=(-j); u <= j; u++)
2201       kernel[i++]=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
2202         (2.0*MagickPI*MagickSigma*MagickSigma);
2203   }
2204   blur_image=ConvolveImageChannel(image,channel,width,kernel,exception);
2205   kernel=(double *) RelinquishMagickMemory(kernel);
2206   return(blur_image);
2207 }
2208 \f
2209 /*
2210 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2211 %                                                                             %
2212 %                                                                             %
2213 %                                                                             %
2214 %     M e d i a n F i l t e r I m a g e                                       %
2215 %                                                                             %
2216 %                                                                             %
2217 %                                                                             %
2218 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2219 %
2220 %  MedianFilterImage() applies a digital filter that improves the quality
2221 %  of a noisy image.  Each pixel is replaced by the median in a set of
2222 %  neighboring pixels as defined by radius.
2223 %
2224 %  The algorithm was contributed by Mike Edmonds and implements an insertion
2225 %  sort for selecting median color-channel values.  For more on this algorithm
2226 %  see "Skip Lists: A probabilistic Alternative to Balanced Trees" by William
2227 %  Pugh in the June 1990 of Communications of the ACM.
2228 %
2229 %  The format of the MedianFilterImage method is:
2230 %
2231 %      Image *MedianFilterImage(const Image *image,const double radius,
2232 %        ExceptionInfo *exception)
2233 %
2234 %  A description of each parameter follows:
2235 %
2236 %    o image: the image.
2237 %
2238 %    o radius: the radius of the pixel neighborhood.
2239 %
2240 %    o exception: return any errors or warnings in this structure.
2241 %
2242 */
2243
2244 #define MedianListChannels  5
2245
2246 typedef struct _MedianListNode
2247 {
2248   unsigned long
2249     next[9],
2250     count,
2251     signature;
2252 } MedianListNode;
2253
2254 typedef struct _MedianSkipList
2255 {
2256   long
2257     level;
2258
2259   MedianListNode
2260     *nodes;
2261 } MedianSkipList;
2262
2263 typedef struct _MedianPixelList
2264 {
2265   unsigned long
2266     center,
2267     seed,
2268     signature;
2269
2270   MedianSkipList
2271     lists[MedianListChannels];
2272 } MedianPixelList;
2273
2274 static MedianPixelList *DestroyMedianPixelList(MedianPixelList *pixel_list)
2275 {
2276   register long
2277     i;
2278
2279   if (pixel_list == (MedianPixelList *) NULL)
2280     return((MedianPixelList *) NULL);
2281   for (i=0; i < MedianListChannels; i++)
2282     if (pixel_list->lists[i].nodes != (MedianListNode *) NULL)
2283       pixel_list->lists[i].nodes=(MedianListNode *) RelinquishMagickMemory(
2284         pixel_list->lists[i].nodes);
2285   pixel_list=(MedianPixelList *) RelinquishAlignedMemory(pixel_list);
2286   return(pixel_list);
2287 }
2288
2289 static MedianPixelList **DestroyMedianPixelListThreadSet(
2290   MedianPixelList **pixel_list)
2291 {
2292   register long
2293     i;
2294
2295   assert(pixel_list != (MedianPixelList **) NULL);
2296   for (i=0; i < (long) GetOpenMPMaximumThreads(); i++)
2297     if (pixel_list[i] != (MedianPixelList *) NULL)
2298       pixel_list[i]=DestroyMedianPixelList(pixel_list[i]);
2299   pixel_list=(MedianPixelList **) RelinquishAlignedMemory(pixel_list);
2300   return(pixel_list);
2301 }
2302
2303 static MedianPixelList *AcquireMedianPixelList(const unsigned long width)
2304 {
2305   MedianPixelList
2306     *pixel_list;
2307
2308   register long
2309     i;
2310
2311   pixel_list=(MedianPixelList *) AcquireAlignedMemory(1,sizeof(*pixel_list));
2312   if (pixel_list == (MedianPixelList *) NULL)
2313     return(pixel_list);
2314   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
2315   pixel_list->center=width*width/2;
2316   for (i=0; i < MedianListChannels; i++)
2317   {
2318     pixel_list->lists[i].nodes=(MedianListNode *) AcquireQuantumMemory(65537UL,
2319       sizeof(*pixel_list->lists[i].nodes));
2320     if (pixel_list->lists[i].nodes == (MedianListNode *) NULL)
2321       return(DestroyMedianPixelList(pixel_list));
2322     (void) ResetMagickMemory(pixel_list->lists[i].nodes,0,65537UL*
2323       sizeof(*pixel_list->lists[i].nodes));
2324   }
2325   pixel_list->signature=MagickSignature;
2326   return(pixel_list);
2327 }
2328
2329 static MedianPixelList **AcquireMedianPixelListThreadSet(
2330   const unsigned long width)
2331 {
2332   register long
2333     i;
2334
2335   MedianPixelList
2336     **pixel_list;
2337
2338   unsigned long
2339     number_threads;
2340
2341   number_threads=GetOpenMPMaximumThreads();
2342   pixel_list=(MedianPixelList **) AcquireAlignedMemory(number_threads,
2343     sizeof(*pixel_list));
2344   if (pixel_list == (MedianPixelList **) NULL)
2345     return((MedianPixelList **) NULL);
2346   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
2347   for (i=0; i < (long) number_threads; i++)
2348   {
2349     pixel_list[i]=AcquireMedianPixelList(width);
2350     if (pixel_list[i] == (MedianPixelList *) NULL)
2351       return(DestroyMedianPixelListThreadSet(pixel_list));
2352   }
2353   return(pixel_list);
2354 }
2355
2356 static void AddNodeMedianPixelList(MedianPixelList *pixel_list,
2357   const long channel,const unsigned long color)
2358 {
2359   register long
2360     level;
2361
2362   register MedianSkipList
2363     *list;
2364
2365   unsigned long
2366     search,
2367     update[9];
2368
2369   /*
2370     Initialize the node.
2371   */
2372   list=pixel_list->lists+channel;
2373   list->nodes[color].signature=pixel_list->signature;
2374   list->nodes[color].count=1;
2375   /*
2376     Determine where it belongs in the list.
2377   */
2378   search=65536UL;
2379   for (level=list->level; level >= 0; level--)
2380   {
2381     while (list->nodes[search].next[level] < color)
2382       search=list->nodes[search].next[level];
2383     update[level]=search;
2384   }
2385   /*
2386     Generate a pseudo-random level for this node.
2387   */
2388   for (level=0; ; level++)
2389   {
2390     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2391     if ((pixel_list->seed & 0x300) != 0x300)
2392       break;
2393   }
2394   if (level > 8)
2395     level=8;
2396   if (level > (list->level+2))
2397     level=list->level+2;
2398   /*
2399     If we're raising the list's level, link back to the root node.
2400   */
2401   while (level > list->level)
2402   {
2403     list->level++;
2404     update[list->level]=65536UL;
2405   }
2406   /*
2407     Link the node into the skip-list.
2408   */
2409   do
2410   {
2411     list->nodes[color].next[level]=list->nodes[update[level]].next[level];
2412     list->nodes[update[level]].next[level]=color;
2413   }
2414   while (level-- > 0);
2415 }
2416
2417 static MagickPixelPacket GetMedianPixelList(MedianPixelList *pixel_list)
2418 {
2419   MagickPixelPacket
2420     pixel;
2421
2422   register long
2423     channel;
2424
2425   register MedianSkipList
2426     *list;
2427
2428   unsigned long
2429     center,
2430     color,
2431     count;
2432
2433   unsigned short
2434     channels[MedianListChannels];
2435
2436   /*
2437     Find the median value for each of the color.
2438   */
2439   center=pixel_list->center;
2440   for (channel=0; channel < 5; channel++)
2441   {
2442     list=pixel_list->lists+channel;
2443     color=65536UL;
2444     count=0;
2445     do
2446     {
2447       color=list->nodes[color].next[0];
2448       count+=list->nodes[color].count;
2449     }
2450     while (count <= center);
2451     channels[channel]=(unsigned short) color;
2452   }
2453   GetMagickPixelPacket((const Image *) NULL,&pixel);
2454   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
2455   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
2456   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
2457   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
2458   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
2459   return(pixel);
2460 }
2461
2462 static inline void InsertMedianPixelList(const Image *image,
2463   const PixelPacket *pixel,const IndexPacket *indexes,
2464   MedianPixelList *pixel_list)
2465 {
2466   unsigned long
2467     signature;
2468
2469   unsigned short
2470     index;
2471
2472   index=ScaleQuantumToShort(pixel->red);
2473   signature=pixel_list->lists[0].nodes[index].signature;
2474   if (signature == pixel_list->signature)
2475     pixel_list->lists[0].nodes[index].count++;
2476   else
2477     AddNodeMedianPixelList(pixel_list,0,index);
2478   index=ScaleQuantumToShort(pixel->green);
2479   signature=pixel_list->lists[1].nodes[index].signature;
2480   if (signature == pixel_list->signature)
2481     pixel_list->lists[1].nodes[index].count++;
2482   else
2483     AddNodeMedianPixelList(pixel_list,1,index);
2484   index=ScaleQuantumToShort(pixel->blue);
2485   signature=pixel_list->lists[2].nodes[index].signature;
2486   if (signature == pixel_list->signature)
2487     pixel_list->lists[2].nodes[index].count++;
2488   else
2489     AddNodeMedianPixelList(pixel_list,2,index);
2490   index=ScaleQuantumToShort(pixel->opacity);
2491   signature=pixel_list->lists[3].nodes[index].signature;
2492   if (signature == pixel_list->signature)
2493     pixel_list->lists[3].nodes[index].count++;
2494   else
2495     AddNodeMedianPixelList(pixel_list,3,index);
2496   if (image->colorspace == CMYKColorspace)
2497     index=ScaleQuantumToShort(*indexes);
2498   signature=pixel_list->lists[4].nodes[index].signature;
2499   if (signature == pixel_list->signature)
2500     pixel_list->lists[4].nodes[index].count++;
2501   else
2502     AddNodeMedianPixelList(pixel_list,4,index);
2503 }
2504
2505 static void ResetMedianPixelList(MedianPixelList *pixel_list)
2506 {
2507   int
2508     level;
2509
2510   register long
2511     channel;
2512
2513   register MedianListNode
2514     *root;
2515
2516   register MedianSkipList
2517     *list;
2518
2519   /*
2520     Reset the skip-list.
2521   */
2522   for (channel=0; channel < 5; channel++)
2523   {
2524     list=pixel_list->lists+channel;
2525     root=list->nodes+65536UL;
2526     list->level=0;
2527     for (level=0; level < 9; level++)
2528       root->next[level]=65536UL;
2529   }
2530   pixel_list->seed=pixel_list->signature++;
2531 }
2532
2533 MagickExport Image *MedianFilterImage(const Image *image,const double radius,
2534   ExceptionInfo *exception)
2535 {
2536 #define MedianFilterImageTag  "MedianFilter/Image"
2537
2538   Image
2539     *median_image;
2540
2541   long
2542     progress,
2543     y;
2544
2545   MagickBooleanType
2546     status;
2547
2548   MedianPixelList
2549     **pixel_list;
2550
2551   unsigned long
2552     width;
2553
2554   CacheView
2555     *image_view,
2556     *median_view;
2557
2558   /*
2559     Initialize median image attributes.
2560   */
2561   assert(image != (Image *) NULL);
2562   assert(image->signature == MagickSignature);
2563   if (image->debug != MagickFalse)
2564     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2565   assert(exception != (ExceptionInfo *) NULL);
2566   assert(exception->signature == MagickSignature);
2567   width=GetOptimalKernelWidth2D(radius,0.5);
2568   if ((image->columns < width) || (image->rows < width))
2569     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
2570   median_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2571     exception);
2572   if (median_image == (Image *) NULL)
2573     return((Image *) NULL);
2574   if (SetImageStorageClass(median_image,DirectClass) == MagickFalse)
2575     {
2576       InheritException(exception,&median_image->exception);
2577       median_image=DestroyImage(median_image);
2578       return((Image *) NULL);
2579     }
2580   pixel_list=AcquireMedianPixelListThreadSet(width);
2581   if (pixel_list == (MedianPixelList **) NULL)
2582     {
2583       median_image=DestroyImage(median_image);
2584       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2585     }
2586   /*
2587     Median filter each image row.
2588   */
2589   status=MagickTrue;
2590   progress=0;
2591   image_view=AcquireCacheView(image);
2592   median_view=AcquireCacheView(median_image);
2593 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2594   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2595 #endif
2596   for (y=0; y < (long) median_image->rows; y++)
2597   {
2598     register const IndexPacket
2599       *restrict indexes;
2600
2601     register const PixelPacket
2602       *restrict p;
2603
2604     register IndexPacket
2605       *restrict median_indexes;
2606
2607     register long
2608       id,
2609       x;
2610
2611     register PixelPacket
2612       *restrict q;
2613
2614     if (status == MagickFalse)
2615       continue;
2616     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
2617       2L),image->columns+width,width,exception);
2618     q=QueueCacheViewAuthenticPixels(median_view,0,y,median_image->columns,1,
2619       exception);
2620     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2621       {
2622         status=MagickFalse;
2623         continue;
2624       }
2625     indexes=GetCacheViewVirtualIndexQueue(image_view);
2626     median_indexes=GetCacheViewAuthenticIndexQueue(median_view);
2627     id=GetOpenMPThreadId();
2628     for (x=0; x < (long) median_image->columns; x++)
2629     {
2630       MagickPixelPacket
2631         pixel;
2632
2633       register const PixelPacket
2634         *restrict r;
2635
2636       register const IndexPacket
2637         *restrict s;
2638
2639       register long
2640         u,
2641         v;
2642
2643       r=p;
2644       s=indexes+x;
2645       ResetMedianPixelList(pixel_list[id]);
2646       for (v=0; v < (long) width; v++)
2647       {
2648         for (u=0; u < (long) width; u++)
2649           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
2650         r+=image->columns+width;
2651         s+=image->columns+width;
2652       }
2653       pixel=GetMedianPixelList(pixel_list[id]);
2654       SetPixelPacket(median_image,&pixel,q,median_indexes+x);
2655       p++;
2656       q++;
2657     }
2658     if (SyncCacheViewAuthenticPixels(median_view,exception) == MagickFalse)
2659       status=MagickFalse;
2660     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2661       {
2662         MagickBooleanType
2663           proceed;
2664
2665 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2666   #pragma omp critical (MagickCore_MedianFilterImage)
2667 #endif
2668         proceed=SetImageProgress(image,MedianFilterImageTag,progress++,
2669           image->rows);
2670         if (proceed == MagickFalse)
2671           status=MagickFalse;
2672       }
2673   }
2674   median_view=DestroyCacheView(median_view);
2675   image_view=DestroyCacheView(image_view);
2676   pixel_list=DestroyMedianPixelListThreadSet(pixel_list);
2677   return(median_image);
2678 }
2679 \f
2680 /*
2681 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2682 %                                                                             %
2683 %                                                                             %
2684 %                                                                             %
2685 %     M o t i o n B l u r I m a g e                                           %
2686 %                                                                             %
2687 %                                                                             %
2688 %                                                                             %
2689 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2690 %
2691 %  MotionBlurImage() simulates motion blur.  We convolve the image with a
2692 %  Gaussian operator of the given radius and standard deviation (sigma).
2693 %  For reasonable results, radius should be larger than sigma.  Use a
2694 %  radius of 0 and MotionBlurImage() selects a suitable radius for you.
2695 %  Angle gives the angle of the blurring motion.
2696 %
2697 %  Andrew Protano contributed this effect.
2698 %
2699 %  The format of the MotionBlurImage method is:
2700 %
2701 %    Image *MotionBlurImage(const Image *image,const double radius,
2702 %      const double sigma,const double angle,ExceptionInfo *exception)
2703 %    Image *MotionBlurImageChannel(const Image *image,const ChannelType channel,
2704 %      const double radius,const double sigma,const double angle,
2705 %      ExceptionInfo *exception)
2706 %
2707 %  A description of each parameter follows:
2708 %
2709 %    o image: the image.
2710 %
2711 %    o channel: the channel type.
2712 %
2713 %    o radius: the radius of the Gaussian, in pixels, not counting the center
2714 %    o radius: the radius of the Gaussian, in pixels, not counting
2715 %      the center pixel.
2716 %
2717 %    o sigma: the standard deviation of the Gaussian, in pixels.
2718 %
2719 %    o angle: Apply the effect along this angle.
2720 %
2721 %    o exception: return any errors or warnings in this structure.
2722 %
2723 */
2724
2725 static double *GetMotionBlurKernel(const unsigned long width,const double sigma)
2726 {
2727   double
2728     *kernel,
2729     normalize;
2730
2731   register long
2732     i;
2733
2734   /*
2735    Generate a 1-D convolution kernel.
2736   */
2737   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
2738   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
2739   if (kernel == (double *) NULL)
2740     return(kernel);
2741   normalize=0.0;
2742   for (i=0; i < (long) width; i++)
2743   {
2744     kernel[i]=exp((-((double) i*i)/(double) (2.0*MagickSigma*MagickSigma)))/
2745       (MagickSQ2PI*MagickSigma);
2746     normalize+=kernel[i];
2747   }
2748   for (i=0; i < (long) width; i++)
2749     kernel[i]/=normalize;
2750   return(kernel);
2751 }
2752
2753 MagickExport Image *MotionBlurImage(const Image *image,const double radius,
2754   const double sigma,const double angle,ExceptionInfo *exception)
2755 {
2756   Image
2757     *motion_blur;
2758
2759   motion_blur=MotionBlurImageChannel(image,DefaultChannels,radius,sigma,angle,
2760     exception);
2761   return(motion_blur);
2762 }
2763
2764 MagickExport Image *MotionBlurImageChannel(const Image *image,
2765   const ChannelType channel,const double radius,const double sigma,
2766   const double angle,ExceptionInfo *exception)
2767 {
2768   typedef struct _OffsetInfo
2769   {
2770     long
2771       x,
2772       y;
2773   } OffsetInfo;
2774
2775   double
2776     *kernel;
2777
2778   Image
2779     *blur_image;
2780
2781   long
2782     progress,
2783     y;
2784
2785   MagickBooleanType
2786     status;
2787
2788   MagickPixelPacket
2789     bias;
2790
2791   OffsetInfo
2792     *offset;
2793
2794   PointInfo
2795     point;
2796
2797   register long
2798     i;
2799
2800   unsigned long
2801     width;
2802
2803   CacheView
2804     *blur_view,
2805     *image_view;
2806
2807   assert(image != (Image *) NULL);
2808   assert(image->signature == MagickSignature);
2809   if (image->debug != MagickFalse)
2810     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2811   assert(exception != (ExceptionInfo *) NULL);
2812   width=GetOptimalKernelWidth1D(radius,sigma);
2813   kernel=GetMotionBlurKernel(width,sigma);
2814   if (kernel == (double *) NULL)
2815     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2816   offset=(OffsetInfo *) AcquireQuantumMemory(width,sizeof(*offset));
2817   if (offset == (OffsetInfo *) NULL)
2818     {
2819       kernel=(double *) RelinquishMagickMemory(kernel);
2820       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2821     }
2822   blur_image=CloneImage(image,0,0,MagickTrue,exception);
2823   if (blur_image == (Image *) NULL)
2824     {
2825       kernel=(double *) RelinquishMagickMemory(kernel);
2826       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2827       return((Image *) NULL);
2828     }
2829   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
2830     {
2831       kernel=(double *) RelinquishMagickMemory(kernel);
2832       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2833       InheritException(exception,&blur_image->exception);
2834       blur_image=DestroyImage(blur_image);
2835       return((Image *) NULL);
2836     }
2837   point.x=(double) width*sin(DegreesToRadians(angle));
2838   point.y=(double) width*cos(DegreesToRadians(angle));
2839   for (i=0; i < (long) width; i++)
2840   {
2841     offset[i].x=(long) ((i*point.y)/hypot(point.x,point.y)+0.5);
2842     offset[i].y=(long) ((i*point.x)/hypot(point.x,point.y)+0.5);
2843   }
2844   /*
2845     Motion blur image.
2846   */
2847   status=MagickTrue;
2848   progress=0;
2849   GetMagickPixelPacket(image,&bias);
2850   image_view=AcquireCacheView(image);
2851   blur_view=AcquireCacheView(blur_image);
2852 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2853   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2854 #endif
2855   for (y=0; y < (long) image->rows; y++)
2856   {
2857     register IndexPacket
2858       *restrict blur_indexes;
2859
2860     register long
2861       x;
2862
2863     register PixelPacket
2864       *restrict q;
2865
2866     if (status == MagickFalse)
2867       continue;
2868     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
2869       exception);
2870     if (q == (PixelPacket *) NULL)
2871       {
2872         status=MagickFalse;
2873         continue;
2874       }
2875     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
2876     for (x=0; x < (long) image->columns; x++)
2877     {
2878       MagickPixelPacket
2879         qixel;
2880
2881       PixelPacket
2882         pixel;
2883
2884       register double
2885         *restrict k;
2886
2887       register long
2888         i;
2889
2890       register const IndexPacket
2891         *restrict indexes;
2892
2893       k=kernel;
2894       qixel=bias;
2895       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
2896         {
2897           for (i=0; i < (long) width; i++)
2898           {
2899             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
2900               offset[i].y,&pixel,exception);
2901             qixel.red+=(*k)*pixel.red;
2902             qixel.green+=(*k)*pixel.green;
2903             qixel.blue+=(*k)*pixel.blue;
2904             qixel.opacity+=(*k)*pixel.opacity;
2905             if (image->colorspace == CMYKColorspace)
2906               {
2907                 indexes=GetCacheViewVirtualIndexQueue(image_view);
2908                 qixel.index+=(*k)*(*indexes);
2909               }
2910             k++;
2911           }
2912           if ((channel & RedChannel) != 0)
2913             q->red=RoundToQuantum(qixel.red);
2914           if ((channel & GreenChannel) != 0)
2915             q->green=RoundToQuantum(qixel.green);
2916           if ((channel & BlueChannel) != 0)
2917             q->blue=RoundToQuantum(qixel.blue);
2918           if ((channel & OpacityChannel) != 0)
2919             q->opacity=RoundToQuantum(qixel.opacity);
2920           if (((channel & IndexChannel) != 0) &&
2921               (image->colorspace == CMYKColorspace))
2922             blur_indexes[x]=(IndexPacket) RoundToQuantum(qixel.index);
2923         }
2924       else
2925         {
2926           MagickRealType
2927             alpha,
2928             gamma;
2929
2930           alpha=0.0;
2931           gamma=0.0;
2932           for (i=0; i < (long) width; i++)
2933           {
2934             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
2935               offset[i].y,&pixel,exception);
2936             alpha=(MagickRealType) (QuantumScale*(QuantumRange-pixel.opacity));
2937             qixel.red+=(*k)*alpha*pixel.red;
2938             qixel.green+=(*k)*alpha*pixel.green;
2939             qixel.blue+=(*k)*alpha*pixel.blue;
2940             qixel.opacity+=(*k)*pixel.opacity;
2941             if (image->colorspace == CMYKColorspace)
2942               {
2943                 indexes=GetCacheViewVirtualIndexQueue(image_view);
2944                 qixel.index+=(*k)*alpha*(*indexes);
2945               }
2946             gamma+=(*k)*alpha;
2947             k++;
2948           }
2949           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2950           if ((channel & RedChannel) != 0)
2951             q->red=RoundToQuantum(gamma*qixel.red);
2952           if ((channel & GreenChannel) != 0)
2953             q->green=RoundToQuantum(gamma*qixel.green);
2954           if ((channel & BlueChannel) != 0)
2955             q->blue=RoundToQuantum(gamma*qixel.blue);
2956           if ((channel & OpacityChannel) != 0)
2957             q->opacity=RoundToQuantum(qixel.opacity);
2958           if (((channel & IndexChannel) != 0) &&
2959               (image->colorspace == CMYKColorspace))
2960             blur_indexes[x]=(IndexPacket) RoundToQuantum(gamma*qixel.index);
2961         }
2962       q++;
2963     }
2964     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
2965       status=MagickFalse;
2966     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2967       {
2968         MagickBooleanType
2969           proceed;
2970
2971 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2972   #pragma omp critical (MagickCore_MotionBlurImageChannel)
2973 #endif
2974         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
2975         if (proceed == MagickFalse)
2976           status=MagickFalse;
2977       }
2978   }
2979   blur_view=DestroyCacheView(blur_view);
2980   image_view=DestroyCacheView(image_view);
2981   kernel=(double *) RelinquishMagickMemory(kernel);
2982   offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2983   if (status == MagickFalse)
2984     blur_image=DestroyImage(blur_image);
2985   return(blur_image);
2986 }
2987 \f
2988 /*
2989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2990 %                                                                             %
2991 %                                                                             %
2992 %                                                                             %
2993 %     P r e v i e w I m a g e                                                 %
2994 %                                                                             %
2995 %                                                                             %
2996 %                                                                             %
2997 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2998 %
2999 %  PreviewImage() tiles 9 thumbnails of the specified image with an image
3000 %  processing operation applied with varying parameters.  This may be helpful
3001 %  pin-pointing an appropriate parameter for a particular image processing
3002 %  operation.
3003 %
3004 %  The format of the PreviewImages method is:
3005 %
3006 %      Image *PreviewImages(const Image *image,const PreviewType preview,
3007 %        ExceptionInfo *exception)
3008 %
3009 %  A description of each parameter follows:
3010 %
3011 %    o image: the image.
3012 %
3013 %    o preview: the image processing operation.
3014 %
3015 %    o exception: return any errors or warnings in this structure.
3016 %
3017 */
3018 MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
3019   ExceptionInfo *exception)
3020 {
3021 #define NumberTiles  9
3022 #define PreviewImageTag  "Preview/Image"
3023 #define DefaultPreviewGeometry  "204x204+10+10"
3024
3025   char
3026     factor[MaxTextExtent],
3027     label[MaxTextExtent];
3028
3029   double
3030     degrees,
3031     gamma,
3032     percentage,
3033     radius,
3034     sigma,
3035     threshold;
3036
3037   Image
3038     *images,
3039     *montage_image,
3040     *preview_image,
3041     *thumbnail;
3042
3043   ImageInfo
3044     *preview_info;
3045
3046   long
3047     y;
3048
3049   MagickBooleanType
3050     proceed;
3051
3052   MontageInfo
3053     *montage_info;
3054
3055   QuantizeInfo
3056     quantize_info;
3057
3058   RectangleInfo
3059     geometry;
3060
3061   register long
3062     i,
3063     x;
3064
3065   unsigned long
3066     colors;
3067
3068   /*
3069     Open output image file.
3070   */
3071   assert(image != (Image *) NULL);
3072   assert(image->signature == MagickSignature);
3073   if (image->debug != MagickFalse)
3074     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3075   colors=2;
3076   degrees=0.0;
3077   gamma=(-0.2f);
3078   preview_info=AcquireImageInfo();
3079   SetGeometry(image,&geometry);
3080   (void) ParseMetaGeometry(DefaultPreviewGeometry,&geometry.x,&geometry.y,
3081     &geometry.width,&geometry.height);
3082   images=NewImageList();
3083   percentage=12.5;
3084   GetQuantizeInfo(&quantize_info);
3085   radius=0.0;
3086   sigma=1.0;
3087   threshold=0.0;
3088   x=0;
3089   y=0;
3090   for (i=0; i < NumberTiles; i++)
3091   {
3092     thumbnail=ThumbnailImage(image,geometry.width,geometry.height,exception);
3093     if (thumbnail == (Image *) NULL)
3094       break;
3095     (void) SetImageProgressMonitor(thumbnail,(MagickProgressMonitor) NULL,
3096       (void *) NULL);
3097     (void) SetImageProperty(thumbnail,"label",DefaultTileLabel);
3098     if (i == (NumberTiles/2))
3099       {
3100         (void) QueryColorDatabase("#dfdfdf",&thumbnail->matte_color,exception);
3101         AppendImageToList(&images,thumbnail);
3102         continue;
3103       }
3104     switch (preview)
3105     {
3106       case RotatePreview:
3107       {
3108         degrees+=45.0;
3109         preview_image=RotateImage(thumbnail,degrees,exception);
3110         (void) FormatMagickString(label,MaxTextExtent,"rotate %g",degrees);
3111         break;
3112       }
3113       case ShearPreview:
3114       {
3115         degrees+=5.0;
3116         preview_image=ShearImage(thumbnail,degrees,degrees,exception);
3117         (void) FormatMagickString(label,MaxTextExtent,"shear %gx%g",
3118           degrees,2.0*degrees);
3119         break;
3120       }
3121       case RollPreview:
3122       {
3123         x=(long) ((i+1)*thumbnail->columns)/NumberTiles;
3124         y=(long) ((i+1)*thumbnail->rows)/NumberTiles;
3125         preview_image=RollImage(thumbnail,x,y,exception);
3126         (void) FormatMagickString(label,MaxTextExtent,"roll %ldx%ld",x,y);
3127         break;
3128       }
3129       case HuePreview:
3130       {
3131         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3132         if (preview_image == (Image *) NULL)
3133           break;
3134         (void) FormatMagickString(factor,MaxTextExtent,"100,100,%g",
3135           2.0*percentage);
3136         (void) ModulateImage(preview_image,factor);
3137         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3138         break;
3139       }
3140       case SaturationPreview:
3141       {
3142         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3143         if (preview_image == (Image *) NULL)
3144           break;
3145         (void) FormatMagickString(factor,MaxTextExtent,"100,%g",2.0*percentage);
3146         (void) ModulateImage(preview_image,factor);
3147         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3148         break;
3149       }
3150       case BrightnessPreview:
3151       {
3152         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3153         if (preview_image == (Image *) NULL)
3154           break;
3155         (void) FormatMagickString(factor,MaxTextExtent,"%g",2.0*percentage);
3156         (void) ModulateImage(preview_image,factor);
3157         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3158         break;
3159       }
3160       case GammaPreview:
3161       default:
3162       {
3163         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3164         if (preview_image == (Image *) NULL)
3165           break;
3166         gamma+=0.4f;
3167         (void) GammaImageChannel(preview_image,DefaultChannels,gamma);
3168         (void) FormatMagickString(label,MaxTextExtent,"gamma %g",gamma);
3169         break;
3170       }
3171       case SpiffPreview:
3172       {
3173         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3174         if (preview_image != (Image *) NULL)
3175           for (x=0; x < i; x++)
3176             (void) ContrastImage(preview_image,MagickTrue);
3177         (void) FormatMagickString(label,MaxTextExtent,"contrast (%ld)",i+1);
3178         break;
3179       }
3180       case DullPreview:
3181       {
3182         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3183         if (preview_image == (Image *) NULL)
3184           break;
3185         for (x=0; x < i; x++)
3186           (void) ContrastImage(preview_image,MagickFalse);
3187         (void) FormatMagickString(label,MaxTextExtent,"+contrast (%ld)",i+1);
3188         break;
3189       }
3190       case GrayscalePreview:
3191       {
3192         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3193         if (preview_image == (Image *) NULL)
3194           break;
3195         colors<<=1;
3196         quantize_info.number_colors=colors;
3197         quantize_info.colorspace=GRAYColorspace;
3198         (void) QuantizeImage(&quantize_info,preview_image);
3199         (void) FormatMagickString(label,MaxTextExtent,
3200           "-colorspace gray -colors %ld",colors);
3201         break;
3202       }
3203       case QuantizePreview:
3204       {
3205         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3206         if (preview_image == (Image *) NULL)
3207           break;
3208         colors<<=1;
3209         quantize_info.number_colors=colors;
3210         (void) QuantizeImage(&quantize_info,preview_image);
3211         (void) FormatMagickString(label,MaxTextExtent,"colors %ld",colors);
3212         break;
3213       }
3214       case DespecklePreview:
3215       {
3216         for (x=0; x < (i-1); x++)
3217         {
3218           preview_image=DespeckleImage(thumbnail,exception);
3219           if (preview_image == (Image *) NULL)
3220             break;
3221           thumbnail=DestroyImage(thumbnail);
3222           thumbnail=preview_image;
3223         }
3224         preview_image=DespeckleImage(thumbnail,exception);
3225         if (preview_image == (Image *) NULL)
3226           break;
3227         (void) FormatMagickString(label,MaxTextExtent,"despeckle (%ld)",i+1);
3228         break;
3229       }
3230       case ReduceNoisePreview:
3231       {
3232         preview_image=ReduceNoiseImage(thumbnail,radius,exception);
3233         (void) FormatMagickString(label,MaxTextExtent,"noise %g",radius);
3234         break;
3235       }
3236       case AddNoisePreview:
3237       {
3238         switch ((int) i)
3239         {
3240           case 0:
3241           {
3242             (void) CopyMagickString(factor,"uniform",MaxTextExtent);
3243             break;
3244           }
3245           case 1:
3246           {
3247             (void) CopyMagickString(factor,"gaussian",MaxTextExtent);
3248             break;
3249           }
3250           case 2:
3251           {
3252             (void) CopyMagickString(factor,"multiplicative",MaxTextExtent);
3253             break;
3254           }
3255           case 3:
3256           {
3257             (void) CopyMagickString(factor,"impulse",MaxTextExtent);
3258             break;
3259           }
3260           case 4:
3261           {
3262             (void) CopyMagickString(factor,"laplacian",MaxTextExtent);
3263             break;
3264           }
3265           case 5:
3266           {
3267             (void) CopyMagickString(factor,"Poisson",MaxTextExtent);
3268             break;
3269           }
3270           default:
3271           {
3272             (void) CopyMagickString(thumbnail->magick,"NULL",MaxTextExtent);
3273             break;
3274           }
3275         }
3276         preview_image=ReduceNoiseImage(thumbnail,(double) i,exception);
3277         (void) FormatMagickString(label,MaxTextExtent,"+noise %s",factor);
3278         break;
3279       }
3280       case SharpenPreview:
3281       {
3282         preview_image=SharpenImage(thumbnail,radius,sigma,exception);
3283         (void) FormatMagickString(label,MaxTextExtent,"sharpen %gx%g",radius,
3284           sigma);
3285         break;
3286       }
3287       case BlurPreview:
3288       {
3289         preview_image=BlurImage(thumbnail,radius,sigma,exception);
3290         (void) FormatMagickString(label,MaxTextExtent,"blur %gx%g",radius,
3291           sigma);
3292         break;
3293       }
3294       case ThresholdPreview:
3295       {
3296         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3297         if (preview_image == (Image *) NULL)
3298           break;
3299         (void) BilevelImage(thumbnail,
3300           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
3301         (void) FormatMagickString(label,MaxTextExtent,"threshold %g",
3302           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
3303         break;
3304       }
3305       case EdgeDetectPreview:
3306       {
3307         preview_image=EdgeImage(thumbnail,radius,exception);
3308         (void) FormatMagickString(label,MaxTextExtent,"edge %g",radius);
3309         break;
3310       }
3311       case SpreadPreview:
3312       {
3313         preview_image=SpreadImage(thumbnail,radius,exception);
3314         (void) FormatMagickString(label,MaxTextExtent,"spread %g",radius+0.5);
3315         break;
3316       }
3317       case SolarizePreview:
3318       {
3319         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3320         if (preview_image == (Image *) NULL)
3321           break;
3322         (void) SolarizeImage(preview_image,(double) QuantumRange*
3323           percentage/100.0);
3324         (void) FormatMagickString(label,MaxTextExtent,"solarize %g",
3325           (QuantumRange*percentage)/100.0);
3326         break;
3327       }
3328       case ShadePreview:
3329       {
3330         degrees+=10.0;
3331         preview_image=ShadeImage(thumbnail,MagickTrue,degrees,degrees,
3332           exception);
3333         (void) FormatMagickString(label,MaxTextExtent,"shade %gx%g",degrees,
3334           degrees);
3335         break;
3336       }
3337       case RaisePreview:
3338       {
3339         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3340         if (preview_image == (Image *) NULL)
3341           break;
3342         geometry.width=(unsigned long) (2*i+2);
3343         geometry.height=(unsigned long) (2*i+2);
3344         geometry.x=i/2;
3345         geometry.y=i/2;
3346         (void) RaiseImage(preview_image,&geometry,MagickTrue);
3347         (void) FormatMagickString(label,MaxTextExtent,"raise %lux%lu%+ld%+ld",
3348           geometry.width,geometry.height,geometry.x,geometry.y);
3349         break;
3350       }
3351       case SegmentPreview:
3352       {
3353         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3354         if (preview_image == (Image *) NULL)
3355           break;
3356         threshold+=0.4f;
3357         (void) SegmentImage(preview_image,RGBColorspace,MagickFalse,threshold,
3358           threshold);
3359         (void) FormatMagickString(label,MaxTextExtent,"segment %gx%g",
3360           threshold,threshold);
3361         break;
3362       }
3363       case SwirlPreview:
3364       {
3365         preview_image=SwirlImage(thumbnail,degrees,exception);
3366         (void) FormatMagickString(label,MaxTextExtent,"swirl %g",degrees);
3367         degrees+=45.0;
3368         break;
3369       }
3370       case ImplodePreview:
3371       {
3372         degrees+=0.1f;
3373         preview_image=ImplodeImage(thumbnail,degrees,exception);
3374         (void) FormatMagickString(label,MaxTextExtent,"implode %g",degrees);
3375         break;
3376       }
3377       case WavePreview:
3378       {
3379         degrees+=5.0f;
3380         preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,exception);
3381         (void) FormatMagickString(label,MaxTextExtent,"wave %gx%g",0.5*degrees,
3382           2.0*degrees);
3383         break;
3384       }
3385       case OilPaintPreview:
3386       {
3387         preview_image=OilPaintImage(thumbnail,(double) radius,exception);
3388         (void) FormatMagickString(label,MaxTextExtent,"paint %g",radius);
3389         break;
3390       }
3391       case CharcoalDrawingPreview:
3392       {
3393         preview_image=CharcoalImage(thumbnail,(double) radius,(double) sigma,
3394           exception);
3395         (void) FormatMagickString(label,MaxTextExtent,"charcoal %gx%g",radius,
3396           sigma);
3397         break;
3398       }
3399       case JPEGPreview:
3400       {
3401         char
3402           filename[MaxTextExtent];
3403
3404         int
3405           file;
3406
3407         MagickBooleanType
3408           status;
3409
3410         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3411         if (preview_image == (Image *) NULL)
3412           break;
3413         preview_info->quality=(unsigned long) percentage;
3414         (void) FormatMagickString(factor,MaxTextExtent,"%lu",
3415           preview_info->quality);
3416         file=AcquireUniqueFileResource(filename);
3417         if (file != -1)
3418           file=close(file)-1;
3419         (void) FormatMagickString(preview_image->filename,MaxTextExtent,
3420           "jpeg:%s",filename);
3421         status=WriteImage(preview_info,preview_image);
3422         if (status != MagickFalse)
3423           {
3424             Image
3425               *quality_image;
3426
3427             (void) CopyMagickString(preview_info->filename,
3428               preview_image->filename,MaxTextExtent);
3429             quality_image=ReadImage(preview_info,exception);
3430             if (quality_image != (Image *) NULL)
3431               {
3432                 preview_image=DestroyImage(preview_image);
3433                 preview_image=quality_image;
3434               }
3435           }
3436         (void) RelinquishUniqueFileResource(preview_image->filename);
3437         if ((GetBlobSize(preview_image)/1024) >= 1024)
3438           (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%gmb ",
3439             factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
3440             1024.0/1024.0);
3441         else
3442           if (GetBlobSize(preview_image) >= 1024)
3443             (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%gkb ",
3444               factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
3445               1024.0);
3446           else
3447             (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%lub ",
3448               factor,(unsigned long) GetBlobSize(thumbnail));
3449         break;
3450       }
3451     }
3452     thumbnail=DestroyImage(thumbnail);
3453     percentage+=12.5;
3454     radius+=0.5;
3455     sigma+=0.25;
3456     if (preview_image == (Image *) NULL)
3457       break;
3458     (void) DeleteImageProperty(preview_image,"label");
3459     (void) SetImageProperty(preview_image,"label",label);
3460     AppendImageToList(&images,preview_image);
3461     proceed=SetImageProgress(image,PreviewImageTag,i,NumberTiles);
3462     if (proceed == MagickFalse)
3463       break;
3464   }
3465   if (images == (Image *) NULL)
3466     {
3467       preview_info=DestroyImageInfo(preview_info);
3468       return((Image *) NULL);
3469     }
3470   /*
3471     Create the montage.
3472   */
3473   montage_info=CloneMontageInfo(preview_info,(MontageInfo *) NULL);
3474   (void) CopyMagickString(montage_info->filename,image->filename,MaxTextExtent);
3475   montage_info->shadow=MagickTrue;
3476   (void) CloneString(&montage_info->tile,"3x3");
3477   (void) CloneString(&montage_info->geometry,DefaultPreviewGeometry);
3478   (void) CloneString(&montage_info->frame,DefaultTileFrame);
3479   montage_image=MontageImages(images,montage_info,exception);
3480   montage_info=DestroyMontageInfo(montage_info);
3481   images=DestroyImageList(images);
3482   if (montage_image == (Image *) NULL)
3483     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3484   if (montage_image->montage != (char *) NULL)
3485     {
3486       /*
3487         Free image directory.
3488       */
3489       montage_image->montage=(char *) RelinquishMagickMemory(
3490         montage_image->montage);
3491       if (image->directory != (char *) NULL)
3492         montage_image->directory=(char *) RelinquishMagickMemory(
3493           montage_image->directory);
3494     }
3495   preview_info=DestroyImageInfo(preview_info);
3496   return(montage_image);
3497 }
3498 \f
3499 /*
3500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3501 %                                                                             %
3502 %                                                                             %
3503 %                                                                             %
3504 %     R a d i a l B l u r I m a g e                                           %
3505 %                                                                             %
3506 %                                                                             %
3507 %                                                                             %
3508 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3509 %
3510 %  RadialBlurImage() applies a radial blur to the image.
3511 %
3512 %  Andrew Protano contributed this effect.
3513 %
3514 %  The format of the RadialBlurImage method is:
3515 %
3516 %    Image *RadialBlurImage(const Image *image,const double angle,
3517 %      ExceptionInfo *exception)
3518 %    Image *RadialBlurImageChannel(const Image *image,const ChannelType channel,
3519 %      const double angle,ExceptionInfo *exception)
3520 %
3521 %  A description of each parameter follows:
3522 %
3523 %    o image: the image.
3524 %
3525 %    o channel: the channel type.
3526 %
3527 %    o angle: the angle of the radial blur.
3528 %
3529 %    o exception: return any errors or warnings in this structure.
3530 %
3531 */
3532
3533 MagickExport Image *RadialBlurImage(const Image *image,const double angle,
3534   ExceptionInfo *exception)
3535 {
3536   Image
3537     *blur_image;
3538
3539   blur_image=RadialBlurImageChannel(image,DefaultChannels,angle,exception);
3540   return(blur_image);
3541 }
3542
3543 MagickExport Image *RadialBlurImageChannel(const Image *image,
3544   const ChannelType channel,const double angle,ExceptionInfo *exception)
3545 {
3546   Image
3547     *blur_image;
3548
3549   long
3550     progress,
3551     y;
3552
3553   MagickBooleanType
3554     status;
3555
3556   MagickPixelPacket
3557     bias;
3558
3559   MagickRealType
3560     blur_radius,
3561     *cos_theta,
3562     offset,
3563     *sin_theta,
3564     theta;
3565
3566   PointInfo
3567     blur_center;
3568
3569   register long
3570     i;
3571
3572   unsigned long
3573     n;
3574
3575   CacheView
3576     *blur_view,
3577     *image_view;
3578
3579   /*
3580     Allocate blur image.
3581   */
3582   assert(image != (Image *) NULL);
3583   assert(image->signature == MagickSignature);
3584   if (image->debug != MagickFalse)
3585     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3586   assert(exception != (ExceptionInfo *) NULL);
3587   assert(exception->signature == MagickSignature);
3588   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3589   if (blur_image == (Image *) NULL)
3590     return((Image *) NULL);
3591   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
3592     {
3593       InheritException(exception,&blur_image->exception);
3594       blur_image=DestroyImage(blur_image);
3595       return((Image *) NULL);
3596     }
3597   blur_center.x=(double) image->columns/2.0;
3598   blur_center.y=(double) image->rows/2.0;
3599   blur_radius=hypot(blur_center.x,blur_center.y);
3600   n=(unsigned long) fabs(4.0*DegreesToRadians(angle)*sqrt((double) blur_radius)+
3601     2UL);
3602   theta=DegreesToRadians(angle)/(MagickRealType) (n-1);
3603   cos_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
3604     sizeof(*cos_theta));
3605   sin_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
3606     sizeof(*sin_theta));
3607   if ((cos_theta == (MagickRealType *) NULL) ||
3608       (sin_theta == (MagickRealType *) NULL))
3609     {
3610       blur_image=DestroyImage(blur_image);
3611       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3612     }
3613   offset=theta*(MagickRealType) (n-1)/2.0;
3614   for (i=0; i < (long) n; i++)
3615   {
3616     cos_theta[i]=cos((double) (theta*i-offset));
3617     sin_theta[i]=sin((double) (theta*i-offset));
3618   }
3619   /*
3620     Radial blur image.
3621   */
3622   status=MagickTrue;
3623   progress=0;
3624   GetMagickPixelPacket(image,&bias);
3625   image_view=AcquireCacheView(image);
3626   blur_view=AcquireCacheView(blur_image);
3627 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3628   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
3629 #endif
3630   for (y=0; y < (long) blur_image->rows; y++)
3631   {
3632     register const IndexPacket
3633       *restrict indexes;
3634
3635     register IndexPacket
3636       *restrict blur_indexes;
3637
3638     register long
3639       x;
3640
3641     register PixelPacket
3642       *restrict q;
3643
3644     if (status == MagickFalse)
3645       continue;
3646     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
3647       exception);
3648     if (q == (PixelPacket *) NULL)
3649       {
3650         status=MagickFalse;
3651         continue;
3652       }
3653     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
3654     for (x=0; x < (long) blur_image->columns; x++)
3655     {
3656       MagickPixelPacket
3657         qixel;
3658
3659       MagickRealType
3660         normalize,
3661         radius;
3662
3663       PixelPacket
3664         pixel;
3665
3666       PointInfo
3667         center;
3668
3669       register long
3670         i;
3671
3672       unsigned long
3673         step;
3674
3675       center.x=(double) x-blur_center.x;
3676       center.y=(double) y-blur_center.y;
3677       radius=hypot((double) center.x,center.y);
3678       if (radius == 0)
3679         step=1;
3680       else
3681         {
3682           step=(unsigned long) (blur_radius/radius);
3683           if (step == 0)
3684             step=1;
3685           else
3686             if (step >= n)
3687               step=n-1;
3688         }
3689       normalize=0.0;
3690       qixel=bias;
3691       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
3692         {
3693           for (i=0; i < (long) n; i+=step)
3694           {
3695             (void) GetOneCacheViewVirtualPixel(image_view,(long) (blur_center.x+
3696               center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),(long) (
3697               blur_center.y+center.x*sin_theta[i]+center.y*cos_theta[i]+0.5),
3698               &pixel,exception);
3699             qixel.red+=pixel.red;
3700             qixel.green+=pixel.green;
3701             qixel.blue+=pixel.blue;
3702             qixel.opacity+=pixel.opacity;
3703             if (image->colorspace == CMYKColorspace)
3704               {
3705                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3706                 qixel.index+=(*indexes);
3707               }
3708             normalize+=1.0;
3709           }
3710           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
3711             normalize);
3712           if ((channel & RedChannel) != 0)
3713             q->red=RoundToQuantum(normalize*qixel.red);
3714           if ((channel & GreenChannel) != 0)
3715             q->green=RoundToQuantum(normalize*qixel.green);
3716           if ((channel & BlueChannel) != 0)
3717             q->blue=RoundToQuantum(normalize*qixel.blue);
3718           if ((channel & OpacityChannel) != 0)
3719             q->opacity=RoundToQuantum(normalize*qixel.opacity);
3720           if (((channel & IndexChannel) != 0) &&
3721               (image->colorspace == CMYKColorspace))
3722             blur_indexes[x]=(IndexPacket) RoundToQuantum(normalize*qixel.index);
3723         }
3724       else
3725         {
3726           MagickRealType
3727             alpha,
3728             gamma;
3729
3730           alpha=1.0;
3731           gamma=0.0;
3732           for (i=0; i < (long) n; i+=step)
3733           {
3734             (void) GetOneCacheViewVirtualPixel(image_view,(long) (blur_center.x+
3735               center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),(long) (
3736               blur_center.y+center.x*sin_theta[i]+center.y*cos_theta[i]+0.5),
3737               &pixel,exception);
3738             alpha=(MagickRealType) (QuantumScale*(QuantumRange-pixel.opacity));
3739             qixel.red+=alpha*pixel.red;
3740             qixel.green+=alpha*pixel.green;
3741             qixel.blue+=alpha*pixel.blue;
3742             qixel.opacity+=pixel.opacity;
3743             if (image->colorspace == CMYKColorspace)
3744               {
3745                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3746                 qixel.index+=alpha*(*indexes);
3747               }
3748             gamma+=alpha;
3749             normalize+=1.0;
3750           }
3751           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
3752           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
3753             normalize);
3754           if ((channel & RedChannel) != 0)
3755             q->red=RoundToQuantum(gamma*qixel.red);
3756           if ((channel & GreenChannel) != 0)
3757             q->green=RoundToQuantum(gamma*qixel.green);
3758           if ((channel & BlueChannel) != 0)
3759             q->blue=RoundToQuantum(gamma*qixel.blue);
3760           if ((channel & OpacityChannel) != 0)
3761             q->opacity=RoundToQuantum(normalize*qixel.opacity);
3762           if (((channel & IndexChannel) != 0) &&
3763               (image->colorspace == CMYKColorspace))
3764             blur_indexes[x]=(IndexPacket) RoundToQuantum(gamma*qixel.index);
3765         }
3766       q++;
3767     }
3768     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
3769       status=MagickFalse;
3770     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3771       {
3772         MagickBooleanType
3773           proceed;
3774
3775 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3776   #pragma omp critical (MagickCore_RadialBlurImageChannel)
3777 #endif
3778         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
3779         if (proceed == MagickFalse)
3780           status=MagickFalse;
3781       }
3782   }
3783   blur_view=DestroyCacheView(blur_view);
3784   image_view=DestroyCacheView(image_view);
3785   cos_theta=(MagickRealType *) RelinquishMagickMemory(cos_theta);
3786   sin_theta=(MagickRealType *) RelinquishMagickMemory(sin_theta);
3787   if (status == MagickFalse)
3788     blur_image=DestroyImage(blur_image);
3789   return(blur_image);
3790 }
3791 \f
3792 /*
3793 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3794 %                                                                             %
3795 %                                                                             %
3796 %                                                                             %
3797 %     R e d u c e N o i s e I m a g e                                         %
3798 %                                                                             %
3799 %                                                                             %
3800 %                                                                             %
3801 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3802 %
3803 %  ReduceNoiseImage() smooths the contours of an image while still preserving
3804 %  edge information.  The algorithm works by replacing each pixel with its
3805 %  neighbor closest in value.  A neighbor is defined by radius.  Use a radius
3806 %  of 0 and ReduceNoise() selects a suitable radius for you.
3807 %
3808 %  The format of the ReduceNoiseImage method is:
3809 %
3810 %      Image *ReduceNoiseImage(const Image *image,const double radius,
3811 %        ExceptionInfo *exception)
3812 %
3813 %  A description of each parameter follows:
3814 %
3815 %    o image: the image.
3816 %
3817 %    o radius: the radius of the pixel neighborhood.
3818 %
3819 %    o exception: return any errors or warnings in this structure.
3820 %
3821 */
3822
3823 static MagickPixelPacket GetNonpeakMedianPixelList(MedianPixelList *pixel_list)
3824 {
3825   MagickPixelPacket
3826     pixel;
3827
3828   register long
3829     channel;
3830
3831   register MedianSkipList
3832     *list;
3833
3834   unsigned long
3835     center,
3836     color,
3837     count,
3838     previous,
3839     next;
3840
3841   unsigned short
3842     channels[5];
3843
3844   /*
3845     Finds the median value for each of the color.
3846   */
3847   center=pixel_list->center;
3848   for (channel=0; channel < 5; channel++)
3849   {
3850     list=pixel_list->lists+channel;
3851     color=65536UL;
3852     next=list->nodes[color].next[0];
3853     count=0;
3854     do
3855     {
3856       previous=color;
3857       color=next;
3858       next=list->nodes[color].next[0];
3859       count+=list->nodes[color].count;
3860     }
3861     while (count <= center);
3862     if ((previous == 65536UL) && (next != 65536UL))
3863       color=next;
3864     else
3865       if ((previous != 65536UL) && (next == 65536UL))
3866         color=previous;
3867     channels[channel]=(unsigned short) color;
3868   }
3869   GetMagickPixelPacket((const Image *) NULL,&pixel);
3870   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3871   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3872   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3873   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3874   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3875   return(pixel);
3876 }
3877
3878 MagickExport Image *ReduceNoiseImage(const Image *image,const double radius,
3879   ExceptionInfo *exception)
3880 {
3881 #define ReduceNoiseImageTag  "ReduceNoise/Image"
3882
3883   Image
3884     *noise_image;
3885
3886   long
3887     progress,
3888     y;
3889
3890   MagickBooleanType
3891     status;
3892
3893   MedianPixelList
3894     **pixel_list;
3895
3896   unsigned long
3897     width;
3898
3899   CacheView
3900     *image_view,
3901     *noise_view;
3902
3903   /*
3904     Initialize noise image attributes.
3905   */
3906   assert(image != (Image *) NULL);
3907   assert(image->signature == MagickSignature);
3908   if (image->debug != MagickFalse)
3909     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3910   assert(exception != (ExceptionInfo *) NULL);
3911   assert(exception->signature == MagickSignature);
3912   width=GetOptimalKernelWidth2D(radius,0.5);
3913   if ((image->columns < width) || (image->rows < width))
3914     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
3915   noise_image=CloneImage(image,image->columns,image->rows,MagickTrue,
3916     exception);
3917   if (noise_image == (Image *) NULL)
3918     return((Image *) NULL);
3919   if (SetImageStorageClass(noise_image,DirectClass) == MagickFalse)
3920     {
3921       InheritException(exception,&noise_image->exception);
3922       noise_image=DestroyImage(noise_image);
3923       return((Image *) NULL);
3924     }
3925   pixel_list=AcquireMedianPixelListThreadSet(width);
3926   if (pixel_list == (MedianPixelList **) NULL)
3927     {
3928       noise_image=DestroyImage(noise_image);
3929       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3930     }
3931   /*
3932     Reduce noise image.
3933   */
3934   status=MagickTrue;
3935   progress=0;
3936   image_view=AcquireCacheView(image);
3937   noise_view=AcquireCacheView(noise_image);
3938 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3939   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
3940 #endif
3941   for (y=0; y < (long) noise_image->rows; y++)
3942   {
3943     register const IndexPacket
3944       *restrict indexes;
3945
3946     register const PixelPacket
3947       *restrict p;
3948
3949     register IndexPacket
3950       *restrict noise_indexes;
3951
3952     register long
3953       id,
3954       x;
3955
3956     register PixelPacket
3957       *restrict q;
3958
3959     if (status == MagickFalse)
3960       continue;
3961     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
3962       2L),image->columns+width,width,exception);
3963     q=QueueCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
3964       exception);
3965     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3966       {
3967         status=MagickFalse;
3968         continue;
3969       }
3970     indexes=GetCacheViewVirtualIndexQueue(image_view);
3971     noise_indexes=GetCacheViewAuthenticIndexQueue(noise_view);
3972     id=GetOpenMPThreadId();
3973     for (x=0; x < (long) noise_image->columns; x++)
3974     {
3975       MagickPixelPacket
3976         pixel;
3977
3978       register const PixelPacket
3979         *restrict r;
3980
3981       register const IndexPacket
3982         *restrict s;
3983
3984       register long
3985         u,
3986         v;
3987
3988       r=p;
3989       s=indexes+x;
3990       ResetMedianPixelList(pixel_list[id]);
3991       for (v=0; v < (long) width; v++)
3992       {
3993         for (u=0; u < (long) width; u++)
3994           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
3995         r+=image->columns+width;
3996         s+=image->columns+width;
3997       }
3998       pixel=GetNonpeakMedianPixelList(pixel_list[id]);
3999       SetPixelPacket(noise_image,&pixel,q,noise_indexes+x);
4000       p++;
4001       q++;
4002     }
4003     if (SyncCacheViewAuthenticPixels(noise_view,exception) == MagickFalse)
4004       status=MagickFalse;
4005     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4006       {
4007         MagickBooleanType
4008           proceed;
4009
4010 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4011   #pragma omp critical (MagickCore_ReduceNoiseImage)
4012 #endif
4013         proceed=SetImageProgress(image,ReduceNoiseImageTag,progress++,
4014           image->rows);
4015         if (proceed == MagickFalse)
4016           status=MagickFalse;
4017       }
4018   }
4019   noise_view=DestroyCacheView(noise_view);
4020   image_view=DestroyCacheView(image_view);
4021   pixel_list=DestroyMedianPixelListThreadSet(pixel_list);
4022   return(noise_image);
4023 }
4024 \f
4025 /*
4026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4027 %                                                                             %
4028 %                                                                             %
4029 %                                                                             %
4030 %     S e l e c t i v e B l u r I m a g e                                     %
4031 %                                                                             %
4032 %                                                                             %
4033 %                                                                             %
4034 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4035 %
4036 %  SelectiveBlurImage() selectively blur pixels within a contrast threshold.
4037 %  It is similar to the unsharpen mask that sharpens everything with contrast
4038 %  above a certain threshold.
4039 %
4040 %  The format of the SelectiveBlurImage method is:
4041 %
4042 %      Image *SelectiveBlurImage(const Image *image,const double radius,
4043 %        const double sigma,const double threshold,ExceptionInfo *exception)
4044 %      Image *SelectiveBlurImageChannel(const Image *image,
4045 %        const ChannelType channel,const double radius,const double sigma,
4046 %        const double threshold,ExceptionInfo *exception)
4047 %
4048 %  A description of each parameter follows:
4049 %
4050 %    o image: the image.
4051 %
4052 %    o channel: the channel type.
4053 %
4054 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4055 %      pixel.
4056 %
4057 %    o sigma: the standard deviation of the Gaussian, in pixels.
4058 %
4059 %    o threshold: only pixels within this contrast threshold are included
4060 %      in the blur operation.
4061 %
4062 %    o exception: return any errors or warnings in this structure.
4063 %
4064 */
4065
4066 static inline MagickBooleanType SelectiveContrast(const PixelPacket *p,
4067   const PixelPacket *q,const double threshold)
4068 {
4069   if (fabs(PixelIntensity(p)-PixelIntensity(q)) < threshold)
4070     return(MagickTrue);
4071   return(MagickFalse);
4072 }
4073
4074 MagickExport Image *SelectiveBlurImage(const Image *image,const double radius,
4075   const double sigma,const double threshold,ExceptionInfo *exception)
4076 {
4077   Image
4078     *blur_image;
4079
4080   blur_image=SelectiveBlurImageChannel(image,DefaultChannels,radius,sigma,
4081     threshold,exception);
4082   return(blur_image);
4083 }
4084
4085 MagickExport Image *SelectiveBlurImageChannel(const Image *image,
4086   const ChannelType channel,const double radius,const double sigma,
4087   const double threshold,ExceptionInfo *exception)
4088 {
4089 #define SelectiveBlurImageTag  "SelectiveBlur/Image"
4090
4091   CacheView
4092     *blur_view,
4093     *image_view;
4094
4095   double
4096     *kernel;
4097
4098   Image
4099     *blur_image;
4100
4101   long
4102     j,
4103     progress,
4104     u,
4105     v,
4106     y;
4107
4108   MagickBooleanType
4109     status;
4110
4111   MagickPixelPacket
4112     bias;
4113
4114   register long
4115     i;
4116
4117   unsigned long
4118     width;
4119
4120   /*
4121     Initialize blur image attributes.
4122   */
4123   assert(image != (Image *) NULL);
4124   assert(image->signature == MagickSignature);
4125   if (image->debug != MagickFalse)
4126     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4127   assert(exception != (ExceptionInfo *) NULL);
4128   assert(exception->signature == MagickSignature);
4129   width=GetOptimalKernelWidth1D(radius,sigma);
4130   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
4131   if (kernel == (double *) NULL)
4132     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4133   j=(long) width/2;
4134   i=0;
4135   for (v=(-j); v <= j; v++)
4136   {
4137     for (u=(-j); u <= j; u++)
4138       kernel[i++]=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
4139         (2.0*MagickPI*MagickSigma*MagickSigma);
4140   }
4141   if (image->debug != MagickFalse)
4142     {
4143       char
4144         format[MaxTextExtent],
4145         *message;
4146
4147       long
4148         u,
4149         v;
4150
4151       register const double
4152         *k;
4153
4154       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
4155         "  SelectiveBlurImage with %ldx%ld kernel:",width,width);
4156       message=AcquireString("");
4157       k=kernel;
4158       for (v=0; v < (long) width; v++)
4159       {
4160         *message='\0';
4161         (void) FormatMagickString(format,MaxTextExtent,"%ld: ",v);
4162         (void) ConcatenateString(&message,format);
4163         for (u=0; u < (long) width; u++)
4164         {
4165           (void) FormatMagickString(format,MaxTextExtent,"%+f ",*k++);
4166           (void) ConcatenateString(&message,format);
4167         }
4168         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
4169       }
4170       message=DestroyString(message);
4171     }
4172   blur_image=CloneImage(image,0,0,MagickTrue,exception);
4173   if (blur_image == (Image *) NULL)
4174     return((Image *) NULL);
4175   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
4176     {
4177       InheritException(exception,&blur_image->exception);
4178       blur_image=DestroyImage(blur_image);
4179       return((Image *) NULL);
4180     }
4181   /*
4182     Threshold blur image.
4183   */
4184   status=MagickTrue;
4185   progress=0;
4186   GetMagickPixelPacket(image,&bias);
4187   SetMagickPixelPacketBias(image,&bias);
4188   image_view=AcquireCacheView(image);
4189   blur_view=AcquireCacheView(blur_image);
4190 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4191   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4192 #endif
4193   for (y=0; y < (long) image->rows; y++)
4194   {
4195     MagickBooleanType
4196       sync;
4197
4198     MagickRealType
4199       gamma;
4200
4201     register const IndexPacket
4202       *restrict indexes;
4203
4204     register const PixelPacket
4205       *restrict p;
4206
4207     register IndexPacket
4208       *restrict blur_indexes;
4209
4210     register long
4211       x;
4212
4213     register PixelPacket
4214       *restrict q;
4215
4216     if (status == MagickFalse)
4217       continue;
4218     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
4219       2L),image->columns+width,width,exception);
4220     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
4221       exception);
4222     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4223       {
4224         status=MagickFalse;
4225         continue;
4226       }
4227     indexes=GetCacheViewVirtualIndexQueue(image_view);
4228     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
4229     for (x=0; x < (long) image->columns; x++)
4230     {
4231       long
4232         j,
4233         v;
4234
4235       MagickPixelPacket
4236         pixel;
4237
4238       register const double
4239         *restrict k;
4240
4241       register long
4242         u;
4243
4244       pixel=bias;
4245       k=kernel;
4246       gamma=0.0;
4247       j=0;
4248       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
4249         {
4250           for (v=0; v < (long) width; v++)
4251           {
4252             for (u=0; u < (long) width; u++)
4253             {
4254               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4255                 {
4256                   pixel.red+=(*k)*(p+u+j)->red;
4257                   pixel.green+=(*k)*(p+u+j)->green;
4258                   pixel.blue+=(*k)*(p+u+j)->blue;
4259                   gamma+=(*k);
4260                   k++;
4261                 }
4262             }
4263             j+=image->columns+width;
4264           }
4265           if (gamma != 0.0)
4266             {
4267               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4268               if ((channel & RedChannel) != 0)
4269                 q->red=RoundToQuantum(gamma*pixel.red);
4270               if ((channel & GreenChannel) != 0)
4271                 q->green=RoundToQuantum(gamma*pixel.green);
4272               if ((channel & BlueChannel) != 0)
4273                 q->blue=RoundToQuantum(gamma*pixel.blue);
4274             }
4275           if ((channel & OpacityChannel) != 0)
4276             {
4277               gamma=0.0;
4278               j=0;
4279               for (v=0; v < (long) width; v++)
4280               {
4281                 for (u=0; u < (long) width; u++)
4282                 {
4283                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4284                     {
4285                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4286                       gamma+=(*k);
4287                       k++;
4288                     }
4289                 }
4290                 j+=image->columns+width;
4291               }
4292               if (gamma != 0.0)
4293                 {
4294                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4295                     gamma);
4296                   q->opacity=RoundToQuantum(gamma*pixel.opacity);
4297                 }
4298             }
4299           if (((channel & IndexChannel) != 0) &&
4300               (image->colorspace == CMYKColorspace))
4301             {
4302               gamma=0.0;
4303               j=0;
4304               for (v=0; v < (long) width; v++)
4305               {
4306                 for (u=0; u < (long) width; u++)
4307                 {
4308                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4309                     {
4310                       pixel.index+=(*k)*indexes[x+u+j];
4311                       gamma+=(*k);
4312                       k++;
4313                     }
4314                 }
4315                 j+=image->columns+width;
4316               }
4317               if (gamma != 0.0)
4318                 {
4319                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4320                     gamma);
4321                   blur_indexes[x]=RoundToQuantum(gamma*pixel.index);
4322                 }
4323             }
4324         }
4325       else
4326         {
4327           MagickRealType
4328             alpha;
4329
4330           for (v=0; v < (long) width; v++)
4331           {
4332             for (u=0; u < (long) width; u++)
4333             {
4334               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4335                 {
4336                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
4337                     (p+u+j)->opacity));
4338                   pixel.red+=(*k)*alpha*(p+u+j)->red;
4339                   pixel.green+=(*k)*alpha*(p+u+j)->green;
4340                   pixel.blue+=(*k)*alpha*(p+u+j)->blue;
4341                   pixel.opacity+=(*k)*(p+u+j)->opacity;
4342                   gamma+=(*k)*alpha;
4343                   k++;
4344                 }
4345             }
4346             j+=image->columns+width;
4347           }
4348           if (gamma != 0.0)
4349             {
4350               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4351               if ((channel & RedChannel) != 0)
4352                 q->red=RoundToQuantum(gamma*pixel.red);
4353               if ((channel & GreenChannel) != 0)
4354                 q->green=RoundToQuantum(gamma*pixel.green);
4355               if ((channel & BlueChannel) != 0)
4356                 q->blue=RoundToQuantum(gamma*pixel.blue);
4357             }
4358           if ((channel & OpacityChannel) != 0)
4359             {
4360               gamma=0.0;
4361               j=0;
4362               for (v=0; v < (long) width; v++)
4363               {
4364                 for (u=0; u < (long) width; u++)
4365                 {
4366                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4367                     {
4368                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4369                       gamma+=(*k);
4370                       k++;
4371                     }
4372                 }
4373                 j+=image->columns+width;
4374               }
4375               if (gamma != 0.0)
4376                 {
4377                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4378                     gamma);
4379                   q->opacity=RoundToQuantum(pixel.opacity);
4380                 }
4381             }
4382           if (((channel & IndexChannel) != 0) &&
4383               (image->colorspace == CMYKColorspace))
4384             {
4385               gamma=0.0;
4386               j=0;
4387               for (v=0; v < (long) width; v++)
4388               {
4389                 for (u=0; u < (long) width; u++)
4390                 {
4391                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4392                     {
4393                       alpha=(MagickRealType) (QuantumScale*(QuantumRange-
4394                         (p+u+j)->opacity));
4395                       pixel.index+=(*k)*alpha*indexes[x+u+j];
4396                       gamma+=(*k);
4397                       k++;
4398                     }
4399                 }
4400                 j+=image->columns+width;
4401               }
4402               if (gamma != 0.0)
4403                 {
4404                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4405                     gamma);
4406                   blur_indexes[x]=RoundToQuantum(gamma*pixel.index);
4407                 }
4408             }
4409         }
4410       p++;
4411       q++;
4412     }
4413     sync=SyncCacheViewAuthenticPixels(blur_view,exception);
4414     if (sync == MagickFalse)
4415       status=MagickFalse;
4416     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4417       {
4418         MagickBooleanType
4419           proceed;
4420
4421 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4422   #pragma omp critical (MagickCore_SelectiveBlurImageChannel)
4423 #endif
4424         proceed=SetImageProgress(image,SelectiveBlurImageTag,progress++,
4425           image->rows);
4426         if (proceed == MagickFalse)
4427           status=MagickFalse;
4428       }
4429   }
4430   blur_image->type=image->type;
4431   blur_view=DestroyCacheView(blur_view);
4432   image_view=DestroyCacheView(image_view);
4433   kernel=(double *) RelinquishMagickMemory(kernel);
4434   if (status == MagickFalse)
4435     blur_image=DestroyImage(blur_image);
4436   return(blur_image);
4437 }
4438 \f
4439 /*
4440 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4441 %                                                                             %
4442 %                                                                             %
4443 %                                                                             %
4444 %     S h a d e I m a g e                                                     %
4445 %                                                                             %
4446 %                                                                             %
4447 %                                                                             %
4448 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4449 %
4450 %  ShadeImage() shines a distant light on an image to create a
4451 %  three-dimensional effect. You control the positioning of the light with
4452 %  azimuth and elevation; azimuth is measured in degrees off the x axis
4453 %  and elevation is measured in pixels above the Z axis.
4454 %
4455 %  The format of the ShadeImage method is:
4456 %
4457 %      Image *ShadeImage(const Image *image,const MagickBooleanType gray,
4458 %        const double azimuth,const double elevation,ExceptionInfo *exception)
4459 %
4460 %  A description of each parameter follows:
4461 %
4462 %    o image: the image.
4463 %
4464 %    o gray: A value other than zero shades the intensity of each pixel.
4465 %
4466 %    o azimuth, elevation:  Define the light source direction.
4467 %
4468 %    o exception: return any errors or warnings in this structure.
4469 %
4470 */
4471 MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
4472   const double azimuth,const double elevation,ExceptionInfo *exception)
4473 {
4474 #define ShadeImageTag  "Shade/Image"
4475
4476   Image
4477     *shade_image;
4478
4479   long
4480     progress,
4481     y;
4482
4483   MagickBooleanType
4484     status;
4485
4486   PrimaryInfo
4487     light;
4488
4489   CacheView
4490     *image_view,
4491     *shade_view;
4492
4493   /*
4494     Initialize shaded image attributes.
4495   */
4496   assert(image != (const Image *) NULL);
4497   assert(image->signature == MagickSignature);
4498   if (image->debug != MagickFalse)
4499     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4500   assert(exception != (ExceptionInfo *) NULL);
4501   assert(exception->signature == MagickSignature);
4502   shade_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
4503   if (shade_image == (Image *) NULL)
4504     return((Image *) NULL);
4505   if (SetImageStorageClass(shade_image,DirectClass) == MagickFalse)
4506     {
4507       InheritException(exception,&shade_image->exception);
4508       shade_image=DestroyImage(shade_image);
4509       return((Image *) NULL);
4510     }
4511   /*
4512     Compute the light vector.
4513   */
4514   light.x=(double) QuantumRange*cos(DegreesToRadians(azimuth))*
4515     cos(DegreesToRadians(elevation));
4516   light.y=(double) QuantumRange*sin(DegreesToRadians(azimuth))*
4517     cos(DegreesToRadians(elevation));
4518   light.z=(double) QuantumRange*sin(DegreesToRadians(elevation));
4519   /*
4520     Shade image.
4521   */
4522   status=MagickTrue;
4523   progress=0;
4524   image_view=AcquireCacheView(image);
4525   shade_view=AcquireCacheView(shade_image);
4526 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4527   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4528 #endif
4529   for (y=0; y < (long) image->rows; y++)
4530   {
4531     MagickRealType
4532       distance,
4533       normal_distance,
4534       shade;
4535
4536     PrimaryInfo
4537       normal;
4538
4539     register const PixelPacket
4540       *restrict p,
4541       *restrict s0,
4542       *restrict s1,
4543       *restrict s2;
4544
4545     register long
4546       x;
4547
4548     register PixelPacket
4549       *restrict q;
4550
4551     if (status == MagickFalse)
4552       continue;
4553     p=GetCacheViewVirtualPixels(image_view,-1,y-1,image->columns+2,3,exception);
4554     q=QueueCacheViewAuthenticPixels(shade_view,0,y,shade_image->columns,1,
4555       exception);
4556     if ((p == (PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4557       {
4558         status=MagickFalse;
4559         continue;
4560       }
4561     /*
4562       Shade this row of pixels.
4563     */
4564     normal.z=2.0*(double) QuantumRange;  /* constant Z of surface normal */
4565     s0=p+1;
4566     s1=s0+image->columns+2;
4567     s2=s1+image->columns+2;
4568     for (x=0; x < (long) image->columns; x++)
4569     {
4570       /*
4571         Determine the surface normal and compute shading.
4572       */
4573       normal.x=(double) (PixelIntensity(s0-1)+PixelIntensity(s1-1)+
4574         PixelIntensity(s2-1)-PixelIntensity(s0+1)-PixelIntensity(s1+1)-
4575         PixelIntensity(s2+1));
4576       normal.y=(double) (PixelIntensity(s2-1)+PixelIntensity(s2)+
4577         PixelIntensity(s2+1)-PixelIntensity(s0-1)-PixelIntensity(s0)-
4578         PixelIntensity(s0+1));
4579       if ((normal.x == 0.0) && (normal.y == 0.0))
4580         shade=light.z;
4581       else
4582         {
4583           shade=0.0;
4584           distance=normal.x*light.x+normal.y*light.y+normal.z*light.z;
4585           if (distance > MagickEpsilon)
4586             {
4587               normal_distance=
4588                 normal.x*normal.x+normal.y*normal.y+normal.z*normal.z;
4589               if (normal_distance > (MagickEpsilon*MagickEpsilon))
4590                 shade=distance/sqrt((double) normal_distance);
4591             }
4592         }
4593       if (gray != MagickFalse)
4594         {
4595           q->red=(Quantum) shade;
4596           q->green=(Quantum) shade;
4597           q->blue=(Quantum) shade;
4598         }
4599       else
4600         {
4601           q->red=RoundToQuantum(QuantumScale*shade*s1->red);
4602           q->green=RoundToQuantum(QuantumScale*shade*s1->green);
4603           q->blue=RoundToQuantum(QuantumScale*shade*s1->blue);
4604         }
4605       q->opacity=s1->opacity;
4606       s0++;
4607       s1++;
4608       s2++;
4609       q++;
4610     }
4611     if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
4612       status=MagickFalse;
4613     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4614       {
4615         MagickBooleanType
4616           proceed;
4617
4618 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4619   #pragma omp critical (MagickCore_ShadeImage)
4620 #endif
4621         proceed=SetImageProgress(image,ShadeImageTag,progress++,image->rows);
4622         if (proceed == MagickFalse)
4623           status=MagickFalse;
4624       }
4625   }
4626   shade_view=DestroyCacheView(shade_view);
4627   image_view=DestroyCacheView(image_view);
4628   if (status == MagickFalse)
4629     shade_image=DestroyImage(shade_image);
4630   return(shade_image);
4631 }
4632 \f
4633 /*
4634 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4635 %                                                                             %
4636 %                                                                             %
4637 %                                                                             %
4638 %     S h a r p e n I m a g e                                                 %
4639 %                                                                             %
4640 %                                                                             %
4641 %                                                                             %
4642 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4643 %
4644 %  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
4645 %  operator of the given radius and standard deviation (sigma).  For
4646 %  reasonable results, radius should be larger than sigma.  Use a radius of 0
4647 %  and SharpenImage() selects a suitable radius for you.
4648 %
4649 %  Using a separable kernel would be faster, but the negative weights cancel
4650 %  out on the corners of the kernel producing often undesirable ringing in the
4651 %  filtered result; this can be avoided by using a 2D gaussian shaped image
4652 %  sharpening kernel instead.
4653 %
4654 %  The format of the SharpenImage method is:
4655 %
4656 %    Image *SharpenImage(const Image *image,const double radius,
4657 %      const double sigma,ExceptionInfo *exception)
4658 %    Image *SharpenImageChannel(const Image *image,const ChannelType channel,
4659 %      const double radius,const double sigma,ExceptionInfo *exception)
4660 %
4661 %  A description of each parameter follows:
4662 %
4663 %    o image: the image.
4664 %
4665 %    o channel: the channel type.
4666 %
4667 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4668 %      pixel.
4669 %
4670 %    o sigma: the standard deviation of the Laplacian, in pixels.
4671 %
4672 %    o exception: return any errors or warnings in this structure.
4673 %
4674 */
4675
4676 MagickExport Image *SharpenImage(const Image *image,const double radius,
4677   const double sigma,ExceptionInfo *exception)
4678 {
4679   Image
4680     *sharp_image;
4681
4682   sharp_image=SharpenImageChannel(image,DefaultChannels,radius,sigma,exception);
4683   return(sharp_image);
4684 }
4685
4686 MagickExport Image *SharpenImageChannel(const Image *image,
4687   const ChannelType channel,const double radius,const double sigma,
4688   ExceptionInfo *exception)
4689 {
4690   double
4691     *kernel,
4692     normalize;
4693
4694   Image
4695     *sharp_image;
4696
4697   long
4698     j,
4699     u,
4700     v;
4701
4702   register long
4703     i;
4704
4705   unsigned long
4706     width;
4707
4708   assert(image != (const Image *) NULL);
4709   assert(image->signature == MagickSignature);
4710   if (image->debug != MagickFalse)
4711     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4712   assert(exception != (ExceptionInfo *) NULL);
4713   assert(exception->signature == MagickSignature);
4714   width=GetOptimalKernelWidth2D(radius,sigma);
4715   kernel=(double *) AcquireQuantumMemory((size_t) width*width,sizeof(*kernel));
4716   if (kernel == (double *) NULL)
4717     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4718   normalize=0.0;
4719   j=(long) width/2;
4720   i=0;
4721   for (v=(-j); v <= j; v++)
4722   {
4723     for (u=(-j); u <= j; u++)
4724     {
4725       kernel[i]=(-exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
4726         (2.0*MagickPI*MagickSigma*MagickSigma));
4727       normalize+=kernel[i];
4728       i++;
4729     }
4730   }
4731   kernel[i/2]=(double) ((-2.0)*normalize);
4732   sharp_image=ConvolveImageChannel(image,channel,width,kernel,exception);
4733   kernel=(double *) RelinquishMagickMemory(kernel);
4734   return(sharp_image);
4735 }
4736 \f
4737 /*
4738 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4739 %                                                                             %
4740 %                                                                             %
4741 %                                                                             %
4742 %     S p r e a d I m a g e                                                   %
4743 %                                                                             %
4744 %                                                                             %
4745 %                                                                             %
4746 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4747 %
4748 %  SpreadImage() is a special effects method that randomly displaces each
4749 %  pixel in a block defined by the radius parameter.
4750 %
4751 %  The format of the SpreadImage method is:
4752 %
4753 %      Image *SpreadImage(const Image *image,const double radius,
4754 %        ExceptionInfo *exception)
4755 %
4756 %  A description of each parameter follows:
4757 %
4758 %    o image: the image.
4759 %
4760 %    o radius:  Choose a random pixel in a neighborhood of this extent.
4761 %
4762 %    o exception: return any errors or warnings in this structure.
4763 %
4764 */
4765 MagickExport Image *SpreadImage(const Image *image,const double radius,
4766   ExceptionInfo *exception)
4767 {
4768 #define SpreadImageTag  "Spread/Image"
4769
4770   Image
4771     *spread_image;
4772
4773   long
4774     progress,
4775     y;
4776
4777   MagickBooleanType
4778     status;
4779
4780   MagickPixelPacket
4781     bias;
4782
4783   RandomInfo
4784     **random_info;
4785
4786   ResampleFilter
4787     **resample_filter;
4788
4789   unsigned long
4790     width;
4791
4792   CacheView
4793     *image_view;
4794
4795   /*
4796     Initialize spread image attributes.
4797   */
4798   assert(image != (Image *) NULL);
4799   assert(image->signature == MagickSignature);
4800   if (image->debug != MagickFalse)
4801     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4802   assert(exception != (ExceptionInfo *) NULL);
4803   assert(exception->signature == MagickSignature);
4804   spread_image=CloneImage(image,image->columns,image->rows,MagickTrue,
4805     exception);
4806   if (spread_image == (Image *) NULL)
4807     return((Image *) NULL);
4808   if (SetImageStorageClass(spread_image,DirectClass) == MagickFalse)
4809     {
4810       InheritException(exception,&spread_image->exception);
4811       spread_image=DestroyImage(spread_image);
4812       return((Image *) NULL);
4813     }
4814   /*
4815     Spread image.
4816   */
4817   status=MagickTrue;
4818   progress=0;
4819   GetMagickPixelPacket(spread_image,&bias);
4820   width=GetOptimalKernelWidth1D(radius,0.5);
4821   resample_filter=AcquireResampleFilterThreadSet(image,MagickTrue,exception);
4822   random_info=AcquireRandomInfoThreadSet();
4823   image_view=AcquireCacheView(spread_image);
4824 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4825   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4826 #endif
4827   for (y=0; y < (long) spread_image->rows; y++)
4828   {
4829     MagickPixelPacket
4830       pixel;
4831
4832     register IndexPacket
4833       *restrict indexes;
4834
4835     register long
4836       id,
4837       x;
4838
4839     register PixelPacket
4840       *restrict q;
4841
4842     if (status == MagickFalse)
4843       continue;
4844     q=QueueCacheViewAuthenticPixels(image_view,0,y,spread_image->columns,1,
4845       exception);
4846     if (q == (PixelPacket *) NULL)
4847       {
4848         status=MagickFalse;
4849         continue;
4850       }
4851     indexes=GetCacheViewAuthenticIndexQueue(image_view);
4852     pixel=bias;
4853     id=GetOpenMPThreadId();
4854     for (x=0; x < (long) spread_image->columns; x++)
4855     {
4856       (void) ResamplePixelColor(resample_filter[id],(double) x+width*
4857         (GetPseudoRandomValue(random_info[id])-0.5),(double) y+width*
4858         (GetPseudoRandomValue(random_info[id])-0.5),&pixel);
4859       SetPixelPacket(spread_image,&pixel,q,indexes+x);
4860       q++;
4861     }
4862     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
4863       status=MagickFalse;
4864     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4865       {
4866         MagickBooleanType
4867           proceed;
4868
4869 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4870   #pragma omp critical (MagickCore_SpreadImage)
4871 #endif
4872         proceed=SetImageProgress(image,SpreadImageTag,progress++,image->rows);
4873         if (proceed == MagickFalse)
4874           status=MagickFalse;
4875       }
4876   }
4877   image_view=DestroyCacheView(image_view);
4878   random_info=DestroyRandomInfoThreadSet(random_info);
4879   resample_filter=DestroyResampleFilterThreadSet(resample_filter);
4880   return(spread_image);
4881 }
4882 \f
4883 /*
4884 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4885 %                                                                             %
4886 %                                                                             %
4887 %                                                                             %
4888 %     U n s h a r p M a s k I m a g e                                         %
4889 %                                                                             %
4890 %                                                                             %
4891 %                                                                             %
4892 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4893 %
4894 %  UnsharpMaskImage() sharpens one or more image channels.  We convolve the
4895 %  image with a Gaussian operator of the given radius and standard deviation
4896 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
4897 %  radius of 0 and UnsharpMaskImage() selects a suitable radius for you.
4898 %
4899 %  The format of the UnsharpMaskImage method is:
4900 %
4901 %    Image *UnsharpMaskImage(const Image *image,const double radius,
4902 %      const double sigma,const double amount,const double threshold,
4903 %      ExceptionInfo *exception)
4904 %    Image *UnsharpMaskImageChannel(const Image *image,
4905 %      const ChannelType channel,const double radius,const double sigma,
4906 %      const double amount,const double threshold,ExceptionInfo *exception)
4907 %
4908 %  A description of each parameter follows:
4909 %
4910 %    o image: the image.
4911 %
4912 %    o channel: the channel type.
4913 %
4914 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4915 %      pixel.
4916 %
4917 %    o sigma: the standard deviation of the Gaussian, in pixels.
4918 %
4919 %    o amount: the percentage of the difference between the original and the
4920 %      blur image that is added back into the original.
4921 %
4922 %    o threshold: the threshold in pixels needed to apply the diffence amount.
4923 %
4924 %    o exception: return any errors or warnings in this structure.
4925 %
4926 */
4927
4928 MagickExport Image *UnsharpMaskImage(const Image *image,const double radius,
4929   const double sigma,const double amount,const double threshold,
4930   ExceptionInfo *exception)
4931 {
4932   Image
4933     *sharp_image;
4934
4935   sharp_image=UnsharpMaskImageChannel(image,DefaultChannels,radius,sigma,amount,
4936     threshold,exception);
4937   return(sharp_image);
4938 }
4939
4940 MagickExport Image *UnsharpMaskImageChannel(const Image *image,
4941   const ChannelType channel,const double radius,const double sigma,
4942   const double amount,const double threshold,ExceptionInfo *exception)
4943 {
4944 #define SharpenImageTag  "Sharpen/Image"
4945
4946   Image
4947     *unsharp_image;
4948
4949   long
4950     progress,
4951     y;
4952
4953   MagickBooleanType
4954     status;
4955
4956   MagickPixelPacket
4957     bias;
4958
4959   MagickRealType
4960     quantum_threshold;
4961
4962   CacheView
4963     *image_view,
4964     *unsharp_view;
4965
4966   assert(image != (const Image *) NULL);
4967   assert(image->signature == MagickSignature);
4968   if (image->debug != MagickFalse)
4969     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4970   assert(exception != (ExceptionInfo *) NULL);
4971   unsharp_image=BlurImageChannel(image,channel,radius,sigma,exception);
4972   if (unsharp_image == (Image *) NULL)
4973     return((Image *) NULL);
4974   quantum_threshold=(MagickRealType) QuantumRange*threshold;
4975   /*
4976     Unsharp-mask image.
4977   */
4978   status=MagickTrue;
4979   progress=0;
4980   GetMagickPixelPacket(image,&bias);
4981   image_view=AcquireCacheView(image);
4982   unsharp_view=AcquireCacheView(unsharp_image);
4983 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4984   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4985 #endif
4986   for (y=0; y < (long) image->rows; y++)
4987   {
4988     MagickPixelPacket
4989       pixel;
4990
4991     register const IndexPacket
4992       *restrict indexes;
4993
4994     register const PixelPacket
4995       *restrict p;
4996
4997     register IndexPacket
4998       *restrict unsharp_indexes;
4999
5000     register long
5001       x;
5002
5003     register PixelPacket
5004       *restrict q;
5005
5006     if (status == MagickFalse)
5007       continue;
5008     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
5009     q=GetCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
5010       exception);
5011     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
5012       {
5013         status=MagickFalse;
5014         continue;
5015       }
5016     indexes=GetCacheViewVirtualIndexQueue(image_view);
5017     unsharp_indexes=GetCacheViewAuthenticIndexQueue(unsharp_view);
5018     pixel=bias;
5019     for (x=0; x < (long) image->columns; x++)
5020     {
5021       if ((channel & RedChannel) != 0)
5022         {
5023           pixel.red=p->red-(MagickRealType) q->red;
5024           if (fabs(2.0*pixel.red) < quantum_threshold)
5025             pixel.red=(MagickRealType) p->red;
5026           else
5027             pixel.red=(MagickRealType) p->red+(pixel.red*amount);
5028           q->red=RoundToQuantum(pixel.red);
5029         }
5030       if ((channel & GreenChannel) != 0)
5031         {
5032           pixel.green=p->green-(MagickRealType) q->green;
5033           if (fabs(2.0*pixel.green) < quantum_threshold)
5034             pixel.green=(MagickRealType) p->green;
5035           else
5036             pixel.green=(MagickRealType) p->green+(pixel.green*amount);
5037           q->green=RoundToQuantum(pixel.green);
5038         }
5039       if ((channel & BlueChannel) != 0)
5040         {
5041           pixel.blue=p->blue-(MagickRealType) q->blue;
5042           if (fabs(2.0*pixel.blue) < quantum_threshold)
5043             pixel.blue=(MagickRealType) p->blue;
5044           else
5045             pixel.blue=(MagickRealType) p->blue+(pixel.blue*amount);
5046           q->blue=RoundToQuantum(pixel.blue);
5047         }
5048       if ((channel & OpacityChannel) != 0)
5049         {
5050           pixel.opacity=p->opacity-(MagickRealType) q->opacity;
5051           if (fabs(2.0*pixel.opacity) < quantum_threshold)
5052             pixel.opacity=(MagickRealType) p->opacity;
5053           else
5054             pixel.opacity=p->opacity+(pixel.opacity*amount);
5055           q->opacity=RoundToQuantum(pixel.opacity);
5056         }
5057       if (((channel & IndexChannel) != 0) &&
5058           (image->colorspace == CMYKColorspace))
5059         {
5060           pixel.index=unsharp_indexes[x]-(MagickRealType) indexes[x];
5061           if (fabs(2.0*pixel.index) < quantum_threshold)
5062             pixel.index=(MagickRealType) unsharp_indexes[x];
5063           else
5064             pixel.index=(MagickRealType) unsharp_indexes[x]+(pixel.index*
5065               amount);
5066           unsharp_indexes[x]=RoundToQuantum(pixel.index);
5067         }
5068       p++;
5069       q++;
5070     }
5071     if (SyncCacheViewAuthenticPixels(unsharp_view,exception) == MagickFalse)
5072       status=MagickFalse;
5073     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5074       {
5075         MagickBooleanType
5076           proceed;
5077
5078 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5079   #pragma omp critical (MagickCore_UnsharpMaskImageChannel)
5080 #endif
5081         proceed=SetImageProgress(image,SharpenImageTag,progress++,image->rows);
5082         if (proceed == MagickFalse)
5083           status=MagickFalse;
5084       }
5085   }
5086   unsharp_image->type=image->type;
5087   unsharp_view=DestroyCacheView(unsharp_view);
5088   image_view=DestroyCacheView(image_view);
5089   if (status == MagickFalse)
5090     unsharp_image=DestroyImage(unsharp_image);
5091   return(unsharp_image);
5092 }