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