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