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

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