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