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

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %                   EEEEE  FFFFF  FFFFF  EEEEE  CCCC  TTTTT                   %
   7 %                   E      F      F      E     C        T                     %
   8 %                   EEE    FFF    FFF    EEE   C        T                     %
   9 %                   E      F      F      E     C        T                     %
  10 %                   EEEEE  F      F      EEEEE  CCCC    T                     %
  11 %                                                                             %
  12 %                                                                             %
  13 %                       MagickCore Image Effects Methods                      %
  14 %                                                                             %
  15 %                               Software Design                               %
  16 %                                 John Cristy                                 %
  17 %                                 October 1996                                %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "magick/studio.h"
  44 #include "magick/accelerate.h"
  45 #include "magick/blob.h"
  46 #include "magick/cache-view.h"
  47 #include "magick/color.h"
  48 #include "magick/color-private.h"
  49 #include "magick/colorspace.h"
  50 #include "magick/constitute.h"
  51 #include "magick/decorate.h"
  52 #include "magick/draw.h"
  53 #include "magick/enhance.h"
  54 #include "magick/exception.h"
  55 #include "magick/exception-private.h"
  56 #include "magick/effect.h"
  57 #include "magick/fx.h"
  58 #include "magick/gem.h"
  59 #include "magick/geometry.h"
  60 #include "magick/image-private.h"
  61 #include "magick/list.h"
  62 #include "magick/log.h"
  63 #include "magick/memory_.h"
  64 #include "magick/monitor.h"
  65 #include "magick/monitor-private.h"
  66 #include "magick/montage.h"
  67 #include "magick/morphology.h"
  68 #include "magick/paint.h"
  69 #include "magick/pixel-private.h"
  70 #include "magick/property.h"
  71 #include "magick/quantize.h"
  72 #include "magick/quantum.h"
  73 #include "magick/random_.h"
  74 #include "magick/random-private.h"
  75 #include "magick/resample.h"
  76 #include "magick/resample-private.h"
  77 #include "magick/resize.h"
  78 #include "magick/resource_.h"
  79 #include "magick/segment.h"
  80 #include "magick/shear.h"
  81 #include "magick/signature-private.h"
  82 #include "magick/string_.h"
  83 #include "magick/thread-private.h"
  84 #include "magick/transform.h"
  85 #include "magick/threshold.h"
  86 \f
  87 /*
  88 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  89 %                                                                             %
  90 %                                                                             %
  91 %                                                                             %
  92 %     A d a p t i v e B l u r I m a g e                                       %
  93 %                                                                             %
  94 %                                                                             %
  95 %                                                                             %
  96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  97 %
  98 %  AdaptiveBlurImage() adaptively blurs the image by blurring less
  99 %  intensely near image edges and more intensely far from edges.  We blur the
 100 %  image with a Gaussian operator of the given radius and standard deviation
 101 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
 102 %  radius of 0 and AdaptiveBlurImage() selects a suitable radius for you.
 103 %
 104 %  The format of the AdaptiveBlurImage method is:
 105 %
 106 %      Image *AdaptiveBlurImage(const Image *image,const double radius,
 107 %        const double sigma,ExceptionInfo *exception)
 108 %      Image *AdaptiveBlurImageChannel(const Image *image,
 109 %        const ChannelType channel,double radius,const double sigma,
 110 %        ExceptionInfo *exception)
 111 %
 112 %  A description of each parameter follows:
 113 %
 114 %    o image: the image.
 115 %
 116 %    o channel: the channel type.
 117 %
 118 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 119 %      pixel.
 120 %
 121 %    o sigma: the standard deviation of the Laplacian, in pixels.
 122 %
 123 %    o exception: return any errors or warnings in this structure.
 124 %
 125 */
 126
 127 MagickExport Image *AdaptiveBlurImage(const Image *image,const double radius,
 128   const double sigma,ExceptionInfo *exception)
 129 {
 130   Image
 131     *blur_image;
 132
 133   blur_image=AdaptiveBlurImageChannel(image,DefaultChannels,radius,sigma,
 134     exception);
 135   return(blur_image);
 136 }
 137
 138 MagickExport Image *AdaptiveBlurImageChannel(const Image *image,
 139   const ChannelType channel,const double radius,const double sigma,
 140   ExceptionInfo *exception)
 141 {
 142 #define AdaptiveBlurImageTag  "Convolve/Image"
 143 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
 144
 145   CacheView
 146     *blur_view,
 147     *edge_view,
 148     *image_view;
 149
 150   double
 151     **kernel,
 152     normalize;
 153
 154   Image
 155     *blur_image,
 156     *edge_image,
 157     *gaussian_image;
 158
 159   MagickBooleanType
 160     status;
 161
 162   MagickOffsetType
 163     progress;
 164
 165   MagickPixelPacket
 166     bias;
 167
 168   register ssize_t
 169     i;
 170
 171   size_t
 172     width;
 173
 174   ssize_t
 175     j,
 176     k,
 177     u,
 178     v,
 179     y;
 180
 181   assert(image != (const Image *) NULL);
 182   assert(image->signature == MagickSignature);
 183   if (image->debug != MagickFalse)
 184     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 185   assert(exception != (ExceptionInfo *) NULL);
 186   assert(exception->signature == MagickSignature);
 187   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
 188   if (blur_image == (Image *) NULL)
 189     return((Image *) NULL);
 190   if (fabs(sigma) <= MagickEpsilon)
 191     return(blur_image);
 192   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
 193     {
 194       InheritException(exception,&blur_image->exception);
 195       blur_image=DestroyImage(blur_image);
 196       return((Image *) NULL);
 197     }
 198   /*
 199     Edge detect the image brighness channel, level, blur, and level again.
 200   */
 201   edge_image=EdgeImage(image,radius,exception);
 202   if (edge_image == (Image *) NULL)
 203     {
 204       blur_image=DestroyImage(blur_image);
 205       return((Image *) NULL);
 206     }
 207   (void) LevelImage(edge_image,"20%,95%");
 208   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
 209   if (gaussian_image != (Image *) NULL)
 210     {
 211       edge_image=DestroyImage(edge_image);
 212       edge_image=gaussian_image;
 213     }
 214   (void) LevelImage(edge_image,"10%,95%");
 215   /*
 216     Create a set of kernels from maximum (radius,sigma) to minimum.
 217   */
 218   width=GetOptimalKernelWidth2D(radius,sigma);
 219   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 220   if (kernel == (double **) NULL)
 221     {
 222       edge_image=DestroyImage(edge_image);
 223       blur_image=DestroyImage(blur_image);
 224       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 225     }
 226   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
 227   for (i=0; i < (ssize_t) width; i+=2)
 228   {
 229     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
 230       sizeof(**kernel));
 231     if (kernel[i] == (double *) NULL)
 232       break;
 233     normalize=0.0;
 234     j=(ssize_t) (width-i)/2;
 235     k=0;
 236     for (v=(-j); v <= j; v++)
 237     {
 238       for (u=(-j); u <= j; u++)
 239       {
 240         kernel[i][k]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
 241           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
 242         normalize+=kernel[i][k];
 243         k++;
 244       }
 245     }
 246     if (fabs(normalize) <= MagickEpsilon)
 247       normalize=1.0;
 248     normalize=1.0/normalize;
 249     for (k=0; k < (j*j); k++)
 250       kernel[i][k]=normalize*kernel[i][k];
 251   }
 252   if (i < (ssize_t) width)
 253     {
 254       for (i-=2; i >= 0; i-=2)
 255         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 256       kernel=(double **) RelinquishMagickMemory(kernel);
 257       edge_image=DestroyImage(edge_image);
 258       blur_image=DestroyImage(blur_image);
 259       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 260     }
 261   /*
 262     Adaptively blur image.
 263   */
 264   status=MagickTrue;
 265   progress=0;
 266   GetMagickPixelPacket(image,&bias);
 267   SetMagickPixelPacketBias(image,&bias);
 268   image_view=AcquireCacheView(image);
 269   edge_view=AcquireCacheView(edge_image);
 270   blur_view=AcquireCacheView(blur_image);
 271 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 272   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 273 #endif
 274   for (y=0; y < (ssize_t) blur_image->rows; y++)
 275   {
 276     register const IndexPacket
 277       *restrict indexes;
 278
 279     register const PixelPacket
 280       *restrict p,
 281       *restrict r;
 282
 283     register IndexPacket
 284       *restrict blur_indexes;
 285
 286     register PixelPacket
 287       *restrict q;
 288
 289     register ssize_t
 290       x;
 291
 292     if (status == MagickFalse)
 293       continue;
 294     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
 295     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
 296       exception);
 297     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 298       {
 299         status=MagickFalse;
 300         continue;
 301       }
 302     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
 303     for (x=0; x < (ssize_t) blur_image->columns; x++)
 304     {
 305       MagickPixelPacket
 306         pixel;
 307
 308       MagickRealType
 309         alpha,
 310         gamma;
 311
 312       register const double
 313         *restrict k;
 314
 315       register ssize_t
 316         i,
 317         u,
 318         v;
 319
 320       gamma=0.0;
 321       i=(ssize_t) ceil((double) width*QuantumScale*PixelIntensity(r)-0.5);
 322       if (i < 0)
 323         i=0;
 324       else
 325         if (i > (ssize_t) width)
 326           i=(ssize_t) width;
 327       if ((i & 0x01) != 0)
 328         i--;
 329       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-i)/2L),y-
 330         (ssize_t) ((width-i)/2L),width-i,width-i,exception);
 331       if (p == (const PixelPacket *) NULL)
 332         break;
 333       indexes=GetCacheViewVirtualIndexQueue(image_view);
 334       pixel=bias;
 335       k=kernel[i];
 336       for (v=0; v < (ssize_t) (width-i); v++)
 337       {
 338         for (u=0; u < (ssize_t) (width-i); u++)
 339         {
 340           alpha=1.0;
 341           if (((channel & OpacityChannel) != 0) &&
 342               (image->matte != MagickFalse))
 343             alpha=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(p));
 344           if ((channel & RedChannel) != 0)
 345             pixel.red+=(*k)*alpha*GetRedPixelComponent(p);
 346           if ((channel & GreenChannel) != 0)
 347             pixel.green+=(*k)*alpha*GetGreenPixelComponent(p);
 348           if ((channel & BlueChannel) != 0)
 349             pixel.blue+=(*k)*alpha*GetBluePixelComponent(p);
 350           if ((channel & OpacityChannel) != 0)
 351             pixel.opacity+=(*k)*GetOpacityPixelComponent(p);
 352           if (((channel & IndexChannel) != 0) &&
 353               (image->colorspace == CMYKColorspace))
 354             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
 355           gamma+=(*k)*alpha;
 356           k++;
 357           p++;
 358         }
 359       }
 360       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 361       if ((channel & RedChannel) != 0)
 362         q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
 363       if ((channel & GreenChannel) != 0)
 364         q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
 365       if ((channel & BlueChannel) != 0)
 366         q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
 367       if ((channel & OpacityChannel) != 0)
 368         SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
 369       if (((channel & IndexChannel) != 0) &&
 370           (image->colorspace == CMYKColorspace))
 371         blur_indexes[x]=ClampToQuantum(gamma*GetIndexPixelComponent(&pixel));
 372       q++;
 373       r++;
 374     }
 375     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
 376       status=MagickFalse;
 377     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 378       {
 379         MagickBooleanType
 380           proceed;
 381
 382 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 383   #pragma omp critical (MagickCore_AdaptiveBlurImageChannel)
 384 #endif
 385         proceed=SetImageProgress(image,AdaptiveBlurImageTag,progress++,
 386           image->rows);
 387         if (proceed == MagickFalse)
 388           status=MagickFalse;
 389       }
 390   }
 391   blur_image->type=image->type;
 392   blur_view=DestroyCacheView(blur_view);
 393   edge_view=DestroyCacheView(edge_view);
 394   image_view=DestroyCacheView(image_view);
 395   edge_image=DestroyImage(edge_image);
 396   for (i=0; i < (ssize_t) width;  i+=2)
 397     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 398   kernel=(double **) RelinquishMagickMemory(kernel);
 399   if (status == MagickFalse)
 400     blur_image=DestroyImage(blur_image);
 401   return(blur_image);
 402 }
 403 \f
 404 /*
 405 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 406 %                                                                             %
 407 %                                                                             %
 408 %                                                                             %
 409 %     A d a p t i v e S h a r p e n I m a g e                                 %
 410 %                                                                             %
 411 %                                                                             %
 412 %                                                                             %
 413 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 414 %
 415 %  AdaptiveSharpenImage() adaptively sharpens the image by sharpening more
 416 %  intensely near image edges and less intensely far from edges. We sharpen the
 417 %  image with a Gaussian operator of the given radius and standard deviation
 418 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
 419 %  radius of 0 and AdaptiveSharpenImage() selects a suitable radius for you.
 420 %
 421 %  The format of the AdaptiveSharpenImage method is:
 422 %
 423 %      Image *AdaptiveSharpenImage(const Image *image,const double radius,
 424 %        const double sigma,ExceptionInfo *exception)
 425 %      Image *AdaptiveSharpenImageChannel(const Image *image,
 426 %        const ChannelType channel,double radius,const double sigma,
 427 %        ExceptionInfo *exception)
 428 %
 429 %  A description of each parameter follows:
 430 %
 431 %    o image: the image.
 432 %
 433 %    o channel: the channel type.
 434 %
 435 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 436 %      pixel.
 437 %
 438 %    o sigma: the standard deviation of the Laplacian, in pixels.
 439 %
 440 %    o exception: return any errors or warnings in this structure.
 441 %
 442 */
 443
 444 MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
 445   const double sigma,ExceptionInfo *exception)
 446 {
 447   Image
 448     *sharp_image;
 449
 450   sharp_image=AdaptiveSharpenImageChannel(image,DefaultChannels,radius,sigma,
 451     exception);
 452   return(sharp_image);
 453 }
 454
 455 MagickExport Image *AdaptiveSharpenImageChannel(const Image *image,
 456   const ChannelType channel,const double radius,const double sigma,
 457   ExceptionInfo *exception)
 458 {
 459 #define AdaptiveSharpenImageTag  "Convolve/Image"
 460 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
 461
 462   CacheView
 463     *sharp_view,
 464     *edge_view,
 465     *image_view;
 466
 467   double
 468     **kernel,
 469     normalize;
 470
 471   Image
 472     *sharp_image,
 473     *edge_image,
 474     *gaussian_image;
 475
 476   MagickBooleanType
 477     status;
 478
 479   MagickOffsetType
 480     progress;
 481
 482   MagickPixelPacket
 483     bias;
 484
 485   register ssize_t
 486     i;
 487
 488   size_t
 489     width;
 490
 491   ssize_t
 492     j,
 493     k,
 494     u,
 495     v,
 496     y;
 497
 498   assert(image != (const Image *) NULL);
 499   assert(image->signature == MagickSignature);
 500   if (image->debug != MagickFalse)
 501     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 502   assert(exception != (ExceptionInfo *) NULL);
 503   assert(exception->signature == MagickSignature);
 504   sharp_image=CloneImage(image,0,0,MagickTrue,exception);
 505   if (sharp_image == (Image *) NULL)
 506     return((Image *) NULL);
 507   if (fabs(sigma) <= MagickEpsilon)
 508     return(sharp_image);
 509   if (SetImageStorageClass(sharp_image,DirectClass) == MagickFalse)
 510     {
 511       InheritException(exception,&sharp_image->exception);
 512       sharp_image=DestroyImage(sharp_image);
 513       return((Image *) NULL);
 514     }
 515   /*
 516     Edge detect the image brighness channel, level, sharp, and level again.
 517   */
 518   edge_image=EdgeImage(image,radius,exception);
 519   if (edge_image == (Image *) NULL)
 520     {
 521       sharp_image=DestroyImage(sharp_image);
 522       return((Image *) NULL);
 523     }
 524   (void) LevelImage(edge_image,"20%,95%");
 525   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
 526   if (gaussian_image != (Image *) NULL)
 527     {
 528       edge_image=DestroyImage(edge_image);
 529       edge_image=gaussian_image;
 530     }
 531   (void) LevelImage(edge_image,"10%,95%");
 532   /*
 533     Create a set of kernels from maximum (radius,sigma) to minimum.
 534   */
 535   width=GetOptimalKernelWidth2D(radius,sigma);
 536   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 537   if (kernel == (double **) NULL)
 538     {
 539       edge_image=DestroyImage(edge_image);
 540       sharp_image=DestroyImage(sharp_image);
 541       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 542     }
 543   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
 544   for (i=0; i < (ssize_t) width; i+=2)
 545   {
 546     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
 547       sizeof(**kernel));
 548     if (kernel[i] == (double *) NULL)
 549       break;
 550     normalize=0.0;
 551     j=(ssize_t) (width-i)/2;
 552     k=0;
 553     for (v=(-j); v <= j; v++)
 554     {
 555       for (u=(-j); u <= j; u++)
 556       {
 557         kernel[i][k]=(double) (-exp(-((double) u*u+v*v)/(2.0*MagickSigma*
 558           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
 559         normalize+=kernel[i][k];
 560         k++;
 561       }
 562     }
 563     if (fabs(normalize) <= MagickEpsilon)
 564       normalize=1.0;
 565     normalize=1.0/normalize;
 566     for (k=0; k < (j*j); k++)
 567       kernel[i][k]=normalize*kernel[i][k];
 568   }
 569   if (i < (ssize_t) width)
 570     {
 571       for (i-=2; i >= 0; i-=2)
 572         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 573       kernel=(double **) RelinquishMagickMemory(kernel);
 574       edge_image=DestroyImage(edge_image);
 575       sharp_image=DestroyImage(sharp_image);
 576       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 577     }
 578   /*
 579     Adaptively sharpen image.
 580   */
 581   status=MagickTrue;
 582   progress=0;
 583   GetMagickPixelPacket(image,&bias);
 584   SetMagickPixelPacketBias(image,&bias);
 585   image_view=AcquireCacheView(image);
 586   edge_view=AcquireCacheView(edge_image);
 587   sharp_view=AcquireCacheView(sharp_image);
 588 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 589   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 590 #endif
 591   for (y=0; y < (ssize_t) sharp_image->rows; y++)
 592   {
 593     register const IndexPacket
 594       *restrict indexes;
 595
 596     register const PixelPacket
 597       *restrict p,
 598       *restrict r;
 599
 600     register IndexPacket
 601       *restrict sharp_indexes;
 602
 603     register PixelPacket
 604       *restrict q;
 605
 606     register ssize_t
 607       x;
 608
 609     if (status == MagickFalse)
 610       continue;
 611     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
 612     q=QueueCacheViewAuthenticPixels(sharp_view,0,y,sharp_image->columns,1,
 613       exception);
 614     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 615       {
 616         status=MagickFalse;
 617         continue;
 618       }
 619     sharp_indexes=GetCacheViewAuthenticIndexQueue(sharp_view);
 620     for (x=0; x < (ssize_t) sharp_image->columns; x++)
 621     {
 622       MagickPixelPacket
 623         pixel;
 624
 625       MagickRealType
 626         alpha,
 627         gamma;
 628
 629       register const double
 630         *restrict k;
 631
 632       register ssize_t
 633         i,
 634         u,
 635         v;
 636
 637       gamma=0.0;
 638       i=(ssize_t) ceil((double) width*(QuantumRange-QuantumScale*
 639         PixelIntensity(r))-0.5);
 640       if (i < 0)
 641         i=0;
 642       else
 643         if (i > (ssize_t) width)
 644           i=(ssize_t) width;
 645       if ((i & 0x01) != 0)
 646         i--;
 647       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-i)/2L),y-
 648         (ssize_t) ((width-i)/2L),width-i,width-i,exception);
 649       if (p == (const PixelPacket *) NULL)
 650         break;
 651       indexes=GetCacheViewVirtualIndexQueue(image_view);
 652       k=kernel[i];
 653       pixel=bias;
 654       for (v=0; v < (ssize_t) (width-i); v++)
 655       {
 656         for (u=0; u < (ssize_t) (width-i); u++)
 657         {
 658           alpha=1.0;
 659           if (((channel & OpacityChannel) != 0) &&
 660               (image->matte != MagickFalse))
 661             alpha=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(p));
 662           if ((channel & RedChannel) != 0)
 663             pixel.red+=(*k)*alpha*GetRedPixelComponent(p);
 664           if ((channel & GreenChannel) != 0)
 665             pixel.green+=(*k)*alpha*GetGreenPixelComponent(p);
 666           if ((channel & BlueChannel) != 0)
 667             pixel.blue+=(*k)*alpha*GetBluePixelComponent(p);
 668           if ((channel & OpacityChannel) != 0)
 669             pixel.opacity+=(*k)*GetOpacityPixelComponent(p);
 670           if (((channel & IndexChannel) != 0) &&
 671               (image->colorspace == CMYKColorspace))
 672             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
 673           gamma+=(*k)*alpha;
 674           k++;
 675           p++;
 676         }
 677       }
 678       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 679       if ((channel & RedChannel) != 0)
 680         q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
 681       if ((channel & GreenChannel) != 0)
 682         q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
 683       if ((channel & BlueChannel) != 0)
 684         q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
 685       if ((channel & OpacityChannel) != 0)
 686         SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
 687       if (((channel & IndexChannel) != 0) &&
 688           (image->colorspace == CMYKColorspace))
 689         sharp_indexes[x]=ClampToQuantum(gamma*GetIndexPixelComponent(&pixel));
 690       q++;
 691       r++;
 692     }
 693     if (SyncCacheViewAuthenticPixels(sharp_view,exception) == MagickFalse)
 694       status=MagickFalse;
 695     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 696       {
 697         MagickBooleanType
 698           proceed;
 699
 700 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 701   #pragma omp critical (MagickCore_AdaptiveSharpenImageChannel)
 702 #endif
 703         proceed=SetImageProgress(image,AdaptiveSharpenImageTag,progress++,
 704           image->rows);
 705         if (proceed == MagickFalse)
 706           status=MagickFalse;
 707       }
 708   }
 709   sharp_image->type=image->type;
 710   sharp_view=DestroyCacheView(sharp_view);
 711   edge_view=DestroyCacheView(edge_view);
 712   image_view=DestroyCacheView(image_view);
 713   edge_image=DestroyImage(edge_image);
 714   for (i=0; i < (ssize_t) width;  i+=2)
 715     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
 716   kernel=(double **) RelinquishMagickMemory(kernel);
 717   if (status == MagickFalse)
 718     sharp_image=DestroyImage(sharp_image);
 719   return(sharp_image);
 720 }
 721 \f
 722 /*
 723 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 724 %                                                                             %
 725 %                                                                             %
 726 %                                                                             %
 727 %     B l u r I m a g e                                                       %
 728 %                                                                             %
 729 %                                                                             %
 730 %                                                                             %
 731 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 732 %
 733 %  BlurImage() blurs an image.  We convolve the image with a Gaussian operator
 734 %  of the given radius and standard deviation (sigma).  For reasonable results,
 735 %  the radius should be larger than sigma.  Use a radius of 0 and BlurImage()
 736 %  selects a suitable radius for you.
 737 %
 738 %  BlurImage() differs from GaussianBlurImage() in that it uses a separable
 739 %  kernel which is faster but mathematically equivalent to the non-separable
 740 %  kernel.
 741 %
 742 %  The format of the BlurImage method is:
 743 %
 744 %      Image *BlurImage(const Image *image,const double radius,
 745 %        const double sigma,ExceptionInfo *exception)
 746 %      Image *BlurImageChannel(const Image *image,const ChannelType channel,
 747 %        const double radius,const double sigma,ExceptionInfo *exception)
 748 %
 749 %  A description of each parameter follows:
 750 %
 751 %    o image: the image.
 752 %
 753 %    o channel: the channel type.
 754 %
 755 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 756 %      pixel.
 757 %
 758 %    o sigma: the standard deviation of the Gaussian, in pixels.
 759 %
 760 %    o exception: return any errors or warnings in this structure.
 761 %
 762 */
 763
 764 MagickExport Image *BlurImage(const Image *image,const double radius,
 765   const double sigma,ExceptionInfo *exception)
 766 {
 767   Image
 768     *blur_image;
 769
 770   blur_image=BlurImageChannel(image,DefaultChannels,radius,sigma,exception);
 771   return(blur_image);
 772 }
 773
 774 static double *GetBlurKernel(const size_t width,const double sigma)
 775 {
 776   double
 777     *kernel,
 778     normalize;
 779
 780   register ssize_t
 781     i;
 782
 783   ssize_t
 784     j,
 785     k;
 786
 787   /*
 788     Generate a 1-D convolution kernel.
 789   */
 790   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
 791   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
 792   if (kernel == (double *) NULL)
 793     return(0);
 794   normalize=0.0;
 795   j=(ssize_t) width/2;
 796   i=0;
 797   for (k=(-j); k <= j; k++)
 798   {
 799     kernel[i]=(double) (exp(-((double) k*k)/(2.0*MagickSigma*MagickSigma))/
 800       (MagickSQ2PI*MagickSigma));
 801     normalize+=kernel[i];
 802     i++;
 803   }
 804   for (i=0; i < (ssize_t) width; i++)
 805     kernel[i]/=normalize;
 806   return(kernel);
 807 }
 808
 809 MagickExport Image *BlurImageChannel(const Image *image,
 810   const ChannelType channel,const double radius,const double sigma,
 811   ExceptionInfo *exception)
 812 {
 813 #define BlurImageTag  "Blur/Image"
 814
 815   CacheView
 816     *blur_view,
 817     *image_view;
 818
 819   double
 820     *kernel;
 821
 822   Image
 823     *blur_image;
 824
 825   MagickBooleanType
 826     status;
 827
 828   MagickOffsetType
 829     progress;
 830
 831   MagickPixelPacket
 832     bias;
 833
 834   register ssize_t
 835     i;
 836
 837   size_t
 838     width;
 839
 840   ssize_t
 841     x,
 842     y;
 843
 844   /*
 845     Initialize blur image attributes.
 846   */
 847   assert(image != (Image *) NULL);
 848   assert(image->signature == MagickSignature);
 849   if (image->debug != MagickFalse)
 850     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 851   assert(exception != (ExceptionInfo *) NULL);
 852   assert(exception->signature == MagickSignature);
 853   blur_image=CloneImage(image,0,0,MagickTrue,exception);
 854   if (blur_image == (Image *) NULL)
 855     return((Image *) NULL);
 856   if (fabs(sigma) <= MagickEpsilon)
 857     return(blur_image);
 858   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
 859     {
 860       InheritException(exception,&blur_image->exception);
 861       blur_image=DestroyImage(blur_image);
 862       return((Image *) NULL);
 863     }
 864   width=GetOptimalKernelWidth1D(radius,sigma);
 865   kernel=GetBlurKernel(width,sigma);
 866   if (kernel == (double *) NULL)
 867     {
 868       blur_image=DestroyImage(blur_image);
 869       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 870     }
 871   if (image->debug != MagickFalse)
 872     {
 873       char
 874         format[MaxTextExtent],
 875         *message;
 876
 877       register const double
 878         *k;
 879
 880       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
 881         "  BlurImage with %.20g kernel:",(double) width);
 882       message=AcquireString("");
 883       k=kernel;
 884       for (i=0; i < (ssize_t) width; i++)
 885       {
 886         *message='\0';
 887         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) i);
 888         (void) ConcatenateString(&message,format);
 889         (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
 890         (void) ConcatenateString(&message,format);
 891         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
 892       }
 893       message=DestroyString(message);
 894     }
 895   /*
 896     Blur rows.
 897   */
 898   status=MagickTrue;
 899   progress=0;
 900   GetMagickPixelPacket(image,&bias);
 901   SetMagickPixelPacketBias(image,&bias);
 902   image_view=AcquireCacheView(image);
 903   blur_view=AcquireCacheView(blur_image);
 904 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 905   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 906 #endif
 907   for (y=0; y < (ssize_t) blur_image->rows; y++)
 908   {
 909     register const IndexPacket
 910       *restrict indexes;
 911
 912     register const PixelPacket
 913       *restrict p;
 914
 915     register IndexPacket
 916       *restrict blur_indexes;
 917
 918     register PixelPacket
 919       *restrict q;
 920
 921     register ssize_t
 922       x;
 923
 924     if (status == MagickFalse)
 925       continue;
 926     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y,
 927       image->columns+width,1,exception);
 928     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
 929       exception);
 930     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
 931       {
 932         status=MagickFalse;
 933         continue;
 934       }
 935     indexes=GetCacheViewVirtualIndexQueue(image_view);
 936     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
 937     for (x=0; x < (ssize_t) blur_image->columns; x++)
 938     {
 939       MagickPixelPacket
 940         pixel;
 941
 942       register const double
 943         *restrict k;
 944
 945       register const PixelPacket
 946         *restrict kernel_pixels;
 947
 948       register ssize_t
 949         i;
 950
 951       pixel=bias;
 952       k=kernel;
 953       kernel_pixels=p;
 954       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
 955         {
 956           for (i=0; i < (ssize_t) width; i++)
 957           {
 958             pixel.red+=(*k)*kernel_pixels->red;
 959             pixel.green+=(*k)*kernel_pixels->green;
 960             pixel.blue+=(*k)*kernel_pixels->blue;
 961             k++;
 962             kernel_pixels++;
 963           }
 964           if ((channel & RedChannel) != 0)
 965             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
 966           if ((channel & GreenChannel) != 0)
 967             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
 968           if ((channel & BlueChannel) != 0)
 969             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
 970           if ((channel & OpacityChannel) != 0)
 971             {
 972               k=kernel;
 973               kernel_pixels=p;
 974               for (i=0; i < (ssize_t) width; i++)
 975               {
 976                 pixel.opacity+=(*k)*kernel_pixels->opacity;
 977                 k++;
 978                 kernel_pixels++;
 979               }
 980               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
 981             }
 982           if (((channel & IndexChannel) != 0) &&
 983               (image->colorspace == CMYKColorspace))
 984             {
 985               register const IndexPacket
 986                 *restrict kernel_indexes;
 987
 988               k=kernel;
 989               kernel_indexes=indexes;
 990               for (i=0; i < (ssize_t) width; i++)
 991               {
 992                 pixel.index+=(*k)*(*kernel_indexes);
 993                 k++;
 994                 kernel_indexes++;
 995               }
 996               blur_indexes[x]=ClampToQuantum(pixel.index);
 997             }
 998         }
 999       else
1000         {
1001           MagickRealType
1002             alpha,
1003             gamma;
1004
1005           gamma=0.0;
1006           for (i=0; i < (ssize_t) width; i++)
1007           {
1008             alpha=(MagickRealType) (QuantumScale*
1009               GetAlphaPixelComponent(kernel_pixels));
1010             pixel.red+=(*k)*alpha*kernel_pixels->red;
1011             pixel.green+=(*k)*alpha*kernel_pixels->green;
1012             pixel.blue+=(*k)*alpha*kernel_pixels->blue;
1013             gamma+=(*k)*alpha;
1014             k++;
1015             kernel_pixels++;
1016           }
1017           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1018           if ((channel & RedChannel) != 0)
1019             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
1020           if ((channel & GreenChannel) != 0)
1021             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
1022           if ((channel & BlueChannel) != 0)
1023             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
1024           if ((channel & OpacityChannel) != 0)
1025             {
1026               k=kernel;
1027               kernel_pixels=p;
1028               for (i=0; i < (ssize_t) width; i++)
1029               {
1030                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1031                 k++;
1032                 kernel_pixels++;
1033               }
1034               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1035             }
1036           if (((channel & IndexChannel) != 0) &&
1037               (image->colorspace == CMYKColorspace))
1038             {
1039               register const IndexPacket
1040                 *restrict kernel_indexes;
1041
1042               k=kernel;
1043               kernel_pixels=p;
1044               kernel_indexes=indexes;
1045               for (i=0; i < (ssize_t) width; i++)
1046               {
1047                 alpha=(MagickRealType) (QuantumScale*
1048                   GetAlphaPixelComponent(kernel_pixels));
1049                 pixel.index+=(*k)*alpha*(*kernel_indexes);
1050                 k++;
1051                 kernel_pixels++;
1052                 kernel_indexes++;
1053               }
1054               blur_indexes[x]=ClampToQuantum(gamma*
1055                 GetIndexPixelComponent(&pixel));
1056             }
1057         }
1058       indexes++;
1059       p++;
1060       q++;
1061     }
1062     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
1063       status=MagickFalse;
1064     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1065       {
1066         MagickBooleanType
1067           proceed;
1068
1069 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1070   #pragma omp critical (MagickCore_BlurImageChannel)
1071 #endif
1072         proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
1073           blur_image->columns);
1074         if (proceed == MagickFalse)
1075           status=MagickFalse;
1076       }
1077   }
1078   blur_view=DestroyCacheView(blur_view);
1079   image_view=DestroyCacheView(image_view);
1080   /*
1081     Blur columns.
1082   */
1083   image_view=AcquireCacheView(blur_image);
1084   blur_view=AcquireCacheView(blur_image);
1085 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1086   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1087 #endif
1088   for (x=0; x < (ssize_t) blur_image->columns; x++)
1089   {
1090     register const IndexPacket
1091       *restrict indexes;
1092
1093     register const PixelPacket
1094       *restrict p;
1095
1096     register IndexPacket
1097       *restrict blur_indexes;
1098
1099     register PixelPacket
1100       *restrict q;
1101
1102     register ssize_t
1103       y;
1104
1105     if (status == MagickFalse)
1106       continue;
1107     p=GetCacheViewVirtualPixels(image_view,x,-((ssize_t) width/2L),1,
1108       image->rows+width,exception);
1109     q=GetCacheViewAuthenticPixels(blur_view,x,0,1,blur_image->rows,exception);
1110     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1111       {
1112         status=MagickFalse;
1113         continue;
1114       }
1115     indexes=GetCacheViewVirtualIndexQueue(image_view);
1116     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
1117     for (y=0; y < (ssize_t) blur_image->rows; y++)
1118     {
1119       MagickPixelPacket
1120         pixel;
1121
1122       register const double
1123         *restrict k;
1124
1125       register const PixelPacket
1126         *restrict kernel_pixels;
1127
1128       register ssize_t
1129         i;
1130
1131       pixel=bias;
1132       k=kernel;
1133       kernel_pixels=p;
1134       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
1135         {
1136           for (i=0; i < (ssize_t) width; i++)
1137           {
1138             pixel.red+=(*k)*kernel_pixels->red;
1139             pixel.green+=(*k)*kernel_pixels->green;
1140             pixel.blue+=(*k)*kernel_pixels->blue;
1141             k++;
1142             kernel_pixels++;
1143           }
1144           if ((channel & RedChannel) != 0)
1145             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
1146           if ((channel & GreenChannel) != 0)
1147             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
1148           if ((channel & BlueChannel) != 0)
1149             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
1150           if ((channel & OpacityChannel) != 0)
1151             {
1152               k=kernel;
1153               kernel_pixels=p;
1154               for (i=0; i < (ssize_t) width; i++)
1155               {
1156                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1157                 k++;
1158                 kernel_pixels++;
1159               }
1160               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1161             }
1162           if (((channel & IndexChannel) != 0) &&
1163               (image->colorspace == CMYKColorspace))
1164             {
1165               register const IndexPacket
1166                 *restrict kernel_indexes;
1167
1168               k=kernel;
1169               kernel_indexes=indexes;
1170               for (i=0; i < (ssize_t) width; i++)
1171               {
1172                 pixel.index+=(*k)*(*kernel_indexes);
1173                 k++;
1174                 kernel_indexes++;
1175               }
1176               blur_indexes[y]=ClampToQuantum(pixel.index);
1177             }
1178         }
1179       else
1180         {
1181           MagickRealType
1182             alpha,
1183             gamma;
1184
1185           gamma=0.0;
1186           for (i=0; i < (ssize_t) width; i++)
1187           {
1188             alpha=(MagickRealType) (QuantumScale*
1189               GetAlphaPixelComponent(kernel_pixels));
1190             pixel.red+=(*k)*alpha*kernel_pixels->red;
1191             pixel.green+=(*k)*alpha*kernel_pixels->green;
1192             pixel.blue+=(*k)*alpha*kernel_pixels->blue;
1193             gamma+=(*k)*alpha;
1194             k++;
1195             kernel_pixels++;
1196           }
1197           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1198           if ((channel & RedChannel) != 0)
1199             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
1200           if ((channel & GreenChannel) != 0)
1201             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
1202           if ((channel & BlueChannel) != 0)
1203             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
1204           if ((channel & OpacityChannel) != 0)
1205             {
1206               k=kernel;
1207               kernel_pixels=p;
1208               for (i=0; i < (ssize_t) width; i++)
1209               {
1210                 pixel.opacity+=(*k)*kernel_pixels->opacity;
1211                 k++;
1212                 kernel_pixels++;
1213               }
1214               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1215             }
1216           if (((channel & IndexChannel) != 0) &&
1217               (image->colorspace == CMYKColorspace))
1218             {
1219               register const IndexPacket
1220                 *restrict kernel_indexes;
1221
1222               k=kernel;
1223               kernel_pixels=p;
1224               kernel_indexes=indexes;
1225               for (i=0; i < (ssize_t) width; i++)
1226               {
1227                 alpha=(MagickRealType) (QuantumScale*
1228                   GetAlphaPixelComponent(kernel_pixels));
1229                 pixel.index+=(*k)*alpha*(*kernel_indexes);
1230                 k++;
1231                 kernel_pixels++;
1232                 kernel_indexes++;
1233               }
1234               blur_indexes[y]=ClampToQuantum(gamma*
1235                 GetIndexPixelComponent(&pixel));
1236             }
1237         }
1238       indexes++;
1239       p++;
1240       q++;
1241     }
1242     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
1243       status=MagickFalse;
1244     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1245       {
1246         MagickBooleanType
1247           proceed;
1248
1249 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1250   #pragma omp critical (MagickCore_BlurImageChannel)
1251 #endif
1252         proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
1253           blur_image->columns);
1254         if (proceed == MagickFalse)
1255           status=MagickFalse;
1256       }
1257   }
1258   blur_view=DestroyCacheView(blur_view);
1259   image_view=DestroyCacheView(image_view);
1260   kernel=(double *) RelinquishMagickMemory(kernel);
1261   if (status == MagickFalse)
1262     blur_image=DestroyImage(blur_image);
1263   blur_image->type=image->type;
1264   return(blur_image);
1265 }
1266 \f
1267 /*
1268 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1269 %                                                                             %
1270 %                                                                             %
1271 %                                                                             %
1272 %     C o n v o l v e I m a g e                                               %
1273 %                                                                             %
1274 %                                                                             %
1275 %                                                                             %
1276 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1277 %
1278 %  ConvolveImage() applies a custom convolution kernel to the image.
1279 %
1280 %  The format of the ConvolveImage method is:
1281 %
1282 %      Image *ConvolveImage(const Image *image,const size_t order,
1283 %        const double *kernel,ExceptionInfo *exception)
1284 %      Image *ConvolveImageChannel(const Image *image,const ChannelType channel,
1285 %        const size_t order,const double *kernel,ExceptionInfo *exception)
1286 %
1287 %  A description of each parameter follows:
1288 %
1289 %    o image: the image.
1290 %
1291 %    o channel: the channel type.
1292 %
1293 %    o order: the number of columns and rows in the filter kernel.
1294 %
1295 %    o kernel: An array of double representing the convolution kernel.
1296 %
1297 %    o exception: return any errors or warnings in this structure.
1298 %
1299 */
1300
1301 MagickExport Image *ConvolveImage(const Image *image,const size_t order,
1302   const double *kernel,ExceptionInfo *exception)
1303 {
1304   Image
1305     *convolve_image;
1306
1307   convolve_image=ConvolveImageChannel(image,DefaultChannels,order,kernel,
1308     exception);
1309   return(convolve_image);
1310 }
1311
1312 MagickExport Image *ConvolveImageChannel(const Image *image,
1313   const ChannelType channel,const size_t order,const double *kernel,
1314   ExceptionInfo *exception)
1315 {
1316 #define ConvolveImageTag  "Convolve/Image"
1317
1318   CacheView
1319     *convolve_view,
1320     *image_view;
1321
1322   double
1323     *normal_kernel;
1324
1325   Image
1326     *convolve_image;
1327
1328   MagickBooleanType
1329     status;
1330
1331   MagickOffsetType
1332     progress;
1333
1334   MagickPixelPacket
1335     bias;
1336
1337   MagickRealType
1338     gamma;
1339
1340   register ssize_t
1341     i;
1342
1343   size_t
1344     width;
1345
1346   ssize_t
1347     y;
1348
1349   /*
1350     Initialize convolve image attributes.
1351   */
1352   assert(image != (Image *) NULL);
1353   assert(image->signature == MagickSignature);
1354   if (image->debug != MagickFalse)
1355     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1356   assert(exception != (ExceptionInfo *) NULL);
1357   assert(exception->signature == MagickSignature);
1358   width=order;
1359   if ((width % 2) == 0)
1360     ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
1361   convolve_image=CloneImage(image,0,0,MagickTrue,exception);
1362   if (convolve_image == (Image *) NULL)
1363     return((Image *) NULL);
1364   if (SetImageStorageClass(convolve_image,DirectClass) == MagickFalse)
1365     {
1366       InheritException(exception,&convolve_image->exception);
1367       convolve_image=DestroyImage(convolve_image);
1368       return((Image *) NULL);
1369     }
1370   if (image->debug != MagickFalse)
1371     {
1372       char
1373         format[MaxTextExtent],
1374         *message;
1375
1376       register const double
1377         *k;
1378
1379       ssize_t
1380         u,
1381         v;
1382
1383       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1384         "  ConvolveImage with %.20gx%.20g kernel:",(double) width,(double)
1385         width);
1386       message=AcquireString("");
1387       k=kernel;
1388       for (v=0; v < (ssize_t) width; v++)
1389       {
1390         *message='\0';
1391         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) v);
1392         (void) ConcatenateString(&message,format);
1393         for (u=0; u < (ssize_t) width; u++)
1394         {
1395           (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
1396           (void) ConcatenateString(&message,format);
1397         }
1398         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
1399       }
1400       message=DestroyString(message);
1401     }
1402   /*
1403     Normalize kernel.
1404   */
1405   normal_kernel=(double *) AcquireQuantumMemory(width*width,
1406     sizeof(*normal_kernel));
1407   if (normal_kernel == (double *) NULL)
1408     {
1409       convolve_image=DestroyImage(convolve_image);
1410       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1411     }
1412   gamma=0.0;
1413   for (i=0; i < (ssize_t) (width*width); i++)
1414     gamma+=kernel[i];
1415   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1416   for (i=0; i < (ssize_t) (width*width); i++)
1417     normal_kernel[i]=gamma*kernel[i];
1418   /*
1419     Convolve image.
1420   */
1421   status=MagickTrue;
1422   progress=0;
1423   GetMagickPixelPacket(image,&bias);
1424   SetMagickPixelPacketBias(image,&bias);
1425   image_view=AcquireCacheView(image);
1426   convolve_view=AcquireCacheView(convolve_image);
1427 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1428   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1429 #endif
1430   for (y=0; y < (ssize_t) image->rows; y++)
1431   {
1432     MagickBooleanType
1433       sync;
1434
1435     register const IndexPacket
1436       *restrict indexes;
1437
1438     register const PixelPacket
1439       *restrict p;
1440
1441     register IndexPacket
1442       *restrict convolve_indexes;
1443
1444     register PixelPacket
1445       *restrict q;
1446
1447     register ssize_t
1448       x;
1449
1450     if (status == MagickFalse)
1451       continue;
1452     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
1453       (width/2L),image->columns+width,width,exception);
1454     q=GetCacheViewAuthenticPixels(convolve_view,0,y,convolve_image->columns,1,
1455       exception);
1456     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1457       {
1458         status=MagickFalse;
1459         continue;
1460       }
1461     indexes=GetCacheViewVirtualIndexQueue(image_view);
1462     convolve_indexes=GetCacheViewAuthenticIndexQueue(convolve_view);
1463     for (x=0; x < (ssize_t) image->columns; x++)
1464     {
1465       MagickPixelPacket
1466         pixel;
1467
1468       register const double
1469         *restrict k;
1470
1471       register const PixelPacket
1472         *restrict kernel_pixels;
1473
1474       register ssize_t
1475         u;
1476
1477       ssize_t
1478         v;
1479
1480       pixel=bias;
1481       k=normal_kernel;
1482       kernel_pixels=p;
1483       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
1484         {
1485           for (v=0; v < (ssize_t) width; v++)
1486           {
1487             for (u=0; u < (ssize_t) width; u++)
1488             {
1489               pixel.red+=(*k)*kernel_pixels[u].red;
1490               pixel.green+=(*k)*kernel_pixels[u].green;
1491               pixel.blue+=(*k)*kernel_pixels[u].blue;
1492               k++;
1493             }
1494             kernel_pixels+=image->columns+width;
1495           }
1496           if ((channel & RedChannel) != 0)
1497             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
1498           if ((channel & GreenChannel) != 0)
1499             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
1500           if ((channel & BlueChannel) != 0)
1501             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
1502           if ((channel & OpacityChannel) != 0)
1503             {
1504               k=normal_kernel;
1505               kernel_pixels=p;
1506               for (v=0; v < (ssize_t) width; v++)
1507               {
1508                 for (u=0; u < (ssize_t) width; u++)
1509                 {
1510                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
1511                   k++;
1512                 }
1513                 kernel_pixels+=image->columns+width;
1514               }
1515               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1516             }
1517           if (((channel & IndexChannel) != 0) &&
1518               (image->colorspace == CMYKColorspace))
1519             {
1520               register const IndexPacket
1521                 *restrict kernel_indexes;
1522
1523               k=normal_kernel;
1524               kernel_indexes=indexes;
1525               for (v=0; v < (ssize_t) width; v++)
1526               {
1527                 for (u=0; u < (ssize_t) width; u++)
1528                 {
1529                   pixel.index+=(*k)*kernel_indexes[u];
1530                   k++;
1531                 }
1532                 kernel_indexes+=image->columns+width;
1533               }
1534               convolve_indexes[x]=ClampToQuantum(pixel.index);
1535             }
1536         }
1537       else
1538         {
1539           MagickRealType
1540             alpha,
1541             gamma;
1542
1543           gamma=0.0;
1544           for (v=0; v < (ssize_t) width; v++)
1545           {
1546             for (u=0; u < (ssize_t) width; u++)
1547             {
1548               alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1549                 kernel_pixels[u].opacity));
1550               pixel.red+=(*k)*alpha*kernel_pixels[u].red;
1551               pixel.green+=(*k)*alpha*kernel_pixels[u].green;
1552               pixel.blue+=(*k)*alpha*kernel_pixels[u].blue;
1553               gamma+=(*k)*alpha;
1554               k++;
1555             }
1556             kernel_pixels+=image->columns+width;
1557           }
1558           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1559           if ((channel & RedChannel) != 0)
1560             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
1561           if ((channel & GreenChannel) != 0)
1562             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
1563           if ((channel & BlueChannel) != 0)
1564             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
1565           if ((channel & OpacityChannel) != 0)
1566             {
1567               k=normal_kernel;
1568               kernel_pixels=p;
1569               for (v=0; v < (ssize_t) width; v++)
1570               {
1571                 for (u=0; u < (ssize_t) width; u++)
1572                 {
1573                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
1574                   k++;
1575                 }
1576                 kernel_pixels+=image->columns+width;
1577               }
1578               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
1579             }
1580           if (((channel & IndexChannel) != 0) &&
1581               (image->colorspace == CMYKColorspace))
1582             {
1583               register const IndexPacket
1584                 *restrict kernel_indexes;
1585
1586               k=normal_kernel;
1587               kernel_pixels=p;
1588               kernel_indexes=indexes;
1589               for (v=0; v < (ssize_t) width; v++)
1590               {
1591                 for (u=0; u < (ssize_t) width; u++)
1592                 {
1593                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
1594                     kernel_pixels[u].opacity));
1595                   pixel.index+=(*k)*alpha*kernel_indexes[u];
1596                   k++;
1597                 }
1598                 kernel_pixels+=image->columns+width;
1599                 kernel_indexes+=image->columns+width;
1600               }
1601               convolve_indexes[x]=ClampToQuantum(gamma*
1602                 GetIndexPixelComponent(&pixel));
1603             }
1604         }
1605       indexes++;
1606       p++;
1607       q++;
1608     }
1609     sync=SyncCacheViewAuthenticPixels(convolve_view,exception);
1610     if (sync == MagickFalse)
1611       status=MagickFalse;
1612     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1613       {
1614         MagickBooleanType
1615           proceed;
1616
1617 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1618   #pragma omp critical (MagickCore_ConvolveImageChannel)
1619 #endif
1620         proceed=SetImageProgress(image,ConvolveImageTag,progress++,image->rows);
1621         if (proceed == MagickFalse)
1622           status=MagickFalse;
1623       }
1624   }
1625   convolve_image->type=image->type;
1626   convolve_view=DestroyCacheView(convolve_view);
1627   image_view=DestroyCacheView(image_view);
1628   normal_kernel=(double *) RelinquishMagickMemory(normal_kernel);
1629   if (status == MagickFalse)
1630     convolve_image=DestroyImage(convolve_image);
1631   return(convolve_image);
1632 }
1633 \f
1634 /*
1635 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1636 %                                                                             %
1637 %                                                                             %
1638 %                                                                             %
1639 %     D e s p e c k l e I m a g e                                             %
1640 %                                                                             %
1641 %                                                                             %
1642 %                                                                             %
1643 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1644 %
1645 %  DespeckleImage() reduces the speckle noise in an image while perserving the
1646 %  edges of the original image.
1647 %
1648 %  The format of the DespeckleImage method is:
1649 %
1650 %      Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
1651 %
1652 %  A description of each parameter follows:
1653 %
1654 %    o image: the image.
1655 %
1656 %    o exception: return any errors or warnings in this structure.
1657 %
1658 */
1659
1660 static void Hull(const ssize_t x_offset,const ssize_t y_offset,
1661   const size_t columns,const size_t rows,Quantum *f,Quantum *g,
1662   const int polarity)
1663 {
1664   MagickRealType
1665     v;
1666
1667   register Quantum
1668     *p,
1669     *q,
1670     *r,
1671     *s;
1672
1673   register ssize_t
1674     x;
1675
1676   ssize_t
1677     y;
1678
1679   assert(f != (Quantum *) NULL);
1680   assert(g != (Quantum *) NULL);
1681   p=f+(columns+2);
1682   q=g+(columns+2);
1683   r=p+(y_offset*((ssize_t) columns+2)+x_offset);
1684   for (y=0; y < (ssize_t) rows; y++)
1685   {
1686     p++;
1687     q++;
1688     r++;
1689     if (polarity > 0)
1690       for (x=(ssize_t) columns; x != 0; x--)
1691       {
1692         v=(MagickRealType) (*p);
1693         if ((MagickRealType) *r >= (v+(MagickRealType) ScaleCharToQuantum(2)))
1694           v+=ScaleCharToQuantum(1);
1695         *q=(Quantum) v;
1696         p++;
1697         q++;
1698         r++;
1699       }
1700     else
1701       for (x=(ssize_t) columns; x != 0; x--)
1702       {
1703         v=(MagickRealType) (*p);
1704         if ((MagickRealType) *r <= (v-(MagickRealType) ScaleCharToQuantum(2)))
1705           v-=(ssize_t) ScaleCharToQuantum(1);
1706         *q=(Quantum) v;
1707         p++;
1708         q++;
1709         r++;
1710       }
1711     p++;
1712     q++;
1713     r++;
1714   }
1715   p=f+(columns+2);
1716   q=g+(columns+2);
1717   r=q+(y_offset*((ssize_t) columns+2)+x_offset);
1718   s=q-(y_offset*((ssize_t) columns+2)+x_offset);
1719   for (y=0; y < (ssize_t) rows; y++)
1720   {
1721     p++;
1722     q++;
1723     r++;
1724     s++;
1725     if (polarity > 0)
1726       for (x=(ssize_t) columns; x != 0; x--)
1727       {
1728         v=(MagickRealType) (*q);
1729         if (((MagickRealType) *s >=
1730              (v+(MagickRealType) ScaleCharToQuantum(2))) &&
1731             ((MagickRealType) *r > v))
1732           v+=ScaleCharToQuantum(1);
1733         *p=(Quantum) v;
1734         p++;
1735         q++;
1736         r++;
1737         s++;
1738       }
1739     else
1740       for (x=(ssize_t) columns; x != 0; x--)
1741       {
1742         v=(MagickRealType) (*q);
1743         if (((MagickRealType) *s <=
1744              (v-(MagickRealType) ScaleCharToQuantum(2))) &&
1745             ((MagickRealType) *r < v))
1746           v-=(MagickRealType) ScaleCharToQuantum(1);
1747         *p=(Quantum) v;
1748         p++;
1749         q++;
1750         r++;
1751         s++;
1752       }
1753     p++;
1754     q++;
1755     r++;
1756     s++;
1757   }
1758 }
1759
1760 MagickExport Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
1761 {
1762 #define DespeckleImageTag  "Despeckle/Image"
1763
1764   CacheView
1765     *despeckle_view,
1766     *image_view;
1767
1768   Image
1769     *despeckle_image;
1770
1771   MagickBooleanType
1772     status;
1773
1774   register ssize_t
1775     i;
1776
1777   Quantum
1778     *restrict buffers,
1779     *restrict pixels;
1780
1781   size_t
1782     length,
1783     number_channels;
1784
1785   static const ssize_t
1786     X[4] = {0, 1, 1,-1},
1787     Y[4] = {1, 0, 1, 1};
1788
1789   /*
1790     Allocate despeckled image.
1791   */
1792   assert(image != (const Image *) NULL);
1793   assert(image->signature == MagickSignature);
1794   if (image->debug != MagickFalse)
1795     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1796   assert(exception != (ExceptionInfo *) NULL);
1797   assert(exception->signature == MagickSignature);
1798   despeckle_image=CloneImage(image,image->columns,image->rows,MagickTrue,
1799     exception);
1800   if (despeckle_image == (Image *) NULL)
1801     return((Image *) NULL);
1802   if (SetImageStorageClass(despeckle_image,DirectClass) == MagickFalse)
1803     {
1804       InheritException(exception,&despeckle_image->exception);
1805       despeckle_image=DestroyImage(despeckle_image);
1806       return((Image *) NULL);
1807     }
1808   /*
1809     Allocate image buffers.
1810   */
1811   length=(size_t) ((image->columns+2)*(image->rows+2));
1812   pixels=(Quantum *) AcquireQuantumMemory(length,2*sizeof(*pixels));
1813   buffers=(Quantum *) AcquireQuantumMemory(length,2*sizeof(*pixels));
1814   if ((pixels == (Quantum *) NULL) || (buffers == (Quantum *) NULL))
1815     {
1816       if (buffers != (Quantum *) NULL)
1817         buffers=(Quantum *) RelinquishMagickMemory(buffers);
1818       if (pixels != (Quantum *) NULL)
1819         pixels=(Quantum *) RelinquishMagickMemory(pixels);
1820       despeckle_image=DestroyImage(despeckle_image);
1821       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1822     }
1823   /*
1824     Reduce speckle in the image.
1825   */
1826   status=MagickTrue;
1827   number_channels=(size_t) (image->colorspace == CMYKColorspace ? 5 : 4);
1828   image_view=AcquireCacheView(image);
1829   despeckle_view=AcquireCacheView(despeckle_image);
1830   for (i=0; i < (ssize_t) number_channels; i++)
1831   {
1832     register Quantum
1833       *buffer,
1834       *pixel;
1835
1836     register ssize_t
1837       k,
1838       x;
1839
1840     ssize_t
1841       j,
1842       y;
1843
1844     if (status == MagickFalse)
1845       continue;
1846     pixel=pixels;
1847     (void) ResetMagickMemory(pixel,0,length*sizeof(*pixel));
1848     buffer=buffers;
1849     j=(ssize_t) image->columns+2;
1850     for (y=0; y < (ssize_t) image->rows; y++)
1851     {
1852       register const IndexPacket
1853         *restrict indexes;
1854
1855       register const PixelPacket
1856         *restrict p;
1857
1858       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1859       if (p == (const PixelPacket *) NULL)
1860         break;
1861       indexes=GetCacheViewVirtualIndexQueue(image_view);
1862       j++;
1863       for (x=0; x < (ssize_t) image->columns; x++)
1864       {
1865         switch (i)
1866         {
1867           case 0: pixel[j]=GetRedPixelComponent(p); break;
1868           case 1: pixel[j]=GetGreenPixelComponent(p); break;
1869           case 2: pixel[j]=GetBluePixelComponent(p); break;
1870           case 3: pixel[j]=GetOpacityPixelComponent(p); break;
1871           case 4: pixel[j]=GetBlackPixelComponent(indexes,x); break;
1872           default: break;
1873         }
1874         p++;
1875         j++;
1876       }
1877       j++;
1878     }
1879     (void) ResetMagickMemory(buffer,0,length*sizeof(*buffer));
1880     for (k=0; k < 4; k++)
1881     {
1882       Hull(X[k],Y[k],image->columns,image->rows,pixel,buffer,1);
1883       Hull(-X[k],-Y[k],image->columns,image->rows,pixel,buffer,1);
1884       Hull(-X[k],-Y[k],image->columns,image->rows,pixel,buffer,-1);
1885       Hull(X[k],Y[k],image->columns,image->rows,pixel,buffer,-1);
1886     }
1887     j=(ssize_t) image->columns+2;
1888     for (y=0; y < (ssize_t) image->rows; y++)
1889     {
1890       MagickBooleanType
1891         sync;
1892
1893       register IndexPacket
1894         *restrict indexes;
1895
1896       register PixelPacket
1897         *restrict q;
1898
1899       q=GetCacheViewAuthenticPixels(despeckle_view,0,y,despeckle_image->columns,
1900         1,exception);
1901       if (q == (PixelPacket *) NULL)
1902         break;
1903       indexes=GetCacheViewAuthenticIndexQueue(image_view);
1904       j++;
1905       for (x=0; x < (ssize_t) image->columns; x++)
1906       {
1907         switch (i)
1908         {
1909           case 0: q->red=pixel[j]; break;
1910           case 1: q->green=pixel[j]; break;
1911           case 2: q->blue=pixel[j]; break;
1912           case 3: q->opacity=pixel[j]; break;
1913           case 4: indexes[x]=pixel[j]; break;
1914           default: break;
1915         }
1916         q++;
1917         j++;
1918       }
1919       sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
1920       if (sync == MagickFalse)
1921         {
1922           status=MagickFalse;
1923           break;
1924         }
1925       j++;
1926     }
1927     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1928       {
1929         MagickBooleanType
1930           proceed;
1931
1932         proceed=SetImageProgress(image,DespeckleImageTag,(MagickOffsetType) i,
1933           number_channels);
1934         if (proceed == MagickFalse)
1935           status=MagickFalse;
1936       }
1937   }
1938   despeckle_view=DestroyCacheView(despeckle_view);
1939   image_view=DestroyCacheView(image_view);
1940   buffers=(Quantum *) RelinquishMagickMemory(buffers);
1941   pixels=(Quantum *) RelinquishMagickMemory(pixels);
1942   despeckle_image->type=image->type;
1943   if (status == MagickFalse)
1944     despeckle_image=DestroyImage(despeckle_image);
1945   return(despeckle_image);
1946 }
1947 \f
1948 /*
1949 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1950 %                                                                             %
1951 %                                                                             %
1952 %                                                                             %
1953 %     E d g e I m a g e                                                       %
1954 %                                                                             %
1955 %                                                                             %
1956 %                                                                             %
1957 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1958 %
1959 %  EdgeImage() finds edges in an image.  Radius defines the radius of the
1960 %  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
1961 %  radius for you.
1962 %
1963 %  The format of the EdgeImage method is:
1964 %
1965 %      Image *EdgeImage(const Image *image,const double radius,
1966 %        ExceptionInfo *exception)
1967 %
1968 %  A description of each parameter follows:
1969 %
1970 %    o image: the image.
1971 %
1972 %    o radius: the radius of the pixel neighborhood.
1973 %
1974 %    o exception: return any errors or warnings in this structure.
1975 %
1976 */
1977 MagickExport Image *EdgeImage(const Image *image,const double radius,
1978   ExceptionInfo *exception)
1979 {
1980   Image
1981     *edge_image;
1982
1983   double
1984     *kernel;
1985
1986   register ssize_t
1987     i;
1988
1989   size_t
1990     width;
1991
1992   assert(image != (const Image *) NULL);
1993   assert(image->signature == MagickSignature);
1994   if (image->debug != MagickFalse)
1995     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1996   assert(exception != (ExceptionInfo *) NULL);
1997   assert(exception->signature == MagickSignature);
1998   width=GetOptimalKernelWidth1D(radius,0.5);
1999   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2000   if (kernel == (double *) NULL)
2001     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2002   for (i=0; i < (ssize_t) (width*width); i++)
2003     kernel[i]=(-1.0);
2004   kernel[i/2]=(double) (width*width-1.0);
2005   edge_image=ConvolveImage(image,width,kernel,exception);
2006   kernel=(double *) RelinquishMagickMemory(kernel);
2007   return(edge_image);
2008 }
2009 \f
2010 /*
2011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2012 %                                                                             %
2013 %                                                                             %
2014 %                                                                             %
2015 %     E m b o s s I m a g e                                                   %
2016 %                                                                             %
2017 %                                                                             %
2018 %                                                                             %
2019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2020 %
2021 %  EmbossImage() returns a grayscale image with a three-dimensional effect.
2022 %  We convolve the image with a Gaussian operator of the given radius and
2023 %  standard deviation (sigma).  For reasonable results, radius should be
2024 %  larger than sigma.  Use a radius of 0 and Emboss() selects a suitable
2025 %  radius for you.
2026 %
2027 %  The format of the EmbossImage method is:
2028 %
2029 %      Image *EmbossImage(const Image *image,const double radius,
2030 %        const double sigma,ExceptionInfo *exception)
2031 %
2032 %  A description of each parameter follows:
2033 %
2034 %    o image: the image.
2035 %
2036 %    o radius: the radius of the pixel neighborhood.
2037 %
2038 %    o sigma: the standard deviation of the Gaussian, in pixels.
2039 %
2040 %    o exception: return any errors or warnings in this structure.
2041 %
2042 */
2043 MagickExport Image *EmbossImage(const Image *image,const double radius,
2044   const double sigma,ExceptionInfo *exception)
2045 {
2046   double
2047     *kernel;
2048
2049   Image
2050     *emboss_image;
2051
2052   register ssize_t
2053     i;
2054
2055   size_t
2056     width;
2057
2058   ssize_t
2059     j,
2060     k,
2061     u,
2062     v;
2063
2064   assert(image != (Image *) NULL);
2065   assert(image->signature == MagickSignature);
2066   if (image->debug != MagickFalse)
2067     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2068   assert(exception != (ExceptionInfo *) NULL);
2069   assert(exception->signature == MagickSignature);
2070   width=GetOptimalKernelWidth2D(radius,sigma);
2071   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2072   if (kernel == (double *) NULL)
2073     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2074   j=(ssize_t) width/2;
2075   k=j;
2076   i=0;
2077   for (v=(-j); v <= j; v++)
2078   {
2079     for (u=(-j); u <= j; u++)
2080     {
2081       kernel[i]=(double) (((u < 0) || (v < 0) ? -8.0 : 8.0)*
2082         exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
2083         (2.0*MagickPI*MagickSigma*MagickSigma));
2084       if (u != k)
2085         kernel[i]=0.0;
2086       i++;
2087     }
2088     k--;
2089   }
2090   emboss_image=ConvolveImage(image,width,kernel,exception);
2091   if (emboss_image != (Image *) NULL)
2092     (void) EqualizeImage(emboss_image);
2093   kernel=(double *) RelinquishMagickMemory(kernel);
2094   return(emboss_image);
2095 }
2096 \f
2097 /*
2098 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2099 %                                                                             %
2100 %                                                                             %
2101 %                                                                             %
2102 %     F i l t e r I m a g e                                                   %
2103 %                                                                             %
2104 %                                                                             %
2105 %                                                                             %
2106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2107 %
2108 %  FilterImage() applies a custom convolution kernel to the image.
2109 %
2110 %  The format of the FilterImage method is:
2111 %
2112 %      Image *FilterImage(const Image *image,const KernelInfo *kernel,
2113 %        ExceptionInfo *exception)
2114 %      Image *FilterImageChannel(const Image *image,const ChannelType channel,
2115 %        const KernelInfo *kernel,ExceptionInfo *exception)
2116 %
2117 %  A description of each parameter follows:
2118 %
2119 %    o image: the image.
2120 %
2121 %    o channel: the channel type.
2122 %
2123 %    o kernel: the filtering kernel.
2124 %
2125 %    o exception: return any errors or warnings in this structure.
2126 %
2127 */
2128
2129 MagickExport Image *FilterImage(const Image *image,const KernelInfo *kernel,
2130   ExceptionInfo *exception)
2131 {
2132   Image
2133     *filter_image;
2134
2135   filter_image=FilterImageChannel(image,DefaultChannels,kernel,exception);
2136   return(filter_image);
2137 }
2138
2139 MagickExport Image *FilterImageChannel(const Image *image,
2140   const ChannelType channel,const KernelInfo *kernel,ExceptionInfo *exception)
2141 {
2142 #define FilterImageTag  "Filter/Image"
2143
2144   CacheView
2145     *filter_view,
2146     *image_view;
2147
2148   Image
2149     *filter_image;
2150
2151   MagickBooleanType
2152     status;
2153
2154   MagickOffsetType
2155     progress;
2156
2157   MagickPixelPacket
2158     bias;
2159
2160   ssize_t
2161     y;
2162
2163   /*
2164     Initialize filter image attributes.
2165   */
2166   assert(image != (Image *) NULL);
2167   assert(image->signature == MagickSignature);
2168   if (image->debug != MagickFalse)
2169     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2170   assert(exception != (ExceptionInfo *) NULL);
2171   assert(exception->signature == MagickSignature);
2172   if ((kernel->width % 2) == 0)
2173     ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
2174   filter_image=CloneImage(image,0,0,MagickTrue,exception);
2175   if (filter_image == (Image *) NULL)
2176     return((Image *) NULL);
2177   if (SetImageStorageClass(filter_image,DirectClass) == MagickFalse)
2178     {
2179       InheritException(exception,&filter_image->exception);
2180       filter_image=DestroyImage(filter_image);
2181       return((Image *) NULL);
2182     }
2183   if (image->debug != MagickFalse)
2184     {
2185       char
2186         format[MaxTextExtent],
2187         *message;
2188
2189       register const double
2190         *k;
2191
2192       ssize_t
2193         u,
2194         v;
2195
2196       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
2197         "  FilterImage with %.20gx%.20g kernel:",(double) kernel->width,(double)
2198         kernel->height);
2199       message=AcquireString("");
2200       k=kernel->values;
2201       for (v=0; v < (ssize_t) kernel->height; v++)
2202       {
2203         *message='\0';
2204         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) v);
2205         (void) ConcatenateString(&message,format);
2206         for (u=0; u < (ssize_t) kernel->width; u++)
2207         {
2208           (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
2209           (void) ConcatenateString(&message,format);
2210         }
2211         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
2212       }
2213       message=DestroyString(message);
2214     }
2215   status=AccelerateConvolveImage(image,kernel,filter_image,exception);
2216   if (status == MagickTrue)
2217     return(filter_image);
2218   /*
2219     Filter image.
2220   */
2221   status=MagickTrue;
2222   progress=0;
2223   GetMagickPixelPacket(image,&bias);
2224   SetMagickPixelPacketBias(image,&bias);
2225   image_view=AcquireCacheView(image);
2226   filter_view=AcquireCacheView(filter_image);
2227 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2228   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2229 #endif
2230   for (y=0; y < (ssize_t) image->rows; y++)
2231   {
2232     MagickBooleanType
2233       sync;
2234
2235     register const IndexPacket
2236       *restrict indexes;
2237
2238     register const PixelPacket
2239       *restrict p;
2240
2241     register IndexPacket
2242       *restrict filter_indexes;
2243
2244     register PixelPacket
2245       *restrict q;
2246
2247     register ssize_t
2248       x;
2249
2250     if (status == MagickFalse)
2251       continue;
2252     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel->width/2L),
2253       y-(ssize_t) (kernel->height/2L),image->columns+kernel->width,
2254       kernel->height,exception);
2255     q=GetCacheViewAuthenticPixels(filter_view,0,y,filter_image->columns,1,
2256       exception);
2257     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2258       {
2259         status=MagickFalse;
2260         continue;
2261       }
2262     indexes=GetCacheViewVirtualIndexQueue(image_view);
2263     filter_indexes=GetCacheViewAuthenticIndexQueue(filter_view);
2264     for (x=0; x < (ssize_t) image->columns; x++)
2265     {
2266       MagickPixelPacket
2267         pixel;
2268
2269       register const double
2270         *restrict k;
2271
2272       register const PixelPacket
2273         *restrict kernel_pixels;
2274
2275       register ssize_t
2276         u;
2277
2278       ssize_t
2279         v;
2280
2281       pixel=bias;
2282       k=kernel->values;
2283       kernel_pixels=p;
2284       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
2285         {
2286           for (v=0; v < (ssize_t) kernel->width; v++)
2287           {
2288             for (u=0; u < (ssize_t) kernel->height; u++)
2289             {
2290               pixel.red+=(*k)*kernel_pixels[u].red;
2291               pixel.green+=(*k)*kernel_pixels[u].green;
2292               pixel.blue+=(*k)*kernel_pixels[u].blue;
2293               k++;
2294             }
2295             kernel_pixels+=image->columns+kernel->width;
2296           }
2297           if ((channel & RedChannel) != 0)
2298             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
2299           if ((channel & GreenChannel) != 0)
2300             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
2301           if ((channel & BlueChannel) != 0)
2302             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
2303           if ((channel & OpacityChannel) != 0)
2304             {
2305               k=kernel->values;
2306               kernel_pixels=p;
2307               for (v=0; v < (ssize_t) kernel->width; v++)
2308               {
2309                 for (u=0; u < (ssize_t) kernel->height; u++)
2310                 {
2311                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
2312                   k++;
2313                 }
2314                 kernel_pixels+=image->columns+kernel->width;
2315               }
2316               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2317             }
2318           if (((channel & IndexChannel) != 0) &&
2319               (image->colorspace == CMYKColorspace))
2320             {
2321               register const IndexPacket
2322                 *restrict kernel_indexes;
2323
2324               k=kernel->values;
2325               kernel_indexes=indexes;
2326               for (v=0; v < (ssize_t) kernel->width; v++)
2327               {
2328                 for (u=0; u < (ssize_t) kernel->height; u++)
2329                 {
2330                   pixel.index+=(*k)*kernel_indexes[u];
2331                   k++;
2332                 }
2333                 kernel_indexes+=image->columns+kernel->width;
2334               }
2335               filter_indexes[x]=ClampToQuantum(pixel.index);
2336             }
2337         }
2338       else
2339         {
2340           MagickRealType
2341             alpha,
2342             gamma;
2343
2344           gamma=0.0;
2345           for (v=0; v < (ssize_t) kernel->width; v++)
2346           {
2347             for (u=0; u < (ssize_t) kernel->height; u++)
2348             {
2349               alpha=(MagickRealType) (QuantumScale*(QuantumRange-
2350                 kernel_pixels[u].opacity));
2351               pixel.red+=(*k)*alpha*kernel_pixels[u].red;
2352               pixel.green+=(*k)*alpha*kernel_pixels[u].green;
2353               pixel.blue+=(*k)*alpha*kernel_pixels[u].blue;
2354               gamma+=(*k)*alpha;
2355               k++;
2356             }
2357             kernel_pixels+=image->columns+kernel->width;
2358           }
2359           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2360           if ((channel & RedChannel) != 0)
2361             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
2362           if ((channel & GreenChannel) != 0)
2363             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
2364           if ((channel & BlueChannel) != 0)
2365             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
2366           if ((channel & OpacityChannel) != 0)
2367             {
2368               k=kernel->values;
2369               kernel_pixels=p;
2370               for (v=0; v < (ssize_t) kernel->width; v++)
2371               {
2372                 for (u=0; u < (ssize_t) kernel->height; u++)
2373                 {
2374                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
2375                   k++;
2376                 }
2377                 kernel_pixels+=image->columns+kernel->width;
2378               }
2379               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
2380             }
2381           if (((channel & IndexChannel) != 0) &&
2382               (image->colorspace == CMYKColorspace))
2383             {
2384               register const IndexPacket
2385                 *restrict kernel_indexes;
2386
2387               k=kernel->values;
2388               kernel_pixels=p;
2389               kernel_indexes=indexes;
2390               for (v=0; v < (ssize_t) kernel->width; v++)
2391               {
2392                 for (u=0; u < (ssize_t) kernel->height; u++)
2393                 {
2394                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
2395                     kernel_pixels[u].opacity));
2396                   pixel.index+=(*k)*alpha*kernel_indexes[u];
2397                   k++;
2398                 }
2399                 kernel_pixels+=image->columns+kernel->width;
2400                 kernel_indexes+=image->columns+kernel->width;
2401               }
2402               filter_indexes[x]=ClampToQuantum(gamma*
2403                 GetIndexPixelComponent(&pixel));
2404             }
2405         }
2406       indexes++;
2407       p++;
2408       q++;
2409     }
2410     sync=SyncCacheViewAuthenticPixels(filter_view,exception);
2411     if (sync == MagickFalse)
2412       status=MagickFalse;
2413     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2414       {
2415         MagickBooleanType
2416           proceed;
2417
2418 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2419   #pragma omp critical (MagickCore_FilterImageChannel)
2420 #endif
2421         proceed=SetImageProgress(image,FilterImageTag,progress++,image->rows);
2422         if (proceed == MagickFalse)
2423           status=MagickFalse;
2424       }
2425   }
2426   filter_image->type=image->type;
2427   filter_view=DestroyCacheView(filter_view);
2428   image_view=DestroyCacheView(image_view);
2429   if (status == MagickFalse)
2430     filter_image=DestroyImage(filter_image);
2431   return(filter_image);
2432 }
2433 \f
2434 /*
2435 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2436 %                                                                             %
2437 %                                                                             %
2438 %                                                                             %
2439 %     G a u s s i a n B l u r I m a g e                                       %
2440 %                                                                             %
2441 %                                                                             %
2442 %                                                                             %
2443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2444 %
2445 %  GaussianBlurImage() blurs an image.  We convolve the image with a
2446 %  Gaussian operator of the given radius and standard deviation (sigma).
2447 %  For reasonable results, the radius should be larger than sigma.  Use a
2448 %  radius of 0 and GaussianBlurImage() selects a suitable radius for you
2449 %
2450 %  The format of the GaussianBlurImage method is:
2451 %
2452 %      Image *GaussianBlurImage(const Image *image,onst double radius,
2453 %        const double sigma,ExceptionInfo *exception)
2454 %      Image *GaussianBlurImageChannel(const Image *image,
2455 %        const ChannelType channel,const double radius,const double sigma,
2456 %        ExceptionInfo *exception)
2457 %
2458 %  A description of each parameter follows:
2459 %
2460 %    o image: the image.
2461 %
2462 %    o channel: the channel type.
2463 %
2464 %    o radius: the radius of the Gaussian, in pixels, not counting the center
2465 %      pixel.
2466 %
2467 %    o sigma: the standard deviation of the Gaussian, in pixels.
2468 %
2469 %    o exception: return any errors or warnings in this structure.
2470 %
2471 */
2472
2473 MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
2474   const double sigma,ExceptionInfo *exception)
2475 {
2476   Image
2477     *blur_image;
2478
2479   blur_image=GaussianBlurImageChannel(image,DefaultChannels,radius,sigma,
2480     exception);
2481   return(blur_image);
2482 }
2483
2484 MagickExport Image *GaussianBlurImageChannel(const Image *image,
2485   const ChannelType channel,const double radius,const double sigma,
2486   ExceptionInfo *exception)
2487 {
2488   double
2489     *kernel;
2490
2491   Image
2492     *blur_image;
2493
2494   register ssize_t
2495     i;
2496
2497   size_t
2498     width;
2499
2500   ssize_t
2501     j,
2502     u,
2503     v;
2504
2505   assert(image != (const Image *) NULL);
2506   assert(image->signature == MagickSignature);
2507   if (image->debug != MagickFalse)
2508     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2509   assert(exception != (ExceptionInfo *) NULL);
2510   assert(exception->signature == MagickSignature);
2511   width=GetOptimalKernelWidth2D(radius,sigma);
2512   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
2513   if (kernel == (double *) NULL)
2514     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2515   j=(ssize_t) width/2;
2516   i=0;
2517   for (v=(-j); v <= j; v++)
2518   {
2519     for (u=(-j); u <= j; u++)
2520       kernel[i++]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
2521         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
2522   }
2523   blur_image=ConvolveImageChannel(image,channel,width,kernel,exception);
2524   kernel=(double *) RelinquishMagickMemory(kernel);
2525   return(blur_image);
2526 }
2527 \f
2528 /*
2529 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2530 %                                                                             %
2531 %                                                                             %
2532 %                                                                             %
2533 %     M e d i a n F i l t e r I m a g e                                       %
2534 %                                                                             %
2535 %                                                                             %
2536 %                                                                             %
2537 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2538 %
2539 %  MedianFilterImage() applies a digital filter that improves the quality
2540 %  of a noisy image.  Each pixel is replaced by the median in a set of
2541 %  neighboring pixels as defined by radius.
2542 %
2543 %  The algorithm was contributed by Mike Edmonds and implements an insertion
2544 %  sort for selecting median color-channel values.  For more on this algorithm
2545 %  see "Skip Lists: A probabilistic Alternative to Balanced Trees" by William
2546 %  Pugh in the June 1990 of Communications of the ACM.
2547 %
2548 %  The format of the MedianFilterImage method is:
2549 %
2550 %      Image *MedianFilterImage(const Image *image,const double radius,
2551 %        ExceptionInfo *exception)
2552 %
2553 %  A description of each parameter follows:
2554 %
2555 %    o image: the image.
2556 %
2557 %    o radius: the radius of the pixel neighborhood.
2558 %
2559 %    o exception: return any errors or warnings in this structure.
2560 %
2561 */
2562
2563 #define ListChannels  5
2564
2565 typedef struct _ListNode
2566 {
2567   size_t
2568     next[9],
2569     count,
2570     signature;
2571 } ListNode;
2572
2573 typedef struct _SkipList
2574 {
2575   ssize_t
2576     level;
2577
2578   ListNode
2579     *nodes;
2580 } SkipList;
2581
2582 typedef struct _PixelList
2583 {
2584   size_t
2585     center,
2586     seed,
2587     signature;
2588
2589   SkipList
2590     lists[ListChannels];
2591 } PixelList;
2592
2593 static PixelList *DestroyPixelList(PixelList *pixel_list)
2594 {
2595   register ssize_t
2596     i;
2597
2598   if (pixel_list == (PixelList *) NULL)
2599     return((PixelList *) NULL);
2600   for (i=0; i < ListChannels; i++)
2601     if (pixel_list->lists[i].nodes != (ListNode *) NULL)
2602       pixel_list->lists[i].nodes=(ListNode *) RelinquishMagickMemory(
2603         pixel_list->lists[i].nodes);
2604   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
2605   return(pixel_list);
2606 }
2607
2608 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
2609 {
2610   register ssize_t
2611     i;
2612
2613   assert(pixel_list != (PixelList **) NULL);
2614   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2615     if (pixel_list[i] != (PixelList *) NULL)
2616       pixel_list[i]=DestroyPixelList(pixel_list[i]);
2617   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
2618   return(pixel_list);
2619 }
2620
2621 static PixelList *AcquirePixelList(const size_t width)
2622 {
2623   PixelList
2624     *pixel_list;
2625
2626   register ssize_t
2627     i;
2628
2629   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2630   if (pixel_list == (PixelList *) NULL)
2631     return(pixel_list);
2632   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
2633   pixel_list->center=width*width/2;
2634   for (i=0; i < ListChannels; i++)
2635   {
2636     pixel_list->lists[i].nodes=(ListNode *) AcquireQuantumMemory(65537UL,
2637       sizeof(*pixel_list->lists[i].nodes));
2638     if (pixel_list->lists[i].nodes == (ListNode *) NULL)
2639       return(DestroyPixelList(pixel_list));
2640     (void) ResetMagickMemory(pixel_list->lists[i].nodes,0,65537UL*
2641       sizeof(*pixel_list->lists[i].nodes));
2642   }
2643   pixel_list->signature=MagickSignature;
2644   return(pixel_list);
2645 }
2646
2647 static PixelList **AcquirePixelListThreadSet(const size_t width)
2648 {
2649   PixelList
2650     **pixel_list;
2651
2652   register ssize_t
2653     i;
2654
2655   size_t
2656     number_threads;
2657
2658   number_threads=GetOpenMPMaximumThreads();
2659   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2660     sizeof(*pixel_list));
2661   if (pixel_list == (PixelList **) NULL)
2662     return((PixelList **) NULL);
2663   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
2664   for (i=0; i < (ssize_t) number_threads; i++)
2665   {
2666     pixel_list[i]=AcquirePixelList(width);
2667     if (pixel_list[i] == (PixelList *) NULL)
2668       return(DestroyPixelListThreadSet(pixel_list));
2669   }
2670   return(pixel_list);
2671 }
2672
2673 static void AddNodePixelList(PixelList *pixel_list,const ssize_t channel,
2674   const size_t color)
2675 {
2676   register SkipList
2677     *list;
2678
2679   register ssize_t
2680     level;
2681
2682   size_t
2683     search,
2684     update[9];
2685
2686   /*
2687     Initialize the node.
2688   */
2689   list=pixel_list->lists+channel;
2690   list->nodes[color].signature=pixel_list->signature;
2691   list->nodes[color].count=1;
2692   /*
2693     Determine where it belongs in the list.
2694   */
2695   search=65536UL;
2696   for (level=list->level; level >= 0; level--)
2697   {
2698     while (list->nodes[search].next[level] < color)
2699       search=list->nodes[search].next[level];
2700     update[level]=search;
2701   }
2702   /*
2703     Generate a pseudo-random level for this node.
2704   */
2705   for (level=0; ; level++)
2706   {
2707     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2708     if ((pixel_list->seed & 0x300) != 0x300)
2709       break;
2710   }
2711   if (level > 8)
2712     level=8;
2713   if (level > (list->level+2))
2714     level=list->level+2;
2715   /*
2716     If we're raising the list's level, link back to the root node.
2717   */
2718   while (level > list->level)
2719   {
2720     list->level++;
2721     update[list->level]=65536UL;
2722   }
2723   /*
2724     Link the node into the skip-list.
2725   */
2726   do
2727   {
2728     list->nodes[color].next[level]=list->nodes[update[level]].next[level];
2729     list->nodes[update[level]].next[level]=color;
2730   }
2731   while (level-- > 0);
2732 }
2733
2734 static MagickPixelPacket GetPixelList(PixelList *pixel_list)
2735 {
2736   MagickPixelPacket
2737     pixel;
2738
2739   register SkipList
2740     *list;
2741
2742   register ssize_t
2743     channel;
2744
2745   size_t
2746     center,
2747     color,
2748     count;
2749
2750   unsigned short
2751     channels[ListChannels];
2752
2753   /*
2754     Find the median value for each of the color.
2755   */
2756   center=pixel_list->center;
2757   for (channel=0; channel < 5; channel++)
2758   {
2759     list=pixel_list->lists+channel;
2760     color=65536UL;
2761     count=0;
2762     do
2763     {
2764       color=list->nodes[color].next[0];
2765       count+=list->nodes[color].count;
2766     }
2767     while (count <= center);
2768     channels[channel]=(unsigned short) color;
2769   }
2770   GetMagickPixelPacket((const Image *) NULL,&pixel);
2771   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
2772   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
2773   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
2774   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
2775   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
2776   return(pixel);
2777 }
2778
2779 static inline void InsertPixelList(const Image *image,const PixelPacket *pixel,
2780   const IndexPacket *indexes,PixelList *pixel_list)
2781 {
2782   size_t
2783     signature;
2784
2785   unsigned short
2786     index;
2787
2788   index=ScaleQuantumToShort(pixel->red);
2789   signature=pixel_list->lists[0].nodes[index].signature;
2790   if (signature == pixel_list->signature)
2791     pixel_list->lists[0].nodes[index].count++;
2792   else
2793     AddNodePixelList(pixel_list,0,index);
2794   index=ScaleQuantumToShort(pixel->green);
2795   signature=pixel_list->lists[1].nodes[index].signature;
2796   if (signature == pixel_list->signature)
2797     pixel_list->lists[1].nodes[index].count++;
2798   else
2799     AddNodePixelList(pixel_list,1,index);
2800   index=ScaleQuantumToShort(pixel->blue);
2801   signature=pixel_list->lists[2].nodes[index].signature;
2802   if (signature == pixel_list->signature)
2803     pixel_list->lists[2].nodes[index].count++;
2804   else
2805     AddNodePixelList(pixel_list,2,index);
2806   index=ScaleQuantumToShort(pixel->opacity);
2807   signature=pixel_list->lists[3].nodes[index].signature;
2808   if (signature == pixel_list->signature)
2809     pixel_list->lists[3].nodes[index].count++;
2810   else
2811     AddNodePixelList(pixel_list,3,index);
2812   if (image->colorspace == CMYKColorspace)
2813     index=ScaleQuantumToShort(*indexes);
2814   signature=pixel_list->lists[4].nodes[index].signature;
2815   if (signature == pixel_list->signature)
2816     pixel_list->lists[4].nodes[index].count++;
2817   else
2818     AddNodePixelList(pixel_list,4,index);
2819 }
2820
2821 static void ResetPixelList(PixelList *pixel_list)
2822 {
2823   int
2824     level;
2825
2826   register ListNode
2827     *root;
2828
2829   register SkipList
2830     *list;
2831
2832   register ssize_t
2833     channel;
2834
2835   /*
2836     Reset the skip-list.
2837   */
2838   for (channel=0; channel < 5; channel++)
2839   {
2840     list=pixel_list->lists+channel;
2841     root=list->nodes+65536UL;
2842     list->level=0;
2843     for (level=0; level < 9; level++)
2844       root->next[level]=65536UL;
2845   }
2846   pixel_list->seed=pixel_list->signature++;
2847 }
2848
2849 MagickExport Image *MedianFilterImage(const Image *image,const double radius,
2850   ExceptionInfo *exception)
2851 {
2852 #define MedianFilterImageTag  "MedianFilter/Image"
2853
2854   CacheView
2855     *image_view,
2856     *median_view;
2857
2858   Image
2859     *median_image;
2860
2861   MagickBooleanType
2862     status;
2863
2864   MagickOffsetType
2865     progress;
2866
2867   PixelList
2868     **restrict pixel_list;
2869
2870   size_t
2871     width;
2872
2873   ssize_t
2874     y;
2875
2876   /*
2877     Initialize median image attributes.
2878   */
2879   assert(image != (Image *) NULL);
2880   assert(image->signature == MagickSignature);
2881   if (image->debug != MagickFalse)
2882     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2883   assert(exception != (ExceptionInfo *) NULL);
2884   assert(exception->signature == MagickSignature);
2885   width=GetOptimalKernelWidth2D(radius,0.5);
2886   median_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2887     exception);
2888   if (median_image == (Image *) NULL)
2889     return((Image *) NULL);
2890   if (SetImageStorageClass(median_image,DirectClass) == MagickFalse)
2891     {
2892       InheritException(exception,&median_image->exception);
2893       median_image=DestroyImage(median_image);
2894       return((Image *) NULL);
2895     }
2896   pixel_list=AcquirePixelListThreadSet(width);
2897   if (pixel_list == (PixelList **) NULL)
2898     {
2899       median_image=DestroyImage(median_image);
2900       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2901     }
2902   /*
2903     Median filter each image row.
2904   */
2905   status=MagickTrue;
2906   progress=0;
2907   image_view=AcquireCacheView(image);
2908   median_view=AcquireCacheView(median_image);
2909 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2910   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2911 #endif
2912   for (y=0; y < (ssize_t) median_image->rows; y++)
2913   {
2914     const int
2915       id = GetOpenMPThreadId();
2916
2917     register const IndexPacket
2918       *restrict indexes;
2919
2920     register const PixelPacket
2921       *restrict p;
2922
2923     register IndexPacket
2924       *restrict median_indexes;
2925
2926     register PixelPacket
2927       *restrict q;
2928
2929     register ssize_t
2930       x;
2931
2932     if (status == MagickFalse)
2933       continue;
2934     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
2935       (width/2L),image->columns+width,width,exception);
2936     q=QueueCacheViewAuthenticPixels(median_view,0,y,median_image->columns,1,
2937       exception);
2938     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2939       {
2940         status=MagickFalse;
2941         continue;
2942       }
2943     indexes=GetCacheViewVirtualIndexQueue(image_view);
2944     median_indexes=GetCacheViewAuthenticIndexQueue(median_view);
2945     for (x=0; x < (ssize_t) median_image->columns; x++)
2946     {
2947       MagickPixelPacket
2948         pixel;
2949
2950       register const IndexPacket
2951         *restrict s;
2952
2953       register const PixelPacket
2954         *restrict r;
2955
2956       register ssize_t
2957         u,
2958         v;
2959
2960       r=p;
2961       s=indexes+x;
2962       ResetPixelList(pixel_list[id]);
2963       for (v=0; v < (ssize_t) width; v++)
2964       {
2965         for (u=0; u < (ssize_t) width; u++)
2966           InsertPixelList(image,r+u,s+u,pixel_list[id]);
2967         r+=image->columns+width;
2968         s+=image->columns+width;
2969       }
2970       pixel=GetPixelList(pixel_list[id]);
2971       SetPixelPacket(median_image,&pixel,q,median_indexes+x);
2972       p++;
2973       q++;
2974     }
2975     if (SyncCacheViewAuthenticPixels(median_view,exception) == MagickFalse)
2976       status=MagickFalse;
2977     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2978       {
2979         MagickBooleanType
2980           proceed;
2981
2982 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2983   #pragma omp critical (MagickCore_MedianFilterImage)
2984 #endif
2985         proceed=SetImageProgress(image,MedianFilterImageTag,progress++,
2986           image->rows);
2987         if (proceed == MagickFalse)
2988           status=MagickFalse;
2989       }
2990   }
2991   median_view=DestroyCacheView(median_view);
2992   image_view=DestroyCacheView(image_view);
2993   pixel_list=DestroyPixelListThreadSet(pixel_list);
2994   return(median_image);
2995 }
2996 \f
2997 /*
2998 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2999 %                                                                             %
3000 %                                                                             %
3001 %                                                                             %
3002 %     M o d e I m a g e                                                       %
3003 %                                                                             %
3004 %                                                                             %
3005 %                                                                             %
3006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3007 %
3008 %  ModeImage() makes each pixel the 'predominate color' of the neighborhood
3009 %  of the specified radius.
3010 %
3011 %  The format of the ModeImage method is:
3012 %
3013 %      Image *ModeImage(const Image *image,const double radius,
3014 %        ExceptionInfo *exception)
3015 %
3016 %  A description of each parameter follows:
3017 %
3018 %    o image: the image.
3019 %
3020 %    o radius: the radius of the pixel neighborhood.
3021 %
3022 %    o exception: return any errors or warnings in this structure.
3023 %
3024 */
3025
3026 static MagickPixelPacket GetModePixelList(PixelList *pixel_list)
3027 {
3028   MagickPixelPacket
3029     pixel;
3030
3031   register SkipList
3032     *list;
3033
3034   register ssize_t
3035     channel;
3036
3037   size_t
3038     color,
3039     count,
3040     max_count,
3041     mode,
3042     width;
3043
3044   unsigned short
3045     channels[5];
3046
3047   /*
3048     Make each pixel the 'predominate color' of the specified neighborhood.
3049   */
3050   width=pixel_list->center << 1;
3051   for (channel=0; channel < 5; channel++)
3052   {
3053     list=pixel_list->lists+channel;
3054     color=65536UL;
3055     mode=color;
3056     max_count=list->nodes[mode].count;
3057     count=0;
3058     do
3059     {
3060       color=list->nodes[color].next[0];
3061       if (list->nodes[color].count > max_count)
3062         {
3063           mode=color;
3064           max_count=list->nodes[mode].count;
3065         }
3066       count+=list->nodes[color].count;
3067     }
3068     while (count <= width);
3069     channels[channel]=(unsigned short) mode;
3070   }
3071   GetMagickPixelPacket((const Image *) NULL,&pixel);
3072   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3073   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3074   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3075   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3076   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3077   return(pixel);
3078 }
3079
3080 MagickExport Image *ModeImage(const Image *image,const double radius,
3081   ExceptionInfo *exception)
3082 {
3083 #define ModeImageTag  "Mode/Image"
3084
3085   CacheView
3086     *image_view,
3087     *mode_view;
3088
3089   Image
3090     *mode_image;
3091
3092   MagickBooleanType
3093     status;
3094
3095   MagickOffsetType
3096     progress;
3097
3098   PixelList
3099     **restrict pixel_list;
3100
3101   size_t
3102     width;
3103
3104   ssize_t
3105     y;
3106
3107   /*
3108     Initialize mode image attributes.
3109   */
3110   assert(image != (Image *) NULL);
3111   assert(image->signature == MagickSignature);
3112   if (image->debug != MagickFalse)
3113     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3114   assert(exception != (ExceptionInfo *) NULL);
3115   assert(exception->signature == MagickSignature);
3116   width=GetOptimalKernelWidth2D(radius,0.5);
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   noise_image=CloneImage(image,image->columns,image->rows,MagickTrue,
4470     exception);
4471   if (noise_image == (Image *) NULL)
4472     return((Image *) NULL);
4473   if (SetImageStorageClass(noise_image,DirectClass) == MagickFalse)
4474     {
4475       InheritException(exception,&noise_image->exception);
4476       noise_image=DestroyImage(noise_image);
4477       return((Image *) NULL);
4478     }
4479   pixel_list=AcquirePixelListThreadSet(width);
4480   if (pixel_list == (PixelList **) NULL)
4481     {
4482       noise_image=DestroyImage(noise_image);
4483       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4484     }
4485   /*
4486     Reduce noise image.
4487   */
4488   status=MagickTrue;
4489   progress=0;
4490   image_view=AcquireCacheView(image);
4491   noise_view=AcquireCacheView(noise_image);
4492 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4493   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4494 #endif
4495   for (y=0; y < (ssize_t) noise_image->rows; y++)
4496   {
4497     const int
4498       id = GetOpenMPThreadId();
4499
4500     register const IndexPacket
4501       *restrict indexes;
4502
4503     register const PixelPacket
4504       *restrict p;
4505
4506     register IndexPacket
4507       *restrict noise_indexes;
4508
4509     register PixelPacket
4510       *restrict q;
4511
4512     register ssize_t
4513       x;
4514
4515     if (status == MagickFalse)
4516       continue;
4517     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
4518       (width/2L),image->columns+width,width,exception);
4519     q=QueueCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
4520       exception);
4521     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4522       {
4523         status=MagickFalse;
4524         continue;
4525       }
4526     indexes=GetCacheViewVirtualIndexQueue(image_view);
4527     noise_indexes=GetCacheViewAuthenticIndexQueue(noise_view);
4528     for (x=0; x < (ssize_t) noise_image->columns; x++)
4529     {
4530       MagickPixelPacket
4531         pixel;
4532
4533       register const PixelPacket
4534         *restrict r;
4535
4536       register const IndexPacket
4537         *restrict s;
4538
4539       register ssize_t
4540         u,
4541         v;
4542
4543       r=p;
4544       s=indexes+x;
4545       ResetPixelList(pixel_list[id]);
4546       for (v=0; v < (ssize_t) width; v++)
4547       {
4548         for (u=0; u < (ssize_t) width; u++)
4549           InsertPixelList(image,r+u,s+u,pixel_list[id]);
4550         r+=image->columns+width;
4551         s+=image->columns+width;
4552       }
4553       pixel=GetNonpeakPixelList(pixel_list[id]);
4554       SetPixelPacket(noise_image,&pixel,q,noise_indexes+x);
4555       p++;
4556       q++;
4557     }
4558     if (SyncCacheViewAuthenticPixels(noise_view,exception) == MagickFalse)
4559       status=MagickFalse;
4560     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4561       {
4562         MagickBooleanType
4563           proceed;
4564
4565 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4566   #pragma omp critical (MagickCore_ReduceNoiseImage)
4567 #endif
4568         proceed=SetImageProgress(image,ReduceNoiseImageTag,progress++,
4569           image->rows);
4570         if (proceed == MagickFalse)
4571           status=MagickFalse;
4572       }
4573   }
4574   noise_view=DestroyCacheView(noise_view);
4575   image_view=DestroyCacheView(image_view);
4576   pixel_list=DestroyPixelListThreadSet(pixel_list);
4577   return(noise_image);
4578 }
4579 \f
4580 /*
4581 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4582 %                                                                             %
4583 %                                                                             %
4584 %                                                                             %
4585 %     S e l e c t i v e B l u r I m a g e                                     %
4586 %                                                                             %
4587 %                                                                             %
4588 %                                                                             %
4589 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4590 %
4591 %  SelectiveBlurImage() selectively blur pixels within a contrast threshold.
4592 %  It is similar to the unsharpen mask that sharpens everything with contrast
4593 %  above a certain threshold.
4594 %
4595 %  The format of the SelectiveBlurImage method is:
4596 %
4597 %      Image *SelectiveBlurImage(const Image *image,const double radius,
4598 %        const double sigma,const double threshold,ExceptionInfo *exception)
4599 %      Image *SelectiveBlurImageChannel(const Image *image,
4600 %        const ChannelType channel,const double radius,const double sigma,
4601 %        const double threshold,ExceptionInfo *exception)
4602 %
4603 %  A description of each parameter follows:
4604 %
4605 %    o image: the image.
4606 %
4607 %    o channel: the channel type.
4608 %
4609 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4610 %      pixel.
4611 %
4612 %    o sigma: the standard deviation of the Gaussian, in pixels.
4613 %
4614 %    o threshold: only pixels within this contrast threshold are included
4615 %      in the blur operation.
4616 %
4617 %    o exception: return any errors or warnings in this structure.
4618 %
4619 */
4620
4621 static inline MagickBooleanType SelectiveContrast(const PixelPacket *p,
4622   const PixelPacket *q,const double threshold)
4623 {
4624   if (fabs(PixelIntensity(p)-PixelIntensity(q)) < threshold)
4625     return(MagickTrue);
4626   return(MagickFalse);
4627 }
4628
4629 MagickExport Image *SelectiveBlurImage(const Image *image,const double radius,
4630   const double sigma,const double threshold,ExceptionInfo *exception)
4631 {
4632   Image
4633     *blur_image;
4634
4635   blur_image=SelectiveBlurImageChannel(image,DefaultChannels,radius,sigma,
4636     threshold,exception);
4637   return(blur_image);
4638 }
4639
4640 MagickExport Image *SelectiveBlurImageChannel(const Image *image,
4641   const ChannelType channel,const double radius,const double sigma,
4642   const double threshold,ExceptionInfo *exception)
4643 {
4644 #define SelectiveBlurImageTag  "SelectiveBlur/Image"
4645
4646   CacheView
4647     *blur_view,
4648     *image_view;
4649
4650   double
4651     *kernel;
4652
4653   Image
4654     *blur_image;
4655
4656   MagickBooleanType
4657     status;
4658
4659   MagickOffsetType
4660     progress;
4661
4662   MagickPixelPacket
4663     bias;
4664
4665   register ssize_t
4666     i;
4667
4668   size_t
4669     width;
4670
4671   ssize_t
4672     j,
4673     u,
4674     v,
4675     y;
4676
4677   /*
4678     Initialize blur image attributes.
4679   */
4680   assert(image != (Image *) NULL);
4681   assert(image->signature == MagickSignature);
4682   if (image->debug != MagickFalse)
4683     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4684   assert(exception != (ExceptionInfo *) NULL);
4685   assert(exception->signature == MagickSignature);
4686   width=GetOptimalKernelWidth1D(radius,sigma);
4687   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
4688   if (kernel == (double *) NULL)
4689     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4690   j=(ssize_t) width/2;
4691   i=0;
4692   for (v=(-j); v <= j; v++)
4693   {
4694     for (u=(-j); u <= j; u++)
4695       kernel[i++]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
4696         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
4697   }
4698   if (image->debug != MagickFalse)
4699     {
4700       char
4701         format[MaxTextExtent],
4702         *message;
4703
4704       register const double
4705         *k;
4706
4707       ssize_t
4708         u,
4709         v;
4710
4711       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
4712         "  SelectiveBlurImage with %.20gx%.20g kernel:",(double) width,(double)
4713         width);
4714       message=AcquireString("");
4715       k=kernel;
4716       for (v=0; v < (ssize_t) width; v++)
4717       {
4718         *message='\0';
4719         (void) FormatMagickString(format,MaxTextExtent,"%.20g: ",(double) v);
4720         (void) ConcatenateString(&message,format);
4721         for (u=0; u < (ssize_t) width; u++)
4722         {
4723           (void) FormatMagickString(format,MaxTextExtent,"%+f ",*k++);
4724           (void) ConcatenateString(&message,format);
4725         }
4726         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
4727       }
4728       message=DestroyString(message);
4729     }
4730   blur_image=CloneImage(image,0,0,MagickTrue,exception);
4731   if (blur_image == (Image *) NULL)
4732     return((Image *) NULL);
4733   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
4734     {
4735       InheritException(exception,&blur_image->exception);
4736       blur_image=DestroyImage(blur_image);
4737       return((Image *) NULL);
4738     }
4739   /*
4740     Threshold blur image.
4741   */
4742   status=MagickTrue;
4743   progress=0;
4744   GetMagickPixelPacket(image,&bias);
4745   SetMagickPixelPacketBias(image,&bias);
4746   image_view=AcquireCacheView(image);
4747   blur_view=AcquireCacheView(blur_image);
4748 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4749   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
4750 #endif
4751   for (y=0; y < (ssize_t) image->rows; y++)
4752   {
4753     MagickBooleanType
4754       sync;
4755
4756     MagickRealType
4757       gamma;
4758
4759     register const IndexPacket
4760       *restrict indexes;
4761
4762     register const PixelPacket
4763       *restrict p;
4764
4765     register IndexPacket
4766       *restrict blur_indexes;
4767
4768     register PixelPacket
4769       *restrict q;
4770
4771     register ssize_t
4772       x;
4773
4774     if (status == MagickFalse)
4775       continue;
4776     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
4777       (width/2L),image->columns+width,width,exception);
4778     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
4779       exception);
4780     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
4781       {
4782         status=MagickFalse;
4783         continue;
4784       }
4785     indexes=GetCacheViewVirtualIndexQueue(image_view);
4786     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
4787     for (x=0; x < (ssize_t) image->columns; x++)
4788     {
4789       MagickPixelPacket
4790         pixel;
4791
4792       register const double
4793         *restrict k;
4794
4795       register ssize_t
4796         u;
4797
4798       ssize_t
4799         j,
4800         v;
4801
4802       pixel=bias;
4803       k=kernel;
4804       gamma=0.0;
4805       j=0;
4806       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
4807         {
4808           for (v=0; v < (ssize_t) width; v++)
4809           {
4810             for (u=0; u < (ssize_t) width; u++)
4811             {
4812               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4813                 {
4814                   pixel.red+=(*k)*(p+u+j)->red;
4815                   pixel.green+=(*k)*(p+u+j)->green;
4816                   pixel.blue+=(*k)*(p+u+j)->blue;
4817                   gamma+=(*k);
4818                   k++;
4819                 }
4820             }
4821             j+=(ssize_t) (image->columns+width);
4822           }
4823           if (gamma != 0.0)
4824             {
4825               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4826               if ((channel & RedChannel) != 0)
4827                 q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
4828               if ((channel & GreenChannel) != 0)
4829                 q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
4830               if ((channel & BlueChannel) != 0)
4831                 q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
4832             }
4833           if ((channel & OpacityChannel) != 0)
4834             {
4835               gamma=0.0;
4836               j=0;
4837               for (v=0; v < (ssize_t) width; v++)
4838               {
4839                 for (u=0; u < (ssize_t) width; u++)
4840                 {
4841                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4842                     {
4843                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4844                       gamma+=(*k);
4845                       k++;
4846                     }
4847                 }
4848                 j+=(ssize_t) (image->columns+width);
4849               }
4850               if (gamma != 0.0)
4851                 {
4852                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4853                     gamma);
4854                   SetOpacityPixelComponent(q,ClampToQuantum(gamma*
4855                     GetOpacityPixelComponent(&pixel)));
4856                 }
4857             }
4858           if (((channel & IndexChannel) != 0) &&
4859               (image->colorspace == CMYKColorspace))
4860             {
4861               gamma=0.0;
4862               j=0;
4863               for (v=0; v < (ssize_t) width; v++)
4864               {
4865                 for (u=0; u < (ssize_t) width; u++)
4866                 {
4867                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4868                     {
4869                       pixel.index+=(*k)*indexes[x+u+j];
4870                       gamma+=(*k);
4871                       k++;
4872                     }
4873                 }
4874                 j+=(ssize_t) (image->columns+width);
4875               }
4876               if (gamma != 0.0)
4877                 {
4878                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4879                     gamma);
4880                   blur_indexes[x]=ClampToQuantum(gamma*
4881                     GetIndexPixelComponent(&pixel));
4882                 }
4883             }
4884         }
4885       else
4886         {
4887           MagickRealType
4888             alpha;
4889
4890           for (v=0; v < (ssize_t) width; v++)
4891           {
4892             for (u=0; u < (ssize_t) width; u++)
4893             {
4894               if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4895                 {
4896                   alpha=(MagickRealType) (QuantumScale*
4897                     GetAlphaPixelComponent(p+u+j));
4898                   pixel.red+=(*k)*alpha*(p+u+j)->red;
4899                   pixel.green+=(*k)*alpha*(p+u+j)->green;
4900                   pixel.blue+=(*k)*alpha*(p+u+j)->blue;
4901                   pixel.opacity+=(*k)*(p+u+j)->opacity;
4902                   gamma+=(*k)*alpha;
4903                   k++;
4904                 }
4905             }
4906             j+=(ssize_t) (image->columns+width);
4907           }
4908           if (gamma != 0.0)
4909             {
4910               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
4911               if ((channel & RedChannel) != 0)
4912                 q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
4913               if ((channel & GreenChannel) != 0)
4914                 q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
4915               if ((channel & BlueChannel) != 0)
4916                 q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
4917             }
4918           if ((channel & OpacityChannel) != 0)
4919             {
4920               gamma=0.0;
4921               j=0;
4922               for (v=0; v < (ssize_t) width; v++)
4923               {
4924                 for (u=0; u < (ssize_t) width; u++)
4925                 {
4926                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4927                     {
4928                       pixel.opacity+=(*k)*(p+u+j)->opacity;
4929                       gamma+=(*k);
4930                       k++;
4931                     }
4932                 }
4933                 j+=(ssize_t) (image->columns+width);
4934               }
4935               if (gamma != 0.0)
4936                 {
4937                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4938                     gamma);
4939                   SetOpacityPixelComponent(q,
4940                     ClampOpacityPixelComponent(&pixel));
4941                 }
4942             }
4943           if (((channel & IndexChannel) != 0) &&
4944               (image->colorspace == CMYKColorspace))
4945             {
4946               gamma=0.0;
4947               j=0;
4948               for (v=0; v < (ssize_t) width; v++)
4949               {
4950                 for (u=0; u < (ssize_t) width; u++)
4951                 {
4952                   if (SelectiveContrast(p+u+j,q,threshold) != MagickFalse)
4953                     {
4954                       alpha=(MagickRealType) (QuantumScale*
4955                         GetAlphaPixelComponent(p+u+j));
4956                       pixel.index+=(*k)*alpha*indexes[x+u+j];
4957                       gamma+=(*k);
4958                       k++;
4959                     }
4960                 }
4961                 j+=(ssize_t) (image->columns+width);
4962               }
4963               if (gamma != 0.0)
4964                 {
4965                   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
4966                     gamma);
4967                   blur_indexes[x]=ClampToQuantum(gamma*
4968                     GetIndexPixelComponent(&pixel));
4969                 }
4970             }
4971         }
4972       p++;
4973       q++;
4974     }
4975     sync=SyncCacheViewAuthenticPixels(blur_view,exception);
4976     if (sync == MagickFalse)
4977       status=MagickFalse;
4978     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4979       {
4980         MagickBooleanType
4981           proceed;
4982
4983 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4984   #pragma omp critical (MagickCore_SelectiveBlurImageChannel)
4985 #endif
4986         proceed=SetImageProgress(image,SelectiveBlurImageTag,progress++,
4987           image->rows);
4988         if (proceed == MagickFalse)
4989           status=MagickFalse;
4990       }
4991   }
4992   blur_image->type=image->type;
4993   blur_view=DestroyCacheView(blur_view);
4994   image_view=DestroyCacheView(image_view);
4995   kernel=(double *) RelinquishMagickMemory(kernel);
4996   if (status == MagickFalse)
4997     blur_image=DestroyImage(blur_image);
4998   return(blur_image);
4999 }
5000 \f
5001 /*
5002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5003 %                                                                             %
5004 %                                                                             %
5005 %                                                                             %
5006 %     S h a d e I m a g e                                                     %
5007 %                                                                             %
5008 %                                                                             %
5009 %                                                                             %
5010 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5011 %
5012 %  ShadeImage() shines a distant light on an image to create a
5013 %  three-dimensional effect. You control the positioning of the light with
5014 %  azimuth and elevation; azimuth is measured in degrees off the x axis
5015 %  and elevation is measured in pixels above the Z axis.
5016 %
5017 %  The format of the ShadeImage method is:
5018 %
5019 %      Image *ShadeImage(const Image *image,const MagickBooleanType gray,
5020 %        const double azimuth,const double elevation,ExceptionInfo *exception)
5021 %
5022 %  A description of each parameter follows:
5023 %
5024 %    o image: the image.
5025 %
5026 %    o gray: A value other than zero shades the intensity of each pixel.
5027 %
5028 %    o azimuth, elevation:  Define the light source direction.
5029 %
5030 %    o exception: return any errors or warnings in this structure.
5031 %
5032 */
5033 MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
5034   const double azimuth,const double elevation,ExceptionInfo *exception)
5035 {
5036 #define ShadeImageTag  "Shade/Image"
5037
5038   CacheView
5039     *image_view,
5040     *shade_view;
5041
5042   Image
5043     *shade_image;
5044
5045   MagickBooleanType
5046     status;
5047
5048   MagickOffsetType
5049     progress;
5050
5051   PrimaryInfo
5052     light;
5053
5054   ssize_t
5055     y;
5056
5057   /*
5058     Initialize shaded image attributes.
5059   */
5060   assert(image != (const Image *) NULL);
5061   assert(image->signature == MagickSignature);
5062   if (image->debug != MagickFalse)
5063     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5064   assert(exception != (ExceptionInfo *) NULL);
5065   assert(exception->signature == MagickSignature);
5066   shade_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
5067   if (shade_image == (Image *) NULL)
5068     return((Image *) NULL);
5069   if (SetImageStorageClass(shade_image,DirectClass) == MagickFalse)
5070     {
5071       InheritException(exception,&shade_image->exception);
5072       shade_image=DestroyImage(shade_image);
5073       return((Image *) NULL);
5074     }
5075   /*
5076     Compute the light vector.
5077   */
5078   light.x=(double) QuantumRange*cos(DegreesToRadians(azimuth))*
5079     cos(DegreesToRadians(elevation));
5080   light.y=(double) QuantumRange*sin(DegreesToRadians(azimuth))*
5081     cos(DegreesToRadians(elevation));
5082   light.z=(double) QuantumRange*sin(DegreesToRadians(elevation));
5083   /*
5084     Shade image.
5085   */
5086   status=MagickTrue;
5087   progress=0;
5088   image_view=AcquireCacheView(image);
5089   shade_view=AcquireCacheView(shade_image);
5090 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5091   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
5092 #endif
5093   for (y=0; y < (ssize_t) image->rows; y++)
5094   {
5095     MagickRealType
5096       distance,
5097       normal_distance,
5098       shade;
5099
5100     PrimaryInfo
5101       normal;
5102
5103     register const PixelPacket
5104       *restrict p,
5105       *restrict s0,
5106       *restrict s1,
5107       *restrict s2;
5108
5109     register PixelPacket
5110       *restrict q;
5111
5112     register ssize_t
5113       x;
5114
5115     if (status == MagickFalse)
5116       continue;
5117     p=GetCacheViewVirtualPixels(image_view,-1,y-1,image->columns+2,3,exception);
5118     q=QueueCacheViewAuthenticPixels(shade_view,0,y,shade_image->columns,1,
5119       exception);
5120     if ((p == (PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
5121       {
5122         status=MagickFalse;
5123         continue;
5124       }
5125     /*
5126       Shade this row of pixels.
5127     */
5128     normal.z=2.0*(double) QuantumRange;  /* constant Z of surface normal */
5129     s0=p+1;
5130     s1=s0+image->columns+2;
5131     s2=s1+image->columns+2;
5132     for (x=0; x < (ssize_t) image->columns; x++)
5133     {
5134       /*
5135         Determine the surface normal and compute shading.
5136       */
5137       normal.x=(double) (PixelIntensity(s0-1)+PixelIntensity(s1-1)+
5138         PixelIntensity(s2-1)-PixelIntensity(s0+1)-PixelIntensity(s1+1)-
5139         PixelIntensity(s2+1));
5140       normal.y=(double) (PixelIntensity(s2-1)+PixelIntensity(s2)+
5141         PixelIntensity(s2+1)-PixelIntensity(s0-1)-PixelIntensity(s0)-
5142         PixelIntensity(s0+1));
5143       if ((normal.x == 0.0) && (normal.y == 0.0))
5144         shade=light.z;
5145       else
5146         {
5147           shade=0.0;
5148           distance=normal.x*light.x+normal.y*light.y+normal.z*light.z;
5149           if (distance > MagickEpsilon)
5150             {
5151               normal_distance=
5152                 normal.x*normal.x+normal.y*normal.y+normal.z*normal.z;
5153               if (normal_distance > (MagickEpsilon*MagickEpsilon))
5154                 shade=distance/sqrt((double) normal_distance);
5155             }
5156         }
5157       if (gray != MagickFalse)
5158         {
5159           q->red=(Quantum) shade;
5160           q->green=(Quantum) shade;
5161           q->blue=(Quantum) shade;
5162         }
5163       else
5164         {
5165           q->red=ClampToQuantum(QuantumScale*shade*s1->red);
5166           q->green=ClampToQuantum(QuantumScale*shade*s1->green);
5167           q->blue=ClampToQuantum(QuantumScale*shade*s1->blue);
5168         }
5169       q->opacity=s1->opacity;
5170       s0++;
5171       s1++;
5172       s2++;
5173       q++;
5174     }
5175     if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
5176       status=MagickFalse;
5177     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5178       {
5179         MagickBooleanType
5180           proceed;
5181
5182 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5183   #pragma omp critical (MagickCore_ShadeImage)
5184 #endif
5185         proceed=SetImageProgress(image,ShadeImageTag,progress++,image->rows);
5186         if (proceed == MagickFalse)
5187           status=MagickFalse;
5188       }
5189   }
5190   shade_view=DestroyCacheView(shade_view);
5191   image_view=DestroyCacheView(image_view);
5192   if (status == MagickFalse)
5193     shade_image=DestroyImage(shade_image);
5194   return(shade_image);
5195 }
5196 \f
5197 /*
5198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5199 %                                                                             %
5200 %                                                                             %
5201 %                                                                             %
5202 %     S h a r p e n I m a g e                                                 %
5203 %                                                                             %
5204 %                                                                             %
5205 %                                                                             %
5206 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5207 %
5208 %  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
5209 %  operator of the given radius and standard deviation (sigma).  For
5210 %  reasonable results, radius should be larger than sigma.  Use a radius of 0
5211 %  and SharpenImage() selects a suitable radius for you.
5212 %
5213 %  Using a separable kernel would be faster, but the negative weights cancel
5214 %  out on the corners of the kernel producing often undesirable ringing in the
5215 %  filtered result; this can be avoided by using a 2D gaussian shaped image
5216 %  sharpening kernel instead.
5217 %
5218 %  The format of the SharpenImage method is:
5219 %
5220 %    Image *SharpenImage(const Image *image,const double radius,
5221 %      const double sigma,ExceptionInfo *exception)
5222 %    Image *SharpenImageChannel(const Image *image,const ChannelType channel,
5223 %      const double radius,const double sigma,ExceptionInfo *exception)
5224 %
5225 %  A description of each parameter follows:
5226 %
5227 %    o image: the image.
5228 %
5229 %    o channel: the channel type.
5230 %
5231 %    o radius: the radius of the Gaussian, in pixels, not counting the center
5232 %      pixel.
5233 %
5234 %    o sigma: the standard deviation of the Laplacian, in pixels.
5235 %
5236 %    o exception: return any errors or warnings in this structure.
5237 %
5238 */
5239
5240 MagickExport Image *SharpenImage(const Image *image,const double radius,
5241   const double sigma,ExceptionInfo *exception)
5242 {
5243   Image
5244     *sharp_image;
5245
5246   sharp_image=SharpenImageChannel(image,DefaultChannels,radius,sigma,exception);
5247   return(sharp_image);
5248 }
5249
5250 MagickExport Image *SharpenImageChannel(const Image *image,
5251   const ChannelType channel,const double radius,const double sigma,
5252   ExceptionInfo *exception)
5253 {
5254   double
5255     *kernel,
5256     normalize;
5257
5258   Image
5259     *sharp_image;
5260
5261   register ssize_t
5262     i;
5263
5264   size_t
5265     width;
5266
5267   ssize_t
5268     j,
5269     u,
5270     v;
5271
5272   assert(image != (const Image *) NULL);
5273   assert(image->signature == MagickSignature);
5274   if (image->debug != MagickFalse)
5275     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5276   assert(exception != (ExceptionInfo *) NULL);
5277   assert(exception->signature == MagickSignature);
5278   width=GetOptimalKernelWidth2D(radius,sigma);
5279   kernel=(double *) AcquireQuantumMemory((size_t) width*width,sizeof(*kernel));
5280   if (kernel == (double *) NULL)
5281     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
5282   normalize=0.0;
5283   j=(ssize_t) width/2;
5284   i=0;
5285   for (v=(-j); v <= j; v++)
5286   {
5287     for (u=(-j); u <= j; u++)
5288     {
5289       kernel[i]=(double) (-exp(-((double) u*u+v*v)/(2.0*MagickSigma*
5290         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
5291       normalize+=kernel[i];
5292       i++;
5293     }
5294   }
5295   kernel[i/2]=(double) ((-2.0)*normalize);
5296   sharp_image=ConvolveImageChannel(image,channel,width,kernel,exception);
5297   kernel=(double *) RelinquishMagickMemory(kernel);
5298   return(sharp_image);
5299 }
5300 \f
5301 /*
5302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5303 %                                                                             %
5304 %                                                                             %
5305 %                                                                             %
5306 %     S p r e a d I m a g e                                                   %
5307 %                                                                             %
5308 %                                                                             %
5309 %                                                                             %
5310 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5311 %
5312 %  SpreadImage() is a special effects method that randomly displaces each
5313 %  pixel in a block defined by the radius parameter.
5314 %
5315 %  The format of the SpreadImage method is:
5316 %
5317 %      Image *SpreadImage(const Image *image,const double radius,
5318 %        ExceptionInfo *exception)
5319 %
5320 %  A description of each parameter follows:
5321 %
5322 %    o image: the image.
5323 %
5324 %    o radius:  Choose a random pixel in a neighborhood of this extent.
5325 %
5326 %    o exception: return any errors or warnings in this structure.
5327 %
5328 */
5329 MagickExport Image *SpreadImage(const Image *image,const double radius,
5330   ExceptionInfo *exception)
5331 {
5332 #define SpreadImageTag  "Spread/Image"
5333
5334   CacheView
5335     *image_view;
5336
5337   Image
5338     *spread_image;
5339
5340   MagickBooleanType
5341     status;
5342
5343   MagickOffsetType
5344     progress;
5345
5346   MagickPixelPacket
5347     bias;
5348
5349   RandomInfo
5350     **restrict random_info;
5351
5352   ResampleFilter
5353     **restrict resample_filter;
5354
5355   size_t
5356     width;
5357
5358   ssize_t
5359     y;
5360
5361   /*
5362     Initialize spread image attributes.
5363   */
5364   assert(image != (Image *) NULL);
5365   assert(image->signature == MagickSignature);
5366   if (image->debug != MagickFalse)
5367     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5368   assert(exception != (ExceptionInfo *) NULL);
5369   assert(exception->signature == MagickSignature);
5370   spread_image=CloneImage(image,image->columns,image->rows,MagickTrue,
5371     exception);
5372   if (spread_image == (Image *) NULL)
5373     return((Image *) NULL);
5374   if (SetImageStorageClass(spread_image,DirectClass) == MagickFalse)
5375     {
5376       InheritException(exception,&spread_image->exception);
5377       spread_image=DestroyImage(spread_image);
5378       return((Image *) NULL);
5379     }
5380   /*
5381     Spread image.
5382   */
5383   status=MagickTrue;
5384   progress=0;
5385   GetMagickPixelPacket(spread_image,&bias);
5386   width=GetOptimalKernelWidth1D(radius,0.5);
5387   resample_filter=AcquireResampleFilterThreadSet(image,
5388     UndefinedVirtualPixelMethod,MagickTrue,exception);
5389   random_info=AcquireRandomInfoThreadSet();
5390   image_view=AcquireCacheView(spread_image);
5391 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5392   #pragma omp parallel for schedule(dynamic,4) shared(progress,status) omp_throttle(1)
5393 #endif
5394   for (y=0; y < (ssize_t) spread_image->rows; y++)
5395   {
5396     const int
5397       id = GetOpenMPThreadId();
5398
5399     MagickPixelPacket
5400       pixel;
5401
5402     register IndexPacket
5403       *restrict indexes;
5404
5405     register PixelPacket
5406       *restrict q;
5407
5408     register ssize_t
5409       x;
5410
5411     if (status == MagickFalse)
5412       continue;
5413     q=QueueCacheViewAuthenticPixels(image_view,0,y,spread_image->columns,1,
5414       exception);
5415     if (q == (PixelPacket *) NULL)
5416       {
5417         status=MagickFalse;
5418         continue;
5419       }
5420     indexes=GetCacheViewAuthenticIndexQueue(image_view);
5421     pixel=bias;
5422     for (x=0; x < (ssize_t) spread_image->columns; x++)
5423     {
5424       (void) ResamplePixelColor(resample_filter[id],(double) x+width*
5425         (GetPseudoRandomValue(random_info[id])-0.5),(double) y+width*
5426         (GetPseudoRandomValue(random_info[id])-0.5),&pixel);
5427       SetPixelPacket(spread_image,&pixel,q,indexes+x);
5428       q++;
5429     }
5430     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
5431       status=MagickFalse;
5432     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5433       {
5434         MagickBooleanType
5435           proceed;
5436
5437 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5438   #pragma omp critical (MagickCore_SpreadImage)
5439 #endif
5440         proceed=SetImageProgress(image,SpreadImageTag,progress++,image->rows);
5441         if (proceed == MagickFalse)
5442           status=MagickFalse;
5443       }
5444   }
5445   image_view=DestroyCacheView(image_view);
5446   random_info=DestroyRandomInfoThreadSet(random_info);
5447   resample_filter=DestroyResampleFilterThreadSet(resample_filter);
5448   return(spread_image);
5449 }
5450 \f
5451 /*
5452 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5453 %                                                                             %
5454 %                                                                             %
5455 %                                                                             %
5456 %     S t a t i s t i c I m a g e                                             %
5457 %                                                                             %
5458 %                                                                             %
5459 %                                                                             %
5460 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5461 %
5462 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
5463 %  the neighborhood of the specified radius.
5464 %
5465 %  The format of the StatisticImage method is:
5466 %
5467 %      Image *StatisticImage(const Image *image,const StatisticType type,
5468 %        const double radius,ExceptionInfo *exception)
5469 %      Image *StatisticImageChannel(const Image *image,
5470 %        const ChannelType channel,const StatisticType type,
5471 %        const double radius,ExceptionInfo *exception)
5472 %
5473 %  A description of each parameter follows:
5474 %
5475 %    o image: the image.
5476 %
5477 %    o channel: the image channel.
5478 %
5479 %    o type: the statistic type (median, mode, etc.).
5480 %
5481 %    o radius: the radius of the pixel neighborhood.
5482 %
5483 %    o exception: return any errors or warnings in this structure.
5484 %
5485 */
5486
5487 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
5488   const double radius,ExceptionInfo *exception)
5489 {
5490   return(StatisticImageChannel(image,DefaultChannels,type,radius,exception));
5491 }
5492
5493 MagickExport Image *StatisticImageChannel(const Image *image,
5494   const ChannelType channel,const StatisticType type,const double radius,
5495   ExceptionInfo *exception)
5496 {
5497 #define StatisticImageTag  "Statistic/Image"
5498
5499   CacheView
5500     *image_view,
5501     *statistic_view;
5502
5503   Image
5504     *statistic_image;
5505
5506   MagickBooleanType
5507     status;
5508
5509   MagickOffsetType
5510     progress;
5511
5512   PixelList
5513     **restrict pixel_list;
5514
5515   size_t
5516     width;
5517
5518   ssize_t
5519     y;
5520
5521   /*
5522     Initialize statistics image attributes.
5523   */
5524   assert(image != (Image *) NULL);
5525   assert(image->signature == MagickSignature);
5526   if (image->debug != MagickFalse)
5527     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5528   assert(exception != (ExceptionInfo *) NULL);
5529   assert(exception->signature == MagickSignature);
5530   width=GetOptimalKernelWidth2D(radius,0.5);
5531   statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
5532     exception);
5533   if (statistic_image == (Image *) NULL)
5534     return((Image *) NULL);
5535   if (SetImageStorageClass(statistic_image,DirectClass) == MagickFalse)
5536     {
5537       InheritException(exception,&statistic_image->exception);
5538       statistic_image=DestroyImage(statistic_image);
5539       return((Image *) NULL);
5540     }
5541   pixel_list=AcquirePixelListThreadSet(width);
5542   if (pixel_list == (PixelList **) NULL)
5543     {
5544       statistic_image=DestroyImage(statistic_image);
5545       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
5546     }
5547   /*
5548     Reduce statistics image.
5549   */
5550   status=MagickTrue;
5551   progress=0;
5552   image_view=AcquireCacheView(image);
5553   statistic_view=AcquireCacheView(statistic_image);
5554 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5555   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
5556 #endif
5557   for (y=0; y < (ssize_t) statistic_image->rows; y++)
5558   {
5559     const int
5560       id = GetOpenMPThreadId();
5561
5562     register const IndexPacket
5563       *restrict indexes;
5564
5565     register const PixelPacket
5566       *restrict p;
5567
5568     register IndexPacket
5569       *restrict statistic_indexes;
5570
5571     register PixelPacket
5572       *restrict q;
5573
5574     register ssize_t
5575       x;
5576
5577     if (status == MagickFalse)
5578       continue;
5579     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
5580       (width/2L),image->columns+width,width,exception);
5581     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,
5582       statistic_image->columns,1,exception);
5583     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
5584       {
5585         status=MagickFalse;
5586         continue;
5587       }
5588     indexes=GetCacheViewVirtualIndexQueue(image_view);
5589     statistic_indexes=GetCacheViewAuthenticIndexQueue(statistic_view);
5590     for (x=0; x < (ssize_t) statistic_image->columns; x++)
5591     {
5592       MagickPixelPacket
5593         pixel;
5594
5595       register const PixelPacket
5596         *restrict r;
5597
5598       register const IndexPacket
5599         *restrict s;
5600
5601       register ssize_t
5602         u,
5603         v;
5604
5605       r=p;
5606       s=indexes+x;
5607       ResetPixelList(pixel_list[id]);
5608       for (v=0; v < (ssize_t) width; v++)
5609       {
5610         for (u=0; u < (ssize_t) width; u++)
5611           InsertPixelList(image,r+u,s+u,pixel_list[id]);
5612         r+=image->columns+width;
5613         s+=image->columns+width;
5614       }
5615       switch (type)
5616       {
5617         case ModeStatistic: pixel=GetModePixelList(pixel_list[id]);
5618         default: break;
5619       }
5620       if ((channel & RedChannel) != 0)
5621         q->red=ClampToQuantum(pixel.red);
5622       if ((channel & GreenChannel) != 0)
5623         q->green=ClampToQuantum(pixel.green);
5624       if ((channel & BlueChannel) != 0)
5625         q->blue=ClampToQuantum(pixel.blue);
5626       if ((image->matte != MagickFalse) && ((channel & OpacityChannel) != 0))
5627         q->opacity=ClampToQuantum(pixel.opacity);
5628       if (((channel & IndexChannel) != 0) &&
5629           (image->colorspace == CMYKColorspace))
5630         statistic_indexes[x]=(IndexPacket) ClampToQuantum(pixel.index);
5631       p++;
5632       q++;
5633     }
5634     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
5635       status=MagickFalse;
5636     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5637       {
5638         MagickBooleanType
5639           proceed;
5640
5641 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5642   #pragma omp critical (MagickCore_StatisticImage)
5643 #endif
5644         proceed=SetImageProgress(image,StatisticImageTag,progress++,
5645           image->rows);
5646         if (proceed == MagickFalse)
5647           status=MagickFalse;
5648       }
5649   }
5650   statistic_view=DestroyCacheView(statistic_view);
5651   image_view=DestroyCacheView(image_view);
5652   pixel_list=DestroyPixelListThreadSet(pixel_list);
5653   return(statistic_image);
5654 }
5655 \f
5656 /*
5657 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5658 %                                                                             %
5659 %                                                                             %
5660 %                                                                             %
5661 %     U n s h a r p M a s k I m a g e                                         %
5662 %                                                                             %
5663 %                                                                             %
5664 %                                                                             %
5665 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5666 %
5667 %  UnsharpMaskImage() sharpens one or more image channels.  We convolve the
5668 %  image with a Gaussian operator of the given radius and standard deviation
5669 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
5670 %  radius of 0 and UnsharpMaskImage() selects a suitable radius for you.
5671 %
5672 %  The format of the UnsharpMaskImage method is:
5673 %
5674 %    Image *UnsharpMaskImage(const Image *image,const double radius,
5675 %      const double sigma,const double amount,const double threshold,
5676 %      ExceptionInfo *exception)
5677 %    Image *UnsharpMaskImageChannel(const Image *image,
5678 %      const ChannelType channel,const double radius,const double sigma,
5679 %      const double amount,const double threshold,ExceptionInfo *exception)
5680 %
5681 %  A description of each parameter follows:
5682 %
5683 %    o image: the image.
5684 %
5685 %    o channel: the channel type.
5686 %
5687 %    o radius: the radius of the Gaussian, in pixels, not counting the center
5688 %      pixel.
5689 %
5690 %    o sigma: the standard deviation of the Gaussian, in pixels.
5691 %
5692 %    o amount: the percentage of the difference between the original and the
5693 %      blur image that is added back into the original.
5694 %
5695 %    o threshold: the threshold in pixels needed to apply the diffence amount.
5696 %
5697 %    o exception: return any errors or warnings in this structure.
5698 %
5699 */
5700
5701 MagickExport Image *UnsharpMaskImage(const Image *image,const double radius,
5702   const double sigma,const double amount,const double threshold,
5703   ExceptionInfo *exception)
5704 {
5705   Image
5706     *sharp_image;
5707
5708   sharp_image=UnsharpMaskImageChannel(image,DefaultChannels,radius,sigma,amount,
5709     threshold,exception);
5710   return(sharp_image);
5711 }
5712
5713 MagickExport Image *UnsharpMaskImageChannel(const Image *image,
5714   const ChannelType channel,const double radius,const double sigma,
5715   const double amount,const double threshold,ExceptionInfo *exception)
5716 {
5717 #define SharpenImageTag  "Sharpen/Image"
5718
5719   CacheView
5720     *image_view,
5721     *unsharp_view;
5722
5723   Image
5724     *unsharp_image;
5725
5726   MagickBooleanType
5727     status;
5728
5729   MagickOffsetType
5730     progress;
5731
5732   MagickPixelPacket
5733     bias;
5734
5735   MagickRealType
5736     quantum_threshold;
5737
5738   ssize_t
5739     y;
5740
5741   assert(image != (const Image *) NULL);
5742   assert(image->signature == MagickSignature);
5743   if (image->debug != MagickFalse)
5744     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5745   assert(exception != (ExceptionInfo *) NULL);
5746   unsharp_image=BlurImageChannel(image,channel,radius,sigma,exception);
5747   if (unsharp_image == (Image *) NULL)
5748     return((Image *) NULL);
5749   quantum_threshold=(MagickRealType) QuantumRange*threshold;
5750   /*
5751     Unsharp-mask image.
5752   */
5753   status=MagickTrue;
5754   progress=0;
5755   GetMagickPixelPacket(image,&bias);
5756   image_view=AcquireCacheView(image);
5757   unsharp_view=AcquireCacheView(unsharp_image);
5758 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5759   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
5760 #endif
5761   for (y=0; y < (ssize_t) image->rows; y++)
5762   {
5763     MagickPixelPacket
5764       pixel;
5765
5766     register const IndexPacket
5767       *restrict indexes;
5768
5769     register const PixelPacket
5770       *restrict p;
5771
5772     register IndexPacket
5773       *restrict unsharp_indexes;
5774
5775     register PixelPacket
5776       *restrict q;
5777
5778     register ssize_t
5779       x;
5780
5781     if (status == MagickFalse)
5782       continue;
5783     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
5784     q=GetCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
5785       exception);
5786     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
5787       {
5788         status=MagickFalse;
5789         continue;
5790       }
5791     indexes=GetCacheViewVirtualIndexQueue(image_view);
5792     unsharp_indexes=GetCacheViewAuthenticIndexQueue(unsharp_view);
5793     pixel=bias;
5794     for (x=0; x < (ssize_t) image->columns; x++)
5795     {
5796       if ((channel & RedChannel) != 0)
5797         {
5798           pixel.red=p->red-(MagickRealType) q->red;
5799           if (fabs(2.0*pixel.red) < quantum_threshold)
5800             pixel.red=(MagickRealType) GetRedPixelComponent(p);
5801           else
5802             pixel.red=(MagickRealType) p->red+(pixel.red*amount);
5803           SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
5804         }
5805       if ((channel & GreenChannel) != 0)
5806         {
5807           pixel.green=p->green-(MagickRealType) q->green;
5808           if (fabs(2.0*pixel.green) < quantum_threshold)
5809             pixel.green=(MagickRealType) GetGreenPixelComponent(p);
5810           else
5811             pixel.green=(MagickRealType) p->green+(pixel.green*amount);
5812           SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
5813         }
5814       if ((channel & BlueChannel) != 0)
5815         {
5816           pixel.blue=p->blue-(MagickRealType) q->blue;
5817           if (fabs(2.0*pixel.blue) < quantum_threshold)
5818             pixel.blue=(MagickRealType) GetBluePixelComponent(p);
5819           else
5820             pixel.blue=(MagickRealType) p->blue+(pixel.blue*amount);
5821           SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
5822         }
5823       if ((channel & OpacityChannel) != 0)
5824         {
5825           pixel.opacity=p->opacity-(MagickRealType) q->opacity;
5826           if (fabs(2.0*pixel.opacity) < quantum_threshold)
5827             pixel.opacity=(MagickRealType) GetOpacityPixelComponent(p);
5828           else
5829             pixel.opacity=p->opacity+(pixel.opacity*amount);
5830           SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
5831         }
5832       if (((channel & IndexChannel) != 0) &&
5833           (image->colorspace == CMYKColorspace))
5834         {
5835           pixel.index=indexes[x]-(MagickRealType) unsharp_indexes[x];
5836           if (fabs(2.0*pixel.index) < quantum_threshold)
5837             pixel.index=(MagickRealType) indexes[x];
5838           else
5839             pixel.index=(MagickRealType) indexes[x]+(pixel.index*amount);
5840           unsharp_indexes[x]=ClampToQuantum(pixel.index);
5841         }
5842       p++;
5843       q++;
5844     }
5845     if (SyncCacheViewAuthenticPixels(unsharp_view,exception) == MagickFalse)
5846       status=MagickFalse;
5847     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5848       {
5849         MagickBooleanType
5850           proceed;
5851
5852 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5853   #pragma omp critical (MagickCore_UnsharpMaskImageChannel)
5854 #endif
5855         proceed=SetImageProgress(image,SharpenImageTag,progress++,image->rows);
5856         if (proceed == MagickFalse)
5857           status=MagickFalse;
5858       }
5859   }
5860   unsharp_image->type=image->type;
5861   unsharp_view=DestroyCacheView(unsharp_view);
5862   image_view=DestroyCacheView(image_view);
5863   if (status == MagickFalse)
5864     unsharp_image=DestroyImage(unsharp_image);
5865   return(unsharp_image);
5866 }