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

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %                   EEEEE  FFFFF  FFFFF  EEEEE  CCCC  TTTTT                   %
   7 %                   E      F      F      E     C        T                     %
   8 %                   EEE    FFF    FFF    EEE   C        T                     %
   9 %                   E      F      F      E     C        T                     %
  10 %                   EEEEE  F      F      EEEEE  CCCC    T                     %
  11 %                                                                             %
  12 %                                                                             %
  13 %                       MagickCore Image Effects Methods                      %
  14 %                                                                             %
  15 %                               Software Design                               %
  16 %                                 John Cristy                                 %
  17 %                                 October 1996                                %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "magick/studio.h"
  44 #include "magick/accelerate.h"
  45 #include "magick/blob.h"
  46 #include "magick/cache-view.h"
  47 #include "magick/color.h"
  48 #include "magick/color-private.h"
  49 #include "magick/colorspace.h"
  50 #include "magick/constitute.h"
  51 #include "magick/decorate.h"
  52 #include "magick/draw.h"
  53 #include "magick/enhance.h"
  54 #include "magick/exception.h"
  55 #include "magick/exception-private.h"
  56 #include "magick/effect.h"
  57 #include "magick/fx.h"
  58 #include "magick/gem.h"
  59 #include "magick/geometry.h"
  60 #include "magick/image-private.h"
  61 #include "magick/list.h"
  62 #include "magick/log.h"
  63 #include "magick/memory_.h"
  64 #include "magick/monitor.h"
  65 #include "magick/monitor-private.h"
  66 #include "magick/montage.h"
  67 #include "magick/morphology.h"
  68 #include "magick/paint.h"
  69 #include "magick/pixel-private.h"
  70 #include "magick/property.h"
  71 #include "magick/quantize.h"
  72 #include "magick/quantum.h"
  73 #include "magick/random_.h"
  74 #include "magick/random-private.h"
  75 #include "magick/resample.h"
  76 #include "magick/resample-private.h"
  77 #include "magick/resize.h"
  78 #include "magick/resource_.h"
  79 #include "magick/segment.h"
  80 #include "magick/shear.h"
  81 #include "magick/signature-private.h"
  82 #include "magick/string_.h"
  83 #include "magick/thread-private.h"
  84 #include "magick/transform.h"
  85 #include "magick/threshold.h"
  86 \f
  87 /*
  88 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  89 %                                                                             %
  90 %                                                                             %
  91 %                                                                             %
  92 %     A d a p t i v e B l u r I m a g e                                       %
  93 %                                                                             %
  94 %                                                                             %
  95 %                                                                             %
  96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  97 %
  98 %  AdaptiveBlurImage() adaptively blurs the image by blurring less
  99 %  intensely near image edges and more intensely far from edges.  We blur the
 100 %  image with a Gaussian operator of the given radius and standard deviation
 101 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
 102 %  radius of 0 and AdaptiveBlurImage() selects a suitable radius for you.
 103 %
 104 %  The format of the AdaptiveBlurImage method is:
 105 %
 106 %      Image *AdaptiveBlurImage(const Image *image,const double radius,
 107 %        const double sigma,ExceptionInfo *exception)
 108 %      Image *AdaptiveBlurImageChannel(const Image *image,
 109 %        const ChannelType channel,double radius,const double sigma,
 110 %        ExceptionInfo *exception)
 111 %
 112 %  A description of each parameter follows:
 113 %
 114 %    o image: the image.
 115 %
 116 %    o channel: the channel type.
 117 %
 118 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 119 %      pixel.
 120 %
 121 %    o sigma: the standard deviation of the Laplacian, in pixels.
 122 %
 123 %    o exception: return any errors or warnings in this structure.
 124 %
 125 */
 126
 127 MagickExport Image *AdaptiveBlurImage(const Image *image,const double radius,
 128   const double sigma,ExceptionInfo *exception)
 129 {
 130   Image
 131     *blur_image;
 132
 133   blur_image=AdaptiveBlurImageChannel(image,DefaultChannels,radius,sigma,
 134     exception);
 135   return(blur_image);
 136 }
 137
 138 MagickExport Image *AdaptiveBlurImageChannel(const Image *image,
 139   const ChannelType channel,const double radius,const double sigma,
 140   ExceptionInfo *exception)
 141 {
 142 #define AdaptiveBlurImageTag  "Convolve/Image"
 143 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
 144
 145   CacheView
 146     *blur_view,
 147     *edge_view,
 148     *image_view;
 149
 150   double
 151     **kernel,
 152     normalize;
 153
 154   Image
 155     *blur_image,
 156     *edge_image,
 157     *gaussian_image;
 158
 159   MagickBooleanType
 160     status;
 161
 162   MagickOffsetType
 163     progress;
 164
 165   MagickPixelPacket
 166     bias;
 167
 168   register ssize_t
 169     i;
 170
 171   size_t
 172     width;
 173
 174   ssize_t
 175     j,
 176     k,
 177     u,
 178     v,
 179     y;
 180
 181   assert(image != (const Image *) NULL);
 182   assert(image->signature == MagickSignature);
 183   if (image->debug != MagickFalse)
 184     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 185   assert(exception != (ExceptionInfo *) NULL);
 186   assert(exception->signature == MagickSignature);
 187   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
 188   if (blur_image == (Image *) NULL)
 189     return((Image *) NULL);
 190   if (fabs(sigma) <= MagickEpsilon)
 191     return(blur_image);
 192   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
 193     {
 194       InheritException(exception,&blur_image->exception);
 195       blur_image=DestroyImage(blur_image);
 196       return((Image *) NULL);
 197     }
 198   /*
 199     Edge detect the image brighness channel, level, blur, and level again.
 200   */
 201   edge_image=EdgeImage(image,radius,exception);
 202   if (edge_image == (Image *) NULL)
 203     {
 204       blur_image=DestroyImage(blur_image);
 205       return((Image *) NULL);
 206     }
 207   (void) LevelImage(edge_image,"20%,95%");
 208   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
 209   if (gaussian_image != (Image *) NULL)
 210     {
 211       edge_image=DestroyImage(edge_image);
 212       edge_image=gaussian_image;
 213     }
 214   (void) LevelImage(edge_image,"10%,95%");
 215   /*
 216     Create a set of kernels from maximum (radius,sigma) to minimum.
 217   */
 218   width=GetOptimalKernelWidth2D(radius,sigma);
 219   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 220   if (kernel == (double **) NULL)
 221     {
 222       edge_image=DestroyImage(edge_image);
 223       blur_image=DestroyImage(blur_image);
 224       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 225     }
 226   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
 227   for (i=0; i < (ssize_t) width; i+=2)
 228   {
 229     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
 230       sizeof(**kernel));
 231     if (kernel[i] == (double *) NULL)
 232       break;
 233     normalize=0.0;
 234     j=(ssize_t) (width-i)/2;
 235     k=0;
 236     for (v=(-j); v <= j; v++)
 237     {
 238       for (u=(-j); u <= j; u++)
 239       {
 240         kernel[i][k]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
 241           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
 242         normalize+=kernel[i][k];
 243         k++;
 244       }
 245     }
 246     if (fabs(normalize) <= MagickEpsilon)
 247       normalize=1.0;
 248     normalize=1.0/normalize;
 249     for (k=0; k < (j*j); k++)
 250       kernel[i][k]=normalize*kernel[i][k];
 251   }
 252   if (i < (ssize_t) width)
 253     {
 254       for (i-=2; i >= 0; i-=2)
 255         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 256       kernel=(double **) RelinquishMagickMemory(kernel);
 257       edge_image=DestroyImage(edge_image);
 258       blur_image=DestroyImage(blur_image);
 259       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 260     }
 261   /*
 262     Adaptively blur image.
 263   */
 264   status=MagickTrue;
 265   progress=0;
 266   GetMagickPixelPacket(image,&bias);
 267   SetMagickPixelPacketBias(image,&bias);
 268   image_view=AcquireCacheView(image);
 269   edge_view=AcquireCacheView(edge_image);
 270   blur_view=AcquireCacheView(blur_image);
 271 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 272   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 273 #endif
 274   for (y=0; y < (ssize_t) blur_image->rows; y++)
 275   {
 276     register const IndexPacket
 277       *restrict indexes;
 278
 279     register const PixelPacket
 280       *restrict p,
 281       *restrict r;
 282
 283     register IndexPacket
 284       *restrict blur_indexes;
 285
 286     register PixelPacket
 287       *restrict q;
 288
 289     register ssize_t
 290       x;
 291
 292     if (status == MagickFalse)
 293       continue;
 294     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
 295     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
 296       exception);
 297     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 298       {
 299         status=MagickFalse;
 300         continue;
 301       }
 302     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
 303     for (x=0; x < (ssize_t) blur_image->columns; x++)
 304     {
 305       MagickPixelPacket
 306         pixel;
 307
 308       MagickRealType
 309         alpha,
 310         gamma;
 311
 312       register const double
 313         *restrict k;
 314
 315       register ssize_t
 316         i,
 317         u,
 318         v;
 319
 320       gamma=0.0;
 321       i=(ssize_t) ceil((double) width*QuantumScale*PixelIntensity(r)-0.5);
 322       if (i < 0)
 323         i=0;
 324       else
 325         if (i > (ssize_t) width)
 326           i=(ssize_t) width;
 327       if ((i & 0x01) != 0)
 328         i--;
 329       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-i)/2L),y-
 330         (ssize_t) ((width-i)/2L),width-i,width-i,exception);
 331       if (p == (const PixelPacket *) NULL)
 332         break;
 333       indexes=GetCacheViewVirtualIndexQueue(image_view);
 334       pixel=bias;
 335       k=kernel[i];
 336       for (v=0; v < (ssize_t) (width-i); v++)
 337       {
 338         for (u=0; u < (ssize_t) (width-i); u++)
 339         {
 340           alpha=1.0;
 341           if (((channel & OpacityChannel) != 0) &&
 342               (image->matte != MagickFalse))
 343             alpha=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(p));
 344           if ((channel & RedChannel) != 0)
 345             pixel.red+=(*k)*alpha*GetRedPixelComponent(p);
 346           if ((channel & GreenChannel) != 0)
 347             pixel.green+=(*k)*alpha*GetGreenPixelComponent(p);
 348           if ((channel & BlueChannel) != 0)
 349             pixel.blue+=(*k)*alpha*GetBluePixelComponent(p);
 350           if ((channel & OpacityChannel) != 0)
 351             pixel.opacity+=(*k)*GetOpacityPixelComponent(p);
 352           if (((channel & IndexChannel) != 0) &&
 353               (image->colorspace == CMYKColorspace))
 354             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
 355           gamma+=(*k)*alpha;
 356           k++;
 357           p++;
 358         }
 359       }
 360       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 361       if ((channel & RedChannel) != 0)
 362         q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
 363       if ((channel & GreenChannel) != 0)
 364         q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
 365       if ((channel & BlueChannel) != 0)
 366         q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
 367       if ((channel & OpacityChannel) != 0)
 368         SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
 369       if (((channel & IndexChannel) != 0) &&
 370           (image->colorspace == CMYKColorspace))
 371         blur_indexes[x]=ClampToQuantum(gamma*GetIndexPixelComponent(&pixel));
 372       q++;
 373       r++;
 374     }
 375     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
 376       status=MagickFalse;
 377     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 378       {
 379         MagickBooleanType
 380           proceed;
 381
 382 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 383   #pragma omp critical (MagickCore_AdaptiveBlurImageChannel)
 384 #endif
 385         proceed=SetImageProgress(image,AdaptiveBlurImageTag,progress++,
 386           image->rows);
 387         if (proceed == MagickFalse)
 388           status=MagickFalse;
 389       }
 390   }
 391   blur_image->type=image->type;
 392   blur_view=DestroyCacheView(blur_view);
 393   edge_view=DestroyCacheView(edge_view);
 394   image_view=DestroyCacheView(image_view);
 395   edge_image=DestroyImage(edge_image);
 396   for (i=0; i < (ssize_t) width;  i+=2)
 397     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 398   kernel=(double **) RelinquishMagickMemory(kernel);
 399   if (status == MagickFalse)
 400     blur_image=DestroyImage(blur_image);
 401   return(blur_image);
 402 }
 403 \f
 404 /*
 405 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 406 %                                                                             %
 407 %                                                                             %
 408 %                                                                             %
 409 %     A d a p t i v e S h a r p e n I m a g e                                 %
 410 %                                                                             %
 411 %                                                                             %
 412 %                                                                             %
 413 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 414 %
 415 %  AdaptiveSharpenImage() adaptively sharpens the image by sharpening more
 416 %  intensely near image edges and less intensely far from edges. We sharpen the
 417 %  image with a Gaussian operator of the given radius and standard deviation
 418 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
 419 %  radius of 0 and AdaptiveSharpenImage() selects a suitable radius for you.
 420 %
 421 %  The format of the AdaptiveSharpenImage method is:
 422 %
 423 %      Image *AdaptiveSharpenImage(const Image *image,const double radius,
 424 %        const double sigma,ExceptionInfo *exception)
 425 %      Image *AdaptiveSharpenImageChannel(const Image *image,
 426 %        const ChannelType channel,double radius,const double sigma,
 427 %        ExceptionInfo *exception)
 428 %
 429 %  A description of each parameter follows:
 430 %
 431 %    o image: the image.
 432 %
 433 %    o channel: the channel type.
 434 %
 435 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 436 %      pixel.
 437 %
 438 %    o sigma: the standard deviation of the Laplacian, in pixels.
 439 %
 440 %    o exception: return any errors or warnings in this structure.
 441 %
 442 */
 443
 444 MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
 445   const double sigma,ExceptionInfo *exception)
 446 {
 447   Image
 448     *sharp_image;
 449
 450   sharp_image=AdaptiveSharpenImageChannel(image,DefaultChannels,radius,sigma,
 451     exception);
 452   return(sharp_image);
 453 }
 454
 455 MagickExport Image *AdaptiveSharpenImageChannel(const Image *image,
 456   const ChannelType channel,const double radius,const double sigma,
 457   ExceptionInfo *exception)
 458 {
 459 #define AdaptiveSharpenImageTag  "Convolve/Image"
 460 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
 461
 462   CacheView
 463     *sharp_view,
 464     *edge_view,
 465     *image_view;
 466
 467   double
 468     **kernel,
 469     normalize;
 470
 471   Image
 472     *sharp_image,
 473     *edge_image,
 474     *gaussian_image;
 475
 476   MagickBooleanType
 477     status;
 478
 479   MagickOffsetType
 480     progress;
 481
 482   MagickPixelPacket
 483     bias;
 484
 485   register ssize_t
 486     i;
 487
 488   size_t
 489     width;
 490
 491   ssize_t
 492     j,
 493     k,
 494     u,
 495     v,
 496     y;
 497
 498   assert(image != (const Image *) NULL);
 499   assert(image->signature == MagickSignature);
 500   if (image->debug != MagickFalse)
 501     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 502   assert(exception != (ExceptionInfo *) NULL);
 503   assert(exception->signature == MagickSignature);
 504   sharp_image=CloneImage(image,0,0,MagickTrue,exception);
 505   if (sharp_image == (Image *) NULL)
 506     return((Image *) NULL);
 507   if (fabs(sigma) <= MagickEpsilon)
 508     return(sharp_image);
 509   if (SetImageStorageClass(sharp_image,DirectClass) == MagickFalse)
 510     {
 511       InheritException(exception,&sharp_image->exception);
 512       sharp_image=DestroyImage(sharp_image);
 513       return((Image *) NULL);
 514     }
 515   /*
 516     Edge detect the image brighness channel, level, sharp, and level again.
 517   */
 518   edge_image=EdgeImage(image,radius,exception);
 519   if (edge_image == (Image *) NULL)
 520     {
 521       sharp_image=DestroyImage(sharp_image);
 522       return((Image *) NULL);
 523     }
 524   (void) LevelImage(edge_image,"20%,95%");
 525   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
 526   if (gaussian_image != (Image *) NULL)
 527     {
 528       edge_image=DestroyImage(edge_image);
 529       edge_image=gaussian_image;
 530     }
 531   (void) LevelImage(edge_image,"10%,95%");
 532   /*
 533     Create a set of kernels from maximum (radius,sigma) to minimum.
 534   */
 535   width=GetOptimalKernelWidth2D(radius,sigma);
 536   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 537   if (kernel == (double **) NULL)
 538     {
 539       edge_image=DestroyImage(edge_image);
 540       sharp_image=DestroyImage(sharp_image);
 541       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 542     }
 543   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
 544   for (i=0; i < (ssize_t) width; i+=2)
 545   {
 546     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
 547       sizeof(**kernel));
 548     if (kernel[i] == (double *) NULL)
 549       break;
 550     normalize=0.0;
 551     j=(ssize_t) (width-i)/2;
 552     k=0;
 553     for (v=(-j); v <= j; v++)
 554     {
 555       for (u=(-j); u <= j; u++)
 556       {
 557         kernel[i][k]=(double) (-exp(-((double) u*u+v*v)/(2.0*MagickSigma*
 558           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
 559         normalize+=kernel[i][k];
 560         k++;
 561       }
 562     }
 563     if (fabs(normalize) <= MagickEpsilon)
 564       normalize=1.0;
 565     normalize=1.0/normalize;
 566     for (k=0; k < (j*j); k++)
 567       kernel[i][k]=normalize*kernel[i][k];
 568   }
 569   if (i < (ssize_t) width)
 570     {
 571       for (i-=2; i >= 0; i-=2)
 572         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 573       kernel=(double **) RelinquishMagickMemory(kernel);
 574       edge_image=DestroyImage(edge_image);
 575       sharp_image=DestroyImage(sharp_image);
 576       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 577     }
 578   /*
 579     Adaptively sharpen image.
 580   */
 581   status=MagickTrue;
 582   progress=0;
 583   GetMagickPixelPacket(image,&bias);
 584   SetMagickPixelPacketBias(image,&bias);
 585   image_view=AcquireCacheView(image);
 586   edge_view=AcquireCacheView(edge_image);
 587   sharp_view=AcquireCacheView(sharp_image);
 588 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 589   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 590 #endif
 591   for (y=0; y < (ssize_t) sharp_image->rows; y++)
 592   {
 593     register const IndexPacket
 594       *restrict indexes;
 595
 596     register const PixelPacket
 597       *restrict p,
 598       *restrict r;
 599
 600     register IndexPacket
 601       *restrict sharp_indexes;
 602
 603     register PixelPacket
 604       *restrict q;
 605
 606     register ssize_t
 607       x;
 608
 609     if (status == MagickFalse)
 610       continue;
 611     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
 612     q=QueueCacheViewAuthenticPixels(sharp_view,0,y,sharp_image->columns,1,
 613       exception);
 614     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 615       {
 616         status=MagickFalse;
 617         continue;
 618       }
 619     sharp_indexes=GetCacheViewAuthenticIndexQueue(sharp_view);
 620     for (x=0; x < (ssize_t) sharp_image->columns; x++)
 621     {
 622       MagickPixelPacket
 623         pixel;
 624
 625       MagickRealType
 626         alpha,
 627         gamma;
 628
 629       register const double
 630         *restrict k;
 631
 632       register ssize_t
 633         i,
 634         u,
 635         v;
 636
 637       gamma=0.0;
 638       i=(ssize_t) ceil((double) width*(QuantumRange-QuantumScale*
 639         PixelIntensity(r))-0.5);
 640       if (i < 0)
 641         i=0;
 642       else
 643         if (i > (ssize_t) width)
 644           i=(ssize_t) width;
 645       if ((i & 0x01) != 0)
 646         i--;
 647       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-i)/2L),y-
 648         (ssize_t) ((width-i)/2L),width-i,width-i,exception);
 649       if (p == (const PixelPacket *) NULL)
 650         break;
 651       indexes=GetCacheViewVirtualIndexQueue(image_view);
 652       k=kernel[i];
 653       pixel=bias;
 654       for (v=0; v < (ssize_t) (width-i); v++)
 655       {
 656         for (u=0; u < (ssize_t) (width-i); u++)
 657         {
 658           alpha=1.0;
 659           if (((channel & OpacityChannel) != 0) &&
 660               (image->matte != MagickFalse))
 661             alpha=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(p));
 662           if ((channel & RedChannel) != 0)
 663             pixel.red+=(*k)*alpha*GetRedPixelComponent(p);
 664           if ((channel & GreenChannel) != 0)
 665             pixel.green+=(*k)*alpha*GetGreenPixelComponent(p);
 666           if ((channel & BlueChannel) != 0)
 667             pixel.blue+=(*k)*alpha*GetBluePixelComponent(p);
 668           if ((channel & OpacityChannel) != 0)
 669             pixel.opacity+=(*k)*GetOpacityPixelComponent(p);
 670           if (((channel & IndexChannel) != 0) &&
 671               (image->colorspace == CMYKColorspace))
 672             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
 673           gamma+=(*k)*alpha;
 674           k++;
 675           p++;
 676         }
 677       }
 678       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 679       if ((channel & RedChannel) != 0)
 680         q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
 681       if ((channel & GreenChannel) != 0)
 682         q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
 683       if ((channel & BlueChannel) != 0)
 684         q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
 685       if ((channel & OpacityChannel) != 0)
 686         SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
 687       if (((channel & IndexChannel) != 0) &&
 688           (image->colorspace == CMYKColorspace))
 689         sharp_indexes[x]=ClampToQuantum(gamma*GetIndexPixelComponent(&pixel));
 690       q++;
 691       r++;
 692     }
 693     if (SyncCacheViewAuthenticPixels(sharp_view,exception) == MagickFalse)
 694       status=MagickFalse;
 695     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 696       {
 697         MagickBooleanType
 698           proceed;
 699
 700 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 701   #pragma omp critical (MagickCore_AdaptiveSharpenImageChannel)
 702 #endif
 703         proceed=SetImageProgress(image,AdaptiveSharpenImageTag,progress++,
 704           image->rows);
 705         if (proceed == MagickFalse)
 706           status=MagickFalse;
 707       }
 708   }
 709   sharp_image->type=image->type;
 710   sharp_view=DestroyCacheView(sharp_view);
 711   edge_view=DestroyCacheView(edge_view);
 712   image_view=DestroyCacheView(image_view);
 713   edge_image=DestroyImage(edge_image);
 714   for (i=0; i < (ssize_t) width;  i+=2)
 715     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 716   kernel=(double **) RelinquishMagickMemory(kernel);
 717   if (status == MagickFalse)
 718     sharp_image=DestroyImage(sharp_image);
 719   return(sharp_image);
 720 }
 721 \f
 722 /*
 723 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 724 %                                                                             %
 725 %                                                                             %
 726 %                                                                             %
 727 %     B l u r I m a g e                                                       %
 728 %                                                                             %
 729 %                                                                             %
 730 %                                                                             %
 731 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 732 %
 733 %  BlurImage() blurs an image.  We convolve the image with a Gaussian operator
 734 %  of the given radius and standard deviation (sigma).  For reasonable results,
 735 %  the radius should be larger than sigma.  Use a radius of 0 and BlurImage()
 736 %  selects a suitable radius for you.
 737 %
 738 %  BlurImage() differs from GaussianBlurImage() in that it uses a separable
 739 %  kernel which is faster but mathematically equivalent to the non-separable
 740 %  kernel.
 741 %
 742 %  The format of the BlurImage method is:
 743 %
 744 %      Image *BlurImage(const Image *image,const double radius,
 745 %        const double sigma,ExceptionInfo *exception)
 746 %      Image *BlurImageChannel(const Image *image,const ChannelType channel,
 747 %        const double radius,const double sigma,ExceptionInfo *exception)
 748 %
 749 %  A description of each parameter follows:
 750 %
 751 %    o image: the image.
 752 %
 753 %    o channel: the channel type.
 754 %
 755 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 756 %      pixel.
 757 %
 758 %    o sigma: the standard deviation of the Gaussian, in pixels.
 759 %
 760 %    o exception: return any errors or warnings in this structure.
 761 %
 762 */
 763
 764 MagickExport Image *BlurImage(const Image *image,const double radius,
 765   const double sigma,ExceptionInfo *exception)
 766 {
 767   Image
 768     *blur_image;
 769
 770   blur_image=BlurImageChannel(image,DefaultChannels,radius,sigma,exception);
 771   return(blur_image);
 772 }
 773
 774 static double *GetBlurKernel(const size_t width,const double sigma)
 775 {
 776   double
 777     *kernel,
 778     normalize;
 779
 780   register ssize_t
 781     i;
 782
 783   ssize_t
 784     j,
 785     k;
 786
 787   /*
 788     Generate a 1-D convolution kernel.
 789   */
 790   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
 791   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 792   if (kernel == (double *) NULL)
 793     return(0);
 794   normalize=0.0;
 795   j=(ssize_t) width/2;
 796   i=0;
 797   for (k=(-j); k <= j; k++)
 798   {
 799     kernel[i]=(double) (exp(-((double) k*k)/(2.0*MagickSigma*MagickSigma))/
 800       (MagickSQ2PI*MagickSigma));
 801     normalize+=kernel[i];
 802     i++;
 803   }
 804   for (i=0; i < (ssize_t) width; i++)
 805     kernel[i]/=normalize;
 806   return(kernel);
 807 }
 808
 809 MagickExport Image *BlurImageChannel(const Image *image,
 810   const ChannelType channel,const double radius,const double sigma,
 811   ExceptionInfo *exception)
 812 {
 813 #define BlurImageTag  "Blur/Image"
 814
 815   CacheView
 816     *blur_view,
 817     *image_view;
 818
 819   double
 820     *kernel;
 821
 822   Image
 823     *blur_image;
 824
 825   MagickBooleanType
 826     status;
 827
 828   MagickOffsetType
 829     progress;
 830
 831   MagickPixelPacket
 832     bias;
 833
 834   register ssize_t
 835     i;
 836
 837   size_t
 838     width;
 839
 840   ssize_t
 841     x,
 842     y;
 843
 844   /*
 845     Initialize blur image attributes.
 846   */
 847   assert(image != (Image *) NULL);
 848   assert(image->signature == MagickSignature);
 849   if (image->debug != MagickFalse)
 850     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 851   assert(exception != (ExceptionInfo *) NULL);
 852   assert(exception->signature == MagickSignature);
 853   blur_image=CloneImage(image,0,0,MagickTrue,exception);
 854   if (blur_image == (Image *) NULL)
 855     return((Image *) NULL);
 856   if (fabs(sigma) <= MagickEpsilon)
 857     return(blur_image);
 858   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
 859     {
 860       InheritException(exception,&blur_image->exception);
 861       blur_image=DestroyImage(blur_image);
 862       return((Image *) NULL);
 863     }
 864   width=GetOptimalKernelWidth1D(radius,sigma);
 865   kernel=GetBlurKernel(width,sigma);
 866   if (kernel == (double *) NULL)
 867     {
 868       blur_image=DestroyImage(blur_image);
 869       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 870     }
 871   if (image->debug != MagickFalse)
 872     {
 873       char
 874         format[MaxTextExtent],
 875         *message;
 876
 877       register const double
 878         *k;
 879
 880       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
 881         "  BlurImage with %.20g kernel:",(double) width);
 882       message=AcquireString("");
 883       k=kernel;
 884       for (i=0; i < (ssize_t) width; i++)
 885       {
 886         *message='\0';
 887         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) i);
 888         (void) ConcatenateString(&message,format);
 889         (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
 890         (void) ConcatenateString(&message,format);
 891         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
 892       }
 893       message=DestroyString(message);
 894     }
 895   /*
 896     Blur rows.
 897   */
 898   status=MagickTrue;
 899   progress=0;
 900   GetMagickPixelPacket(image,&bias);
 901   SetMagickPixelPacketBias(image,&bias);
 902   image_view=AcquireCacheView(image);
 903   blur_view=AcquireCacheView(blur_image);
 904 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 905   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 906 #endif
 907   for (y=0; y < (ssize_t) blur_image->rows; y++)
 908   {
 909     register const IndexPacket
 910       *restrict indexes;
 911
 912     register const PixelPacket
 913       *restrict p;
 914
 915     register IndexPacket
 916       *restrict blur_indexes;
 917
 918     register PixelPacket
 919       *restrict q;
 920
 921     register ssize_t
 922       x;
 923
 924     if (status == MagickFalse)
 925       continue;
 926     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y,
 927       image->columns+width,1,exception);
 928     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
 929       exception);
 930     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 931       {
 932         status=MagickFalse;
 933         continue;
 934       }
 935     indexes=GetCacheViewVirtualIndexQueue(image_view);
 936     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
 937     for (x=0; x < (ssize_t) blur_image->columns; x++)
 938     {
 939       MagickPixelPacket
 940         pixel;
 941
 942       register const double
 943         *restrict k;
 944
 945       register const PixelPacket
 946         *restrict kernel_pixels;
 947
 948       register ssize_t
 949         i;
 950
 951       pixel=bias;
 952       k=kernel;
 953       kernel_pixels=p;
 954       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
 955         {
 956           for (i=0; i < (ssize_t) width; i++)
 957           {
 958             pixel.red+=(*k)*kernel_pixels->red;
 959             pixel.green+=(*k)*kernel_pixels->green;
 960             pixel.blue+=(*k)*kernel_pixels->blue;
 961             k++;
 962             kernel_pixels++;
 963           }
 964           if ((channel & RedChannel) != 0)
 965             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
 966           if ((channel & GreenChannel) != 0)
 967             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
 968           if ((channel & BlueChannel) != 0)
 969             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
 970           if ((channel & OpacityChannel) != 0)
 971             {
 972               k=kernel;
 973               kernel_pixels=p;
 974               for (i=0; i < (ssize_t) width; i++)
 975               {
 976                 pixel.opacity+=(*k)*kernel_pixels->opacity;
 977                 k++;
 978                 kernel_pixels++;
 979               }
 980               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
 981             }
 982           if (((channel & IndexChannel) != 0) &&
 983               (image->colorspace == CMYKColorspace))
 984             {
 985               register const IndexPacket
 986                 *restrict kernel_indexes;
 987
 988               k=kernel;
 989               kernel_indexes=indexes;
 990               for (i=0; i < (ssize_t) width; i++)
 991               {
 992                 pixel.index+=(*k)*(*kernel_indexes);
 993                 k++;
 994                 kernel_indexes++;
 995               }
 996               blur_indexes[x]=ClampToQuantum(pixel.index);
 997             }
 998         }
 999       else
1000         {
1001           MagickRealType
1002             alpha,
1003             gamma;
1004
1005           gamma=0.0;
1006           for (i=0; i < (ssize_t) width; i++)
1007           {
1008             alpha=(MagickRealType) (QuantumScale*
1009               GetAlphaPixelComponent(kernel_pixels));
1010             pixel.red+=(*k)*alpha*kernel_pixels->red;
1011             pixel.green+=(*k)*alpha*kernel_pixels->green;
1012             pixel.blue+=(*k)*alpha*kernel_pixels->blue;
1013             gamma+=(*k)*alpha;
1014             k++;
1015             kernel_pixels++;
1016           }
1017           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1018           if ((channel & RedChannel) != 0)
1019             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
1020           if ((channel & GreenChannel) != 0)
1021             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
1022           if ((channel & BlueChannel) != 0)
1023             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
1024           if ((channel & OpacityChannel) != 0)
1025             {
1026               k=kernel;
1027               kernel_pixels=p;
1028               for (i=0; i < (ssize_t) width; i++)
1029               {
1030                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1031                 k++;
1032                 kernel_pixels++;
1033               }
1034               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1035             }
1036           if (((channel & IndexChannel) != 0) &&
1037               (image->colorspace == CMYKColorspace))
1038             {
1039               register const IndexPacket
1040                 *restrict kernel_indexes;
1041
1042               k=kernel;
1043               kernel_pixels=p;
1044               kernel_indexes=indexes;
1045               for (i=0; i < (ssize_t) width; i++)
1046               {
1047                 alpha=(MagickRealType) (QuantumScale*
1048                   GetAlphaPixelComponent(kernel_pixels));
1049                 pixel.index+=(*k)*alpha*(*kernel_indexes);
1050                 k++;
1051                 kernel_pixels++;
1052                 kernel_indexes++;
1053               }
1054               blur_indexes[x]=ClampToQuantum(gamma*
1055                 GetIndexPixelComponent(&pixel));
1056             }
1057         }
1058       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;
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;
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;
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;
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;
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;
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   if ((image->columns < width) || (image->rows < width))
2883     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
2884   median_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2885     exception);
2886   if (median_image == (Image *) NULL)
2887     return((Image *) NULL);
2888   if (SetImageStorageClass(median_image,DirectClass) == MagickFalse)
2889     {
2890       InheritException(exception,&median_image->exception);
2891       median_image=DestroyImage(median_image);
2892       return((Image *) NULL);
2893     }
2894   pixel_list=AcquirePixelListThreadSet(width);
2895   if (pixel_list == (PixelList **) NULL)
2896     {
2897       median_image=DestroyImage(median_image);
2898       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2899     }
2900   /*
2901     Median filter each image row.
2902   */
2903   status=MagickTrue;
2904   progress=0;
2905   image_view=AcquireCacheView(image);
2906   median_view=AcquireCacheView(median_image);
2907 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2908   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2909 #endif
2910   for (y=0; y < (ssize_t) median_image->rows; y++)
2911   {
2912     const int
2913       id = GetOpenMPThreadId();
2914
2915     register const IndexPacket
2916       *restrict indexes;
2917
2918     register const PixelPacket
2919       *restrict p;
2920
2921     register IndexPacket
2922       *restrict median_indexes;
2923
2924     register PixelPacket
2925       *restrict q;
2926
2927     register ssize_t
2928       x;
2929
2930     if (status == MagickFalse)
2931       continue;
2932     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
2933       (width/2L),image->columns+width,width,exception);
2934     q=QueueCacheViewAuthenticPixels(median_view,0,y,median_image->columns,1,
2935       exception);
2936     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2937       {
2938         status=MagickFalse;
2939         continue;
2940       }
2941     indexes=GetCacheViewVirtualIndexQueue(image_view);
2942     median_indexes=GetCacheViewAuthenticIndexQueue(median_view);
2943     for (x=0; x < (ssize_t) median_image->columns; x++)
2944     {
2945       MagickPixelPacket
2946         pixel;
2947
2948       register const IndexPacket
2949         *restrict s;
2950
2951       register const PixelPacket
2952         *restrict r;
2953
2954       register ssize_t
2955         u,
2956         v;
2957
2958       r=p;
2959       s=indexes+x;
2960       ResetPixelList(pixel_list[id]);
2961       for (v=0; v < (ssize_t) width; v++)
2962       {
2963         for (u=0; u < (ssize_t) width; u++)
2964           InsertPixelList(image,r+u,s+u,pixel_list[id]);
2965         r+=image->columns+width;
2966         s+=image->columns+width;
2967       }
2968       pixel=GetPixelList(pixel_list[id]);
2969       SetPixelPacket(median_image,&pixel,q,median_indexes+x);
2970       p++;
2971       q++;
2972     }
2973     if (SyncCacheViewAuthenticPixels(median_view,exception) == MagickFalse)
2974       status=MagickFalse;
2975     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2976       {
2977         MagickBooleanType
2978           proceed;
2979
2980 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2981   #pragma omp critical (MagickCore_MedianFilterImage)
2982 #endif
2983         proceed=SetImageProgress(image,MedianFilterImageTag,progress++,
2984           image->rows);
2985         if (proceed == MagickFalse)
2986           status=MagickFalse;
2987       }
2988   }
2989   median_view=DestroyCacheView(median_view);
2990   image_view=DestroyCacheView(image_view);
2991   pixel_list=DestroyPixelListThreadSet(pixel_list);
2992   return(median_image);
2993 }
2994 \f
2995 /*
2996 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2997 %                                                                             %
2998 %                                                                             %
2999 %                                                                             %
3000 %     M o d e I m a g e                                                       %
3001 %                                                                             %
3002 %                                                                             %
3003 %                                                                             %
3004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3005 %
3006 %  ModeImage() makes each pixel the 'predominate color' of the neighborhood
3007 %  if the specified radius.
3008 %
3009 %  The format of the ModeImage method is:
3010 %
3011 %      Image *ModeImage(const Image *image,const double radius,
3012 %        ExceptionInfo *exception)
3013 %
3014 %  A description of each parameter follows:
3015 %
3016 %    o image: the image.
3017 %
3018 %    o radius: the radius of the pixel neighborhood.
3019 %
3020 %    o exception: return any errors or warnings in this structure.
3021 %
3022 */
3023
3024 static MagickPixelPacket GetModePixelList(PixelList *pixel_list)
3025 {
3026   MagickPixelPacket
3027     pixel;
3028
3029   register SkipList
3030     *list;
3031
3032   register ssize_t
3033     channel;
3034
3035   size_t
3036     color,
3037     count,
3038     max_count,
3039     mode,
3040     width;
3041
3042   unsigned short
3043     channels[5];
3044
3045   /*
3046     Make each pixel the 'predominate color' of the specified neighborhood.
3047   */
3048   width=pixel_list->center << 1;
3049   for (channel=0; channel < 5; channel++)
3050   {
3051     list=pixel_list->lists+channel;
3052     color=65536UL;
3053     mode=color;
3054     max_count=list->nodes[mode].count;
3055     count=0;
3056     do
3057     {
3058       color=list->nodes[color].next[0];
3059       if (list->nodes[color].count > max_count)
3060         {
3061           mode=color;
3062           max_count=list->nodes[mode].count;
3063         }
3064       count+=list->nodes[color].count;
3065     }
3066     while (count <= width);
3067     channels[channel]=(unsigned short) mode;
3068   }
3069   GetMagickPixelPacket((const Image *) NULL,&pixel);
3070   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3071   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3072   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3073   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3074   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3075   return(pixel);
3076 }
3077
3078 MagickExport Image *ModeImage(const Image *image,const double radius,
3079   ExceptionInfo *exception)
3080 {
3081 #define ModeImageTag  "Mode/Image"
3082
3083   CacheView
3084     *image_view,
3085     *mode_view;
3086
3087   Image
3088     *mode_image;
3089
3090   MagickBooleanType
3091     status;
3092
3093   MagickOffsetType
3094     progress;
3095
3096   PixelList
3097     **restrict pixel_list;
3098
3099   size_t
3100     width;
3101
3102   ssize_t
3103     y;
3104
3105   /*
3106     Initialize mode image attributes.
3107   */
3108   assert(image != (Image *) NULL);
3109   assert(image->signature == MagickSignature);
3110   if (image->debug != MagickFalse)
3111     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3112   assert(exception != (ExceptionInfo *) NULL);
3113   assert(exception->signature == MagickSignature);
3114   width=GetOptimalKernelWidth2D(radius,0.5);
3115   if ((image->columns < width) || (image->rows < width))
3116     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
3117   mode_image=CloneImage(image,image->columns,image->rows,MagickTrue,
3118     exception);
3119   if (mode_image == (Image *) NULL)
3120     return((Image *) NULL);
3121   if (SetImageStorageClass(mode_image,DirectClass) == MagickFalse)
3122     {
3123       InheritException(exception,&mode_image->exception);
3124       mode_image=DestroyImage(mode_image);
3125       return((Image *) NULL);
3126     }
3127   pixel_list=AcquirePixelListThreadSet(width);
3128   if (pixel_list == (PixelList **) NULL)
3129     {
3130       mode_image=DestroyImage(mode_image);
3131       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3132     }
3133   /*
3134     Reduce mode image.
3135   */
3136   status=MagickTrue;
3137   progress=0;
3138   image_view=AcquireCacheView(image);
3139   mode_view=AcquireCacheView(mode_image);
3140 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3141   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
3142 #endif
3143   for (y=0; y < (ssize_t) mode_image->rows; y++)
3144   {
3145     const int
3146       id = GetOpenMPThreadId();
3147
3148     register const IndexPacket
3149       *restrict indexes;
3150
3151     register const PixelPacket
3152       *restrict p;
3153
3154     register IndexPacket
3155       *restrict mode_indexes;
3156
3157     register PixelPacket
3158       *restrict q;
3159
3160     register ssize_t
3161       x;
3162
3163     if (status == MagickFalse)
3164       continue;
3165     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
3166       (width/2L),image->columns+width,width,exception);
3167     q=QueueCacheViewAuthenticPixels(mode_view,0,y,mode_image->columns,1,
3168       exception);
3169     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3170       {
3171         status=MagickFalse;
3172         continue;
3173       }
3174     indexes=GetCacheViewVirtualIndexQueue(image_view);
3175     mode_indexes=GetCacheViewAuthenticIndexQueue(mode_view);
3176     for (x=0; x < (ssize_t) mode_image->columns; x++)
3177     {
3178       MagickPixelPacket
3179         pixel;
3180
3181       register const PixelPacket
3182         *restrict r;
3183
3184       register const IndexPacket
3185         *restrict s;
3186
3187       register ssize_t
3188         u,
3189         v;
3190
3191       r=p;
3192       s=indexes+x;
3193       ResetPixelList(pixel_list[id]);
3194       for (v=0; v < (ssize_t) width; v++)
3195       {
3196         for (u=0; u < (ssize_t) width; u++)
3197           InsertPixelList(image,r+u,s+u,pixel_list[id]);
3198         r+=image->columns+width;
3199         s+=image->columns+width;
3200       }
3201       pixel=GetModePixelList(pixel_list[id]);
3202       SetPixelPacket(mode_image,&pixel,q,mode_indexes+x);
3203       p++;
3204       q++;
3205     }
3206     if (SyncCacheViewAuthenticPixels(mode_view,exception) == MagickFalse)
3207       status=MagickFalse;
3208     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3209       {
3210         MagickBooleanType
3211           proceed;
3212
3213 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3214   #pragma omp critical (MagickCore_ModeImage)
3215 #endif
3216         proceed=SetImageProgress(image,ModeImageTag,progress++,image->rows);
3217         if (proceed == MagickFalse)
3218           status=MagickFalse;
3219       }
3220   }
3221   mode_view=DestroyCacheView(mode_view);
3222   image_view=DestroyCacheView(image_view);
3223   pixel_list=DestroyPixelListThreadSet(pixel_list);
3224   return(mode_image);
3225 }
3226 \f
3227 /*
3228 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3229 %                                                                             %
3230 %                                                                             %
3231 %                                                                             %
3232 %     M o t i o n B l u r I m a g e                                           %
3233 %                                                                             %
3234 %                                                                             %
3235 %                                                                             %
3236 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3237 %
3238 %  MotionBlurImage() simulates motion blur.  We convolve the image with a
3239 %  Gaussian operator of the given radius and standard deviation (sigma).
3240 %  For reasonable results, radius should be larger than sigma.  Use a
3241 %  radius of 0 and MotionBlurImage() selects a suitable radius for you.
3242 %  Angle gives the angle of the blurring motion.
3243 %
3244 %  Andrew Protano contributed this effect.
3245 %
3246 %  The format of the MotionBlurImage method is:
3247 %
3248 %    Image *MotionBlurImage(const Image *image,const double radius,
3249 %      const double sigma,const double angle,ExceptionInfo *exception)
3250 %    Image *MotionBlurImageChannel(const Image *image,const ChannelType channel,
3251 %      const double radius,const double sigma,const double angle,
3252 %      ExceptionInfo *exception)
3253 %
3254 %  A description of each parameter follows:
3255 %
3256 %    o image: the image.
3257 %
3258 %    o channel: the channel type.
3259 %
3260 %    o radius: the radius of the Gaussian, in pixels, not counting the center
3261 %    o radius: the radius of the Gaussian, in pixels, not counting
3262 %      the center pixel.
3263 %
3264 %    o sigma: the standard deviation of the Gaussian, in pixels.
3265 %
3266 %    o angle: Apply the effect along this angle.
3267 %
3268 %    o exception: return any errors or warnings in this structure.
3269 %
3270 */
3271
3272 static double *GetMotionBlurKernel(const size_t width,const double sigma)
3273 {
3274   double
3275     *kernel,
3276     normalize;
3277
3278   register ssize_t
3279     i;
3280
3281   /*
3282    Generate a 1-D convolution kernel.
3283   */
3284   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
3285   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
3286   if (kernel == (double *) NULL)
3287     return(kernel);
3288   normalize=0.0;
3289   for (i=0; i < (ssize_t) width; i++)
3290   {
3291     kernel[i]=(double) (exp((-((double) i*i)/(double) (2.0*MagickSigma*
3292       MagickSigma)))/(MagickSQ2PI*MagickSigma));
3293     normalize+=kernel[i];
3294   }
3295   for (i=0; i < (ssize_t) width; i++)
3296     kernel[i]/=normalize;
3297   return(kernel);
3298 }
3299
3300 MagickExport Image *MotionBlurImage(const Image *image,const double radius,
3301   const double sigma,const double angle,ExceptionInfo *exception)
3302 {
3303   Image
3304     *motion_blur;
3305
3306   motion_blur=MotionBlurImageChannel(image,DefaultChannels,radius,sigma,angle,
3307     exception);
3308   return(motion_blur);
3309 }
3310
3311 MagickExport Image *MotionBlurImageChannel(const Image *image,
3312   const ChannelType channel,const double radius,const double sigma,
3313   const double angle,ExceptionInfo *exception)
3314 {
3315   CacheView
3316     *blur_view,
3317     *image_view;
3318
3319   double
3320     *kernel;
3321
3322   Image
3323     *blur_image;
3324
3325   MagickBooleanType
3326     status;
3327
3328   MagickOffsetType
3329     progress;
3330
3331   MagickPixelPacket
3332     bias;
3333
3334   OffsetInfo
3335     *offset;
3336
3337   PointInfo
3338     point;
3339
3340   register ssize_t
3341     i;
3342
3343   size_t
3344     width;
3345
3346   ssize_t
3347     y;
3348
3349   assert(image != (Image *) NULL);
3350   assert(image->signature == MagickSignature);
3351   if (image->debug != MagickFalse)
3352     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3353   assert(exception != (ExceptionInfo *) NULL);
3354   width=GetOptimalKernelWidth1D(radius,sigma);
3355   kernel=GetMotionBlurKernel(width,sigma);
3356   if (kernel == (double *) NULL)
3357     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3358   offset=(OffsetInfo *) AcquireQuantumMemory(width,sizeof(*offset));
3359   if (offset == (OffsetInfo *) NULL)
3360     {
3361       kernel=(double *) RelinquishMagickMemory(kernel);
3362       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3363     }
3364   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3365   if (blur_image == (Image *) NULL)
3366     {
3367       kernel=(double *) RelinquishMagickMemory(kernel);
3368       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
3369       return((Image *) NULL);
3370     }
3371   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
3372     {
3373       kernel=(double *) RelinquishMagickMemory(kernel);
3374       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
3375       InheritException(exception,&blur_image->exception);
3376       blur_image=DestroyImage(blur_image);
3377       return((Image *) NULL);
3378     }
3379   point.x=(double) width*sin(DegreesToRadians(angle));
3380   point.y=(double) width*cos(DegreesToRadians(angle));
3381   for (i=0; i < (ssize_t) width; i++)
3382   {
3383     offset[i].x=(ssize_t) ceil((double) (i*point.y)/hypot(point.x,point.y)-0.5);
3384     offset[i].y=(ssize_t) ceil((double) (i*point.x)/hypot(point.x,point.y)-0.5);
3385   }
3386   /*
3387     Motion blur image.
3388   */
3389   status=MagickTrue;
3390   progress=0;
3391   GetMagickPixelPacket(image,&bias);
3392   image_view=AcquireCacheView(image);
3393   blur_view=AcquireCacheView(blur_image);
3394 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3395   #pragma omp parallel for schedule(dynamic,4) shared(progress,status) omp_throttle(1)
3396 #endif
3397   for (y=0; y < (ssize_t) image->rows; y++)
3398   {
3399     register IndexPacket
3400       *restrict blur_indexes;
3401
3402     register PixelPacket
3403       *restrict q;
3404
3405     register ssize_t
3406       x;
3407
3408     if (status == MagickFalse)
3409       continue;
3410     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
3411       exception);
3412     if (q == (PixelPacket *) NULL)
3413       {
3414         status=MagickFalse;
3415         continue;
3416       }
3417     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
3418     for (x=0; x < (ssize_t) image->columns; x++)
3419     {
3420       MagickPixelPacket
3421         qixel;
3422
3423       PixelPacket
3424         pixel;
3425
3426       register const IndexPacket
3427         *restrict indexes;
3428
3429       register double
3430         *restrict k;
3431
3432       register ssize_t
3433         i;
3434
3435       k=kernel;
3436       qixel=bias;
3437       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
3438         {
3439           for (i=0; i < (ssize_t) width; i++)
3440           {
3441             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
3442               offset[i].y,&pixel,exception);
3443             qixel.red+=(*k)*pixel.red;
3444             qixel.green+=(*k)*pixel.green;
3445             qixel.blue+=(*k)*pixel.blue;
3446             qixel.opacity+=(*k)*pixel.opacity;
3447             if (image->colorspace == CMYKColorspace)
3448               {
3449                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3450                 qixel.index+=(*k)*(*indexes);
3451               }
3452             k++;
3453           }
3454           if ((channel & RedChannel) != 0)
3455             q->red=ClampToQuantum(qixel.red);
3456           if ((channel & GreenChannel) != 0)
3457             q->green=ClampToQuantum(qixel.green);
3458           if ((channel & BlueChannel) != 0)
3459             q->blue=ClampToQuantum(qixel.blue);
3460           if ((channel & OpacityChannel) != 0)
3461             q->opacity=ClampToQuantum(qixel.opacity);
3462           if (((channel & IndexChannel) != 0) &&
3463               (image->colorspace == CMYKColorspace))
3464             blur_indexes[x]=(IndexPacket) ClampToQuantum(qixel.index);
3465         }
3466       else
3467         {
3468           MagickRealType
3469             alpha,
3470             gamma;
3471
3472           alpha=0.0;
3473           gamma=0.0;
3474           for (i=0; i < (ssize_t) width; i++)
3475           {
3476             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
3477               offset[i].y,&pixel,exception);
3478             alpha=(MagickRealType) (QuantumScale*
3479               GetAlphaPixelComponent(&pixel));
3480             qixel.red+=(*k)*alpha*pixel.red;
3481             qixel.green+=(*k)*alpha*pixel.green;
3482             qixel.blue+=(*k)*alpha*pixel.blue;
3483             qixel.opacity+=(*k)*pixel.opacity;
3484             if (image->colorspace == CMYKColorspace)
3485               {
3486                 indexes=GetCacheViewVirtualIndexQueue(image_view);
3487                 qixel.index+=(*k)*alpha*(*indexes);
3488               }
3489             gamma+=(*k)*alpha;
3490             k++;
3491           }
3492           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
3493           if ((channel & RedChannel) != 0)
3494             q->red=ClampToQuantum(gamma*qixel.red);
3495           if ((channel & GreenChannel) != 0)
3496             q->green=ClampToQuantum(gamma*qixel.green);
3497           if ((channel & BlueChannel) != 0)
3498             q->blue=ClampToQuantum(gamma*qixel.blue);
3499           if ((channel & OpacityChannel) != 0)
3500             q->opacity=ClampToQuantum(qixel.opacity);
3501           if (((channel & IndexChannel) != 0) &&
3502               (image->colorspace == CMYKColorspace))
3503             blur_indexes[x]=(IndexPacket) ClampToQuantum(gamma*qixel.index);
3504         }
3505       q++;
3506     }
3507     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
3508       status=MagickFalse;
3509     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3510       {
3511         MagickBooleanType
3512           proceed;
3513
3514 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3515   #pragma omp critical (MagickCore_MotionBlurImageChannel)
3516 #endif
3517         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
3518         if (proceed == MagickFalse)
3519           status=MagickFalse;
3520       }
3521   }
3522   blur_view=DestroyCacheView(blur_view);
3523   image_view=DestroyCacheView(image_view);
3524   kernel=(double *) RelinquishMagickMemory(kernel);
3525   offset=(OffsetInfo *) RelinquishMagickMemory(offset);
3526   if (status == MagickFalse)
3527     blur_image=DestroyImage(blur_image);
3528   return(blur_image);
3529 }
3530 \f
3531 /*
3532 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3533 %                                                                             %
3534 %                                                                             %
3535 %                                                                             %
3536 %     P r e v i e w I m a g e                                                 %
3537 %                                                                             %
3538 %                                                                             %
3539 %                                                                             %
3540 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3541 %
3542 %  PreviewImage() tiles 9 thumbnails of the specified image with an image
3543 %  processing operation applied with varying parameters.  This may be helpful
3544 %  pin-pointing an appropriate parameter for a particular image processing
3545 %  operation.
3546 %
3547 %  The format of the PreviewImages method is:
3548 %
3549 %      Image *PreviewImages(const Image *image,const PreviewType preview,
3550 %        ExceptionInfo *exception)
3551 %
3552 %  A description of each parameter follows:
3553 %
3554 %    o image: the image.
3555 %
3556 %    o preview: the image processing operation.
3557 %
3558 %    o exception: return any errors or warnings in this structure.
3559 %
3560 */
3561 MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
3562   ExceptionInfo *exception)
3563 {
3564 #define NumberTiles  9
3565 #define PreviewImageTag  "Preview/Image"
3566 #define DefaultPreviewGeometry  "204x204+10+10"
3567
3568   char
3569     factor[MaxTextExtent],
3570     label[MaxTextExtent];
3571
3572   double
3573     degrees,
3574     gamma,
3575     percentage,
3576     radius,
3577     sigma,
3578     threshold;
3579
3580   Image
3581     *images,
3582     *montage_image,
3583     *preview_image,
3584     *thumbnail;
3585
3586   ImageInfo
3587     *preview_info;
3588
3589   MagickBooleanType
3590     proceed;
3591
3592   MontageInfo
3593     *montage_info;
3594
3595   QuantizeInfo
3596     quantize_info;
3597
3598   RectangleInfo
3599     geometry;
3600
3601   register ssize_t
3602     i,
3603     x;
3604
3605   size_t
3606     colors;
3607
3608   ssize_t
3609     y;
3610
3611   /*
3612     Open output image file.
3613   */
3614   assert(image != (Image *) NULL);
3615   assert(image->signature == MagickSignature);
3616   if (image->debug != MagickFalse)
3617     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3618   colors=2;
3619   degrees=0.0;
3620   gamma=(-0.2f);
3621   preview_info=AcquireImageInfo();
3622   SetGeometry(image,&geometry);
3623   (void) ParseMetaGeometry(DefaultPreviewGeometry,&geometry.x,&geometry.y,
3624     &geometry.width,&geometry.height);
3625   images=NewImageList();
3626   percentage=12.5;
3627   GetQuantizeInfo(&quantize_info);
3628   radius=0.0;
3629   sigma=1.0;
3630   threshold=0.0;
3631   x=0;
3632   y=0;
3633   for (i=0; i < NumberTiles; i++)
3634   {
3635     thumbnail=ThumbnailImage(image,geometry.width,geometry.height,exception);
3636     if (thumbnail == (Image *) NULL)
3637       break;
3638     (void) SetImageProgressMonitor(thumbnail,(MagickProgressMonitor) NULL,
3639       (void *) NULL);
3640     (void) SetImageProperty(thumbnail,"label",DefaultTileLabel);
3641     if (i == (NumberTiles/2))
3642       {
3643         (void) QueryColorDatabase("#dfdfdf",&thumbnail->matte_color,exception);
3644         AppendImageToList(&images,thumbnail);
3645         continue;
3646       }
3647     switch (preview)
3648     {
3649       case RotatePreview:
3650       {
3651         degrees+=45.0;
3652         preview_image=RotateImage(thumbnail,degrees,exception);
3653         (void) FormatMagickString(label,MaxTextExtent,"rotate %g",degrees);
3654         break;
3655       }
3656       case ShearPreview:
3657       {
3658         degrees+=5.0;
3659         preview_image=ShearImage(thumbnail,degrees,degrees,exception);
3660         (void) FormatMagickString(label,MaxTextExtent,"shear %gx%g",
3661           degrees,2.0*degrees);
3662         break;
3663       }
3664       case RollPreview:
3665       {
3666         x=(ssize_t) ((i+1)*thumbnail->columns)/NumberTiles;
3667         y=(ssize_t) ((i+1)*thumbnail->rows)/NumberTiles;
3668         preview_image=RollImage(thumbnail,x,y,exception);
3669         (void) FormatMagickString(label,MaxTextExtent,"roll %+.20gx%+.20g",
3670           (double) x,(double) y);
3671         break;
3672       }
3673       case HuePreview:
3674       {
3675         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3676         if (preview_image == (Image *) NULL)
3677           break;
3678         (void) FormatMagickString(factor,MaxTextExtent,"100,100,%g",
3679           2.0*percentage);
3680         (void) ModulateImage(preview_image,factor);
3681         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3682         break;
3683       }
3684       case SaturationPreview:
3685       {
3686         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3687         if (preview_image == (Image *) NULL)
3688           break;
3689         (void) FormatMagickString(factor,MaxTextExtent,"100,%g",
3690           2.0*percentage);
3691         (void) ModulateImage(preview_image,factor);
3692         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3693         break;
3694       }
3695       case BrightnessPreview:
3696       {
3697         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3698         if (preview_image == (Image *) NULL)
3699           break;
3700         (void) FormatMagickString(factor,MaxTextExtent,"%g",2.0*percentage);
3701         (void) ModulateImage(preview_image,factor);
3702         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
3703         break;
3704       }
3705       case GammaPreview:
3706       default:
3707       {
3708         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3709         if (preview_image == (Image *) NULL)
3710           break;
3711         gamma+=0.4f;
3712         (void) GammaImageChannel(preview_image,DefaultChannels,gamma);
3713         (void) FormatMagickString(label,MaxTextExtent,"gamma %g",gamma);
3714         break;
3715       }
3716       case SpiffPreview:
3717       {
3718         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3719         if (preview_image != (Image *) NULL)
3720           for (x=0; x < i; x++)
3721             (void) ContrastImage(preview_image,MagickTrue);
3722         (void) FormatMagickString(label,MaxTextExtent,"contrast (%.20g)",
3723           (double) i+1);
3724         break;
3725       }
3726       case DullPreview:
3727       {
3728         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3729         if (preview_image == (Image *) NULL)
3730           break;
3731         for (x=0; x < i; x++)
3732           (void) ContrastImage(preview_image,MagickFalse);
3733         (void) FormatMagickString(label,MaxTextExtent,"+contrast (%.20g)",
3734           (double) i+1);
3735         break;
3736       }
3737       case GrayscalePreview:
3738       {
3739         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3740         if (preview_image == (Image *) NULL)
3741           break;
3742         colors<<=1;
3743         quantize_info.number_colors=colors;
3744         quantize_info.colorspace=GRAYColorspace;
3745         (void) QuantizeImage(&quantize_info,preview_image);
3746         (void) FormatMagickString(label,MaxTextExtent,
3747           "-colorspace gray -colors %.20g",(double) colors);
3748         break;
3749       }
3750       case QuantizePreview:
3751       {
3752         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3753         if (preview_image == (Image *) NULL)
3754           break;
3755         colors<<=1;
3756         quantize_info.number_colors=colors;
3757         (void) QuantizeImage(&quantize_info,preview_image);
3758         (void) FormatMagickString(label,MaxTextExtent,"colors %.20g",(double)
3759           colors);
3760         break;
3761       }
3762       case DespecklePreview:
3763       {
3764         for (x=0; x < (i-1); x++)
3765         {
3766           preview_image=DespeckleImage(thumbnail,exception);
3767           if (preview_image == (Image *) NULL)
3768             break;
3769           thumbnail=DestroyImage(thumbnail);
3770           thumbnail=preview_image;
3771         }
3772         preview_image=DespeckleImage(thumbnail,exception);
3773         if (preview_image == (Image *) NULL)
3774           break;
3775         (void) FormatMagickString(label,MaxTextExtent,"despeckle (%.20g)",
3776           (double) i+1);
3777         break;
3778       }
3779       case ReduceNoisePreview:
3780       {
3781         preview_image=ReduceNoiseImage(thumbnail,radius,exception);
3782         (void) FormatMagickString(label,MaxTextExtent,"noise %g",radius);
3783         break;
3784       }
3785       case AddNoisePreview:
3786       {
3787         switch ((int) i)
3788         {
3789           case 0:
3790           {
3791             (void) CopyMagickString(factor,"uniform",MaxTextExtent);
3792             break;
3793           }
3794           case 1:
3795           {
3796             (void) CopyMagickString(factor,"gaussian",MaxTextExtent);
3797             break;
3798           }
3799           case 2:
3800           {
3801             (void) CopyMagickString(factor,"multiplicative",MaxTextExtent);
3802             break;
3803           }
3804           case 3:
3805           {
3806             (void) CopyMagickString(factor,"impulse",MaxTextExtent);
3807             break;
3808           }
3809           case 4:
3810           {
3811             (void) CopyMagickString(factor,"laplacian",MaxTextExtent);
3812             break;
3813           }
3814           case 5:
3815           {
3816             (void) CopyMagickString(factor,"Poisson",MaxTextExtent);
3817             break;
3818           }
3819           default:
3820           {
3821             (void) CopyMagickString(thumbnail->magick,"NULL",MaxTextExtent);
3822             break;
3823           }
3824         }
3825         preview_image=ReduceNoiseImage(thumbnail,(double) i,exception);
3826         (void) FormatMagickString(label,MaxTextExtent,"+noise %s",factor);
3827         break;
3828       }
3829       case SharpenPreview:
3830       {
3831         preview_image=SharpenImage(thumbnail,radius,sigma,exception);
3832         (void) FormatMagickString(label,MaxTextExtent,"sharpen %gx%g",
3833           radius,sigma);
3834         break;
3835       }
3836       case BlurPreview:
3837       {
3838         preview_image=BlurImage(thumbnail,radius,sigma,exception);
3839         (void) FormatMagickString(label,MaxTextExtent,"blur %gx%g",radius,
3840           sigma);
3841         break;
3842       }
3843       case ThresholdPreview:
3844       {
3845         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3846         if (preview_image == (Image *) NULL)
3847           break;
3848         (void) BilevelImage(thumbnail,
3849           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
3850         (void) FormatMagickString(label,MaxTextExtent,"threshold %g",
3851           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
3852         break;
3853       }
3854       case EdgeDetectPreview:
3855       {
3856         preview_image=EdgeImage(thumbnail,radius,exception);
3857         (void) FormatMagickString(label,MaxTextExtent,"edge %g",radius);
3858         break;
3859       }
3860       case SpreadPreview:
3861       {
3862         preview_image=SpreadImage(thumbnail,radius,exception);
3863         (void) FormatMagickString(label,MaxTextExtent,"spread %g",
3864           radius+0.5);
3865         break;
3866       }
3867       case SolarizePreview:
3868       {
3869         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3870         if (preview_image == (Image *) NULL)
3871           break;
3872         (void) SolarizeImage(preview_image,(double) QuantumRange*
3873           percentage/100.0);
3874         (void) FormatMagickString(label,MaxTextExtent,"solarize %g",
3875           (QuantumRange*percentage)/100.0);
3876         break;
3877       }
3878       case ShadePreview:
3879       {
3880         degrees+=10.0;
3881         preview_image=ShadeImage(thumbnail,MagickTrue,degrees,degrees,
3882           exception);
3883         (void) FormatMagickString(label,MaxTextExtent,"shade %gx%g",
3884           degrees,degrees);
3885         break;
3886       }
3887       case RaisePreview:
3888       {
3889         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3890         if (preview_image == (Image *) NULL)
3891           break;
3892         geometry.width=(size_t) (2*i+2);
3893         geometry.height=(size_t) (2*i+2);
3894         geometry.x=i/2;
3895         geometry.y=i/2;
3896         (void) RaiseImage(preview_image,&geometry,MagickTrue);
3897         (void) FormatMagickString(label,MaxTextExtent,
3898           "raise %.20gx%.20g%+.20g%+.20g",(double) geometry.width,(double)
3899           geometry.height,(double) geometry.x,(double) geometry.y);
3900         break;
3901       }
3902       case SegmentPreview:
3903       {
3904         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3905         if (preview_image == (Image *) NULL)
3906           break;
3907         threshold+=0.4f;
3908         (void) SegmentImage(preview_image,RGBColorspace,MagickFalse,threshold,
3909           threshold);
3910         (void) FormatMagickString(label,MaxTextExtent,"segment %gx%g",
3911           threshold,threshold);
3912         break;
3913       }
3914       case SwirlPreview:
3915       {
3916         preview_image=SwirlImage(thumbnail,degrees,exception);
3917         (void) FormatMagickString(label,MaxTextExtent,"swirl %g",degrees);
3918         degrees+=45.0;
3919         break;
3920       }
3921       case ImplodePreview:
3922       {
3923         degrees+=0.1f;
3924         preview_image=ImplodeImage(thumbnail,degrees,exception);
3925         (void) FormatMagickString(label,MaxTextExtent,"implode %g",degrees);
3926         break;
3927       }
3928       case WavePreview:
3929       {
3930         degrees+=5.0f;
3931         preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,exception);
3932         (void) FormatMagickString(label,MaxTextExtent,"wave %gx%g",
3933           0.5*degrees,2.0*degrees);
3934         break;
3935       }
3936       case OilPaintPreview:
3937       {
3938         preview_image=OilPaintImage(thumbnail,(double) radius,exception);
3939         (void) FormatMagickString(label,MaxTextExtent,"paint %g",radius);
3940         break;
3941       }
3942       case CharcoalDrawingPreview:
3943       {
3944         preview_image=CharcoalImage(thumbnail,(double) radius,(double) sigma,
3945           exception);
3946         (void) FormatMagickString(label,MaxTextExtent,"charcoal %gx%g",
3947           radius,sigma);
3948         break;
3949       }
3950       case JPEGPreview:
3951       {
3952         char
3953           filename[MaxTextExtent];
3954
3955         int
3956           file;
3957
3958         MagickBooleanType
3959           status;
3960
3961         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3962         if (preview_image == (Image *) NULL)
3963           break;
3964         preview_info->quality=(size_t) percentage;
3965         (void) FormatMagickString(factor,MaxTextExtent,"%.20g",(double)
3966           preview_info->quality);
3967         file=AcquireUniqueFileResource(filename);
3968         if (file != -1)
3969           file=close(file)-1;
3970         (void) FormatMagickString(preview_image->filename,MaxTextExtent,
3971           "jpeg:%s",filename);
3972         status=WriteImage(preview_info,preview_image);
3973         if (status != MagickFalse)
3974           {
3975             Image
3976               *quality_image;
3977
3978             (void) CopyMagickString(preview_info->filename,
3979               preview_image->filename,MaxTextExtent);
3980             quality_image=ReadImage(preview_info,exception);
3981             if (quality_image != (Image *) NULL)
3982               {
3983                 preview_image=DestroyImage(preview_image);
3984                 preview_image=quality_image;
3985               }
3986           }
3987         (void) RelinquishUniqueFileResource(preview_image->filename);
3988         if ((GetBlobSize(preview_image)/1024) >= 1024)
3989           (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%gmb ",
3990             factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
3991             1024.0/1024.0);
3992         else
3993           if (GetBlobSize(preview_image) >= 1024)
3994             (void) FormatMagickString(label,MaxTextExtent,
3995               "quality %s\n%gkb ",factor,(double) ((MagickOffsetType)
3996               GetBlobSize(preview_image))/1024.0);
3997           else
3998             (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%.20gb ",
3999               factor,(double) GetBlobSize(thumbnail));
4000         break;
4001       }
4002     }
4003     thumbnail=DestroyImage(thumbnail);
4004     percentage+=12.5;
4005     radius+=0.5;
4006     sigma+=0.25;
4007     if (preview_image == (Image *) NULL)
4008       break;
4009     (void) DeleteImageProperty(preview_image,"label");
4010     (void) SetImageProperty(preview_image,"label",label);
4011     AppendImageToList(&images,preview_image);
4012     proceed=SetImageProgress(image,PreviewImageTag,(MagickOffsetType) i,
4013       NumberTiles);
4014     if (proceed == MagickFalse)
4015       break;
4016   }
4017   if (images == (Image *) NULL)
4018     {
4019       preview_info=DestroyImageInfo(preview_info);
4020       return((Image *) NULL);
4021     }
4022   /*
4023     Create the montage.
4024   */
4025   montage_info=CloneMontageInfo(preview_info,(MontageInfo *) NULL);
4026   (void) CopyMagickString(montage_info->filename,image->filename,MaxTextExtent);
4027   montage_info->shadow=MagickTrue;
4028   (void) CloneString(&montage_info->tile,"3x3");
4029   (void) CloneString(&montage_info->geometry,DefaultPreviewGeometry);
4030   (void) CloneString(&montage_info->frame,DefaultTileFrame);
4031   montage_image=MontageImages(images,montage_info,exception);
4032   montage_info=DestroyMontageInfo(montage_info);
4033   images=DestroyImageList(images);
4034   if (montage_image == (Image *) NULL)
4035     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4036   if (montage_image->montage != (char *) NULL)
4037     {
4038       /*
4039         Free image directory.
4040       */
4041       montage_image->montage=(char *) RelinquishMagickMemory(
4042         montage_image->montage);
4043       if (image->directory != (char *) NULL)
4044         montage_image->directory=(char *) RelinquishMagickMemory(
4045           montage_image->directory);
4046     }
4047   preview_info=DestroyImageInfo(preview_info);
4048   return(montage_image);
4049 }
4050 \f
4051 /*
4052 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4053 %                                                                             %
4054 %                                                                             %
4055 %                                                                             %
4056 %     R a d i a l B l u r I m a g e                                           %
4057 %                                                                             %
4058 %                                                                             %
4059 %                                                                             %
4060 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4061 %
4062 %  RadialBlurImage() applies a radial blur to the image.
4063 %
4064 %  Andrew Protano contributed this effect.
4065 %
4066 %  The format of the RadialBlurImage method is:
4067 %
4068 %    Image *RadialBlurImage(const Image *image,const double angle,
4069 %      ExceptionInfo *exception)
4070 %    Image *RadialBlurImageChannel(const Image *image,const ChannelType channel,
4071 %      const double angle,ExceptionInfo *exception)
4072 %
4073 %  A description of each parameter follows:
4074 %
4075 %    o image: the image.
4076 %
4077 %    o channel: the channel type.
4078 %
4079 %    o angle: the angle of the radial blur.
4080 %
4081 %    o exception: return any errors or warnings in this structure.
4082 %
4083 */
4084
4085 MagickExport Image *RadialBlurImage(const Image *image,const double angle,
4086   ExceptionInfo *exception)
4087 {
4088   Image
4089     *blur_image;
4090
4091   blur_image=RadialBlurImageChannel(image,DefaultChannels,angle,exception);
4092   return(blur_image);
4093 }
4094
4095 MagickExport Image *RadialBlurImageChannel(const Image *image,
4096   const ChannelType channel,const double angle,ExceptionInfo *exception)
4097 {
4098   CacheView
4099     *blur_view,
4100     *image_view;
4101
4102   Image
4103     *blur_image;
4104
4105   MagickBooleanType
4106     status;
4107
4108   MagickOffsetType
4109     progress;
4110
4111   MagickPixelPacket
4112     bias;
4113
4114   MagickRealType
4115     blur_radius,
4116     *cos_theta,
4117     offset,
4118     *sin_theta,
4119     theta;
4120
4121   PointInfo
4122     blur_center;
4123
4124   register ssize_t
4125     i;
4126
4127   size_t
4128     n;
4129
4130   ssize_t
4131     y;
4132
4133   /*
4134     Allocate blur image.
4135   */
4136   assert(image != (Image *) NULL);
4137   assert(image->signature == MagickSignature);
4138   if (image->debug != MagickFalse)
4139     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4140   assert(exception != (ExceptionInfo *) NULL);
4141   assert(exception->signature == MagickSignature);
4142   blur_image=CloneImage(image,0,0,MagickTrue,exception);
4143   if (blur_image == (Image *) NULL)
4144     return((Image *) NULL);
4145   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
4146     {
4147       InheritException(exception,&blur_image->exception);
4148       blur_image=DestroyImage(blur_image);
4149       return((Image *) NULL);
4150     }
4151   blur_center.x=(double) image->columns/2.0;
4152   blur_center.y=(double) image->rows/2.0;
4153   blur_radius=hypot(blur_center.x,blur_center.y);
4154   n=(size_t) fabs(4.0*DegreesToRadians(angle)*sqrt((double) blur_radius)+2UL);
4155   theta=DegreesToRadians(angle)/(MagickRealType) (n-1);
4156   cos_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
4157     sizeof(*cos_theta));
4158   sin_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
4159     sizeof(*sin_theta));
4160   if ((cos_theta == (MagickRealType *) NULL) ||
4161       (sin_theta == (MagickRealType *) NULL))
4162     {
4163       blur_image=DestroyImage(blur_image);
4164       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4165     }
4166   offset=theta*(MagickRealType) (n-1)/2.0;
4167   for (i=0; i < (ssize_t) n; i++)
4168   {
4169     cos_theta[i]=cos((double) (theta*i-offset));
4170     sin_theta[i]=sin((double) (theta*i-offset));
4171   }
4172   /*
4173     Radial blur image.
4174   */
4175   status=MagickTrue;
4176   progress=0;
4177   GetMagickPixelPacket(image,&bias);
4178   image_view=AcquireCacheView(image);
4179   blur_view=AcquireCacheView(blur_image);
4180 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4181   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4182 #endif
4183   for (y=0; y < (ssize_t) blur_image->rows; y++)
4184   {
4185     register const IndexPacket
4186       *restrict indexes;
4187
4188     register IndexPacket
4189       *restrict blur_indexes;
4190
4191     register PixelPacket
4192       *restrict q;
4193
4194     register ssize_t
4195       x;
4196
4197     if (status == MagickFalse)
4198       continue;
4199     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
4200       exception);
4201     if (q == (PixelPacket *) NULL)
4202       {
4203         status=MagickFalse;
4204         continue;
4205       }
4206     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
4207     for (x=0; x < (ssize_t) blur_image->columns; x++)
4208     {
4209       MagickPixelPacket
4210         qixel;
4211
4212       MagickRealType
4213         normalize,
4214         radius;
4215
4216       PixelPacket
4217         pixel;
4218
4219       PointInfo
4220         center;
4221
4222       register ssize_t
4223         i;
4224
4225       size_t
4226         step;
4227
4228       center.x=(double) x-blur_center.x;
4229       center.y=(double) y-blur_center.y;
4230       radius=hypot((double) center.x,center.y);
4231       if (radius == 0)
4232         step=1;
4233       else
4234         {
4235           step=(size_t) (blur_radius/radius);
4236           if (step == 0)
4237             step=1;
4238           else
4239             if (step >= n)
4240               step=n-1;
4241         }
4242       normalize=0.0;
4243       qixel=bias;
4244       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
4245         {
4246           for (i=0; i < (ssize_t) n; i+=(ssize_t) step)
4247           {
4248             (void) GetOneCacheViewVirtualPixel(image_view,(ssize_t)
4249               (blur_center.x+center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),
4250               (ssize_t) (blur_center.y+center.x*sin_theta[i]+center.y*
4251               cos_theta[i]+0.5),&pixel,exception);
4252             qixel.red+=pixel.red;
4253             qixel.green+=pixel.green;
4254             qixel.blue+=pixel.blue;
4255             qixel.opacity+=pixel.opacity;
4256             if (image->colorspace == CMYKColorspace)
4257               {
4258                 indexes=GetCacheViewVirtualIndexQueue(image_view);
4259                 qixel.index+=(*indexes);
4260               }
4261             normalize+=1.0;
4262           }
4263           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
4264             normalize);
4265           if ((channel & RedChannel) != 0)
4266             q->red=ClampToQuantum(normalize*qixel.red);
4267           if ((channel & GreenChannel) != 0)
4268             q->green=ClampToQuantum(normalize*qixel.green);
4269           if ((channel & BlueChannel) != 0)
4270             q->blue=ClampToQuantum(normalize*qixel.blue);
4271           if ((channel & OpacityChannel) != 0)
4272             q->opacity=ClampToQuantum(normalize*qixel.opacity);
4273           if (((channel & IndexChannel) != 0) &&
4274               (image->colorspace == CMYKColorspace))
4275             blur_indexes[x]=(IndexPacket) ClampToQuantum(normalize*qixel.index);
4276         }
4277       else
4278         {
4279           MagickRealType
4280             alpha,
4281             gamma;
4282
4283           alpha=1.0;
4284           gamma=0.0;
4285           for (i=0; i < (ssize_t) n; i+=(ssize_t) step)
4286           {
4287             (void) GetOneCacheViewVirtualPixel(image_view,(ssize_t)
4288               (blur_center.x+center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),
4289               (ssize_t) (blur_center.y+center.x*sin_theta[i]+center.y*
4290               cos_theta[i]+0.5),&pixel,exception);
4291             alpha=(MagickRealType) (QuantumScale*
4292               GetAlphaPixelComponent(&pixel));
4293             qixel.red+=alpha*pixel.red;
4294             qixel.green+=alpha*pixel.green;
4295             qixel.blue+=alpha*pixel.blue;
4296             qixel.opacity+=pixel.opacity;
4297             if (image->colorspace == CMYKColorspace)
4298               {
4299                 indexes=GetCacheViewVirtualIndexQueue(image_view);
4300                 qixel.index+=alpha*(*indexes);
4301               }
4302             gamma+=alpha;
4303             normalize+=1.0;
4304           }
4305           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4306           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
4307             normalize);
4308           if ((channel & RedChannel) != 0)
4309             q->red=ClampToQuantum(gamma*qixel.red);
4310           if ((channel & GreenChannel) != 0)
4311             q->green=ClampToQuantum(gamma*qixel.green);
4312           if ((channel & BlueChannel) != 0)
4313             q->blue=ClampToQuantum(gamma*qixel.blue);
4314           if ((channel & OpacityChannel) != 0)
4315             q->opacity=ClampToQuantum(normalize*qixel.opacity);
4316           if (((channel & IndexChannel) != 0) &&
4317               (image->colorspace == CMYKColorspace))
4318             blur_indexes[x]=(IndexPacket) ClampToQuantum(gamma*qixel.index);
4319         }
4320       q++;
4321     }
4322     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
4323       status=MagickFalse;
4324     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4325       {
4326         MagickBooleanType
4327           proceed;
4328
4329 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4330   #pragma omp critical (MagickCore_RadialBlurImageChannel)
4331 #endif
4332         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
4333         if (proceed == MagickFalse)
4334           status=MagickFalse;
4335       }
4336   }
4337   blur_view=DestroyCacheView(blur_view);
4338   image_view=DestroyCacheView(image_view);
4339   cos_theta=(MagickRealType *) RelinquishMagickMemory(cos_theta);
4340   sin_theta=(MagickRealType *) RelinquishMagickMemory(sin_theta);
4341   if (status == MagickFalse)
4342     blur_image=DestroyImage(blur_image);
4343   return(blur_image);
4344 }
4345 \f
4346 /*
4347 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4348 %                                                                             %
4349 %                                                                             %
4350 %                                                                             %
4351 %     R e d u c e N o i s e I m a g e                                         %
4352 %                                                                             %
4353 %                                                                             %
4354 %                                                                             %
4355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4356 %
4357 %  ReduceNoiseImage() smooths the contours of an image while still preserving
4358 %  edge information.  The algorithm works by replacing each pixel with its
4359 %  neighbor closest in value.  A neighbor is defined by radius.  Use a radius
4360 %  of 0 and ReduceNoise() selects a suitable radius for you.
4361 %
4362 %  The format of the ReduceNoiseImage method is:
4363 %
4364 %      Image *ReduceNoiseImage(const Image *image,const double radius,
4365 %        ExceptionInfo *exception)
4366 %
4367 %  A description of each parameter follows:
4368 %
4369 %    o image: the image.
4370 %
4371 %    o radius: the radius of the pixel neighborhood.
4372 %
4373 %    o exception: return any errors or warnings in this structure.
4374 %
4375 */
4376
4377 static MagickPixelPacket GetNonpeakPixelList(PixelList *pixel_list)
4378 {
4379   MagickPixelPacket
4380     pixel;
4381
4382   register SkipList
4383     *list;
4384
4385   register ssize_t
4386     channel;
4387
4388   size_t
4389     center,
4390     color,
4391     count,
4392     next,
4393     previous;
4394
4395   unsigned short
4396     channels[5];
4397
4398   /*
4399     Finds the median value for each of the color.
4400   */
4401   center=pixel_list->center;
4402   for (channel=0; channel < 5; channel++)
4403   {
4404     list=pixel_list->lists+channel;
4405     color=65536UL;
4406     next=list->nodes[color].next[0];
4407     count=0;
4408     do
4409     {
4410       previous=color;
4411       color=next;
4412       next=list->nodes[color].next[0];
4413       count+=list->nodes[color].count;
4414     }
4415     while (count <= center);
4416     if ((previous == 65536UL) && (next != 65536UL))
4417       color=next;
4418     else
4419       if ((previous != 65536UL) && (next == 65536UL))
4420         color=previous;
4421     channels[channel]=(unsigned short) color;
4422   }
4423   GetMagickPixelPacket((const Image *) NULL,&pixel);
4424   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
4425   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
4426   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
4427   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
4428   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
4429   return(pixel);
4430 }
4431
4432 MagickExport Image *ReduceNoiseImage(const Image *image,const double radius,
4433   ExceptionInfo *exception)
4434 {
4435 #define ReduceNoiseImageTag  "ReduceNoise/Image"
4436
4437   CacheView
4438     *image_view,
4439     *noise_view;
4440
4441   Image
4442     *noise_image;
4443
4444   MagickBooleanType
4445     status;
4446
4447   MagickOffsetType
4448     progress;
4449
4450   PixelList
4451     **restrict pixel_list;
4452
4453   size_t
4454     width;
4455
4456   ssize_t
4457     y;
4458
4459   /*
4460     Initialize noise image attributes.
4461   */
4462   assert(image != (Image *) NULL);
4463   assert(image->signature == MagickSignature);
4464   if (image->debug != MagickFalse)
4465     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4466   assert(exception != (ExceptionInfo *) NULL);
4467   assert(exception->signature == MagickSignature);
4468   width=GetOptimalKernelWidth2D(radius,0.5);
4469   if ((image->columns < width) || (image->rows < width))
4470     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
4471   noise_image=CloneImage(image,image->columns,image->rows,MagickTrue,
4472     exception);
4473   if (noise_image == (Image *) NULL)
4474     return((Image *) NULL);
4475   if (SetImageStorageClass(noise_image,DirectClass) == MagickFalse)
4476     {
4477       InheritException(exception,&noise_image->exception);
4478       noise_image=DestroyImage(noise_image);
4479       return((Image *) NULL);
4480     }
4481   pixel_list=AcquirePixelListThreadSet(width);
4482   if (pixel_list == (PixelList **) NULL)
4483     {
4484       noise_image=DestroyImage(noise_image);
4485       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4486     }
4487   /*
4488     Reduce noise image.
4489   */
4490   status=MagickTrue;
4491   progress=0;
4492   image_view=AcquireCacheView(image);
4493   noise_view=AcquireCacheView(noise_image);
4494 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4495   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4496 #endif
4497   for (y=0; y < (ssize_t) noise_image->rows; y++)
4498   {
4499     const int
4500       id = GetOpenMPThreadId();
4501
4502     register const IndexPacket
4503       *restrict indexes;
4504
4505     register const PixelPacket
4506       *restrict p;
4507
4508     register IndexPacket
4509       *restrict noise_indexes;
4510
4511     register PixelPacket
4512       *restrict q;
4513
4514     register ssize_t
4515       x;
4516
4517     if (status == MagickFalse)
4518       continue;
4519     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
4520       (width/2L),image->columns+width,width,exception);
4521     q=QueueCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
4522       exception);
4523     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4524       {
4525         status=MagickFalse;
4526         continue;
4527       }
4528     indexes=GetCacheViewVirtualIndexQueue(image_view);
4529     noise_indexes=GetCacheViewAuthenticIndexQueue(noise_view);
4530     for (x=0; x < (ssize_t) noise_image->columns; x++)
4531     {
4532       MagickPixelPacket
4533         pixel;
4534
4535       register const PixelPacket
4536         *restrict r;
4537
4538       register const IndexPacket
4539         *restrict s;
4540
4541       register ssize_t
4542         u,
4543         v;
4544
4545       r=p;
4546       s=indexes+x;
4547       ResetPixelList(pixel_list[id]);
4548       for (v=0; v < (ssize_t) width; v++)
4549       {
4550         for (u=0; u < (ssize_t) width; u++)
4551           InsertPixelList(image,r+u,s+u,pixel_list[id]);
4552         r+=image->columns+width;
4553         s+=image->columns+width;
4554       }
4555       pixel=GetNonpeakPixelList(pixel_list[id]);
4556       SetPixelPacket(noise_image,&pixel,q,noise_indexes+x);
4557       p++;
4558       q++;
4559     }
4560     if (SyncCacheViewAuthenticPixels(noise_view,exception) == MagickFalse)
4561       status=MagickFalse;
4562     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4563       {
4564         MagickBooleanType
4565           proceed;
4566
4567 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4568   #pragma omp critical (MagickCore_ReduceNoiseImage)
4569 #endif
4570         proceed=SetImageProgress(image,ReduceNoiseImageTag,progress++,
4571           image->rows);
4572         if (proceed == MagickFalse)
4573           status=MagickFalse;
4574       }
4575   }
4576   noise_view=DestroyCacheView(noise_view);
4577   image_view=DestroyCacheView(image_view);
4578   pixel_list=DestroyPixelListThreadSet(pixel_list);
4579   return(noise_image);
4580 }
4581 \f
4582 /*
4583 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4584 %                                                                             %
4585 %                                                                             %
4586 %                                                                             %
4587 %     S e l e c t i v e B l u r I m a g e                                     %
4588 %                                                                             %
4589 %                                                                             %
4590 %                                                                             %
4591 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4592 %
4593 %  SelectiveBlurImage() selectively blur pixels within a contrast threshold.
4594 %  It is similar to the unsharpen mask that sharpens everything with contrast
4595 %  above a certain threshold.
4596 %
4597 %  The format of the SelectiveBlurImage method is:
4598 %
4599 %      Image *SelectiveBlurImage(const Image *image,const double radius,
4600 %        const double sigma,const double threshold,ExceptionInfo *exception)
4601 %      Image *SelectiveBlurImageChannel(const Image *image,
4602 %        const ChannelType channel,const double radius,const double sigma,
4603 %        const double threshold,ExceptionInfo *exception)
4604 %
4605 %  A description of each parameter follows:
4606 %
4607 %    o image: the image.
4608 %
4609 %    o channel: the channel type.
4610 %
4611 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4612 %      pixel.
4613 %
4614 %    o sigma: the standard deviation of the Gaussian, in pixels.
4615 %
4616 %    o threshold: only pixels within this contrast threshold are included
4617 %      in the blur operation.
4618 %
4619 %    o exception: return any errors or warnings in this structure.
4620 %
4621 */
4622
4623 static inline MagickBooleanType SelectiveContrast(const PixelPacket *p,
4624   const PixelPacket *q,const double threshold)
4625 {
4626   if (fabs(PixelIntensity(p)-PixelIntensity(q)) < threshold)
4627     return(MagickTrue);
4628   return(MagickFalse);
4629 }
4630
4631 MagickExport Image *SelectiveBlurImage(const Image *image,const double radius,
4632   const double sigma,const double threshold,ExceptionInfo *exception)
4633 {
4634   Image
4635     *blur_image;
4636
4637   blur_image=SelectiveBlurImageChannel(image,DefaultChannels,radius,sigma,
4638     threshold,exception);
4639   return(blur_image);
4640 }
4641
4642 MagickExport Image *SelectiveBlurImageChannel(const Image *image,
4643   const ChannelType channel,const double radius,const double sigma,
4644   const double threshold,ExceptionInfo *exception)
4645 {
4646 #define SelectiveBlurImageTag  "SelectiveBlur/Image"
4647
4648   CacheView
4649     *blur_view,
4650     *image_view;
4651
4652   double
4653     *kernel;
4654
4655   Image
4656     *blur_image;
4657
4658   MagickBooleanType
4659     status;
4660
4661   MagickOffsetType
4662     progress;
4663
4664   MagickPixelPacket
4665     bias;
4666
4667   register ssize_t
4668     i;
4669
4670   size_t
4671     width;
4672
4673   ssize_t
4674     j,
4675     u,
4676     v,
4677     y;
4678
4679   /*
4680     Initialize blur image attributes.
4681   */
4682   assert(image != (Image *) NULL);
4683   assert(image->signature == MagickSignature);
4684   if (image->debug != MagickFalse)
4685     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4686   assert(exception != (ExceptionInfo *) NULL);
4687   assert(exception->signature == MagickSignature);
4688   width=GetOptimalKernelWidth1D(radius,sigma);
4689   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
4690   if (kernel == (double *) NULL)
4691     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4692   j=(ssize_t) width/2;
4693   i=0;
4694   for (v=(-j); v <= j; v++)
4695   {
4696     for (u=(-j); u <= j; u++)
4697       kernel[i++]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
4698         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
4699   }
4700   if (image->debug != MagickFalse)
4701     {
4702       char
4703         format[MaxTextExtent],
4704         *message;
4705
4706       register const double
4707         *k;
4708
4709       ssize_t
4710         u,
4711         v;
4712
4713       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
4714         "  SelectiveBlurImage with %.20gx%.20g kernel:",(double) width,(double)
4715         width);
4716       message=AcquireString("");
4717       k=kernel;
4718       for (v=0; v < (ssize_t) width; v++)
4719       {
4720         *message='\0';
4721         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) v);
4722         (void) ConcatenateString(&message,format);
4723         for (u=0; u < (ssize_t) width; u++)
4724         {
4725           (void) FormatMagickString(format,MaxTextExtent,"%+f ",*k++);
4726           (void) ConcatenateString(&message,format);
4727         }
4728         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
4729       }
4730       message=DestroyString(message);
4731     }
4732   blur_image=CloneImage(image,0,0,MagickTrue,exception);
4733   if (blur_image == (Image *) NULL)
4734     return((Image *) NULL);
4735   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
4736     {
4737       InheritException(exception,&blur_image->exception);
4738       blur_image=DestroyImage(blur_image);
4739       return((Image *) NULL);
4740     }
4741   /*
4742     Threshold blur image.
4743   */
4744   status=MagickTrue;
4745   progress=0;
4746   GetMagickPixelPacket(image,&bias);
4747   SetMagickPixelPacketBias(image,&bias);
4748   image_view=AcquireCacheView(image);
4749   blur_view=AcquireCacheView(blur_image);
4750 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4751   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4752 #endif
4753   for (y=0; y < (ssize_t) image->rows; y++)
4754   {
4755     MagickBooleanType
4756       sync;
4757
4758     MagickRealType
4759       gamma;
4760
4761     register const IndexPacket
4762       *restrict indexes;
4763
4764     register const PixelPacket
4765       *restrict p;
4766
4767     register IndexPacket
4768       *restrict blur_indexes;
4769
4770     register PixelPacket
4771       *restrict q;
4772
4773     register ssize_t
4774       x;
4775
4776     if (status == MagickFalse)
4777       continue;
4778     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
4779       (width/2L),image->columns+width,width,exception);
4780     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
4781       exception);
4782     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4783       {
4784         status=MagickFalse;
4785         continue;
4786       }
4787     indexes=GetCacheViewVirtualIndexQueue(image_view);
4788     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
4789     for (x=0; x < (ssize_t) image->columns; x++)
4790     {
4791       MagickPixelPacket
4792         pixel;
4793
4794       register const double
4795         *restrict k;
4796
4797       register ssize_t
4798         u;
4799
4800       ssize_t
4801         j,
4802         v;
4803
4804       pixel=bias;
4805       k=kernel;
4806       gamma=0.0;
4807       j=0;
4808       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
4809         {
4810           for (v=0; v < (ssize_t) width; v++)
4811           {
4812             for (u=0; u < (ssize_t) width; u++)
4813             {
4814               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4815                 {
4816                   pixel.red+=(*k)*(p+u+j)->red;
4817                   pixel.green+=(*k)*(p+u+j)->green;
4818                   pixel.blue+=(*k)*(p+u+j)->blue;
4819                   gamma+=(*k);
4820                   k++;
4821                 }
4822             }
4823             j+=(ssize_t) (image->columns+width);
4824           }
4825           if (gamma != 0.0)
4826             {
4827               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4828               if ((channel & RedChannel) != 0)
4829                 q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
4830               if ((channel & GreenChannel) != 0)
4831                 q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
4832               if ((channel & BlueChannel) != 0)
4833                 q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
4834             }
4835           if ((channel & OpacityChannel) != 0)
4836             {
4837               gamma=0.0;
4838               j=0;
4839               for (v=0; v < (ssize_t) width; v++)
4840               {
4841                 for (u=0; u < (ssize_t) width; u++)
4842                 {
4843                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4844                     {
4845                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4846                       gamma+=(*k);
4847                       k++;
4848                     }
4849                 }
4850                 j+=(ssize_t) (image->columns+width);
4851               }
4852               if (gamma != 0.0)
4853                 {
4854                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4855                     gamma);
4856                   SetOpacityPixelComponent(q,ClampToQuantum(gamma*
4857                     GetOpacityPixelComponent(&pixel)));
4858                 }
4859             }
4860           if (((channel & IndexChannel) != 0) &&
4861               (image->colorspace == CMYKColorspace))
4862             {
4863               gamma=0.0;
4864               j=0;
4865               for (v=0; v < (ssize_t) width; v++)
4866               {
4867                 for (u=0; u < (ssize_t) width; u++)
4868                 {
4869                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4870                     {
4871                       pixel.index+=(*k)*indexes[x+u+j];
4872                       gamma+=(*k);
4873                       k++;
4874                     }
4875                 }
4876                 j+=(ssize_t) (image->columns+width);
4877               }
4878               if (gamma != 0.0)
4879                 {
4880                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4881                     gamma);
4882                   blur_indexes[x]=ClampToQuantum(gamma*
4883                     GetIndexPixelComponent(&pixel));
4884                 }
4885             }
4886         }
4887       else
4888         {
4889           MagickRealType
4890             alpha;
4891
4892           for (v=0; v < (ssize_t) width; v++)
4893           {
4894             for (u=0; u < (ssize_t) width; u++)
4895             {
4896               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4897                 {
4898                   alpha=(MagickRealType) (QuantumScale*
4899                     GetAlphaPixelComponent(p+u+j));
4900                   pixel.red+=(*k)*alpha*(p+u+j)->red;
4901                   pixel.green+=(*k)*alpha*(p+u+j)->green;
4902                   pixel.blue+=(*k)*alpha*(p+u+j)->blue;
4903                   pixel.opacity+=(*k)*(p+u+j)->opacity;
4904                   gamma+=(*k)*alpha;
4905                   k++;
4906                 }
4907             }
4908             j+=(ssize_t) (image->columns+width);
4909           }
4910           if (gamma != 0.0)
4911             {
4912               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4913               if ((channel & RedChannel) != 0)
4914                 q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
4915               if ((channel & GreenChannel) != 0)
4916                 q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
4917               if ((channel & BlueChannel) != 0)
4918                 q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
4919             }
4920           if ((channel & OpacityChannel) != 0)
4921             {
4922               gamma=0.0;
4923               j=0;
4924               for (v=0; v < (ssize_t) width; v++)
4925               {
4926                 for (u=0; u < (ssize_t) width; u++)
4927                 {
4928                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4929                     {
4930                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4931                       gamma+=(*k);
4932                       k++;
4933                     }
4934                 }
4935                 j+=(ssize_t) (image->columns+width);
4936               }
4937               if (gamma != 0.0)
4938                 {
4939                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4940                     gamma);
4941                   SetOpacityPixelComponent(q,
4942                     ClampOpacityPixelComponent(&pixel));
4943                 }
4944             }
4945           if (((channel & IndexChannel) != 0) &&
4946               (image->colorspace == CMYKColorspace))
4947             {
4948               gamma=0.0;
4949               j=0;
4950               for (v=0; v < (ssize_t) width; v++)
4951               {
4952                 for (u=0; u < (ssize_t) width; u++)
4953                 {
4954                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4955                     {
4956                       alpha=(MagickRealType) (QuantumScale*
4957                         GetAlphaPixelComponent(p+u+j));
4958                       pixel.index+=(*k)*alpha*indexes[x+u+j];
4959                       gamma+=(*k);
4960                       k++;
4961                     }
4962                 }
4963                 j+=(ssize_t) (image->columns+width);
4964               }
4965               if (gamma != 0.0)
4966                 {
4967                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4968                     gamma);
4969                   blur_indexes[x]=ClampToQuantum(gamma*
4970                     GetIndexPixelComponent(&pixel));
4971                 }
4972             }
4973         }
4974       p++;
4975       q++;
4976     }
4977     sync=SyncCacheViewAuthenticPixels(blur_view,exception);
4978     if (sync == MagickFalse)
4979       status=MagickFalse;
4980     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4981       {
4982         MagickBooleanType
4983           proceed;
4984
4985 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4986   #pragma omp critical (MagickCore_SelectiveBlurImageChannel)
4987 #endif
4988         proceed=SetImageProgress(image,SelectiveBlurImageTag,progress++,
4989           image->rows);
4990         if (proceed == MagickFalse)
4991           status=MagickFalse;
4992       }
4993   }
4994   blur_image->type=image->type;
4995   blur_view=DestroyCacheView(blur_view);
4996   image_view=DestroyCacheView(image_view);
4997   kernel=(double *) RelinquishMagickMemory(kernel);
4998   if (status == MagickFalse)
4999     blur_image=DestroyImage(blur_image);
5000   return(blur_image);
5001 }
5002 \f
5003 /*
5004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5005 %                                                                             %
5006 %                                                                             %
5007 %                                                                             %
5008 %     S h a d e I m a g e                                                     %
5009 %                                                                             %
5010 %                                                                             %
5011 %                                                                             %
5012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5013 %
5014 %  ShadeImage() shines a distant light on an image to create a
5015 %  three-dimensional effect. You control the positioning of the light with
5016 %  azimuth and elevation; azimuth is measured in degrees off the x axis
5017 %  and elevation is measured in pixels above the Z axis.
5018 %
5019 %  The format of the ShadeImage method is:
5020 %
5021 %      Image *ShadeImage(const Image *image,const MagickBooleanType gray,
5022 %        const double azimuth,const double elevation,ExceptionInfo *exception)
5023 %
5024 %  A description of each parameter follows:
5025 %
5026 %    o image: the image.
5027 %
5028 %    o gray: A value other than zero shades the intensity of each pixel.
5029 %
5030 %    o azimuth, elevation:  Define the light source direction.
5031 %
5032 %    o exception: return any errors or warnings in this structure.
5033 %
5034 */
5035 MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
5036   const double azimuth,const double elevation,ExceptionInfo *exception)
5037 {
5038 #define ShadeImageTag  "Shade/Image"
5039
5040   CacheView
5041     *image_view,
5042     *shade_view;
5043
5044   Image
5045     *shade_image;
5046
5047   MagickBooleanType
5048     status;
5049
5050   MagickOffsetType
5051     progress;
5052
5053   PrimaryInfo
5054     light;
5055
5056   ssize_t
5057     y;
5058
5059   /*
5060     Initialize shaded image attributes.
5061   */
5062   assert(image != (const Image *) NULL);
5063   assert(image->signature == MagickSignature);
5064   if (image->debug != MagickFalse)
5065     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5066   assert(exception != (ExceptionInfo *) NULL);
5067   assert(exception->signature == MagickSignature);
5068   shade_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
5069   if (shade_image == (Image *) NULL)
5070     return((Image *) NULL);
5071   if (SetImageStorageClass(shade_image,DirectClass) == MagickFalse)
5072     {
5073       InheritException(exception,&shade_image->exception);
5074       shade_image=DestroyImage(shade_image);
5075       return((Image *) NULL);
5076     }
5077   /*
5078     Compute the light vector.
5079   */
5080   light.x=(double) QuantumRange*cos(DegreesToRadians(azimuth))*
5081     cos(DegreesToRadians(elevation));
5082   light.y=(double) QuantumRange*sin(DegreesToRadians(azimuth))*
5083     cos(DegreesToRadians(elevation));
5084   light.z=(double) QuantumRange*sin(DegreesToRadians(elevation));
5085   /*
5086     Shade image.
5087   */
5088   status=MagickTrue;
5089   progress=0;
5090   image_view=AcquireCacheView(image);
5091   shade_view=AcquireCacheView(shade_image);
5092 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5093   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
5094 #endif
5095   for (y=0; y < (ssize_t) image->rows; y++)
5096   {
5097     MagickRealType
5098       distance,
5099       normal_distance,
5100       shade;
5101
5102     PrimaryInfo
5103       normal;
5104
5105     register const PixelPacket
5106       *restrict p,
5107       *restrict s0,
5108       *restrict s1,
5109       *restrict s2;
5110
5111     register PixelPacket
5112       *restrict q;
5113
5114     register ssize_t
5115       x;
5116
5117     if (status == MagickFalse)
5118       continue;
5119     p=GetCacheViewVirtualPixels(image_view,-1,y-1,image->columns+2,3,exception);
5120     q=QueueCacheViewAuthenticPixels(shade_view,0,y,shade_image->columns,1,
5121       exception);
5122     if ((p == (PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
5123       {
5124         status=MagickFalse;
5125         continue;
5126       }
5127     /*
5128       Shade this row of pixels.
5129     */
5130     normal.z=2.0*(double) QuantumRange;  /* constant Z of surface normal */
5131     s0=p+1;
5132     s1=s0+image->columns+2;
5133     s2=s1+image->columns+2;
5134     for (x=0; x < (ssize_t) image->columns; x++)
5135     {
5136       /*
5137         Determine the surface normal and compute shading.
5138       */
5139       normal.x=(double) (PixelIntensity(s0-1)+PixelIntensity(s1-1)+
5140         PixelIntensity(s2-1)-PixelIntensity(s0+1)-PixelIntensity(s1+1)-
5141         PixelIntensity(s2+1));
5142       normal.y=(double) (PixelIntensity(s2-1)+PixelIntensity(s2)+
5143         PixelIntensity(s2+1)-PixelIntensity(s0-1)-PixelIntensity(s0)-
5144         PixelIntensity(s0+1));
5145       if ((normal.x == 0.0) && (normal.y == 0.0))
5146         shade=light.z;
5147       else
5148         {
5149           shade=0.0;
5150           distance=normal.x*light.x+normal.y*light.y+normal.z*light.z;
5151           if (distance > MagickEpsilon)
5152             {
5153               normal_distance=
5154                 normal.x*normal.x+normal.y*normal.y+normal.z*normal.z;
5155               if (normal_distance > (MagickEpsilon*MagickEpsilon))
5156                 shade=distance/sqrt((double) normal_distance);
5157             }
5158         }
5159       if (gray != MagickFalse)
5160         {
5161           q->red=(Quantum) shade;
5162           q->green=(Quantum) shade;
5163           q->blue=(Quantum) shade;
5164         }
5165       else
5166         {
5167           q->red=ClampToQuantum(QuantumScale*shade*s1->red);
5168           q->green=ClampToQuantum(QuantumScale*shade*s1->green);
5169           q->blue=ClampToQuantum(QuantumScale*shade*s1->blue);
5170         }
5171       q->opacity=s1->opacity;
5172       s0++;
5173       s1++;
5174       s2++;
5175       q++;
5176     }
5177     if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
5178       status=MagickFalse;
5179     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5180       {
5181         MagickBooleanType
5182           proceed;
5183
5184 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5185   #pragma omp critical (MagickCore_ShadeImage)
5186 #endif
5187         proceed=SetImageProgress(image,ShadeImageTag,progress++,image->rows);
5188         if (proceed == MagickFalse)
5189           status=MagickFalse;
5190       }
5191   }
5192   shade_view=DestroyCacheView(shade_view);
5193   image_view=DestroyCacheView(image_view);
5194   if (status == MagickFalse)
5195     shade_image=DestroyImage(shade_image);
5196   return(shade_image);
5197 }
5198 \f
5199 /*
5200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5201 %                                                                             %
5202 %                                                                             %
5203 %                                                                             %
5204 %     S h a r p e n I m a g e                                                 %
5205 %                                                                             %
5206 %                                                                             %
5207 %                                                                             %
5208 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5209 %
5210 %  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
5211 %  operator of the given radius and standard deviation (sigma).  For
5212 %  reasonable results, radius should be larger than sigma.  Use a radius of 0
5213 %  and SharpenImage() selects a suitable radius for you.
5214 %
5215 %  Using a separable kernel would be faster, but the negative weights cancel
5216 %  out on the corners of the kernel producing often undesirable ringing in the
5217 %  filtered result; this can be avoided by using a 2D gaussian shaped image
5218 %  sharpening kernel instead.
5219 %
5220 %  The format of the SharpenImage method is:
5221 %
5222 %    Image *SharpenImage(const Image *image,const double radius,
5223 %      const double sigma,ExceptionInfo *exception)
5224 %    Image *SharpenImageChannel(const Image *image,const ChannelType channel,
5225 %      const double radius,const double sigma,ExceptionInfo *exception)
5226 %
5227 %  A description of each parameter follows:
5228 %
5229 %    o image: the image.
5230 %
5231 %    o channel: the channel type.
5232 %
5233 %    o radius: the radius of the Gaussian, in pixels, not counting the center
5234 %      pixel.
5235 %
5236 %    o sigma: the standard deviation of the Laplacian, in pixels.
5237 %
5238 %    o exception: return any errors or warnings in this structure.
5239 %
5240 */
5241
5242 MagickExport Image *SharpenImage(const Image *image,const double radius,
5243   const double sigma,ExceptionInfo *exception)
5244 {
5245   Image
5246     *sharp_image;
5247
5248   sharp_image=SharpenImageChannel(image,DefaultChannels,radius,sigma,exception);
5249   return(sharp_image);
5250 }
5251
5252 MagickExport Image *SharpenImageChannel(const Image *image,
5253   const ChannelType channel,const double radius,const double sigma,
5254   ExceptionInfo *exception)
5255 {
5256   double
5257     *kernel,
5258     normalize;
5259
5260   Image
5261     *sharp_image;
5262
5263   register ssize_t
5264     i;
5265
5266   size_t
5267     width;
5268
5269   ssize_t
5270     j,
5271     u,
5272     v;
5273
5274   assert(image != (const Image *) NULL);
5275   assert(image->signature == MagickSignature);
5276   if (image->debug != MagickFalse)
5277     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5278   assert(exception != (ExceptionInfo *) NULL);
5279   assert(exception->signature == MagickSignature);
5280   width=GetOptimalKernelWidth2D(radius,sigma);
5281   kernel=(double *) AcquireQuantumMemory((size_t) width*width,sizeof(*kernel));
5282   if (kernel == (double *) NULL)
5283     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
5284   normalize=0.0;
5285   j=(ssize_t) width/2;
5286   i=0;
5287   for (v=(-j); v <= j; v++)
5288   {
5289     for (u=(-j); u <= j; u++)
5290     {
5291       kernel[i]=(double) (-exp(-((double) u*u+v*v)/(2.0*MagickSigma*
5292         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
5293       normalize+=kernel[i];
5294       i++;
5295     }
5296   }
5297   kernel[i/2]=(double) ((-2.0)*normalize);
5298   sharp_image=ConvolveImageChannel(image,channel,width,kernel,exception);
5299   kernel=(double *) RelinquishMagickMemory(kernel);
5300   return(sharp_image);
5301 }
5302 \f
5303 /*
5304 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5305 %                                                                             %
5306 %                                                                             %
5307 %                                                                             %
5308 %     S p r e a d I m a g e                                                   %
5309 %                                                                             %
5310 %                                                                             %
5311 %                                                                             %
5312 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5313 %
5314 %  SpreadImage() is a special effects method that randomly displaces each
5315 %  pixel in a block defined by the radius parameter.
5316 %
5317 %  The format of the SpreadImage method is:
5318 %
5319 %      Image *SpreadImage(const Image *image,const double radius,
5320 %        ExceptionInfo *exception)
5321 %
5322 %  A description of each parameter follows:
5323 %
5324 %    o image: the image.
5325 %
5326 %    o radius:  Choose a random pixel in a neighborhood of this extent.
5327 %
5328 %    o exception: return any errors or warnings in this structure.
5329 %
5330 */
5331 MagickExport Image *SpreadImage(const Image *image,const double radius,
5332   ExceptionInfo *exception)
5333 {
5334 #define SpreadImageTag  "Spread/Image"
5335
5336   CacheView
5337     *image_view;
5338
5339   Image
5340     *spread_image;
5341
5342   MagickBooleanType
5343     status;
5344
5345   MagickOffsetType
5346     progress;
5347
5348   MagickPixelPacket
5349     bias;
5350
5351   RandomInfo
5352     **restrict random_info;
5353
5354   ResampleFilter
5355     **restrict resample_filter;
5356
5357   size_t
5358     width;
5359
5360   ssize_t
5361     y;
5362
5363   /*
5364     Initialize spread image attributes.
5365   */
5366   assert(image != (Image *) NULL);
5367   assert(image->signature == MagickSignature);
5368   if (image->debug != MagickFalse)
5369     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5370   assert(exception != (ExceptionInfo *) NULL);
5371   assert(exception->signature == MagickSignature);
5372   spread_image=CloneImage(image,image->columns,image->rows,MagickTrue,
5373     exception);
5374   if (spread_image == (Image *) NULL)
5375     return((Image *) NULL);
5376   if (SetImageStorageClass(spread_image,DirectClass) == MagickFalse)
5377     {
5378       InheritException(exception,&spread_image->exception);
5379       spread_image=DestroyImage(spread_image);
5380       return((Image *) NULL);
5381     }
5382   /*
5383     Spread image.
5384   */
5385   status=MagickTrue;
5386   progress=0;
5387   GetMagickPixelPacket(spread_image,&bias);
5388   width=GetOptimalKernelWidth1D(radius,0.5);
5389   resample_filter=AcquireResampleFilterThreadSet(image,
5390     UndefinedVirtualPixelMethod,MagickTrue,exception);
5391   random_info=AcquireRandomInfoThreadSet();
5392   image_view=AcquireCacheView(spread_image);
5393 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5394   #pragma omp parallel for schedule(dynamic,4) shared(progress,status) omp_throttle(1)
5395 #endif
5396   for (y=0; y < (ssize_t) spread_image->rows; y++)
5397   {
5398     const int
5399       id = GetOpenMPThreadId();
5400
5401     MagickPixelPacket
5402       pixel;
5403
5404     register IndexPacket
5405       *restrict indexes;
5406
5407     register PixelPacket
5408       *restrict q;
5409
5410     register ssize_t
5411       x;
5412
5413     if (status == MagickFalse)
5414       continue;
5415     q=QueueCacheViewAuthenticPixels(image_view,0,y,spread_image->columns,1,
5416       exception);
5417     if (q == (PixelPacket *) NULL)
5418       {
5419         status=MagickFalse;
5420         continue;
5421       }
5422     indexes=GetCacheViewAuthenticIndexQueue(image_view);
5423     pixel=bias;
5424     for (x=0; x < (ssize_t) spread_image->columns; x++)
5425     {
5426       (void) ResamplePixelColor(resample_filter[id],(double) x+width*
5427         (GetPseudoRandomValue(random_info[id])-0.5),(double) y+width*
5428         (GetPseudoRandomValue(random_info[id])-0.5),&pixel);
5429       SetPixelPacket(spread_image,&pixel,q,indexes+x);
5430       q++;
5431     }
5432     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
5433       status=MagickFalse;
5434     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5435       {
5436         MagickBooleanType
5437           proceed;
5438
5439 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5440   #pragma omp critical (MagickCore_SpreadImage)
5441 #endif
5442         proceed=SetImageProgress(image,SpreadImageTag,progress++,image->rows);
5443         if (proceed == MagickFalse)
5444           status=MagickFalse;
5445       }
5446   }
5447   image_view=DestroyCacheView(image_view);
5448   random_info=DestroyRandomInfoThreadSet(random_info);
5449   resample_filter=DestroyResampleFilterThreadSet(resample_filter);
5450   return(spread_image);
5451 }
5452 \f
5453 /*
5454 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5455 %                                                                             %
5456 %                                                                             %
5457 %                                                                             %
5458 %     U n s h a r p M a s k I m a g e                                         %
5459 %                                                                             %
5460 %                                                                             %
5461 %                                                                             %
5462 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5463 %
5464 %  UnsharpMaskImage() sharpens one or more image channels.  We convolve the
5465 %  image with a Gaussian operator of the given radius and standard deviation
5466 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
5467 %  radius of 0 and UnsharpMaskImage() selects a suitable radius for you.
5468 %
5469 %  The format of the UnsharpMaskImage method is:
5470 %
5471 %    Image *UnsharpMaskImage(const Image *image,const double radius,
5472 %      const double sigma,const double amount,const double threshold,
5473 %      ExceptionInfo *exception)
5474 %    Image *UnsharpMaskImageChannel(const Image *image,
5475 %      const ChannelType channel,const double radius,const double sigma,
5476 %      const double amount,const double threshold,ExceptionInfo *exception)
5477 %
5478 %  A description of each parameter follows:
5479 %
5480 %    o image: the image.
5481 %
5482 %    o channel: the channel type.
5483 %
5484 %    o radius: the radius of the Gaussian, in pixels, not counting the center
5485 %      pixel.
5486 %
5487 %    o sigma: the standard deviation of the Gaussian, in pixels.
5488 %
5489 %    o amount: the percentage of the difference between the original and the
5490 %      blur image that is added back into the original.
5491 %
5492 %    o threshold: the threshold in pixels needed to apply the diffence amount.
5493 %
5494 %    o exception: return any errors or warnings in this structure.
5495 %
5496 */
5497
5498 MagickExport Image *UnsharpMaskImage(const Image *image,const double radius,
5499   const double sigma,const double amount,const double threshold,
5500   ExceptionInfo *exception)
5501 {
5502   Image
5503     *sharp_image;
5504
5505   sharp_image=UnsharpMaskImageChannel(image,DefaultChannels,radius,sigma,amount,
5506     threshold,exception);
5507   return(sharp_image);
5508 }
5509
5510 MagickExport Image *UnsharpMaskImageChannel(const Image *image,
5511   const ChannelType channel,const double radius,const double sigma,
5512   const double amount,const double threshold,ExceptionInfo *exception)
5513 {
5514 #define SharpenImageTag  "Sharpen/Image"
5515
5516   CacheView
5517     *image_view,
5518     *unsharp_view;
5519
5520   Image
5521     *unsharp_image;
5522
5523   MagickBooleanType
5524     status;
5525
5526   MagickOffsetType
5527     progress;
5528
5529   MagickPixelPacket
5530     bias;
5531
5532   MagickRealType
5533     quantum_threshold;
5534
5535   ssize_t
5536     y;
5537
5538   assert(image != (const Image *) NULL);
5539   assert(image->signature == MagickSignature);
5540   if (image->debug != MagickFalse)
5541     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5542   assert(exception != (ExceptionInfo *) NULL);
5543   unsharp_image=BlurImageChannel(image,channel,radius,sigma,exception);
5544   if (unsharp_image == (Image *) NULL)
5545     return((Image *) NULL);
5546   quantum_threshold=(MagickRealType) QuantumRange*threshold;
5547   /*
5548     Unsharp-mask image.
5549   */
5550   status=MagickTrue;
5551   progress=0;
5552   GetMagickPixelPacket(image,&bias);
5553   image_view=AcquireCacheView(image);
5554   unsharp_view=AcquireCacheView(unsharp_image);
5555 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5556   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
5557 #endif
5558   for (y=0; y < (ssize_t) image->rows; y++)
5559   {
5560     MagickPixelPacket
5561       pixel;
5562
5563     register const IndexPacket
5564       *restrict indexes;
5565
5566     register const PixelPacket
5567       *restrict p;
5568
5569     register IndexPacket
5570       *restrict unsharp_indexes;
5571
5572     register PixelPacket
5573       *restrict q;
5574
5575     register ssize_t
5576       x;
5577
5578     if (status == MagickFalse)
5579       continue;
5580     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
5581     q=GetCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
5582       exception);
5583     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
5584       {
5585         status=MagickFalse;
5586         continue;
5587       }
5588     indexes=GetCacheViewVirtualIndexQueue(image_view);
5589     unsharp_indexes=GetCacheViewAuthenticIndexQueue(unsharp_view);
5590     pixel=bias;
5591     for (x=0; x < (ssize_t) image->columns; x++)
5592     {
5593       if ((channel & RedChannel) != 0)
5594         {
5595           pixel.red=p->red-(MagickRealType) q->red;
5596           if (fabs(2.0*pixel.red) < quantum_threshold)
5597             pixel.red=(MagickRealType) GetRedPixelComponent(p);
5598           else
5599             pixel.red=(MagickRealType) p->red+(pixel.red*amount);
5600           SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
5601         }
5602       if ((channel & GreenChannel) != 0)
5603         {
5604           pixel.green=p->green-(MagickRealType) q->green;
5605           if (fabs(2.0*pixel.green) < quantum_threshold)
5606             pixel.green=(MagickRealType) GetGreenPixelComponent(p);
5607           else
5608             pixel.green=(MagickRealType) p->green+(pixel.green*amount);
5609           SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
5610         }
5611       if ((channel & BlueChannel) != 0)
5612         {
5613           pixel.blue=p->blue-(MagickRealType) q->blue;
5614           if (fabs(2.0*pixel.blue) < quantum_threshold)
5615             pixel.blue=(MagickRealType) GetBluePixelComponent(p);
5616           else
5617             pixel.blue=(MagickRealType) p->blue+(pixel.blue*amount);
5618           SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
5619         }
5620       if ((channel & OpacityChannel) != 0)
5621         {
5622           pixel.opacity=p->opacity-(MagickRealType) q->opacity;
5623           if (fabs(2.0*pixel.opacity) < quantum_threshold)
5624             pixel.opacity=(MagickRealType) GetOpacityPixelComponent(p);
5625           else
5626             pixel.opacity=p->opacity+(pixel.opacity*amount);
5627           SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
5628         }
5629       if (((channel & IndexChannel) != 0) &&
5630           (image->colorspace == CMYKColorspace))
5631         {
5632           pixel.index=unsharp_indexes[x]-(MagickRealType) indexes[x];
5633           if (fabs(2.0*pixel.index) < quantum_threshold)
5634             pixel.index=(MagickRealType) indexes[x];
5635           else
5636             pixel.index=(MagickRealType) indexes[x]+(pixel.index*amount);
5637           unsharp_indexes[x]=ClampToQuantum(pixel.index);
5638         }
5639       p++;
5640       q++;
5641     }
5642     if (SyncCacheViewAuthenticPixels(unsharp_view,exception) == MagickFalse)
5643       status=MagickFalse;
5644     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5645       {
5646         MagickBooleanType
5647           proceed;
5648
5649 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5650   #pragma omp critical (MagickCore_UnsharpMaskImageChannel)
5651 #endif
5652         proceed=SetImageProgress(image,SharpenImageTag,progress++,image->rows);
5653         if (proceed == MagickFalse)
5654           status=MagickFalse;
5655       }
5656   }
5657   unsharp_image->type=image->type;
5658   unsharp_view=DestroyCacheView(unsharp_view);
5659   image_view=DestroyCacheView(image_view);
5660   if (status == MagickFalse)
5661     unsharp_image=DestroyImage(unsharp_image);
5662   return(unsharp_image);
5663 }