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   Image
1438     *despeckle_image;
1439
1440   long
1441     channel;
1442
1443   MagickBooleanType
1444     status;
1445
1446   Quantum
1447     **buffers,
1448     **pixels;
1449
1450   size_t
1451     length;
1452
1453   static const int
1454     X[4]= {0, 1, 1,-1},
1455     Y[4]= {1, 0, 1, 1};
1456
1457   CacheView
1458     *despeckle_view,
1459     *image_view;
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       x;
1513
1514     register Quantum
1515       *buffer,
1516       *pixel;
1517
1518     if (status == MagickFalse)
1519       continue;
1520     pixel=pixels[GetOpenMPThreadId()];
1521     (void) ResetMagickMemory(pixel,0,length*sizeof(*pixel));
1522     j=(long) image->columns+2;
1523     for (y=0; y < (long) image->rows; y++)
1524     {
1525       register const PixelPacket
1526         *__restrict p;
1527
1528       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1529       if (p == (const PixelPacket *) NULL)
1530         break;
1531       j++;
1532       for (x=0; x < (long) image->columns; x++)
1533       {
1534         switch (channel)
1535         {
1536           case 0: pixel[j]=p->red; break;
1537           case 1: pixel[j]=p->green; break;
1538           case 2: pixel[j]=p->blue; break;
1539           case 3: pixel[j]=p->opacity; break;
1540           default: break;
1541         }
1542         p++;
1543         j++;
1544       }
1545       j++;
1546     }
1547     buffer=buffers[GetOpenMPThreadId()];
1548     (void) ResetMagickMemory(buffer,0,length*sizeof(*buffer));
1549     for (i=0; i < 4; i++)
1550     {
1551       Hull(X[i],Y[i],image->columns,image->rows,pixel,buffer,1);
1552       Hull(-X[i],-Y[i],image->columns,image->rows,pixel,buffer,1);
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     }
1556     j=(long) image->columns+2;
1557     for (y=0; y < (long) image->rows; y++)
1558     {
1559       MagickBooleanType
1560         sync;
1561
1562       register PixelPacket
1563         *__restrict q;
1564
1565       q=GetCacheViewAuthenticPixels(despeckle_view,0,y,despeckle_image->columns,
1566         1,exception);
1567       if (q == (PixelPacket *) NULL)
1568         break;
1569       j++;
1570       for (x=0; x < (long) image->columns; x++)
1571       {
1572         switch (channel)
1573         {
1574           case 0: q->red=pixel[j]; break;
1575           case 1: q->green=pixel[j]; break;
1576           case 2: q->blue=pixel[j]; break;
1577           case 3: q->opacity=pixel[j]; break;
1578           default: break;
1579         }
1580         q++;
1581         j++;
1582       }
1583       sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
1584       if (sync == MagickFalse)
1585         {
1586           status=MagickFalse;
1587           break;
1588         }
1589       j++;
1590     }
1591     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1592       {
1593         MagickBooleanType
1594           proceed;
1595
1596 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1597   #pragma omp critical (MagickCore_DespeckleImage)
1598 #endif
1599         proceed=SetImageProgress(image,DespeckleImageTag,channel,3);
1600         if (proceed == MagickFalse)
1601           status=MagickFalse;
1602       }
1603   }
1604   despeckle_view=DestroyCacheView(despeckle_view);
1605   image_view=DestroyCacheView(image_view);
1606   buffers=DestroyPixelThreadSet(buffers);
1607   pixels=DestroyPixelThreadSet(pixels);
1608   despeckle_image->type=image->type;
1609   if (status == MagickFalse)
1610     despeckle_image=DestroyImage(despeckle_image);
1611   return(despeckle_image);
1612 }
1613 \f
1614 /*
1615 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1616 %                                                                             %
1617 %                                                                             %
1618 %                                                                             %
1619 %     E d g e I m a g e                                                       %
1620 %                                                                             %
1621 %                                                                             %
1622 %                                                                             %
1623 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1624 %
1625 %  EdgeImage() finds edges in an image.  Radius defines the radius of the
1626 %  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
1627 %  radius for you.
1628 %
1629 %  The format of the EdgeImage method is:
1630 %
1631 %      Image *EdgeImage(const Image *image,const double radius,
1632 %        ExceptionInfo *exception)
1633 %
1634 %  A description of each parameter follows:
1635 %
1636 %    o image: the image.
1637 %
1638 %    o radius: the radius of the pixel neighborhood.
1639 %
1640 %    o exception: return any errors or warnings in this structure.
1641 %
1642 */
1643 MagickExport Image *EdgeImage(const Image *image,const double radius,
1644   ExceptionInfo *exception)
1645 {
1646   Image
1647     *edge_image;
1648
1649   double
1650     *kernel;
1651
1652   register long
1653     i;
1654
1655   unsigned long
1656     width;
1657
1658   assert(image != (const Image *) NULL);
1659   assert(image->signature == MagickSignature);
1660   if (image->debug != MagickFalse)
1661     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1662   assert(exception != (ExceptionInfo *) NULL);
1663   assert(exception->signature == MagickSignature);
1664   width=GetOptimalKernelWidth1D(radius,0.5);
1665   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
1666   if (kernel == (double *) NULL)
1667     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1668   for (i=0; i < (long) (width*width); i++)
1669     kernel[i]=(-1.0);
1670   kernel[i/2]=(double) (width*width-1.0);
1671   edge_image=ConvolveImage(image,width,kernel,exception);
1672   kernel=(double *) RelinquishMagickMemory(kernel);
1673   return(edge_image);
1674 }
1675 \f
1676 /*
1677 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1678 %                                                                             %
1679 %                                                                             %
1680 %                                                                             %
1681 %     E m b o s s I m a g e                                                   %
1682 %                                                                             %
1683 %                                                                             %
1684 %                                                                             %
1685 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1686 %
1687 %  EmbossImage() returns a grayscale image with a three-dimensional effect.
1688 %  We convolve the image with a Gaussian operator of the given radius and
1689 %  standard deviation (sigma).  For reasonable results, radius should be
1690 %  larger than sigma.  Use a radius of 0 and Emboss() selects a suitable
1691 %  radius for you.
1692 %
1693 %  The format of the EmbossImage method is:
1694 %
1695 %      Image *EmbossImage(const Image *image,const double radius,
1696 %        const double sigma,ExceptionInfo *exception)
1697 %
1698 %  A description of each parameter follows:
1699 %
1700 %    o image: the image.
1701 %
1702 %    o radius: the radius of the pixel neighborhood.
1703 %
1704 %    o sigma: the standard deviation of the Gaussian, in pixels.
1705 %
1706 %    o exception: return any errors or warnings in this structure.
1707 %
1708 */
1709 MagickExport Image *EmbossImage(const Image *image,const double radius,
1710   const double sigma,ExceptionInfo *exception)
1711 {
1712   double
1713     *kernel;
1714
1715   Image
1716     *emboss_image;
1717
1718   long
1719     j;
1720
1721   MagickRealType
1722     alpha;
1723
1724   register long
1725     i,
1726     u,
1727     v;
1728
1729   unsigned long
1730     width;
1731
1732   assert(image != (Image *) NULL);
1733   assert(image->signature == MagickSignature);
1734   if (image->debug != MagickFalse)
1735     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1736   assert(exception != (ExceptionInfo *) NULL);
1737   assert(exception->signature == MagickSignature);
1738   width=GetOptimalKernelWidth2D(radius,sigma);
1739   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
1740   if (kernel == (double *) NULL)
1741     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1742   i=0;
1743   j=(long) width/2;
1744   for (v=(-((long) width/2)); v <= (long) (width/2); v++)
1745   {
1746     for (u=(-((long) width/2)); u <= (long) (width/2); u++)
1747     {
1748       alpha=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma));
1749       kernel[i]=(double) (((u < 0) || (v < 0) ? -8.0 : 8.0)*alpha/
1750         (2.0*MagickPI*MagickSigma*MagickSigma));
1751       if (u != j)
1752         kernel[i]=0.0;
1753       i++;
1754     }
1755     j--;
1756   }
1757   emboss_image=ConvolveImage(image,width,kernel,exception);
1758   if (emboss_image != (Image *) NULL)
1759     (void) EqualizeImage(emboss_image);
1760   kernel=(double *) RelinquishMagickMemory(kernel);
1761   return(emboss_image);
1762 }
1763 \f
1764 /*
1765 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1766 %                                                                             %
1767 %                                                                             %
1768 %                                                                             %
1769 %     G a u s s i a n B l u r I m a g e                                       %
1770 %                                                                             %
1771 %                                                                             %
1772 %                                                                             %
1773 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1774 %
1775 %  GaussianBlurImage() blurs an image.  We convolve the image with a
1776 %  Gaussian operator of the given radius and standard deviation (sigma).
1777 %  For reasonable results, the radius should be larger than sigma.  Use a
1778 %  radius of 0 and GaussianBlurImage() selects a suitable radius for you
1779 %
1780 %  The format of the GaussianBlurImage method is:
1781 %
1782 %      Image *GaussianBlurImage(const Image *image,onst double radius,
1783 %        const double sigma,ExceptionInfo *exception)
1784 %      Image *GaussianBlurImageChannel(const Image *image,
1785 %        const ChannelType channel,const double radius,const double sigma,
1786 %        ExceptionInfo *exception)
1787 %
1788 %  A description of each parameter follows:
1789 %
1790 %    o image: the image.
1791 %
1792 %    o channel: the channel type.
1793 %
1794 %    o radius: the radius of the Gaussian, in pixels, not counting the center
1795 %      pixel.
1796 %
1797 %    o sigma: the standard deviation of the Gaussian, in pixels.
1798 %
1799 %    o exception: return any errors or warnings in this structure.
1800 %
1801 */
1802
1803 MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
1804   const double sigma,ExceptionInfo *exception)
1805 {
1806   Image
1807     *blur_image;
1808
1809   blur_image=GaussianBlurImageChannel(image,DefaultChannels,radius,sigma,
1810     exception);
1811   return(blur_image);
1812 }
1813
1814 MagickExport Image *GaussianBlurImageChannel(const Image *image,
1815   const ChannelType channel,const double radius,const double sigma,
1816   ExceptionInfo *exception)
1817 {
1818   double
1819     *kernel;
1820
1821   Image
1822     *blur_image;
1823
1824   MagickRealType
1825     alpha;
1826
1827   register long
1828     i,
1829     u,
1830     v;
1831
1832   unsigned long
1833     width;
1834
1835   assert(image != (const Image *) NULL);
1836   assert(image->signature == MagickSignature);
1837   if (image->debug != MagickFalse)
1838     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1839   assert(exception != (ExceptionInfo *) NULL);
1840   assert(exception->signature == MagickSignature);
1841   width=GetOptimalKernelWidth2D(radius,sigma);
1842   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
1843   if (kernel == (double *) NULL)
1844     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1845   i=0;
1846   for (v=(-((long) width/2)); v <= (long) (width/2); v++)
1847   {
1848     for (u=(-((long) width/2)); u <= (long) (width/2); u++)
1849     {
1850       alpha=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma));
1851       kernel[i]=(double) (alpha/(2.0*MagickPI*MagickSigma*MagickSigma));
1852       i++;
1853     }
1854   }
1855   blur_image=ConvolveImageChannel(image,channel,width,kernel,exception);
1856   kernel=(double *) RelinquishMagickMemory(kernel);
1857   return(blur_image);
1858 }
1859 \f
1860 /*
1861 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1862 %                                                                             %
1863 %                                                                             %
1864 %                                                                             %
1865 %     M e d i a n F i l t e r I m a g e                                       %
1866 %                                                                             %
1867 %                                                                             %
1868 %                                                                             %
1869 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1870 %
1871 %  MedianFilterImage() applies a digital filter that improves the quality
1872 %  of a noisy image.  Each pixel is replaced by the median in a set of
1873 %  neighboring pixels as defined by radius.
1874 %
1875 %  The algorithm was contributed by Mike Edmonds and implements an insertion
1876 %  sort for selecting median color-channel values.  For more on this algorithm
1877 %  see "Skip Lists: A probabilistic Alternative to Balanced Trees" by William
1878 %  Pugh in the June 1990 of Communications of the ACM.
1879 %
1880 %  The format of the MedianFilterImage method is:
1881 %
1882 %      Image *MedianFilterImage(const Image *image,const double radius,
1883 %        ExceptionInfo *exception)
1884 %
1885 %  A description of each parameter follows:
1886 %
1887 %    o image: the image.
1888 %
1889 %    o radius: the radius of the pixel neighborhood.
1890 %
1891 %    o exception: return any errors or warnings in this structure.
1892 %
1893 */
1894
1895 #define MedianListChannels  5
1896
1897 typedef struct _MedianListNode
1898 {
1899   unsigned long
1900     next[9],
1901     count,
1902     signature;
1903 } MedianListNode;
1904
1905 typedef struct _MedianSkipList
1906 {
1907   long
1908     level;
1909
1910   MedianListNode
1911     *nodes;
1912 } MedianSkipList;
1913
1914 typedef struct _MedianPixelList
1915 {
1916   unsigned long
1917     center,
1918     seed,
1919     signature;
1920
1921   MedianSkipList
1922     lists[MedianListChannels];
1923 } MedianPixelList;
1924
1925 static MedianPixelList *DestroyMedianPixelList(MedianPixelList *pixel_list)
1926 {
1927   register long
1928     i;
1929
1930   if (pixel_list == (MedianPixelList *) NULL)
1931     return((MedianPixelList *) NULL);
1932   for (i=0; i < MedianListChannels; i++)
1933     if (pixel_list->lists[i].nodes != (MedianListNode *) NULL)
1934       pixel_list->lists[i].nodes=(MedianListNode *) RelinquishMagickMemory(
1935         pixel_list->lists[i].nodes);
1936   pixel_list=(MedianPixelList *) RelinquishAlignedMemory(pixel_list);
1937   return(pixel_list);
1938 }
1939
1940 static MedianPixelList **DestroyMedianPixelListThreadSet(
1941   MedianPixelList **pixel_list)
1942 {
1943   register long
1944     i;
1945
1946   assert(pixel_list != (MedianPixelList **) NULL);
1947   for (i=0; i < (long) GetOpenMPMaximumThreads(); i++)
1948     if (pixel_list[i] != (MedianPixelList *) NULL)
1949       pixel_list[i]=DestroyMedianPixelList(pixel_list[i]);
1950   pixel_list=(MedianPixelList **) RelinquishAlignedMemory(pixel_list);
1951   return(pixel_list);
1952 }
1953
1954 static MedianPixelList *AcquireMedianPixelList(const unsigned long width)
1955 {
1956   MedianPixelList
1957     *pixel_list;
1958
1959   register long
1960     i;
1961
1962   pixel_list=(MedianPixelList *) AcquireAlignedMemory(1,sizeof(*pixel_list));
1963   if (pixel_list == (MedianPixelList *) NULL)
1964     return(pixel_list);
1965   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
1966   pixel_list->center=width*width/2;
1967   for (i=0; i < MedianListChannels; i++)
1968   {
1969     pixel_list->lists[i].nodes=(MedianListNode *) AcquireQuantumMemory(65537UL,
1970       sizeof(*pixel_list->lists[i].nodes));
1971     if (pixel_list->lists[i].nodes == (MedianListNode *) NULL)
1972       return(DestroyMedianPixelList(pixel_list));
1973     (void) ResetMagickMemory(pixel_list->lists[i].nodes,0,65537UL*
1974       sizeof(*pixel_list->lists[i].nodes));
1975   }
1976   pixel_list->signature=MagickSignature;
1977   return(pixel_list);
1978 }
1979
1980 static MedianPixelList **AcquireMedianPixelListThreadSet(
1981   const unsigned long width)
1982 {
1983   register long
1984     i;
1985
1986   MedianPixelList
1987     **pixel_list;
1988
1989   unsigned long
1990     number_threads;
1991
1992   number_threads=GetOpenMPMaximumThreads();
1993   pixel_list=(MedianPixelList **) AcquireAlignedMemory(number_threads,
1994     sizeof(*pixel_list));
1995   if (pixel_list == (MedianPixelList **) NULL)
1996     return((MedianPixelList **) NULL);
1997   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
1998   for (i=0; i < (long) number_threads; i++)
1999   {
2000     pixel_list[i]=AcquireMedianPixelList(width);
2001     if (pixel_list[i] == (MedianPixelList *) NULL)
2002       return(DestroyMedianPixelListThreadSet(pixel_list));
2003   }
2004   return(pixel_list);
2005 }
2006
2007 static void AddNodeMedianPixelList(MedianPixelList *pixel_list,
2008   const long channel,const unsigned long color)
2009 {
2010   register long
2011     level;
2012
2013   register MedianSkipList
2014     *list;
2015
2016   unsigned long
2017     search,
2018     update[9];
2019
2020   /*
2021     Initialize the node.
2022   */
2023   list=pixel_list->lists+channel;
2024   list->nodes[color].signature=pixel_list->signature;
2025   list->nodes[color].count=1;
2026   /*
2027     Determine where it belongs in the list.
2028   */
2029   search=65536UL;
2030   for (level=list->level; level >= 0; level--)
2031   {
2032     while (list->nodes[search].next[level] < color)
2033       search=list->nodes[search].next[level];
2034     update[level]=search;
2035   }
2036   /*
2037     Generate a pseudo-random level for this node.
2038   */
2039   for (level=0; ; level++)
2040   {
2041     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2042     if ((pixel_list->seed & 0x300) != 0x300)
2043       break;
2044   }
2045   if (level > 8)
2046     level=8;
2047   if (level > (list->level+2))
2048     level=list->level+2;
2049   /*
2050     If we're raising the list's level, link back to the root node.
2051   */
2052   while (level > list->level)
2053   {
2054     list->level++;
2055     update[list->level]=65536UL;
2056   }
2057   /*
2058     Link the node into the skip-list.
2059   */
2060   do
2061   {
2062     list->nodes[color].next[level]=list->nodes[update[level]].next[level];
2063     list->nodes[update[level]].next[level]=color;
2064   }
2065   while (level-- > 0);
2066 }
2067
2068 static MagickPixelPacket GetMedianPixelList(MedianPixelList *pixel_list)
2069 {
2070   MagickPixelPacket
2071     pixel;
2072
2073   register long
2074     channel;
2075
2076   register MedianSkipList
2077     *list;
2078
2079   unsigned long
2080     center,
2081     color,
2082     count;
2083
2084   unsigned short
2085     channels[MedianListChannels];
2086
2087   /*
2088     Find the median value for each of the color.
2089   */
2090   center=pixel_list->center;
2091   for (channel=0; channel < 5; channel++)
2092   {
2093     list=pixel_list->lists+channel;
2094     color=65536UL;
2095     count=0;
2096     do
2097     {
2098       color=list->nodes[color].next[0];
2099       count+=list->nodes[color].count;
2100     }
2101     while (count <= center);
2102     channels[channel]=(unsigned short) color;
2103   }
2104   GetMagickPixelPacket((const Image *) NULL,&pixel);
2105   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
2106   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
2107   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
2108   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
2109   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
2110   return(pixel);
2111 }
2112
2113 static inline void InsertMedianPixelList(const Image *image,
2114   const PixelPacket *pixel,const IndexPacket *indexes,
2115   MedianPixelList *pixel_list)
2116 {
2117   unsigned long
2118     signature;
2119
2120   unsigned short
2121     index;
2122
2123   index=ScaleQuantumToShort(pixel->red);
2124   signature=pixel_list->lists[0].nodes[index].signature;
2125   if (signature == pixel_list->signature)
2126     pixel_list->lists[0].nodes[index].count++;
2127   else
2128     AddNodeMedianPixelList(pixel_list,0,index);
2129   index=ScaleQuantumToShort(pixel->green);
2130   signature=pixel_list->lists[1].nodes[index].signature;
2131   if (signature == pixel_list->signature)
2132     pixel_list->lists[1].nodes[index].count++;
2133   else
2134     AddNodeMedianPixelList(pixel_list,1,index);
2135   index=ScaleQuantumToShort(pixel->blue);
2136   signature=pixel_list->lists[2].nodes[index].signature;
2137   if (signature == pixel_list->signature)
2138     pixel_list->lists[2].nodes[index].count++;
2139   else
2140     AddNodeMedianPixelList(pixel_list,2,index);
2141   index=ScaleQuantumToShort(pixel->opacity);
2142   signature=pixel_list->lists[3].nodes[index].signature;
2143   if (signature == pixel_list->signature)
2144     pixel_list->lists[3].nodes[index].count++;
2145   else
2146     AddNodeMedianPixelList(pixel_list,3,index);
2147   if (image->colorspace == CMYKColorspace)
2148     index=ScaleQuantumToShort(*indexes);
2149   signature=pixel_list->lists[4].nodes[index].signature;
2150   if (signature == pixel_list->signature)
2151     pixel_list->lists[4].nodes[index].count++;
2152   else
2153     AddNodeMedianPixelList(pixel_list,4,index);
2154 }
2155
2156 static void ResetMedianPixelList(MedianPixelList *pixel_list)
2157 {
2158   int
2159     level;
2160
2161   register long
2162     channel;
2163
2164   register MedianListNode
2165     *root;
2166
2167   register MedianSkipList
2168     *list;
2169
2170   /*
2171     Reset the skip-list.
2172   */
2173   for (channel=0; channel < 5; channel++)
2174   {
2175     list=pixel_list->lists+channel;
2176     root=list->nodes+65536UL;
2177     list->level=0;
2178     for (level=0; level < 9; level++)
2179       root->next[level]=65536UL;
2180   }
2181   pixel_list->seed=pixel_list->signature++;
2182 }
2183
2184 MagickExport Image *MedianFilterImage(const Image *image,const double radius,
2185   ExceptionInfo *exception)
2186 {
2187 #define MedianFilterImageTag  "MedianFilter/Image"
2188
2189   Image
2190     *median_image;
2191
2192   long
2193     progress,
2194     y;
2195
2196   MagickBooleanType
2197     status;
2198
2199   MedianPixelList
2200     **pixel_list;
2201
2202   unsigned long
2203     width;
2204
2205   CacheView
2206     *image_view,
2207     *median_view;
2208
2209   /*
2210     Initialize median image attributes.
2211   */
2212   assert(image != (Image *) NULL);
2213   assert(image->signature == MagickSignature);
2214   if (image->debug != MagickFalse)
2215     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2216   assert(exception != (ExceptionInfo *) NULL);
2217   assert(exception->signature == MagickSignature);
2218   width=GetOptimalKernelWidth2D(radius,0.5);
2219   if ((image->columns < width) || (image->rows < width))
2220     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
2221   median_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2222     exception);
2223   if (median_image == (Image *) NULL)
2224     return((Image *) NULL);
2225   if (SetImageStorageClass(median_image,DirectClass) == MagickFalse)
2226     {
2227       InheritException(exception,&median_image->exception);
2228       median_image=DestroyImage(median_image);
2229       return((Image *) NULL);
2230     }
2231   pixel_list=AcquireMedianPixelListThreadSet(width);
2232   if (pixel_list == (MedianPixelList **) NULL)
2233     {
2234       median_image=DestroyImage(median_image);
2235       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2236     }
2237   /*
2238     Median filter each image row.
2239   */
2240   status=MagickTrue;
2241   progress=0;
2242   image_view=AcquireCacheView(image);
2243   median_view=AcquireCacheView(median_image);
2244 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2245   #pragma omp parallel for schedule(dynamic) shared(progress,status)
2246 #endif
2247   for (y=0; y < (long) median_image->rows; y++)
2248   {
2249     register const IndexPacket
2250       *__restrict indexes;
2251
2252     register const PixelPacket
2253       *__restrict p;
2254
2255     register IndexPacket
2256       *__restrict median_indexes;
2257
2258     register long
2259       id,
2260       x;
2261
2262     register PixelPacket
2263       *__restrict q;
2264
2265     if (status == MagickFalse)
2266       continue;
2267     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
2268       2L),image->columns+width,width,exception);
2269     q=QueueCacheViewAuthenticPixels(median_view,0,y,median_image->columns,1,
2270       exception);
2271     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2272       {
2273         status=MagickFalse;
2274         continue;
2275       }
2276     indexes=GetCacheViewVirtualIndexQueue(image_view);
2277     median_indexes=GetCacheViewAuthenticIndexQueue(median_view);
2278     id=GetOpenMPThreadId();
2279     for (x=0; x < (long) median_image->columns; x++)
2280     {
2281       MagickPixelPacket
2282         pixel;
2283
2284       register const PixelPacket
2285         *__restrict r;
2286
2287       register const IndexPacket
2288         *__restrict s;
2289
2290       register long
2291         u,
2292         v;
2293
2294       r=p;
2295       s=indexes+x;
2296       ResetMedianPixelList(pixel_list[id]);
2297       for (v=0; v < (long) width; v++)
2298       {
2299         for (u=0; u < (long) width; u++)
2300           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
2301         r+=image->columns+width;
2302         s+=image->columns+width;
2303       }
2304       pixel=GetMedianPixelList(pixel_list[id]);
2305       SetPixelPacket(median_image,&pixel,q,median_indexes+x);
2306       p++;
2307       q++;
2308     }
2309     if (SyncCacheViewAuthenticPixels(median_view,exception) == MagickFalse)
2310       status=MagickFalse;
2311     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2312       {
2313         MagickBooleanType
2314           proceed;
2315
2316 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2317   #pragma omp critical (MagickCore_MedianFilterImage)
2318 #endif
2319         proceed=SetImageProgress(image,MedianFilterImageTag,progress++,
2320           image->rows);
2321         if (proceed == MagickFalse)
2322           status=MagickFalse;
2323       }
2324   }
2325   median_view=DestroyCacheView(median_view);
2326   image_view=DestroyCacheView(image_view);
2327   pixel_list=DestroyMedianPixelListThreadSet(pixel_list);
2328   return(median_image);
2329 }
2330 \f
2331 /*
2332 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2333 %                                                                             %
2334 %                                                                             %
2335 %                                                                             %
2336 %     M o t i o n B l u r I m a g e                                           %
2337 %                                                                             %
2338 %                                                                             %
2339 %                                                                             %
2340 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2341 %
2342 %  MotionBlurImage() simulates motion blur.  We convolve the image with a
2343 %  Gaussian operator of the given radius and standard deviation (sigma).
2344 %  For reasonable results, radius should be larger than sigma.  Use a
2345 %  radius of 0 and MotionBlurImage() selects a suitable radius for you.
2346 %  Angle gives the angle of the blurring motion.
2347 %
2348 %  Andrew Protano contributed this effect.
2349 %
2350 %  The format of the MotionBlurImage method is:
2351 %
2352 %    Image *MotionBlurImage(const Image *image,const double radius,
2353 %      const double sigma,const double angle,ExceptionInfo *exception)
2354 %    Image *MotionBlurImageChannel(const Image *image,const ChannelType channel,
2355 %      const double radius,const double sigma,const double angle,
2356 %      ExceptionInfo *exception)
2357 %
2358 %  A description of each parameter follows:
2359 %
2360 %    o image: the image.
2361 %
2362 %    o channel: the channel type.
2363 %
2364 %    o radius: the radius of the Gaussian, in pixels, not counting the center
2365 %    o radius: the radius of the Gaussian, in pixels, not counting
2366 %      the center pixel.
2367 %
2368 %    o sigma: the standard deviation of the Gaussian, in pixels.
2369 %
2370 %    o angle: Apply the effect along this angle.
2371 %
2372 %    o exception: return any errors or warnings in this structure.
2373 %
2374 */
2375
2376 static double *GetMotionBlurKernel(unsigned long width,
2377   const MagickRealType sigma)
2378 {
2379 #define KernelRank 3
2380
2381   double
2382     *kernel;
2383
2384   long
2385     bias;
2386
2387   MagickRealType
2388     alpha,
2389     normalize;
2390
2391   register long
2392     i;
2393
2394   /*
2395     Generate a 1-D convolution kernel.
2396   */
2397   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
2398   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
2399   if (kernel == (double *) NULL)
2400     return(kernel);
2401   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
2402   bias=(long) (KernelRank*width);
2403   for (i=0; i < (long) bias; i++)
2404   {
2405     alpha=exp((-((double) (i*i))/(double) (2.0*KernelRank*KernelRank*
2406       MagickSigma*MagickSigma)));
2407     kernel[i/KernelRank]+=(double) alpha/(MagickSQ2PI*sigma);
2408   }
2409   normalize=0.0;
2410   for (i=0; i < (long) width; i++)
2411     normalize+=kernel[i];
2412   for (i=0; i < (long) width; i++)
2413     kernel[i]/=normalize;
2414   return(kernel);
2415 }
2416
2417 MagickExport Image *MotionBlurImage(const Image *image,const double radius,
2418   const double sigma,const double angle,ExceptionInfo *exception)
2419 {
2420   Image
2421     *motion_blur;
2422
2423   motion_blur=MotionBlurImageChannel(image,DefaultChannels,radius,sigma,angle,
2424     exception);
2425   return(motion_blur);
2426 }
2427
2428 MagickExport Image *MotionBlurImageChannel(const Image *image,
2429   const ChannelType channel,const double radius,const double sigma,
2430   const double angle,ExceptionInfo *exception)
2431 {
2432   typedef struct _OffsetInfo
2433   {
2434     long
2435       x,
2436       y;
2437   } OffsetInfo;
2438
2439   double
2440     *kernel;
2441
2442   Image
2443     *blur_image;
2444
2445   long
2446     progress,
2447     y;
2448
2449   MagickBooleanType
2450     status;
2451
2452   MagickPixelPacket
2453     zero;
2454
2455   OffsetInfo
2456     *offset;
2457
2458   PointInfo
2459     point;
2460
2461   register long
2462     i;
2463
2464   unsigned long
2465     width;
2466
2467   CacheView
2468     *blur_view,
2469     *image_view;
2470
2471   assert(image != (Image *) NULL);
2472   assert(image->signature == MagickSignature);
2473   if (image->debug != MagickFalse)
2474     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2475   assert(exception != (ExceptionInfo *) NULL);
2476   width=GetOptimalKernelWidth1D(radius,sigma);
2477   kernel=GetMotionBlurKernel(width,sigma);
2478   if (kernel == (double *) NULL)
2479     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2480   offset=(OffsetInfo *) AcquireQuantumMemory(width,sizeof(*offset));
2481   if (offset == (OffsetInfo *) NULL)
2482     {
2483       kernel=(double *) RelinquishMagickMemory(kernel);
2484       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2485     }
2486   blur_image=CloneImage(image,0,0,MagickTrue,exception);
2487   if (blur_image == (Image *) NULL)
2488     {
2489       kernel=(double *) RelinquishMagickMemory(kernel);
2490       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2491       return((Image *) NULL);
2492     }
2493   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
2494     {
2495       kernel=(double *) RelinquishMagickMemory(kernel);
2496       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2497       InheritException(exception,&blur_image->exception);
2498       blur_image=DestroyImage(blur_image);
2499       return((Image *) NULL);
2500     }
2501   point.x=(double) width*sin(DegreesToRadians(angle));
2502   point.y=(double) width*cos(DegreesToRadians(angle));
2503   for (i=0; i < (long) width; i++)
2504   {
2505     offset[i].x=(long) ((i*point.y)/hypot(point.x,point.y)+0.5);
2506     offset[i].y=(long) ((i*point.x)/hypot(point.x,point.y)+0.5);
2507   }
2508   /*
2509     Motion blur image.
2510   */
2511   status=MagickTrue;
2512   progress=0;
2513   GetMagickPixelPacket(image,&zero);
2514   image_view=AcquireCacheView(image);
2515   blur_view=AcquireCacheView(blur_image);
2516 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2517   #pragma omp parallel for schedule(dynamic) shared(progress,status)
2518 #endif
2519   for (y=0; y < (long) image->rows; y++)
2520   {
2521     register IndexPacket
2522       *__restrict blur_indexes;
2523
2524     register long
2525       x;
2526
2527     register PixelPacket
2528       *__restrict q;
2529
2530     if (status == MagickFalse)
2531       continue;
2532     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
2533       exception);
2534     if (q == (PixelPacket *) NULL)
2535       {
2536         status=MagickFalse;
2537         continue;
2538       }
2539     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
2540     for (x=0; x < (long) image->columns; x++)
2541     {
2542       MagickPixelPacket
2543         qixel;
2544
2545       PixelPacket
2546         pixel;
2547
2548       register double
2549         *__restrict k;
2550
2551       register long
2552         i;
2553
2554       register const IndexPacket
2555         *__restrict indexes;
2556
2557       k=kernel;
2558       qixel=zero;
2559       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
2560         {
2561           for (i=0; i < (long) width; i++)
2562           {
2563             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
2564               offset[i].y,&pixel,exception);
2565             qixel.red+=(*k)*pixel.red;
2566             qixel.green+=(*k)*pixel.green;
2567             qixel.blue+=(*k)*pixel.blue;
2568             qixel.opacity+=(*k)*pixel.opacity;
2569             if (image->colorspace == CMYKColorspace)
2570               {
2571                 indexes=GetCacheViewVirtualIndexQueue(image_view);
2572                 qixel.index+=(*k)*(*indexes);
2573               }
2574             k++;
2575           }
2576           if ((channel & RedChannel) != 0)
2577             q->red=RoundToQuantum(qixel.red);
2578           if ((channel & GreenChannel) != 0)
2579             q->green=RoundToQuantum(qixel.green);
2580           if ((channel & BlueChannel) != 0)
2581             q->blue=RoundToQuantum(qixel.blue);
2582           if ((channel & OpacityChannel) != 0)
2583             q->opacity=RoundToQuantum(qixel.opacity);
2584           if (((channel & IndexChannel) != 0) &&
2585               (image->colorspace == CMYKColorspace))
2586             blur_indexes[x]=(IndexPacket) RoundToQuantum(qixel.index);
2587         }
2588       else
2589         {
2590           MagickRealType
2591             alpha,
2592             gamma;
2593
2594           alpha=0.0;
2595           gamma=0.0;
2596           for (i=0; i < (long) width; i++)
2597           {
2598             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
2599               offset[i].y,&pixel,exception);
2600             alpha=(MagickRealType) (QuantumScale*(QuantumRange-pixel.opacity));
2601             qixel.red+=(*k)*alpha*pixel.red;
2602             qixel.green+=(*k)*alpha*pixel.green;
2603             qixel.blue+=(*k)*alpha*pixel.blue;
2604             qixel.opacity+=(*k)*pixel.opacity;
2605             if (image->colorspace == CMYKColorspace)
2606               {
2607                 indexes=GetCacheViewVirtualIndexQueue(image_view);
2608                 qixel.index+=(*k)*alpha*(*indexes);
2609               }
2610             gamma+=(*k)*alpha;
2611             k++;
2612           }
2613           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2614           if ((channel & RedChannel) != 0)
2615             q->red=RoundToQuantum(gamma*qixel.red);
2616           if ((channel & GreenChannel) != 0)
2617             q->green=RoundToQuantum(gamma*qixel.green);
2618           if ((channel & BlueChannel) != 0)
2619             q->blue=RoundToQuantum(gamma*qixel.blue);
2620           if ((channel & OpacityChannel) != 0)
2621             q->opacity=RoundToQuantum(qixel.opacity);
2622           if (((channel & IndexChannel) != 0) &&
2623               (image->colorspace == CMYKColorspace))
2624             blur_indexes[x]=(IndexPacket) RoundToQuantum(gamma*qixel.index);
2625         }
2626       q++;
2627     }
2628     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
2629       status=MagickFalse;
2630     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2631       {
2632         MagickBooleanType
2633           proceed;
2634
2635 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2636   #pragma omp critical (MagickCore_MotionBlurImageChannel)
2637 #endif
2638         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
2639         if (proceed == MagickFalse)
2640           status=MagickFalse;
2641       }
2642   }
2643   blur_view=DestroyCacheView(blur_view);
2644   image_view=DestroyCacheView(image_view);
2645   kernel=(double *) RelinquishMagickMemory(kernel);
2646   offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2647   if (status == MagickFalse)
2648     blur_image=DestroyImage(blur_image);
2649   return(blur_image);
2650 }
2651 \f
2652 /*
2653 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2654 %                                                                             %
2655 %                                                                             %
2656 %                                                                             %
2657 %     P r e v i e w I m a g e                                                 %
2658 %                                                                             %
2659 %                                                                             %
2660 %                                                                             %
2661 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2662 %
2663 %  PreviewImage() tiles 9 thumbnails of the specified image with an image
2664 %  processing operation applied with varying parameters.  This may be helpful
2665 %  pin-pointing an appropriate parameter for a particular image processing
2666 %  operation.
2667 %
2668 %  The format of the PreviewImages method is:
2669 %
2670 %      Image *PreviewImages(const Image *image,const PreviewType preview,
2671 %        ExceptionInfo *exception)
2672 %
2673 %  A description of each parameter follows:
2674 %
2675 %    o image: the image.
2676 %
2677 %    o preview: the image processing operation.
2678 %
2679 %    o exception: return any errors or warnings in this structure.
2680 %
2681 */
2682 MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
2683   ExceptionInfo *exception)
2684 {
2685 #define NumberTiles  9
2686 #define PreviewImageTag  "Preview/Image"
2687 #define DefaultPreviewGeometry  "204x204+10+10"
2688
2689   char
2690     factor[MaxTextExtent],
2691     label[MaxTextExtent];
2692
2693   double
2694     degrees,
2695     gamma,
2696     percentage,
2697     radius,
2698     sigma,
2699     threshold;
2700
2701   Image
2702     *images,
2703     *montage_image,
2704     *preview_image,
2705     *thumbnail;
2706
2707   ImageInfo
2708     *preview_info;
2709
2710   long
2711     y;
2712
2713   MagickBooleanType
2714     proceed;
2715
2716   MontageInfo
2717     *montage_info;
2718
2719   QuantizeInfo
2720     quantize_info;
2721
2722   RectangleInfo
2723     geometry;
2724
2725   register long
2726     i,
2727     x;
2728
2729   unsigned long
2730     colors;
2731
2732   /*
2733     Open output image file.
2734   */
2735   assert(image != (Image *) NULL);
2736   assert(image->signature == MagickSignature);
2737   if (image->debug != MagickFalse)
2738     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2739   colors=2;
2740   degrees=0.0;
2741   gamma=(-0.2f);
2742   preview_info=AcquireImageInfo();
2743   SetGeometry(image,&geometry);
2744   (void) ParseMetaGeometry(DefaultPreviewGeometry,&geometry.x,&geometry.y,
2745     &geometry.width,&geometry.height);
2746   images=NewImageList();
2747   percentage=12.5;
2748   GetQuantizeInfo(&quantize_info);
2749   radius=0.0;
2750   sigma=1.0;
2751   threshold=0.0;
2752   x=0;
2753   y=0;
2754   for (i=0; i < NumberTiles; i++)
2755   {
2756     thumbnail=ThumbnailImage(image,geometry.width,geometry.height,exception);
2757     if (thumbnail == (Image *) NULL)
2758       break;
2759     (void) SetImageProgressMonitor(thumbnail,(MagickProgressMonitor) NULL,
2760       (void *) NULL);
2761     (void) SetImageProperty(thumbnail,"label",DefaultTileLabel);
2762     if (i == (NumberTiles/2))
2763       {
2764         (void) QueryColorDatabase("#dfdfdf",&thumbnail->matte_color,exception);
2765         AppendImageToList(&images,thumbnail);
2766         continue;
2767       }
2768     switch (preview)
2769     {
2770       case RotatePreview:
2771       {
2772         degrees+=45.0;
2773         preview_image=RotateImage(thumbnail,degrees,exception);
2774         (void) FormatMagickString(label,MaxTextExtent,"rotate %g",degrees);
2775         break;
2776       }
2777       case ShearPreview:
2778       {
2779         degrees+=5.0;
2780         preview_image=ShearImage(thumbnail,degrees,degrees,exception);
2781         (void) FormatMagickString(label,MaxTextExtent,"shear %gx%g",
2782           degrees,2.0*degrees);
2783         break;
2784       }
2785       case RollPreview:
2786       {
2787         x=(long) ((i+1)*thumbnail->columns)/NumberTiles;
2788         y=(long) ((i+1)*thumbnail->rows)/NumberTiles;
2789         preview_image=RollImage(thumbnail,x,y,exception);
2790         (void) FormatMagickString(label,MaxTextExtent,"roll %ldx%ld",x,y);
2791         break;
2792       }
2793       case HuePreview:
2794       {
2795         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2796         if (preview_image == (Image *) NULL)
2797           break;
2798         (void) FormatMagickString(factor,MaxTextExtent,"100,100,%g",
2799           2.0*percentage);
2800         (void) ModulateImage(preview_image,factor);
2801         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
2802         break;
2803       }
2804       case SaturationPreview:
2805       {
2806         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2807         if (preview_image == (Image *) NULL)
2808           break;
2809         (void) FormatMagickString(factor,MaxTextExtent,"100,%g",2.0*percentage);
2810         (void) ModulateImage(preview_image,factor);
2811         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
2812         break;
2813       }
2814       case BrightnessPreview:
2815       {
2816         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2817         if (preview_image == (Image *) NULL)
2818           break;
2819         (void) FormatMagickString(factor,MaxTextExtent,"%g",2.0*percentage);
2820         (void) ModulateImage(preview_image,factor);
2821         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
2822         break;
2823       }
2824       case GammaPreview:
2825       default:
2826       {
2827         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2828         if (preview_image == (Image *) NULL)
2829           break;
2830         gamma+=0.4f;
2831         (void) GammaImageChannel(preview_image,DefaultChannels,gamma);
2832         (void) FormatMagickString(label,MaxTextExtent,"gamma %g",gamma);
2833         break;
2834       }
2835       case SpiffPreview:
2836       {
2837         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2838         if (preview_image != (Image *) NULL)
2839           for (x=0; x < i; x++)
2840             (void) ContrastImage(preview_image,MagickTrue);
2841         (void) FormatMagickString(label,MaxTextExtent,"contrast (%ld)",i+1);
2842         break;
2843       }
2844       case DullPreview:
2845       {
2846         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2847         if (preview_image == (Image *) NULL)
2848           break;
2849         for (x=0; x < i; x++)
2850           (void) ContrastImage(preview_image,MagickFalse);
2851         (void) FormatMagickString(label,MaxTextExtent,"+contrast (%ld)",i+1);
2852         break;
2853       }
2854       case GrayscalePreview:
2855       {
2856         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2857         if (preview_image == (Image *) NULL)
2858           break;
2859         colors<<=1;
2860         quantize_info.number_colors=colors;
2861         quantize_info.colorspace=GRAYColorspace;
2862         (void) QuantizeImage(&quantize_info,preview_image);
2863         (void) FormatMagickString(label,MaxTextExtent,
2864           "-colorspace gray -colors %ld",colors);
2865         break;
2866       }
2867       case QuantizePreview:
2868       {
2869         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2870         if (preview_image == (Image *) NULL)
2871           break;
2872         colors<<=1;
2873         quantize_info.number_colors=colors;
2874         (void) QuantizeImage(&quantize_info,preview_image);
2875         (void) FormatMagickString(label,MaxTextExtent,"colors %ld",colors);
2876         break;
2877       }
2878       case DespecklePreview:
2879       {
2880         for (x=0; x < (i-1); x++)
2881         {
2882           preview_image=DespeckleImage(thumbnail,exception);
2883           if (preview_image == (Image *) NULL)
2884             break;
2885           thumbnail=DestroyImage(thumbnail);
2886           thumbnail=preview_image;
2887         }
2888         preview_image=DespeckleImage(thumbnail,exception);
2889         if (preview_image == (Image *) NULL)
2890           break;
2891         (void) FormatMagickString(label,MaxTextExtent,"despeckle (%ld)",i+1);
2892         break;
2893       }
2894       case ReduceNoisePreview:
2895       {
2896         preview_image=ReduceNoiseImage(thumbnail,radius,exception);
2897         (void) FormatMagickString(label,MaxTextExtent,"noise %g",radius);
2898         break;
2899       }
2900       case AddNoisePreview:
2901       {
2902         switch ((int) i)
2903         {
2904           case 0:
2905           {
2906             (void) CopyMagickString(factor,"uniform",MaxTextExtent);
2907             break;
2908           }
2909           case 1:
2910           {
2911             (void) CopyMagickString(factor,"gaussian",MaxTextExtent);
2912             break;
2913           }
2914           case 2:
2915           {
2916             (void) CopyMagickString(factor,"multiplicative",MaxTextExtent);
2917             break;
2918           }
2919           case 3:
2920           {
2921             (void) CopyMagickString(factor,"impulse",MaxTextExtent);
2922             break;
2923           }
2924           case 4:
2925           {
2926             (void) CopyMagickString(factor,"laplacian",MaxTextExtent);
2927             break;
2928           }
2929           case 5:
2930           {
2931             (void) CopyMagickString(factor,"Poisson",MaxTextExtent);
2932             break;
2933           }
2934           default:
2935           {
2936             (void) CopyMagickString(thumbnail->magick,"NULL",MaxTextExtent);
2937             break;
2938           }
2939         }
2940         preview_image=ReduceNoiseImage(thumbnail,(double) i,exception);
2941         (void) FormatMagickString(label,MaxTextExtent,"+noise %s",factor);
2942         break;
2943       }
2944       case SharpenPreview:
2945       {
2946         preview_image=SharpenImage(thumbnail,radius,sigma,exception);
2947         (void) FormatMagickString(label,MaxTextExtent,"sharpen %gx%g",radius,
2948           sigma);
2949         break;
2950       }
2951       case BlurPreview:
2952       {
2953         preview_image=BlurImage(thumbnail,radius,sigma,exception);
2954         (void) FormatMagickString(label,MaxTextExtent,"blur %gx%g",radius,
2955           sigma);
2956         break;
2957       }
2958       case ThresholdPreview:
2959       {
2960         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2961         if (preview_image == (Image *) NULL)
2962           break;
2963         (void) BilevelImage(thumbnail,
2964           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
2965         (void) FormatMagickString(label,MaxTextExtent,"threshold %g",
2966           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
2967         break;
2968       }
2969       case EdgeDetectPreview:
2970       {
2971         preview_image=EdgeImage(thumbnail,radius,exception);
2972         (void) FormatMagickString(label,MaxTextExtent,"edge %g",radius);
2973         break;
2974       }
2975       case SpreadPreview:
2976       {
2977         preview_image=SpreadImage(thumbnail,radius,exception);
2978         (void) FormatMagickString(label,MaxTextExtent,"spread %g",radius+0.5);
2979         break;
2980       }
2981       case SolarizePreview:
2982       {
2983         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2984         if (preview_image == (Image *) NULL)
2985           break;
2986         (void) SolarizeImage(preview_image,(double) QuantumRange*
2987           percentage/100.0);
2988         (void) FormatMagickString(label,MaxTextExtent,"solarize %g",
2989           (QuantumRange*percentage)/100.0);
2990         break;
2991       }
2992       case ShadePreview:
2993       {
2994         degrees+=10.0;
2995         preview_image=ShadeImage(thumbnail,MagickTrue,degrees,degrees,
2996           exception);
2997         (void) FormatMagickString(label,MaxTextExtent,"shade %gx%g",degrees,
2998           degrees);
2999         break;
3000       }
3001       case RaisePreview:
3002       {
3003         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3004         if (preview_image == (Image *) NULL)
3005           break;
3006         geometry.width=(unsigned long) (2*i+2);
3007         geometry.height=(unsigned long) (2*i+2);
3008         geometry.x=i/2;
3009         geometry.y=i/2;
3010         (void) RaiseImage(preview_image,&geometry,MagickTrue);
3011         (void) FormatMagickString(label,MaxTextExtent,"raise %lux%lu%+ld%+ld",
3012           geometry.width,geometry.height,geometry.x,geometry.y);
3013         break;
3014       }
3015       case SegmentPreview:
3016       {
3017         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3018         if (preview_image == (Image *) NULL)
3019           break;
3020         threshold+=0.4f;
3021         (void) SegmentImage(preview_image,RGBColorspace,MagickFalse,threshold,
3022           threshold);
3023         (void) FormatMagickString(label,MaxTextExtent,"segment %gx%g",
3024           threshold,threshold);
3025         break;
3026       }
3027       case SwirlPreview:
3028       {
3029         preview_image=SwirlImage(thumbnail,degrees,exception);
3030         (void) FormatMagickString(label,MaxTextExtent,"swirl %g",degrees);
3031         degrees+=45.0;
3032         break;
3033       }
3034       case ImplodePreview:
3035       {
3036         degrees+=0.1f;
3037         preview_image=ImplodeImage(thumbnail,degrees,exception);
3038         (void) FormatMagickString(label,MaxTextExtent,"implode %g",degrees);
3039         break;
3040       }
3041       case WavePreview:
3042       {
3043         degrees+=5.0f;
3044         preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,exception);
3045         (void) FormatMagickString(label,MaxTextExtent,"wave %gx%g",0.5*degrees,
3046           2.0*degrees);
3047         break;
3048       }
3049       case OilPaintPreview:
3050       {
3051         preview_image=OilPaintImage(thumbnail,(double) radius,exception);
3052         (void) FormatMagickString(label,MaxTextExtent,"paint %g",radius);
3053         break;
3054       }
3055       case CharcoalDrawingPreview:
3056       {
3057         preview_image=CharcoalImage(thumbnail,(double) radius,(double) sigma,
3058           exception);
3059         (void) FormatMagickString(label,MaxTextExtent,"charcoal %gx%g",radius,
3060           sigma);
3061         break;
3062       }
3063       case JPEGPreview:
3064       {
3065         char
3066           filename[MaxTextExtent];
3067
3068         int
3069           file;
3070
3071         MagickBooleanType
3072           status;
3073
3074         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3075         if (preview_image == (Image *) NULL)
3076           break;
3077         preview_info->quality=(unsigned long) percentage;
3078         (void) FormatMagickString(factor,MaxTextExtent,"%lu",
3079           preview_info->quality);
3080         file=AcquireUniqueFileResource(filename);
3081         if (file != -1)
3082           file=close(file)-1;
3083         (void) FormatMagickString(preview_image->filename,MaxTextExtent,
3084           "jpeg:%s",filename);
3085         status=WriteImage(preview_info,preview_image);
3086         if (status != MagickFalse)
3087           {
3088             Image
3089               *quality_image;
3090
3091             (void) CopyMagickString(preview_info->filename,
3092               preview_image->filename,MaxTextExtent);
3093             quality_image=ReadImage(preview_info,exception);
3094             if (quality_image != (Image *) NULL)
3095               {
3096                 preview_image=DestroyImage(preview_image);
3097                 preview_image=quality_image;
3098               }
3099           }
3100         (void) RelinquishUniqueFileResource(preview_image->filename);
3101         if ((GetBlobSize(preview_image)/1024) >= 1024)
3102           (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%gmb ",
3103             factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
3104             1024.0/1024.0);
3105         else
3106           if (GetBlobSize(preview_image) >= 1024)
3107             (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%gkb ",
3108               factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
3109               1024.0);
3110           else
3111             (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%lub ",
3112               factor,(unsigned long) GetBlobSize(thumbnail));
3113         break;
3114       }
3115     }
3116     thumbnail=DestroyImage(thumbnail);
3117     percentage+=12.5;
3118     radius+=0.5;
3119     sigma+=0.25;
3120     if (preview_image == (Image *) NULL)
3121       break;
3122     (void) DeleteImageProperty(preview_image,"label");
3123     (void) SetImageProperty(preview_image,"label",label);
3124     AppendImageToList(&images,preview_image);
3125     proceed=SetImageProgress(image,PreviewImageTag,i,NumberTiles);
3126     if (proceed == MagickFalse)
3127       break;
3128   }
3129   if (images == (Image *) NULL)
3130     {
3131       preview_info=DestroyImageInfo(preview_info);
3132       return((Image *) NULL);
3133     }
3134   /*
3135     Create the montage.
3136   */
3137   montage_info=CloneMontageInfo(preview_info,(MontageInfo *) NULL);
3138   (void) CopyMagickString(montage_info->filename,image->filename,MaxTextExtent);
3139   montage_info->shadow=MagickTrue;
3140   (void) CloneString(&montage_info->tile,"3x3");
3141   (void) CloneString(&montage_info->geometry,DefaultPreviewGeometry);
3142   (void) CloneString(&montage_info->frame,DefaultTileFrame);
3143   montage_image=MontageImages(images,montage_info,exception);
3144   montage_info=DestroyMontageInfo(montage_info);
3145   images=DestroyImageList(images);
3146   if (montage_image == (Image *) NULL)
3147     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3148   if (montage_image->montage != (char *) NULL)
3149     {
3150       /*
3151         Free image directory.
3152       */
3153       montage_image->montage=(char *) RelinquishMagickMemory(
3154         montage_image->montage);
3155       if (image->directory != (char *) NULL)
3156         montage_image->directory=(char *) RelinquishMagickMemory(
3157           montage_image->directory);
3158     }
3159   preview_info=DestroyImageInfo(preview_info);
3160   return(montage_image);
3161 }
3162 \f
3163 /*
3164 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3165 %                                                                             %
3166 %                                                                             %
3167 %                                                                             %
3168 %     R a d i a l B l u r I m a g e                                           %
3169 %                                                                             %
3170 %                                                                             %
3171 %                                                                             %
3172 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3173 %
3174 %  RadialBlurImage() applies a radial blur to the image.
3175 %
3176 %  Andrew Protano contributed this effect.
3177 %
3178 %  The format of the RadialBlurImage method is:
3179 %
3180 %    Image *RadialBlurImage(const Image *image,const double angle,
3181 %      ExceptionInfo *exception)
3182 %    Image *RadialBlurImageChannel(const Image *image,const ChannelType channel,
3183 %      const double angle,ExceptionInfo *exception)
3184 %
3185 %  A description of each parameter follows:
3186 %
3187 %    o image: the image.
3188 %
3189 %    o channel: the channel type.
3190 %
3191 %    o angle: the angle of the radial blur.
3192 %
3193 %    o exception: return any errors or warnings in this structure.
3194 %
3195 */
3196
3197 MagickExport Image *RadialBlurImage(const Image *image,const double angle,
3198   ExceptionInfo *exception)
3199 {
3200   Image
3201     *blur_image;
3202
3203   blur_image=RadialBlurImageChannel(image,DefaultChannels,angle,exception);
3204   return(blur_image);
3205 }
3206
3207 MagickExport Image *RadialBlurImageChannel(const Image *image,
3208   const ChannelType channel,const double angle,ExceptionInfo *exception)
3209 {
3210   Image
3211     *blur_image;
3212
3213   long
3214     progress,
3215     y;
3216
3217   MagickBooleanType
3218     status;
3219
3220   MagickPixelPacket
3221     zero;
3222
3223   MagickRealType
3224     blur_radius,
3225     *cos_theta,
3226     offset,
3227     *sin_theta,
3228     theta;
3229
3230   PointInfo
3231     blur_center;
3232
3233   register long
3234     i;
3235
3236   unsigned long
3237     n;
3238
3239   CacheView
3240     *blur_view,
3241     *image_view;
3242
3243   /*
3244     Allocate blur image.
3245   */
3246   assert(image != (Image *) NULL);
3247   assert(image->signature == MagickSignature);
3248   if (image->debug != MagickFalse)
3249     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3250   assert(exception != (ExceptionInfo *) NULL);
3251   assert(exception->signature == MagickSignature);
3252   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3253   if (blur_image == (Image *) NULL)
3254     return((Image *) NULL);
3255   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
3256     {
3257       InheritException(exception,&blur_image->exception);
3258       blur_image=DestroyImage(blur_image);
3259       return((Image *) NULL);
3260     }
3261   blur_center.x=(double) image->columns/2.0;
3262   blur_center.y=(double) image->rows/2.0;
3263   blur_radius=hypot(blur_center.x,blur_center.y);
3264   n=(unsigned long) fabs(4.0*DegreesToRadians(angle)*sqrt((double) blur_radius)+
3265     2UL);
3266   theta=DegreesToRadians(angle)/(MagickRealType) (n-1);
3267   cos_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
3268     sizeof(*cos_theta));
3269   sin_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
3270     sizeof(*sin_theta));
3271   if ((cos_theta == (MagickRealType *) NULL) ||
3272       (sin_theta == (MagickRealType *) NULL))
3273     {
3274       blur_image=DestroyImage(blur_image);
3275       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3276     }
3277   offset=theta*(MagickRealType) (n-1)/2.0;
3278   for (i=0; i < (long) n; i++)
3279   {
3280     cos_theta[i]=cos((double) (theta*i-offset));
3281     sin_theta[i]=sin((double) (theta*i-offset));
3282   }
3283   /*
3284     Radial blur image.
3285   */
3286   status=MagickTrue;
3287   progress=0;
3288   GetMagickPixelPacket(image,&zero);
3289   image_view=AcquireCacheView(image);
3290   blur_view=AcquireCacheView(blur_image);
3291 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3292   #pragma omp parallel for schedule(dynamic) shared(progress,status)
3293 #endif
3294   for (y=0; y < (long) blur_image->rows; y++)
3295   {
3296     register const IndexPacket
3297       *__restrict indexes;
3298
3299     register IndexPacket
3300       *__restrict blur_indexes;
3301
3302     register long
3303       x;
3304
3305     register PixelPacket
3306       *__restrict q;
3307
3308     if (status == MagickFalse)
3309       continue;
3310     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
3311       exception);
3312     if (q == (PixelPacket *) NULL)
3313       {
3314         status=MagickFalse;
3315         continue;
3316       }
3317     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
3318     for (x=0; x < (long) blur_image->columns; x++)
3319     {
3320       MagickPixelPacket
3321         qixel;
3322
3323       MagickRealType
3324         normalize,
3325         radius;
3326
3327       PixelPacket
3328         pixel;
3329
3330       PointInfo
3331         center;
3332
3333       register long
3334         i;
3335
3336       unsigned long
3337         step;
3338
3339       center.x=(double) x-blur_center.x;
3340       center.y=(double) y-blur_center.y;
3341       radius=hypot((double) center.x,center.y);
3342       if (radius == 0)
3343         step=1;
3344       else
3345         {
3346           step=(unsigned long) (blur_radius/radius);
3347           if (step == 0)
3348             step=1;
3349           else
3350             if (step >= n)
3351               step=n-1;
3352         }
3353       normalize=0.0;
3354       qixel=zero;
3355       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
3356         {
3357           for (i=0; i < (long) n; i+=step)
3358           {
3359             (void) GetOneCacheViewVirtualPixel(image_view,(long) (blur_center.x+
3360               center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),(long) (
3361               blur_center.y+center.x*sin_theta[i]+center.y*cos_theta[i]+0.5),
3362               &pixel,exception);
3363             qixel.red+=pixel.red;
3364             qixel.green+=pixel.green;
3365             qixel.blue+=pixel.blue;
3366             qixel.opacity+=pixel.opacity;
3367             if (image->colorspace == CMYKColorspace)
3368               {
3369                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3370                 qixel.index+=(*indexes);
3371               }
3372             normalize+=1.0;
3373           }
3374           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
3375             normalize);
3376           if ((channel & RedChannel) != 0)
3377             q->red=RoundToQuantum(normalize*qixel.red);
3378           if ((channel & GreenChannel) != 0)
3379             q->green=RoundToQuantum(normalize*qixel.green);
3380           if ((channel & BlueChannel) != 0)
3381             q->blue=RoundToQuantum(normalize*qixel.blue);
3382           if ((channel & OpacityChannel) != 0)
3383             q->opacity=RoundToQuantum(normalize*qixel.opacity);
3384           if (((channel & IndexChannel) != 0) &&
3385               (image->colorspace == CMYKColorspace))
3386             blur_indexes[x]=(IndexPacket) RoundToQuantum(normalize*qixel.index);
3387         }
3388       else
3389         {
3390           MagickRealType
3391             alpha,
3392             gamma;
3393
3394           alpha=1.0;
3395           gamma=0.0;
3396           for (i=0; i < (long) n; i+=step)
3397           {
3398             (void) GetOneCacheViewVirtualPixel(image_view,(long) (blur_center.x+
3399               center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),(long) (
3400               blur_center.y+center.x*sin_theta[i]+center.y*cos_theta[i]+0.5),
3401               &pixel,exception);
3402             alpha=(MagickRealType) (QuantumScale*(QuantumRange-pixel.opacity));
3403             qixel.red+=alpha*pixel.red;
3404             qixel.green+=alpha*pixel.green;
3405             qixel.blue+=alpha*pixel.blue;
3406             qixel.opacity+=pixel.opacity;
3407             if (image->colorspace == CMYKColorspace)
3408               {
3409                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3410                 qixel.index+=alpha*(*indexes);
3411               }
3412             gamma+=alpha;
3413             normalize+=1.0;
3414           }
3415           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
3416           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
3417             normalize);
3418           if ((channel & RedChannel) != 0)
3419             q->red=RoundToQuantum(gamma*qixel.red);
3420           if ((channel & GreenChannel) != 0)
3421             q->green=RoundToQuantum(gamma*qixel.green);
3422           if ((channel & BlueChannel) != 0)
3423             q->blue=RoundToQuantum(gamma*qixel.blue);
3424           if ((channel & OpacityChannel) != 0)
3425             q->opacity=RoundToQuantum(normalize*qixel.opacity);
3426           if (((channel & IndexChannel) != 0) &&
3427               (image->colorspace == CMYKColorspace))
3428             blur_indexes[x]=(IndexPacket) RoundToQuantum(gamma*qixel.index);
3429         }
3430       q++;
3431     }
3432     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
3433       status=MagickFalse;
3434     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3435       {
3436         MagickBooleanType
3437           proceed;
3438
3439 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3440   #pragma omp critical (MagickCore_RadialBlurImageChannel)
3441 #endif
3442         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
3443         if (proceed == MagickFalse)
3444           status=MagickFalse;
3445       }
3446   }
3447   blur_view=DestroyCacheView(blur_view);
3448   image_view=DestroyCacheView(image_view);
3449   cos_theta=(MagickRealType *) RelinquishMagickMemory(cos_theta);
3450   sin_theta=(MagickRealType *) RelinquishMagickMemory(sin_theta);
3451   if (status == MagickFalse)
3452     blur_image=DestroyImage(blur_image);
3453   return(blur_image);
3454 }
3455 \f
3456 /*
3457 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3458 %                                                                             %
3459 %                                                                             %
3460 %                                                                             %
3461 %     R e d u c e N o i s e I m a g e                                         %
3462 %                                                                             %
3463 %                                                                             %
3464 %                                                                             %
3465 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3466 %
3467 %  ReduceNoiseImage() smooths the contours of an image while still preserving
3468 %  edge information.  The algorithm works by replacing each pixel with its
3469 %  neighbor closest in value.  A neighbor is defined by radius.  Use a radius
3470 %  of 0 and ReduceNoise() selects a suitable radius for you.
3471 %
3472 %  The format of the ReduceNoiseImage method is:
3473 %
3474 %      Image *ReduceNoiseImage(const Image *image,const double radius,
3475 %        ExceptionInfo *exception)
3476 %
3477 %  A description of each parameter follows:
3478 %
3479 %    o image: the image.
3480 %
3481 %    o radius: the radius of the pixel neighborhood.
3482 %
3483 %    o exception: return any errors or warnings in this structure.
3484 %
3485 */
3486
3487 static MagickPixelPacket GetNonpeakMedianPixelList(MedianPixelList *pixel_list)
3488 {
3489   MagickPixelPacket
3490     pixel;
3491
3492   register long
3493     channel;
3494
3495   register MedianSkipList
3496     *list;
3497
3498   unsigned long
3499     center,
3500     color,
3501     count,
3502     previous,
3503     next;
3504
3505   unsigned short
3506     channels[5];
3507
3508   /*
3509     Finds the median value for each of the color.
3510   */
3511   center=pixel_list->center;
3512   for (channel=0; channel < 5; channel++)
3513   {
3514     list=pixel_list->lists+channel;
3515     color=65536UL;
3516     next=list->nodes[color].next[0];
3517     count=0;
3518     do
3519     {
3520       previous=color;
3521       color=next;
3522       next=list->nodes[color].next[0];
3523       count+=list->nodes[color].count;
3524     }
3525     while (count <= center);
3526     if ((previous == 65536UL) && (next != 65536UL))
3527       color=next;
3528     else
3529       if ((previous != 65536UL) && (next == 65536UL))
3530         color=previous;
3531     channels[channel]=(unsigned short) color;
3532   }
3533   GetMagickPixelPacket((const Image *) NULL,&pixel);
3534   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3535   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3536   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3537   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3538   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3539   return(pixel);
3540 }
3541
3542 MagickExport Image *ReduceNoiseImage(const Image *image,const double radius,
3543   ExceptionInfo *exception)
3544 {
3545 #define ReduceNoiseImageTag  "ReduceNoise/Image"
3546
3547   Image
3548     *noise_image;
3549
3550   long
3551     progress,
3552     y;
3553
3554   MagickBooleanType
3555     status;
3556
3557   MedianPixelList
3558     **pixel_list;
3559
3560   unsigned long
3561     width;
3562
3563   CacheView
3564     *image_view,
3565     *noise_view;
3566
3567   /*
3568     Initialize noise image attributes.
3569   */
3570   assert(image != (Image *) NULL);
3571   assert(image->signature == MagickSignature);
3572   if (image->debug != MagickFalse)
3573     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3574   assert(exception != (ExceptionInfo *) NULL);
3575   assert(exception->signature == MagickSignature);
3576   width=GetOptimalKernelWidth2D(radius,0.5);
3577   if ((image->columns < width) || (image->rows < width))
3578     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
3579   noise_image=CloneImage(image,image->columns,image->rows,MagickTrue,
3580     exception);
3581   if (noise_image == (Image *) NULL)
3582     return((Image *) NULL);
3583   if (SetImageStorageClass(noise_image,DirectClass) == MagickFalse)
3584     {
3585       InheritException(exception,&noise_image->exception);
3586       noise_image=DestroyImage(noise_image);
3587       return((Image *) NULL);
3588     }
3589   pixel_list=AcquireMedianPixelListThreadSet(width);
3590   if (pixel_list == (MedianPixelList **) NULL)
3591     {
3592       noise_image=DestroyImage(noise_image);
3593       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3594     }
3595   /*
3596     Reduce noise image.
3597   */
3598   status=MagickTrue;
3599   progress=0;
3600   image_view=AcquireCacheView(image);
3601   noise_view=AcquireCacheView(noise_image);
3602 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3603   #pragma omp parallel for schedule(dynamic) shared(progress,status)
3604 #endif
3605   for (y=0; y < (long) noise_image->rows; y++)
3606   {
3607     register const IndexPacket
3608       *__restrict indexes;
3609
3610     register const PixelPacket
3611       *__restrict p;
3612
3613     register IndexPacket
3614       *__restrict noise_indexes;
3615
3616     register long
3617       id,
3618       x;
3619
3620     register PixelPacket
3621       *__restrict q;
3622
3623     if (status == MagickFalse)
3624       continue;
3625     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
3626       2L),image->columns+width,width,exception);
3627     q=QueueCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
3628       exception);
3629     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3630       {
3631         status=MagickFalse;
3632         continue;
3633       }
3634     indexes=GetCacheViewVirtualIndexQueue(image_view);
3635     noise_indexes=GetCacheViewAuthenticIndexQueue(noise_view);
3636     id=GetOpenMPThreadId();
3637     for (x=0; x < (long) noise_image->columns; x++)
3638     {
3639       MagickPixelPacket
3640         pixel;
3641
3642       register const PixelPacket
3643         *__restrict r;
3644
3645       register const IndexPacket
3646         *__restrict s;
3647
3648       register long
3649         u,
3650         v;
3651
3652       r=p;
3653       s=indexes+x;
3654       ResetMedianPixelList(pixel_list[id]);
3655       for (v=0; v < (long) width; v++)
3656       {
3657         for (u=0; u < (long) width; u++)
3658           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
3659         r+=image->columns+width;
3660         s+=image->columns+width;
3661       }
3662       pixel=GetNonpeakMedianPixelList(pixel_list[id]);
3663       SetPixelPacket(noise_image,&pixel,q,noise_indexes+x);
3664       p++;
3665       q++;
3666     }
3667     if (SyncCacheViewAuthenticPixels(noise_view,exception) == MagickFalse)
3668       status=MagickFalse;
3669     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3670       {
3671         MagickBooleanType
3672           proceed;
3673
3674 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3675   #pragma omp critical (MagickCore_ReduceNoiseImage)
3676 #endif
3677         proceed=SetImageProgress(image,ReduceNoiseImageTag,progress++,
3678           image->rows);
3679         if (proceed == MagickFalse)
3680           status=MagickFalse;
3681       }
3682   }
3683   noise_view=DestroyCacheView(noise_view);
3684   image_view=DestroyCacheView(image_view);
3685   pixel_list=DestroyMedianPixelListThreadSet(pixel_list);
3686   return(noise_image);
3687 }
3688 \f
3689 /*
3690 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3691 %                                                                             %
3692 %                                                                             %
3693 %                                                                             %
3694 %     S e l e c t i v e B l u r I m a g e                                     %
3695 %                                                                             %
3696 %                                                                             %
3697 %                                                                             %
3698 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3699 %
3700 %  SelectiveBlurImage() selectively blur pixels within a contrast threshold.
3701 %  It is similar to the unsharpen mask that sharpens everything with contrast
3702 %  above a certain threshold.
3703 %
3704 %  The format of the SelectiveBlurImage method is:
3705 %
3706 %      Image *SelectiveBlurImage(const Image *image,const double radius,
3707 %        const double sigma,const double threshold,ExceptionInfo *exception)
3708 %      Image *SelectiveBlurImageChannel(const Image *image,
3709 %        const ChannelType channel,const double radius,const double sigma,
3710 %        const double threshold,ExceptionInfo *exception)
3711 %
3712 %  A description of each parameter follows:
3713 %
3714 %    o image: the image.
3715 %
3716 %    o channel: the channel type.
3717 %
3718 %    o radius: the radius of the Gaussian, in pixels, not counting the center
3719 %      pixel.
3720 %
3721 %    o sigma: the standard deviation of the Gaussian, in pixels.
3722 %
3723 %    o threshold: only pixels within this contrast threshold are included
3724 %      in the blur operation.
3725 %
3726 %    o exception: return any errors or warnings in this structure.
3727 %
3728 */
3729
3730 static inline MagickBooleanType SelectiveContrast(const PixelPacket *p,
3731   const PixelPacket *q,const double threshold)
3732 {
3733   if (fabs(PixelIntensity(p)-PixelIntensity(q)) < threshold)
3734     return(MagickTrue);
3735   return(MagickFalse);
3736 }
3737
3738 MagickExport Image *SelectiveBlurImage(const Image *image,const double radius,
3739   const double sigma,const double threshold,ExceptionInfo *exception)
3740 {
3741   Image
3742     *blur_image;
3743
3744   blur_image=SelectiveBlurImageChannel(image,DefaultChannels,radius,sigma,
3745     threshold,exception);
3746   return(blur_image);
3747 }
3748
3749 MagickExport Image *SelectiveBlurImageChannel(const Image *image,
3750   const ChannelType channel,const double radius,const double sigma,
3751   const double threshold,ExceptionInfo *exception)
3752 {
3753 #define SelectiveBlurImageTag  "SelectiveBlur/Image"
3754
3755   double
3756     alpha,
3757     *kernel;
3758
3759   Image
3760     *blur_image;
3761
3762   long
3763     progress,
3764     v,
3765     y;
3766
3767   MagickBooleanType
3768     status;
3769
3770   MagickPixelPacket
3771     zero;
3772
3773   MagickRealType
3774     bias;
3775
3776   register long
3777     i,
3778     u;
3779
3780   unsigned long
3781     width;
3782
3783   CacheView
3784     *blur_view,
3785     *image_view;
3786
3787   /*
3788     Initialize blur image attributes.
3789   */
3790   assert(image != (Image *) NULL);
3791   assert(image->signature == MagickSignature);
3792   if (image->debug != MagickFalse)
3793     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3794   assert(exception != (ExceptionInfo *) NULL);
3795   assert(exception->signature == MagickSignature);
3796   width=GetOptimalKernelWidth1D(radius,sigma);
3797   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
3798   if (kernel == (double *) NULL)
3799     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3800   i=0;
3801   for (v=(-((long) width/2)); v <= (long) (width/2); v++)
3802   {
3803     for (u=(-((long) width/2)); u <= (long) (width/2); u++)
3804     {
3805       alpha=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma));
3806       kernel[i]=(double) (alpha/(2.0*MagickPI*MagickSigma*MagickSigma));
3807       i++;
3808     }
3809   }
3810   if (image->debug != MagickFalse)
3811     {
3812       char
3813         format[MaxTextExtent],
3814         *message;
3815
3816       long
3817         u,
3818         v;
3819
3820       register const double
3821         *k;
3822
3823       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
3824         "  SelectiveBlurImage with %ldx%ld kernel:",width,width);
3825       message=AcquireString("");
3826       k=kernel;
3827       for (v=0; v < (long) width; v++)
3828       {
3829         *message='\0';
3830         (void) FormatMagickString(format,MaxTextExtent,"%ld: ",v);
3831         (void) ConcatenateString(&message,format);
3832         for (u=0; u < (long) width; u++)
3833         {
3834           (void) FormatMagickString(format,MaxTextExtent,"%+f ",*k++);
3835           (void) ConcatenateString(&message,format);
3836         }
3837         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
3838       }
3839       message=DestroyString(message);
3840     }
3841   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3842   if (blur_image == (Image *) NULL)
3843     return((Image *) NULL);
3844   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
3845     {
3846       InheritException(exception,&blur_image->exception);
3847       blur_image=DestroyImage(blur_image);
3848       return((Image *) NULL);
3849     }
3850   /*
3851     Threshold blur image.
3852   */
3853   status=MagickTrue;
3854   progress=0;
3855   GetMagickPixelPacket(image,&zero);
3856   bias=image->bias;
3857   image_view=AcquireCacheView(image);
3858   blur_view=AcquireCacheView(blur_image);
3859 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3860   #pragma omp parallel for schedule(dynamic) shared(progress,status)
3861 #endif
3862   for (y=0; y < (long) image->rows; y++)
3863   {
3864     MagickBooleanType
3865       sync;
3866
3867     MagickRealType
3868       gamma;
3869
3870     register const IndexPacket
3871       *__restrict indexes;
3872
3873     register const PixelPacket
3874       *__restrict p;
3875
3876     register IndexPacket
3877       *__restrict blur_indexes;
3878
3879     register long
3880       x;
3881
3882     register PixelPacket
3883       *__restrict q;
3884
3885     if (status == MagickFalse)
3886       continue;
3887     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
3888       2L),image->columns+width,width,exception);
3889     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
3890       exception);
3891     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3892       {
3893         status=MagickFalse;
3894         continue;
3895       }
3896     indexes=GetCacheViewVirtualIndexQueue(image_view);
3897     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
3898     for (x=0; x < (long) image->columns; x++)
3899     {
3900       long
3901         j,
3902         v;
3903
3904       MagickPixelPacket
3905         pixel;
3906
3907       register const double
3908         *__restrict k;
3909
3910       register long
3911         u;
3912
3913       pixel=zero;
3914       k=kernel;
3915       gamma=0.0;
3916       j=0;
3917       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
3918         {
3919           for (v=0; v < (long) width; v++)
3920           {
3921             for (u=0; u < (long) width; u++)
3922             {
3923               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
3924                 {
3925                   pixel.red+=(*k)*(p+u+j)->red;
3926                   pixel.green+=(*k)*(p+u+j)->green;
3927                   pixel.blue+=(*k)*(p+u+j)->blue;
3928                   gamma+=(*k);
3929                   k++;
3930                 }
3931             }
3932             j+=image->columns+width;
3933           }
3934           if (gamma != 0.0)
3935             {
3936               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
3937               if ((channel & RedChannel) != 0)
3938                 q->red=RoundToQuantum(gamma*pixel.red+bias);
3939               if ((channel & GreenChannel) != 0)
3940                 q->green=RoundToQuantum(gamma*pixel.green+bias);
3941               if ((channel & BlueChannel) != 0)
3942                 q->blue=RoundToQuantum(gamma*pixel.blue+bias);
3943             }
3944           if ((channel & OpacityChannel) != 0)
3945             {
3946               gamma=0.0;
3947               j=0;
3948               for (v=0; v < (long) width; v++)
3949               {
3950                 for (u=0; u < (long) width; u++)
3951                 {
3952                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
3953                     {
3954                       pixel.opacity+=(*k)*(p+u+j)->opacity;
3955                       gamma+=(*k);
3956                       k++;
3957                     }
3958                 }
3959                 j+=image->columns+width;
3960               }
3961               if (gamma != 0.0)
3962                 {
3963                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
3964                     gamma);
3965                   q->opacity=RoundToQuantum(gamma*pixel.opacity+bias);
3966                 }
3967             }
3968           if (((channel & IndexChannel) != 0) &&
3969               (image->colorspace == CMYKColorspace))
3970             {
3971               gamma=0.0;
3972               j=0;
3973               for (v=0; v < (long) width; v++)
3974               {
3975                 for (u=0; u < (long) width; u++)
3976                 {
3977                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
3978                     {
3979                       pixel.index+=(*k)*indexes[x+u+j];
3980                       gamma+=(*k);
3981                       k++;
3982                     }
3983                 }
3984                 j+=image->columns+width;
3985               }
3986               if (gamma != 0.0)
3987                 {
3988                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
3989                     gamma);
3990                   blur_indexes[x]=RoundToQuantum(gamma*pixel.index+bias);
3991                 }
3992             }
3993         }
3994       else
3995         {
3996           MagickRealType
3997             alpha;
3998
3999           for (v=0; v < (long) width; v++)
4000           {
4001             for (u=0; u < (long) width; u++)
4002             {
4003               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4004                 {
4005                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
4006                     (p+u+j)->opacity));
4007                   pixel.red+=(*k)*alpha*(p+u+j)->red;
4008                   pixel.green+=(*k)*alpha*(p+u+j)->green;
4009                   pixel.blue+=(*k)*alpha*(p+u+j)->blue;
4010                   pixel.opacity+=(*k)*(p+u+j)->opacity;
4011                   gamma+=(*k)*alpha;
4012                   k++;
4013                 }
4014             }
4015             j+=image->columns+width;
4016           }
4017           if (gamma != 0.0)
4018             {
4019               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4020               if ((channel & RedChannel) != 0)
4021                 q->red=RoundToQuantum(gamma*pixel.red+bias);
4022               if ((channel & GreenChannel) != 0)
4023                 q->green=RoundToQuantum(gamma*pixel.green+bias);
4024               if ((channel & BlueChannel) != 0)
4025                 q->blue=RoundToQuantum(gamma*pixel.blue+bias);
4026             }
4027           if ((channel & OpacityChannel) != 0)
4028             {
4029               gamma=0.0;
4030               j=0;
4031               for (v=0; v < (long) width; v++)
4032               {
4033                 for (u=0; u < (long) width; u++)
4034                 {
4035                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4036                     {
4037                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4038                       gamma+=(*k);
4039                       k++;
4040                     }
4041                 }
4042                 j+=image->columns+width;
4043               }
4044               if (gamma != 0.0)
4045                 {
4046                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4047                     gamma);
4048                   q->opacity=RoundToQuantum(pixel.opacity+bias);
4049                 }
4050             }
4051           if (((channel & IndexChannel) != 0) &&
4052               (image->colorspace == CMYKColorspace))
4053             {
4054               gamma=0.0;
4055               j=0;
4056               for (v=0; v < (long) width; v++)
4057               {
4058                 for (u=0; u < (long) width; u++)
4059                 {
4060                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4061                     {
4062                       alpha=(MagickRealType) (QuantumScale*(QuantumRange-
4063                         (p+u+j)->opacity));
4064                       pixel.index+=(*k)*alpha*indexes[x+u+j];
4065                       gamma+=(*k);
4066                       k++;
4067                     }
4068                 }
4069                 j+=image->columns+width;
4070               }
4071               if (gamma != 0.0)
4072                 {
4073                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4074                     gamma);
4075                   blur_indexes[x]=RoundToQuantum(gamma*pixel.index+bias);
4076                 }
4077             }
4078         }
4079       p++;
4080       q++;
4081     }
4082     sync=SyncCacheViewAuthenticPixels(blur_view,exception);
4083     if (sync == MagickFalse)
4084       status=MagickFalse;
4085     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4086       {
4087         MagickBooleanType
4088           proceed;
4089
4090 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4091   #pragma omp critical (MagickCore_SelectiveBlurImageChannel)
4092 #endif
4093         proceed=SetImageProgress(image,SelectiveBlurImageTag,progress++,
4094           image->rows);
4095         if (proceed == MagickFalse)
4096           status=MagickFalse;
4097       }
4098   }
4099   blur_image->type=image->type;
4100   blur_view=DestroyCacheView(blur_view);
4101   image_view=DestroyCacheView(image_view);
4102   kernel=(double *) RelinquishMagickMemory(kernel);
4103   if (status == MagickFalse)
4104     blur_image=DestroyImage(blur_image);
4105   return(blur_image);
4106 }
4107 \f
4108 /*
4109 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4110 %                                                                             %
4111 %                                                                             %
4112 %                                                                             %
4113 %     S h a d e I m a g e                                                     %
4114 %                                                                             %
4115 %                                                                             %
4116 %                                                                             %
4117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4118 %
4119 %  ShadeImage() shines a distant light on an image to create a
4120 %  three-dimensional effect. You control the positioning of the light with
4121 %  azimuth and elevation; azimuth is measured in degrees off the x axis
4122 %  and elevation is measured in pixels above the Z axis.
4123 %
4124 %  The format of the ShadeImage method is:
4125 %
4126 %      Image *ShadeImage(const Image *image,const MagickBooleanType gray,
4127 %        const double azimuth,const double elevation,ExceptionInfo *exception)
4128 %
4129 %  A description of each parameter follows:
4130 %
4131 %    o image: the image.
4132 %
4133 %    o gray: A value other than zero shades the intensity of each pixel.
4134 %
4135 %    o azimuth, elevation:  Define the light source direction.
4136 %
4137 %    o exception: return any errors or warnings in this structure.
4138 %
4139 */
4140 MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
4141   const double azimuth,const double elevation,ExceptionInfo *exception)
4142 {
4143 #define ShadeImageTag  "Shade/Image"
4144
4145   Image
4146     *shade_image;
4147
4148   long
4149     progress,
4150     y;
4151
4152   MagickBooleanType
4153     status;
4154
4155   PrimaryInfo
4156     light;
4157
4158   CacheView
4159     *image_view,
4160     *shade_view;
4161
4162   /*
4163     Initialize shaded image attributes.
4164   */
4165   assert(image != (const Image *) NULL);
4166   assert(image->signature == MagickSignature);
4167   if (image->debug != MagickFalse)
4168     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4169   assert(exception != (ExceptionInfo *) NULL);
4170   assert(exception->signature == MagickSignature);
4171   shade_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
4172   if (shade_image == (Image *) NULL)
4173     return((Image *) NULL);
4174   if (SetImageStorageClass(shade_image,DirectClass) == MagickFalse)
4175     {
4176       InheritException(exception,&shade_image->exception);
4177       shade_image=DestroyImage(shade_image);
4178       return((Image *) NULL);
4179     }
4180   /*
4181     Compute the light vector.
4182   */
4183   light.x=(double) QuantumRange*cos(DegreesToRadians(azimuth))*
4184     cos(DegreesToRadians(elevation));
4185   light.y=(double) QuantumRange*sin(DegreesToRadians(azimuth))*
4186     cos(DegreesToRadians(elevation));
4187   light.z=(double) QuantumRange*sin(DegreesToRadians(elevation));
4188   /*
4189     Shade image.
4190   */
4191   status=MagickTrue;
4192   progress=0;
4193   image_view=AcquireCacheView(image);
4194   shade_view=AcquireCacheView(shade_image);
4195 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4196   #pragma omp parallel for schedule(dynamic) shared(progress,status)
4197 #endif
4198   for (y=0; y < (long) image->rows; y++)
4199   {
4200     MagickRealType
4201       distance,
4202       normal_distance,
4203       shade;
4204
4205     PrimaryInfo
4206       normal;
4207
4208     register const PixelPacket
4209       *__restrict p,
4210       *__restrict s0,
4211       *__restrict s1,
4212       *__restrict s2;
4213
4214     register long
4215       x;
4216
4217     register PixelPacket
4218       *__restrict q;
4219
4220     if (status == MagickFalse)
4221       continue;
4222     p=GetCacheViewVirtualPixels(image_view,-1,y-1,image->columns+2,3,exception);
4223     q=QueueCacheViewAuthenticPixels(shade_view,0,y,shade_image->columns,1,
4224       exception);
4225     if ((p == (PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4226       {
4227         status=MagickFalse;
4228         continue;
4229       }
4230     /*
4231       Shade this row of pixels.
4232     */
4233     normal.z=2.0*(double) QuantumRange;  /* constant Z of surface normal */
4234     s0=p+1;
4235     s1=s0+image->columns+2;
4236     s2=s1+image->columns+2;
4237     for (x=0; x < (long) image->columns; x++)
4238     {
4239       /*
4240         Determine the surface normal and compute shading.
4241       */
4242       normal.x=(double) (PixelIntensity(s0-1)+PixelIntensity(s1-1)+
4243         PixelIntensity(s2-1)-PixelIntensity(s0+1)-PixelIntensity(s1+1)-
4244         PixelIntensity(s2+1));
4245       normal.y=(double) (PixelIntensity(s2-1)+PixelIntensity(s2)+
4246         PixelIntensity(s2+1)-PixelIntensity(s0-1)-PixelIntensity(s0)-
4247         PixelIntensity(s0+1));
4248       if ((normal.x == 0.0) && (normal.y == 0.0))
4249         shade=light.z;
4250       else
4251         {
4252           shade=0.0;
4253           distance=normal.x*light.x+normal.y*light.y+normal.z*light.z;
4254           if (distance > MagickEpsilon)
4255             {
4256               normal_distance=
4257                 normal.x*normal.x+normal.y*normal.y+normal.z*normal.z;
4258               if (normal_distance > (MagickEpsilon*MagickEpsilon))
4259                 shade=distance/sqrt((double) normal_distance);
4260             }
4261         }
4262       if (gray != MagickFalse)
4263         {
4264           q->red=(Quantum) shade;
4265           q->green=(Quantum) shade;
4266           q->blue=(Quantum) shade;
4267         }
4268       else
4269         {
4270           q->red=RoundToQuantum(QuantumScale*shade*s1->red);
4271           q->green=RoundToQuantum(QuantumScale*shade*s1->green);
4272           q->blue=RoundToQuantum(QuantumScale*shade*s1->blue);
4273         }
4274       q->opacity=s1->opacity;
4275       s0++;
4276       s1++;
4277       s2++;
4278       q++;
4279     }
4280     if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
4281       status=MagickFalse;
4282     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4283       {
4284         MagickBooleanType
4285           proceed;
4286
4287 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4288   #pragma omp critical (MagickCore_ShadeImage)
4289 #endif
4290         proceed=SetImageProgress(image,ShadeImageTag,progress++,image->rows);
4291         if (proceed == MagickFalse)
4292           status=MagickFalse;
4293       }
4294   }
4295   shade_view=DestroyCacheView(shade_view);
4296   image_view=DestroyCacheView(image_view);
4297   if (status == MagickFalse)
4298     shade_image=DestroyImage(shade_image);
4299   return(shade_image);
4300 }
4301 \f
4302 /*
4303 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4304 %                                                                             %
4305 %                                                                             %
4306 %                                                                             %
4307 %     S h a r p e n I m a g e                                                 %
4308 %                                                                             %
4309 %                                                                             %
4310 %                                                                             %
4311 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4312 %
4313 %  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
4314 %  operator of the given radius and standard deviation (sigma).  For
4315 %  reasonable results, radius should be larger than sigma.  Use a radius of 0
4316 %  and SharpenImage() selects a suitable radius for you.
4317 %
4318 %  Using a separable kernel would be faster, but the negative weights cancel
4319 %  out on the corners of the kernel producing often undesirable ringing in the
4320 %  filtered result; this can be avoided by using a 2D gaussian shaped image
4321 %  sharpening kernel instead.
4322 %
4323 %  The format of the SharpenImage method is:
4324 %
4325 %    Image *SharpenImage(const Image *image,const double radius,
4326 %      const double sigma,ExceptionInfo *exception)
4327 %    Image *SharpenImageChannel(const Image *image,const ChannelType channel,
4328 %      const double radius,const double sigma,ExceptionInfo *exception)
4329 %
4330 %  A description of each parameter follows:
4331 %
4332 %    o image: the image.
4333 %
4334 %    o channel: the channel type.
4335 %
4336 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4337 %      pixel.
4338 %
4339 %    o sigma: the standard deviation of the Laplacian, in pixels.
4340 %
4341 %    o exception: return any errors or warnings in this structure.
4342 %
4343 */
4344
4345 MagickExport Image *SharpenImage(const Image *image,const double radius,
4346   const double sigma,ExceptionInfo *exception)
4347 {
4348   Image
4349     *sharp_image;
4350
4351   sharp_image=SharpenImageChannel(image,DefaultChannels,radius,sigma,exception);
4352   return(sharp_image);
4353 }
4354
4355 MagickExport Image *SharpenImageChannel(const Image *image,
4356   const ChannelType channel,const double radius,const double sigma,
4357   ExceptionInfo *exception)
4358 {
4359   double
4360     *kernel;
4361
4362   Image
4363     *sharp_image;
4364
4365   MagickRealType
4366     alpha,
4367     normalize;
4368
4369   register long
4370     i,
4371     u,
4372     v;
4373
4374   unsigned long
4375     width;
4376
4377   assert(image != (const Image *) NULL);
4378   assert(image->signature == MagickSignature);
4379   if (image->debug != MagickFalse)
4380     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4381   assert(exception != (ExceptionInfo *) NULL);
4382   assert(exception->signature == MagickSignature);
4383   width=GetOptimalKernelWidth2D(radius,sigma);
4384   kernel=(double *) AcquireQuantumMemory((size_t) width*width,sizeof(*kernel));
4385   if (kernel == (double *) NULL)
4386     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4387   i=0;
4388   normalize=0.0;
4389   for (v=(-((long) width/2)); v <= (long) (width/2); v++)
4390   {
4391     for (u=(-((long) width/2)); u <= (long) (width/2); u++)
4392     {
4393       alpha=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma));
4394       kernel[i]=(double) (-alpha/(2.0*MagickPI*MagickSigma*MagickSigma));
4395       normalize+=kernel[i];
4396       i++;
4397     }
4398   }
4399   kernel[i/2]=(double) ((-2.0)*normalize);
4400   sharp_image=ConvolveImageChannel(image,channel,width,kernel,exception);
4401   kernel=(double *) RelinquishMagickMemory(kernel);
4402   return(sharp_image);
4403 }
4404 \f
4405 /*
4406 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4407 %                                                                             %
4408 %                                                                             %
4409 %                                                                             %
4410 %     S p r e a d I m a g e                                                   %
4411 %                                                                             %
4412 %                                                                             %
4413 %                                                                             %
4414 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4415 %
4416 %  SpreadImage() is a special effects method that randomly displaces each
4417 %  pixel in a block defined by the radius parameter.
4418 %
4419 %  The format of the SpreadImage method is:
4420 %
4421 %      Image *SpreadImage(const Image *image,const double radius,
4422 %        ExceptionInfo *exception)
4423 %
4424 %  A description of each parameter follows:
4425 %
4426 %    o image: the image.
4427 %
4428 %    o radius:  Choose a random pixel in a neighborhood of this extent.
4429 %
4430 %    o exception: return any errors or warnings in this structure.
4431 %
4432 */
4433 MagickExport Image *SpreadImage(const Image *image,const double radius,
4434   ExceptionInfo *exception)
4435 {
4436 #define SpreadImageTag  "Spread/Image"
4437
4438   Image
4439     *spread_image;
4440
4441   long
4442     progress,
4443     y;
4444
4445   MagickBooleanType
4446     status;
4447
4448   MagickPixelPacket
4449     zero;
4450
4451   RandomInfo
4452     **random_info;
4453
4454   ResampleFilter
4455     **resample_filter;
4456
4457   unsigned long
4458     width;
4459
4460   CacheView
4461     *image_view;
4462
4463   /*
4464     Initialize spread image attributes.
4465   */
4466   assert(image != (Image *) NULL);
4467   assert(image->signature == MagickSignature);
4468   if (image->debug != MagickFalse)
4469     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4470   assert(exception != (ExceptionInfo *) NULL);
4471   assert(exception->signature == MagickSignature);
4472   spread_image=CloneImage(image,image->columns,image->rows,MagickTrue,
4473     exception);
4474   if (spread_image == (Image *) NULL)
4475     return((Image *) NULL);
4476   if (SetImageStorageClass(spread_image,DirectClass) == MagickFalse)
4477     {
4478       InheritException(exception,&spread_image->exception);
4479       spread_image=DestroyImage(spread_image);
4480       return((Image *) NULL);
4481     }
4482   /*
4483     Spread image.
4484   */
4485   status=MagickTrue;
4486   progress=0;
4487   GetMagickPixelPacket(spread_image,&zero);
4488   width=GetOptimalKernelWidth1D(radius,0.5);
4489   resample_filter=AcquireResampleFilterThreadSet(image,MagickTrue,exception);
4490   random_info=AcquireRandomInfoThreadSet();
4491   image_view=AcquireCacheView(spread_image);
4492 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4493   #pragma omp parallel for schedule(dynamic,8) shared(progress,status)
4494 #endif
4495   for (y=0; y < (long) spread_image->rows; y++)
4496   {
4497     MagickPixelPacket
4498       pixel;
4499
4500     register IndexPacket
4501       *__restrict indexes;
4502
4503     register long
4504       id,
4505       x;
4506
4507     register PixelPacket
4508       *__restrict q;
4509
4510     if (status == MagickFalse)
4511       continue;
4512     q=QueueCacheViewAuthenticPixels(image_view,0,y,spread_image->columns,1,
4513       exception);
4514     if (q == (PixelPacket *) NULL)
4515       {
4516         status=MagickFalse;
4517         continue;
4518       }
4519     indexes=GetCacheViewAuthenticIndexQueue(image_view);
4520     pixel=zero;
4521     id=GetOpenMPThreadId();
4522     for (x=0; x < (long) spread_image->columns; x++)
4523     {
4524       (void) ResamplePixelColor(resample_filter[id],(double) x+width*
4525         (GetPseudoRandomValue(random_info[id])-0.5),(double) y+width*
4526         (GetPseudoRandomValue(random_info[id])-0.5),&pixel);
4527       SetPixelPacket(spread_image,&pixel,q,indexes+x);
4528       q++;
4529     }
4530     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
4531       status=MagickFalse;
4532     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4533       {
4534         MagickBooleanType
4535           proceed;
4536
4537 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4538   #pragma omp critical (MagickCore_SpreadImage)
4539 #endif
4540         proceed=SetImageProgress(image,SpreadImageTag,progress++,image->rows);
4541         if (proceed == MagickFalse)
4542           status=MagickFalse;
4543       }
4544   }
4545   image_view=DestroyCacheView(image_view);
4546   random_info=DestroyRandomInfoThreadSet(random_info);
4547   resample_filter=DestroyResampleFilterThreadSet(resample_filter);
4548   return(spread_image);
4549 }
4550 \f
4551 /*
4552 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4553 %                                                                             %
4554 %                                                                             %
4555 %                                                                             %
4556 %     U n s h a r p M a s k I m a g e                                         %
4557 %                                                                             %
4558 %                                                                             %
4559 %                                                                             %
4560 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4561 %
4562 %  UnsharpMaskImage() sharpens one or more image channels.  We convolve the
4563 %  image with a Gaussian operator of the given radius and standard deviation
4564 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
4565 %  radius of 0 and UnsharpMaskImage() selects a suitable radius for you.
4566 %
4567 %  The format of the UnsharpMaskImage method is:
4568 %
4569 %    Image *UnsharpMaskImage(const Image *image,const double radius,
4570 %      const double sigma,const double amount,const double threshold,
4571 %      ExceptionInfo *exception)
4572 %    Image *UnsharpMaskImageChannel(const Image *image,
4573 %      const ChannelType channel,const double radius,const double sigma,
4574 %      const double amount,const double threshold,ExceptionInfo *exception)
4575 %
4576 %  A description of each parameter follows:
4577 %
4578 %    o image: the image.
4579 %
4580 %    o channel: the channel type.
4581 %
4582 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4583 %      pixel.
4584 %
4585 %    o sigma: the standard deviation of the Gaussian, in pixels.
4586 %
4587 %    o amount: the percentage of the difference between the original and the
4588 %      blur image that is added back into the original.
4589 %
4590 %    o threshold: the threshold in pixels needed to apply the diffence amount.
4591 %
4592 %    o exception: return any errors or warnings in this structure.
4593 %
4594 */
4595
4596 MagickExport Image *UnsharpMaskImage(const Image *image,const double radius,
4597   const double sigma,const double amount,const double threshold,
4598   ExceptionInfo *exception)
4599 {
4600   Image
4601     *sharp_image;
4602
4603   sharp_image=UnsharpMaskImageChannel(image,DefaultChannels,radius,sigma,amount,
4604     threshold,exception);
4605   return(sharp_image);
4606 }
4607
4608 MagickExport Image *UnsharpMaskImageChannel(const Image *image,
4609   const ChannelType channel,const double radius,const double sigma,
4610   const double amount,const double threshold,ExceptionInfo *exception)
4611 {
4612 #define SharpenImageTag  "Sharpen/Image"
4613
4614   Image
4615     *unsharp_image;
4616
4617   long
4618     progress,
4619     y;
4620
4621   MagickBooleanType
4622     status;
4623
4624   MagickPixelPacket
4625     zero;
4626
4627   MagickRealType
4628     quantum_threshold;
4629
4630   CacheView
4631     *image_view,
4632     *unsharp_view;
4633
4634   assert(image != (const Image *) NULL);
4635   assert(image->signature == MagickSignature);
4636   if (image->debug != MagickFalse)
4637     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4638   assert(exception != (ExceptionInfo *) NULL);
4639   unsharp_image=BlurImageChannel(image,channel,radius,sigma,exception);
4640   if (unsharp_image == (Image *) NULL)
4641     return((Image *) NULL);
4642   quantum_threshold=(MagickRealType) QuantumRange*threshold;
4643   /*
4644     Unsharp-mask image.
4645   */
4646   status=MagickTrue;
4647   progress=0;
4648   GetMagickPixelPacket(image,&zero);
4649   image_view=AcquireCacheView(image);
4650   unsharp_view=AcquireCacheView(unsharp_image);
4651 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4652   #pragma omp parallel for schedule(dynamic) shared(progress,status)
4653 #endif
4654   for (y=0; y < (long) image->rows; y++)
4655   {
4656     MagickPixelPacket
4657       pixel;
4658
4659     register const IndexPacket
4660       *__restrict indexes;
4661
4662     register const PixelPacket
4663       *__restrict p;
4664
4665     register IndexPacket
4666       *__restrict unsharp_indexes;
4667
4668     register long
4669       x;
4670
4671     register PixelPacket
4672       *__restrict q;
4673
4674     if (status == MagickFalse)
4675       continue;
4676     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
4677     q=GetCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
4678       exception);
4679     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4680       {
4681         status=MagickFalse;
4682         continue;
4683       }
4684     indexes=GetCacheViewVirtualIndexQueue(image_view);
4685     unsharp_indexes=GetCacheViewAuthenticIndexQueue(unsharp_view);
4686     pixel=zero;
4687     for (x=0; x < (long) image->columns; x++)
4688     {
4689       if ((channel & RedChannel) != 0)
4690         {
4691           pixel.red=p->red-(MagickRealType) q->red;
4692           if (fabs(2.0*pixel.red) < quantum_threshold)
4693             pixel.red=(MagickRealType) p->red;
4694           else
4695             pixel.red=(MagickRealType) p->red+(pixel.red*amount);
4696           q->red=RoundToQuantum(pixel.red);
4697         }
4698       if ((channel & GreenChannel) != 0)
4699         {
4700           pixel.green=p->green-(MagickRealType) q->green;
4701           if (fabs(2.0*pixel.green) < quantum_threshold)
4702             pixel.green=(MagickRealType) p->green;
4703           else
4704             pixel.green=(MagickRealType) p->green+(pixel.green*amount);
4705           q->green=RoundToQuantum(pixel.green);
4706         }
4707       if ((channel & BlueChannel) != 0)
4708         {
4709           pixel.blue=p->blue-(MagickRealType) q->blue;
4710           if (fabs(2.0*pixel.blue) < quantum_threshold)
4711             pixel.blue=(MagickRealType) p->blue;
4712           else
4713             pixel.blue=(MagickRealType) p->blue+(pixel.blue*amount);
4714           q->blue=RoundToQuantum(pixel.blue);
4715         }
4716       if ((channel & OpacityChannel) != 0)
4717         {
4718           pixel.opacity=p->opacity-(MagickRealType) q->opacity;
4719           if (fabs(2.0*pixel.opacity) < quantum_threshold)
4720             pixel.opacity=(MagickRealType) p->opacity;
4721           else
4722             pixel.opacity=p->opacity+(pixel.opacity*amount);
4723           q->opacity=RoundToQuantum(pixel.opacity);
4724         }
4725       if (((channel & IndexChannel) != 0) &&
4726           (image->colorspace == CMYKColorspace))
4727         {
4728           pixel.index=unsharp_indexes[x]-(MagickRealType) indexes[x];
4729           if (fabs(2.0*pixel.index) < quantum_threshold)
4730             pixel.index=(MagickRealType) unsharp_indexes[x];
4731           else
4732             pixel.index=(MagickRealType) unsharp_indexes[x]+(pixel.index*
4733               amount);
4734           unsharp_indexes[x]=RoundToQuantum(pixel.index);
4735         }
4736       p++;
4737       q++;
4738     }
4739     if (SyncCacheViewAuthenticPixels(unsharp_view,exception) == MagickFalse)
4740       status=MagickFalse;
4741     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4742       {
4743         MagickBooleanType
4744           proceed;
4745
4746 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4747   #pragma omp critical (MagickCore_UnsharpMaskImageChannel)
4748 #endif
4749         proceed=SetImageProgress(image,SharpenImageTag,progress++,image->rows);
4750         if (proceed == MagickFalse)
4751           status=MagickFalse;
4752       }
4753   }
4754   unsharp_image->type=image->type;
4755   unsharp_view=DestroyCacheView(unsharp_view);
4756   image_view=DestroyCacheView(image_view);
4757   if (status == MagickFalse)
4758     unsharp_image=DestroyImage(unsharp_image);
4759   return(unsharp_image);
4760 }