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