granicus.if.org Git - imagemagick/blob - magick/effect.c

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