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