granicus.if.org Git - imagemagick/blob - MagickCore/fx.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %                                 FFFFF  X   X                                %
   7 %                                 F       X X                                 %
   8 %                                 FFF      X                                  %
   9 %                                 F       X X                                 %
  10 %                                 F      X   X                                %
  11 %                                                                             %
  12 %                                                                             %
  13 %                   MagickCore Image Special Effects Methods                  %
  14 %                                                                             %
  15 %                               Software Design                               %
  16 %                                    Cristy                                   %
  17 %                                 October 1996                                %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2017 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 %    https://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 "MagickCore/studio.h"
  44 #include "MagickCore/accelerate-private.h"
  45 #include "MagickCore/annotate.h"
  46 #include "MagickCore/artifact.h"
  47 #include "MagickCore/attribute.h"
  48 #include "MagickCore/cache.h"
  49 #include "MagickCore/cache-view.h"
  50 #include "MagickCore/channel.h"
  51 #include "MagickCore/color.h"
  52 #include "MagickCore/color-private.h"
  53 #include "MagickCore/colorspace-private.h"
  54 #include "MagickCore/composite.h"
  55 #include "MagickCore/decorate.h"
  56 #include "MagickCore/distort.h"
  57 #include "MagickCore/draw.h"
  58 #include "MagickCore/effect.h"
  59 #include "MagickCore/enhance.h"
  60 #include "MagickCore/exception.h"
  61 #include "MagickCore/exception-private.h"
  62 #include "MagickCore/fx.h"
  63 #include "MagickCore/fx-private.h"
  64 #include "MagickCore/gem.h"
  65 #include "MagickCore/gem-private.h"
  66 #include "MagickCore/geometry.h"
  67 #include "MagickCore/layer.h"
  68 #include "MagickCore/list.h"
  69 #include "MagickCore/log.h"
  70 #include "MagickCore/image.h"
  71 #include "MagickCore/image-private.h"
  72 #include "MagickCore/magick.h"
  73 #include "MagickCore/memory_.h"
  74 #include "MagickCore/monitor.h"
  75 #include "MagickCore/monitor-private.h"
  76 #include "MagickCore/option.h"
  77 #include "MagickCore/pixel.h"
  78 #include "MagickCore/pixel-accessor.h"
  79 #include "MagickCore/property.h"
  80 #include "MagickCore/quantum.h"
  81 #include "MagickCore/quantum-private.h"
  82 #include "MagickCore/random_.h"
  83 #include "MagickCore/random-private.h"
  84 #include "MagickCore/resample.h"
  85 #include "MagickCore/resample-private.h"
  86 #include "MagickCore/resize.h"
  87 #include "MagickCore/resource_.h"
  88 #include "MagickCore/splay-tree.h"
  89 #include "MagickCore/statistic.h"
  90 #include "MagickCore/string_.h"
  91 #include "MagickCore/string-private.h"
  92 #include "MagickCore/thread-private.h"
  93 #include "MagickCore/transform.h"
  94 #include "MagickCore/transform-private.h"
  95 #include "MagickCore/utility.h"
  96 \f
  97 /*
  98   Define declarations.
  99 */
 100 #define LeftShiftOperator  0xf5U
 101 #define RightShiftOperator  0xf6U
 102 #define LessThanEqualOperator  0xf7U
 103 #define GreaterThanEqualOperator  0xf8U
 104 #define EqualOperator  0xf9U
 105 #define NotEqualOperator  0xfaU
 106 #define LogicalAndOperator  0xfbU
 107 #define LogicalOrOperator  0xfcU
 108 #define ExponentialNotation 0xfdU
 109
 110 struct _FxInfo
 111 {
 112   const Image
 113     *images;
 114
 115   char
 116     *expression;
 117
 118   FILE
 119     *file;
 120
 121   SplayTreeInfo
 122     *colors,
 123     *symbols;
 124
 125   CacheView
 126     **view;
 127
 128   RandomInfo
 129     *random_info;
 130
 131   ExceptionInfo
 132     *exception;
 133 };
 134 \f
 135 /*
 136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 137 %                                                                             %
 138 %                                                                             %
 139 %                                                                             %
 140 +   A c q u i r e F x I n f o                                                 %
 141 %                                                                             %
 142 %                                                                             %
 143 %                                                                             %
 144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 145 %
 146 %  AcquireFxInfo() allocates the FxInfo structure.
 147 %
 148 %  The format of the AcquireFxInfo method is:
 149 %
 150 %      FxInfo *AcquireFxInfo(Image *images,const char *expression,
 151 %        ExceptionInfo *exception)
 152 %
 153 %  A description of each parameter follows:
 154 %
 155 %    o images: the image sequence.
 156 %
 157 %    o expression: the expression.
 158 %
 159 %    o exception: return any errors or warnings in this structure.
 160 %
 161 */
 162 MagickPrivate FxInfo *AcquireFxInfo(const Image *images,const char *expression,
 163   ExceptionInfo *exception)
 164 {
 165   char
 166     fx_op[2];
 167
 168   const Image
 169     *next;
 170
 171   FxInfo
 172     *fx_info;
 173
 174   register ssize_t
 175     i;
 176
 177   fx_info=(FxInfo *) AcquireMagickMemory(sizeof(*fx_info));
 178   if (fx_info == (FxInfo *) NULL)
 179     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 180   (void) ResetMagickMemory(fx_info,0,sizeof(*fx_info));
 181   fx_info->exception=AcquireExceptionInfo();
 182   fx_info->images=images;
 183   fx_info->colors=NewSplayTree(CompareSplayTreeString,RelinquishMagickMemory,
 184     RelinquishMagickMemory);
 185   fx_info->symbols=NewSplayTree(CompareSplayTreeString,RelinquishMagickMemory,
 186     RelinquishMagickMemory);
 187   fx_info->view=(CacheView **) AcquireQuantumMemory(GetImageListLength(
 188     fx_info->images),sizeof(*fx_info->view));
 189   if (fx_info->view == (CacheView **) NULL)
 190     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 191   i=0;
 192   next=GetFirstImageInList(fx_info->images);
 193   for ( ; next != (Image *) NULL; next=next->next)
 194   {
 195     fx_info->view[i]=AcquireVirtualCacheView(next,exception);
 196     i++;
 197   }
 198   fx_info->random_info=AcquireRandomInfo();
 199   fx_info->expression=ConstantString(expression);
 200   fx_info->file=stderr;
 201   (void) SubstituteString(&fx_info->expression," ","");  /* compact string */
 202   /*
 203     Force right-to-left associativity for unary negation.
 204   */
 205   (void) SubstituteString(&fx_info->expression,"-","-1.0*");
 206   (void) SubstituteString(&fx_info->expression,"^-1.0*","^-");
 207   (void) SubstituteString(&fx_info->expression,"E-1.0*","E-");
 208   (void) SubstituteString(&fx_info->expression,"e-1.0*","e-");
 209   /*
 210     Convert compound to simple operators.
 211   */
 212   fx_op[1]='\0';
 213   *fx_op=(char) LeftShiftOperator;
 214   (void) SubstituteString(&fx_info->expression,"<<",fx_op);
 215   *fx_op=(char) RightShiftOperator;
 216   (void) SubstituteString(&fx_info->expression,">>",fx_op);
 217   *fx_op=(char) LessThanEqualOperator;
 218   (void) SubstituteString(&fx_info->expression,"<=",fx_op);
 219   *fx_op=(char) GreaterThanEqualOperator;
 220   (void) SubstituteString(&fx_info->expression,">=",fx_op);
 221   *fx_op=(char) EqualOperator;
 222   (void) SubstituteString(&fx_info->expression,"==",fx_op);
 223   *fx_op=(char) NotEqualOperator;
 224   (void) SubstituteString(&fx_info->expression,"!=",fx_op);
 225   *fx_op=(char) LogicalAndOperator;
 226   (void) SubstituteString(&fx_info->expression,"&&",fx_op);
 227   *fx_op=(char) LogicalOrOperator;
 228   (void) SubstituteString(&fx_info->expression,"||",fx_op);
 229   *fx_op=(char) ExponentialNotation;
 230   (void) SubstituteString(&fx_info->expression,"**",fx_op);
 231   return(fx_info);
 232 }
 233 \f
 234 /*
 235 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 236 %                                                                             %
 237 %                                                                             %
 238 %                                                                             %
 239 %     A d d N o i s e I m a g e                                               %
 240 %                                                                             %
 241 %                                                                             %
 242 %                                                                             %
 243 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 244 %
 245 %  AddNoiseImage() adds random noise to the image.
 246 %
 247 %  The format of the AddNoiseImage method is:
 248 %
 249 %      Image *AddNoiseImage(const Image *image,const NoiseType noise_type,
 250 %        const double attenuate,ExceptionInfo *exception)
 251 %
 252 %  A description of each parameter follows:
 253 %
 254 %    o image: the image.
 255 %
 256 %    o channel: the channel type.
 257 %
 258 %    o noise_type:  The type of noise: Uniform, Gaussian, Multiplicative,
 259 %      Impulse, Laplacian, or Poisson.
 260 %
 261 %    o attenuate:  attenuate the random distribution.
 262 %
 263 %    o exception: return any errors or warnings in this structure.
 264 %
 265 */
 266 MagickExport Image *AddNoiseImage(const Image *image,const NoiseType noise_type,
 267   const double attenuate,ExceptionInfo *exception)
 268 {
 269 #define AddNoiseImageTag  "AddNoise/Image"
 270
 271   CacheView
 272     *image_view,
 273     *noise_view;
 274
 275   Image
 276     *noise_image;
 277
 278   MagickBooleanType
 279     status;
 280
 281   MagickOffsetType
 282     progress;
 283
 284   RandomInfo
 285     **magick_restrict random_info;
 286
 287   ssize_t
 288     y;
 289
 290 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 291   unsigned long
 292     key;
 293 #endif
 294
 295   /*
 296     Initialize noise image attributes.
 297   */
 298   assert(image != (const Image *) NULL);
 299   assert(image->signature == MagickCoreSignature);
 300   if (image->debug != MagickFalse)
 301     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 302   assert(exception != (ExceptionInfo *) NULL);
 303   assert(exception->signature == MagickCoreSignature);
 304 #if defined(MAGICKCORE_OPENCL_SUPPORT)
 305   noise_image=AccelerateAddNoiseImage(image,noise_type,exception);
 306   if (noise_image != (Image *) NULL)
 307     return(noise_image);
 308 #endif
 309   noise_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
 310   if (noise_image == (Image *) NULL)
 311     return((Image *) NULL);
 312   if (SetImageStorageClass(noise_image,DirectClass,exception) == MagickFalse)
 313     {
 314       noise_image=DestroyImage(noise_image);
 315       return((Image *) NULL);
 316     }
 317   /*
 318     Add noise in each row.
 319   */
 320   status=MagickTrue;
 321   progress=0;
 322   random_info=AcquireRandomInfoThreadSet();
 323   image_view=AcquireVirtualCacheView(image,exception);
 324   noise_view=AcquireAuthenticCacheView(noise_image,exception);
 325 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 326   key=GetRandomSecretKey(random_info[0]);
 327   #pragma omp parallel for schedule(static,4) shared(progress,status) \
 328     magick_threads(image,noise_image,image->rows,key == ~0UL)
 329 #endif
 330   for (y=0; y < (ssize_t) image->rows; y++)
 331   {
 332     const int
 333       id = GetOpenMPThreadId();
 334
 335     MagickBooleanType
 336       sync;
 337
 338     register const Quantum
 339       *magick_restrict p;
 340
 341     register ssize_t
 342       x;
 343
 344     register Quantum
 345       *magick_restrict q;
 346
 347     if (status == MagickFalse)
 348       continue;
 349     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 350     q=QueueCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
 351       exception);
 352     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
 353       {
 354         status=MagickFalse;
 355         continue;
 356       }
 357     for (x=0; x < (ssize_t) image->columns; x++)
 358     {
 359       register ssize_t
 360         i;
 361
 362       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 363       {
 364         PixelChannel channel=GetPixelChannelChannel(image,i);
 365         PixelTrait traits=GetPixelChannelTraits(image,channel);
 366         PixelTrait noise_traits=GetPixelChannelTraits(noise_image,channel);
 367         if ((traits == UndefinedPixelTrait) ||
 368             (noise_traits == UndefinedPixelTrait))
 369           continue;
 370         if (((noise_traits & CopyPixelTrait) != 0) ||
 371             (GetPixelWriteMask(image,p) == 0))
 372           {
 373             SetPixelChannel(noise_image,channel,p[i],q);
 374             continue;
 375           }
 376         SetPixelChannel(noise_image,channel,ClampToQuantum(
 377           GenerateDifferentialNoise(random_info[id],p[i],noise_type,attenuate)),
 378           q);
 379       }
 380       p+=GetPixelChannels(image);
 381       q+=GetPixelChannels(noise_image);
 382     }
 383     sync=SyncCacheViewAuthenticPixels(noise_view,exception);
 384     if (sync == MagickFalse)
 385       status=MagickFalse;
 386     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 387       {
 388         MagickBooleanType
 389           proceed;
 390
 391 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 392         #pragma omp critical (MagickCore_AddNoiseImage)
 393 #endif
 394         proceed=SetImageProgress(image,AddNoiseImageTag,progress++,
 395           image->rows);
 396         if (proceed == MagickFalse)
 397           status=MagickFalse;
 398       }
 399   }
 400   noise_view=DestroyCacheView(noise_view);
 401   image_view=DestroyCacheView(image_view);
 402   random_info=DestroyRandomInfoThreadSet(random_info);
 403   if (status == MagickFalse)
 404     noise_image=DestroyImage(noise_image);
 405   return(noise_image);
 406 }
 407 \f
 408 /*
 409 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 410 %                                                                             %
 411 %                                                                             %
 412 %                                                                             %
 413 %     B l u e S h i f t I m a g e                                             %
 414 %                                                                             %
 415 %                                                                             %
 416 %                                                                             %
 417 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 418 %
 419 %  BlueShiftImage() mutes the colors of the image to simulate a scene at
 420 %  nighttime in the moonlight.
 421 %
 422 %  The format of the BlueShiftImage method is:
 423 %
 424 %      Image *BlueShiftImage(const Image *image,const double factor,
 425 %        ExceptionInfo *exception)
 426 %
 427 %  A description of each parameter follows:
 428 %
 429 %    o image: the image.
 430 %
 431 %    o factor: the shift factor.
 432 %
 433 %    o exception: return any errors or warnings in this structure.
 434 %
 435 */
 436 MagickExport Image *BlueShiftImage(const Image *image,const double factor,
 437   ExceptionInfo *exception)
 438 {
 439 #define BlueShiftImageTag  "BlueShift/Image"
 440
 441   CacheView
 442     *image_view,
 443     *shift_view;
 444
 445   Image
 446     *shift_image;
 447
 448   MagickBooleanType
 449     status;
 450
 451   MagickOffsetType
 452     progress;
 453
 454   ssize_t
 455     y;
 456
 457   /*
 458     Allocate blue shift image.
 459   */
 460   assert(image != (const Image *) NULL);
 461   assert(image->signature == MagickCoreSignature);
 462   if (image->debug != MagickFalse)
 463     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 464   assert(exception != (ExceptionInfo *) NULL);
 465   assert(exception->signature == MagickCoreSignature);
 466   shift_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
 467   if (shift_image == (Image *) NULL)
 468     return((Image *) NULL);
 469   if (SetImageStorageClass(shift_image,DirectClass,exception) == MagickFalse)
 470     {
 471       shift_image=DestroyImage(shift_image);
 472       return((Image *) NULL);
 473     }
 474   /*
 475     Blue-shift DirectClass image.
 476   */
 477   status=MagickTrue;
 478   progress=0;
 479   image_view=AcquireVirtualCacheView(image,exception);
 480   shift_view=AcquireAuthenticCacheView(shift_image,exception);
 481 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 482   #pragma omp parallel for schedule(static,4) shared(progress,status) \
 483     magick_threads(image,shift_image,image->rows,1)
 484 #endif
 485   for (y=0; y < (ssize_t) image->rows; y++)
 486   {
 487     MagickBooleanType
 488       sync;
 489
 490     PixelInfo
 491       pixel;
 492
 493     Quantum
 494       quantum;
 495
 496     register const Quantum
 497       *magick_restrict p;
 498
 499     register ssize_t
 500       x;
 501
 502     register Quantum
 503       *magick_restrict q;
 504
 505     if (status == MagickFalse)
 506       continue;
 507     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 508     q=QueueCacheViewAuthenticPixels(shift_view,0,y,shift_image->columns,1,
 509       exception);
 510     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
 511       {
 512         status=MagickFalse;
 513         continue;
 514       }
 515     for (x=0; x < (ssize_t) image->columns; x++)
 516     {
 517       quantum=GetPixelRed(image,p);
 518       if (GetPixelGreen(image,p) < quantum)
 519         quantum=GetPixelGreen(image,p);
 520       if (GetPixelBlue(image,p) < quantum)
 521         quantum=GetPixelBlue(image,p);
 522       pixel.red=0.5*(GetPixelRed(image,p)+factor*quantum);
 523       pixel.green=0.5*(GetPixelGreen(image,p)+factor*quantum);
 524       pixel.blue=0.5*(GetPixelBlue(image,p)+factor*quantum);
 525       quantum=GetPixelRed(image,p);
 526       if (GetPixelGreen(image,p) > quantum)
 527         quantum=GetPixelGreen(image,p);
 528       if (GetPixelBlue(image,p) > quantum)
 529         quantum=GetPixelBlue(image,p);
 530       pixel.red=0.5*(pixel.red+factor*quantum);
 531       pixel.green=0.5*(pixel.green+factor*quantum);
 532       pixel.blue=0.5*(pixel.blue+factor*quantum);
 533       SetPixelRed(shift_image,ClampToQuantum(pixel.red),q);
 534       SetPixelGreen(shift_image,ClampToQuantum(pixel.green),q);
 535       SetPixelBlue(shift_image,ClampToQuantum(pixel.blue),q);
 536       p+=GetPixelChannels(image);
 537       q+=GetPixelChannels(shift_image);
 538     }
 539     sync=SyncCacheViewAuthenticPixels(shift_view,exception);
 540     if (sync == MagickFalse)
 541       status=MagickFalse;
 542     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 543       {
 544         MagickBooleanType
 545           proceed;
 546
 547 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 548         #pragma omp critical (MagickCore_BlueShiftImage)
 549 #endif
 550         proceed=SetImageProgress(image,BlueShiftImageTag,progress++,
 551           image->rows);
 552         if (proceed == MagickFalse)
 553           status=MagickFalse;
 554       }
 555   }
 556   image_view=DestroyCacheView(image_view);
 557   shift_view=DestroyCacheView(shift_view);
 558   if (status == MagickFalse)
 559     shift_image=DestroyImage(shift_image);
 560   return(shift_image);
 561 }
 562 \f
 563 /*
 564 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 565 %                                                                             %
 566 %                                                                             %
 567 %                                                                             %
 568 %     C h a r c o a l I m a g e                                               %
 569 %                                                                             %
 570 %                                                                             %
 571 %                                                                             %
 572 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 573 %
 574 %  CharcoalImage() creates a new image that is a copy of an existing one with
 575 %  the edge highlighted.  It allocates the memory necessary for the new Image
 576 %  structure and returns a pointer to the new image.
 577 %
 578 %  The format of the CharcoalImage method is:
 579 %
 580 %      Image *CharcoalImage(const Image *image,const double radius,
 581 %        const double sigma,ExceptionInfo *exception)
 582 %
 583 %  A description of each parameter follows:
 584 %
 585 %    o image: the image.
 586 %
 587 %    o radius: the radius of the pixel neighborhood.
 588 %
 589 %    o sigma: the standard deviation of the Gaussian, in pixels.
 590 %
 591 %    o exception: return any errors or warnings in this structure.
 592 %
 593 */
 594 MagickExport Image *CharcoalImage(const Image *image,const double radius,
 595   const double sigma,ExceptionInfo *exception)
 596 {
 597   Image
 598     *charcoal_image,
 599     *clone_image,
 600     *edge_image;
 601
 602   assert(image != (Image *) NULL);
 603   assert(image->signature == MagickCoreSignature);
 604   if (image->debug != MagickFalse)
 605     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 606   assert(exception != (ExceptionInfo *) NULL);
 607   assert(exception->signature == MagickCoreSignature);
 608   clone_image=CloneImage(image,0,0,MagickTrue,exception);
 609   if (clone_image == (Image *) NULL)
 610     return((Image *) NULL);
 611   edge_image=EdgeImage(clone_image,radius,exception);
 612   clone_image=DestroyImage(clone_image);
 613   if (edge_image == (Image *) NULL)
 614     return((Image *) NULL);
 615   charcoal_image=BlurImage(edge_image,radius,sigma,exception);
 616   edge_image=DestroyImage(edge_image);
 617   if (charcoal_image == (Image *) NULL)
 618     return((Image *) NULL);
 619   (void) NormalizeImage(charcoal_image,exception);
 620   (void) NegateImage(charcoal_image,MagickFalse,exception);
 621   (void) GrayscaleImage(charcoal_image,image->intensity,exception);
 622   return(charcoal_image);
 623 }
 624 \f
 625 /*
 626 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 627 %                                                                             %
 628 %                                                                             %
 629 %                                                                             %
 630 %     C o l o r i z e I m a g e                                               %
 631 %                                                                             %
 632 %                                                                             %
 633 %                                                                             %
 634 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 635 %
 636 %  ColorizeImage() blends the fill color with each pixel in the image.
 637 %  A percentage blend is specified with opacity.  Control the application
 638 %  of different color components by specifying a different percentage for
 639 %  each component (e.g. 90/100/10 is 90% red, 100% green, and 10% blue).
 640 %
 641 %  The format of the ColorizeImage method is:
 642 %
 643 %      Image *ColorizeImage(const Image *image,const char *blend,
 644 %        const PixelInfo *colorize,ExceptionInfo *exception)
 645 %
 646 %  A description of each parameter follows:
 647 %
 648 %    o image: the image.
 649 %
 650 %    o blend:  A character string indicating the level of blending as a
 651 %      percentage.
 652 %
 653 %    o colorize: A color value.
 654 %
 655 %    o exception: return any errors or warnings in this structure.
 656 %
 657 */
 658 MagickExport Image *ColorizeImage(const Image *image,const char *blend,
 659   const PixelInfo *colorize,ExceptionInfo *exception)
 660 {
 661 #define ColorizeImageTag  "Colorize/Image"
 662 #define Colorize(pixel,blend_percentage,colorize)  \
 663   (((pixel)*(100.0-(blend_percentage))+(colorize)*(blend_percentage))/100.0)
 664
 665   CacheView
 666     *image_view;
 667
 668   GeometryInfo
 669     geometry_info;
 670
 671   Image
 672     *colorize_image;
 673
 674   MagickBooleanType
 675     status;
 676
 677   MagickOffsetType
 678     progress;
 679
 680   MagickStatusType
 681     flags;
 682
 683   PixelInfo
 684     blend_percentage;
 685
 686   ssize_t
 687     y;
 688
 689   /*
 690     Allocate colorized image.
 691   */
 692   assert(image != (const Image *) NULL);
 693   assert(image->signature == MagickCoreSignature);
 694   if (image->debug != MagickFalse)
 695     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 696   assert(exception != (ExceptionInfo *) NULL);
 697   assert(exception->signature == MagickCoreSignature);
 698   colorize_image=CloneImage(image,0,0,MagickTrue,exception);
 699   if (colorize_image == (Image *) NULL)
 700     return((Image *) NULL);
 701   if (SetImageStorageClass(colorize_image,DirectClass,exception) == MagickFalse)
 702     {
 703       colorize_image=DestroyImage(colorize_image);
 704       return((Image *) NULL);
 705     }
 706   if ((IsGrayColorspace(colorize_image->colorspace) != MagickFalse) ||
 707       (IsPixelInfoGray(colorize) != MagickFalse))
 708     (void) SetImageColorspace(colorize_image,sRGBColorspace,exception);
 709   if ((colorize_image->alpha_trait == UndefinedPixelTrait) &&
 710       (colorize->alpha_trait != UndefinedPixelTrait))
 711     (void) SetImageAlpha(colorize_image,OpaqueAlpha,exception);
 712   if (blend == (const char *) NULL)
 713     return(colorize_image);
 714   GetPixelInfo(colorize_image,&blend_percentage);
 715   flags=ParseGeometry(blend,&geometry_info);
 716   blend_percentage.red=geometry_info.rho;
 717   blend_percentage.green=geometry_info.rho;
 718   blend_percentage.blue=geometry_info.rho;
 719   blend_percentage.black=geometry_info.rho;
 720   blend_percentage.alpha=(MagickRealType) TransparentAlpha;
 721   if ((flags & SigmaValue) != 0)
 722     blend_percentage.green=geometry_info.sigma;
 723   if ((flags & XiValue) != 0)
 724     blend_percentage.blue=geometry_info.xi;
 725   if ((flags & PsiValue) != 0)
 726     blend_percentage.alpha=geometry_info.psi;
 727   if (blend_percentage.colorspace == CMYKColorspace)
 728     {
 729       if ((flags & PsiValue) != 0)
 730         blend_percentage.black=geometry_info.psi;
 731       if ((flags & ChiValue) != 0)
 732         blend_percentage.alpha=geometry_info.chi;
 733     }
 734   /*
 735     Colorize DirectClass image.
 736   */
 737   status=MagickTrue;
 738   progress=0;
 739   image_view=AcquireVirtualCacheView(colorize_image,exception);
 740 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 741   #pragma omp parallel for schedule(static,4) shared(progress,status) \
 742     magick_threads(colorize_image,colorize_image,colorize_image->rows,1)
 743 #endif
 744   for (y=0; y < (ssize_t) colorize_image->rows; y++)
 745   {
 746     MagickBooleanType
 747       sync;
 748
 749     register Quantum
 750       *magick_restrict q;
 751
 752     register ssize_t
 753       x;
 754
 755     if (status == MagickFalse)
 756       continue;
 757     q=GetCacheViewAuthenticPixels(image_view,0,y,colorize_image->columns,1,
 758       exception);
 759     if (q == (Quantum *) NULL)
 760       {
 761         status=MagickFalse;
 762         continue;
 763       }
 764     for (x=0; x < (ssize_t) colorize_image->columns; x++)
 765     {
 766       register ssize_t
 767         i;
 768
 769       for (i=0; i < (ssize_t) GetPixelChannels(colorize_image); i++)
 770       {
 771         PixelTrait traits=GetPixelChannelTraits(colorize_image,
 772           (PixelChannel) i);
 773         if (traits == UndefinedPixelTrait)
 774           continue;
 775         if (((traits & CopyPixelTrait) != 0) ||
 776             (GetPixelWriteMask(colorize_image,q) == 0))
 777           continue;
 778         SetPixelChannel(colorize_image,(PixelChannel) i,ClampToQuantum(
 779           Colorize(q[i],GetPixelInfoChannel(&blend_percentage,(PixelChannel) i),
 780           GetPixelInfoChannel(colorize,(PixelChannel) i))),q);
 781       }
 782       q+=GetPixelChannels(colorize_image);
 783     }
 784     sync=SyncCacheViewAuthenticPixels(image_view,exception);
 785     if (sync == MagickFalse)
 786       status=MagickFalse;
 787     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 788       {
 789         MagickBooleanType
 790           proceed;
 791
 792 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 793         #pragma omp critical (MagickCore_ColorizeImage)
 794 #endif
 795         proceed=SetImageProgress(image,ColorizeImageTag,progress++,
 796           colorize_image->rows);
 797         if (proceed == MagickFalse)
 798           status=MagickFalse;
 799       }
 800   }
 801   image_view=DestroyCacheView(image_view);
 802   if (status == MagickFalse)
 803     colorize_image=DestroyImage(colorize_image);
 804   return(colorize_image);
 805 }
 806 \f
 807 /*
 808 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 809 %                                                                             %
 810 %                                                                             %
 811 %                                                                             %
 812 %     C o l o r M a t r i x I m a g e                                         %
 813 %                                                                             %
 814 %                                                                             %
 815 %                                                                             %
 816 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 817 %
 818 %  ColorMatrixImage() applies color transformation to an image. This method
 819 %  permits saturation changes, hue rotation, luminance to alpha, and various
 820 %  other effects.  Although variable-sized transformation matrices can be used,
 821 %  typically one uses a 5x5 matrix for an RGBA image and a 6x6 for CMYKA
 822 %  (or RGBA with offsets).  The matrix is similar to those used by Adobe Flash
 823 %  except offsets are in column 6 rather than 5 (in support of CMYKA images)
 824 %  and offsets are normalized (divide Flash offset by 255).
 825 %
 826 %  The format of the ColorMatrixImage method is:
 827 %
 828 %      Image *ColorMatrixImage(const Image *image,
 829 %        const KernelInfo *color_matrix,ExceptionInfo *exception)
 830 %
 831 %  A description of each parameter follows:
 832 %
 833 %    o image: the image.
 834 %
 835 %    o color_matrix:  the color matrix.
 836 %
 837 %    o exception: return any errors or warnings in this structure.
 838 %
 839 */
 840 /* FUTURE: modify to make use of a MagickMatrix Mutliply function
 841    That should be provided in "matrix.c"
 842    (ASIDE: actually distorts should do this too but currently doesn't)
 843 */
 844
 845 MagickExport Image *ColorMatrixImage(const Image *image,
 846   const KernelInfo *color_matrix,ExceptionInfo *exception)
 847 {
 848 #define ColorMatrixImageTag  "ColorMatrix/Image"
 849
 850   CacheView
 851     *color_view,
 852     *image_view;
 853
 854   double
 855     ColorMatrix[6][6] =
 856     {
 857       { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
 858       { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
 859       { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
 860       { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
 861       { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
 862       { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 }
 863     };
 864
 865   Image
 866     *color_image;
 867
 868   MagickBooleanType
 869     status;
 870
 871   MagickOffsetType
 872     progress;
 873
 874   register ssize_t
 875     i;
 876
 877   ssize_t
 878     u,
 879     v,
 880     y;
 881
 882   /*
 883     Map given color_matrix, into a 6x6 matrix   RGBKA and a constant
 884   */
 885   assert(image != (Image *) NULL);
 886   assert(image->signature == MagickCoreSignature);
 887   if (image->debug != MagickFalse)
 888     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 889   assert(exception != (ExceptionInfo *) NULL);
 890   assert(exception->signature == MagickCoreSignature);
 891   i=0;
 892   for (v=0; v < (ssize_t) color_matrix->height; v++)
 893     for (u=0; u < (ssize_t) color_matrix->width; u++)
 894     {
 895       if ((v < 6) && (u < 6))
 896         ColorMatrix[v][u]=color_matrix->values[i];
 897       i++;
 898     }
 899   /*
 900     Initialize color image.
 901   */
 902   color_image=CloneImage(image,0,0,MagickTrue,exception);
 903   if (color_image == (Image *) NULL)
 904     return((Image *) NULL);
 905   if (SetImageStorageClass(color_image,DirectClass,exception) == MagickFalse)
 906     {
 907       color_image=DestroyImage(color_image);
 908       return((Image *) NULL);
 909     }
 910   if (image->debug != MagickFalse)
 911     {
 912       char
 913         format[MagickPathExtent],
 914         *message;
 915
 916       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
 917         "  ColorMatrix image with color matrix:");
 918       message=AcquireString("");
 919       for (v=0; v < 6; v++)
 920       {
 921         *message='\0';
 922         (void) FormatLocaleString(format,MagickPathExtent,"%.20g: ",(double) v);
 923         (void) ConcatenateString(&message,format);
 924         for (u=0; u < 6; u++)
 925         {
 926           (void) FormatLocaleString(format,MagickPathExtent,"%+f ",
 927             ColorMatrix[v][u]);
 928           (void) ConcatenateString(&message,format);
 929         }
 930         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
 931       }
 932       message=DestroyString(message);
 933     }
 934   /*
 935     Apply the ColorMatrix to image.
 936   */
 937   status=MagickTrue;
 938   progress=0;
 939   image_view=AcquireVirtualCacheView(image,exception);
 940   color_view=AcquireAuthenticCacheView(color_image,exception);
 941 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 942   #pragma omp parallel for schedule(static,4) shared(progress,status) \
 943     magick_threads(image,color_image,image->rows,1)
 944 #endif
 945   for (y=0; y < (ssize_t) image->rows; y++)
 946   {
 947     PixelInfo
 948       pixel;
 949
 950     register const Quantum
 951       *magick_restrict p;
 952
 953     register Quantum
 954       *magick_restrict q;
 955
 956     register ssize_t
 957       x;
 958
 959     if (status == MagickFalse)
 960       continue;
 961     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 962     q=GetCacheViewAuthenticPixels(color_view,0,y,color_image->columns,1,
 963       exception);
 964     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
 965       {
 966         status=MagickFalse;
 967         continue;
 968       }
 969     GetPixelInfo(image,&pixel);
 970     for (x=0; x < (ssize_t) image->columns; x++)
 971     {
 972       register ssize_t
 973         v;
 974
 975       size_t
 976         height;
 977
 978       GetPixelInfoPixel(image,p,&pixel);
 979       height=color_matrix->height > 6 ? 6UL : color_matrix->height;
 980       for (v=0; v < (ssize_t) height; v++)
 981       {
 982         double
 983           sum;
 984
 985         sum=ColorMatrix[v][0]*GetPixelRed(image,p)+ColorMatrix[v][1]*
 986           GetPixelGreen(image,p)+ColorMatrix[v][2]*GetPixelBlue(image,p);
 987         if (image->colorspace == CMYKColorspace)
 988           sum+=ColorMatrix[v][3]*GetPixelBlack(image,p);
 989         if (image->alpha_trait != UndefinedPixelTrait)
 990           sum+=ColorMatrix[v][4]*GetPixelAlpha(image,p);
 991         sum+=QuantumRange*ColorMatrix[v][5];
 992         switch (v)
 993         {
 994           case 0: pixel.red=sum; break;
 995           case 1: pixel.green=sum; break;
 996           case 2: pixel.blue=sum; break;
 997           case 3: pixel.black=sum; break;
 998           case 4: pixel.alpha=sum; break;
 999           default: break;
1000         }
1001       }
1002       SetPixelViaPixelInfo(color_image,&pixel,q);
1003       p+=GetPixelChannels(image);
1004       q+=GetPixelChannels(color_image);
1005     }
1006     if (SyncCacheViewAuthenticPixels(color_view,exception) == MagickFalse)
1007       status=MagickFalse;
1008     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1009       {
1010         MagickBooleanType
1011           proceed;
1012
1013 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1014         #pragma omp critical (MagickCore_ColorMatrixImage)
1015 #endif
1016         proceed=SetImageProgress(image,ColorMatrixImageTag,progress++,
1017           image->rows);
1018         if (proceed == MagickFalse)
1019           status=MagickFalse;
1020       }
1021   }
1022   color_view=DestroyCacheView(color_view);
1023   image_view=DestroyCacheView(image_view);
1024   if (status == MagickFalse)
1025     color_image=DestroyImage(color_image);
1026   return(color_image);
1027 }
1028 \f
1029 /*
1030 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1031 %                                                                             %
1032 %                                                                             %
1033 %                                                                             %
1034 +   D e s t r o y F x I n f o                                                 %
1035 %                                                                             %
1036 %                                                                             %
1037 %                                                                             %
1038 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1039 %
1040 %  DestroyFxInfo() deallocates memory associated with an FxInfo structure.
1041 %
1042 %  The format of the DestroyFxInfo method is:
1043 %
1044 %      ImageInfo *DestroyFxInfo(ImageInfo *fx_info)
1045 %
1046 %  A description of each parameter follows:
1047 %
1048 %    o fx_info: the fx info.
1049 %
1050 */
1051 MagickPrivate FxInfo *DestroyFxInfo(FxInfo *fx_info)
1052 {
1053   register ssize_t
1054     i;
1055
1056   fx_info->exception=DestroyExceptionInfo(fx_info->exception);
1057   fx_info->expression=DestroyString(fx_info->expression);
1058   fx_info->symbols=DestroySplayTree(fx_info->symbols);
1059   fx_info->colors=DestroySplayTree(fx_info->colors);
1060   for (i=(ssize_t) GetImageListLength(fx_info->images)-1; i >= 0; i--)
1061     fx_info->view[i]=DestroyCacheView(fx_info->view[i]);
1062   fx_info->view=(CacheView **) RelinquishMagickMemory(fx_info->view);
1063   fx_info->random_info=DestroyRandomInfo(fx_info->random_info);
1064   fx_info=(FxInfo *) RelinquishMagickMemory(fx_info);
1065   return(fx_info);
1066 }
1067 \f
1068 /*
1069 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1070 %                                                                             %
1071 %                                                                             %
1072 %                                                                             %
1073 +     F x E v a l u a t e C h a n n e l E x p r e s s i o n                   %
1074 %                                                                             %
1075 %                                                                             %
1076 %                                                                             %
1077 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1078 %
1079 %  FxEvaluateChannelExpression() evaluates an expression and returns the
1080 %  results.
1081 %
1082 %  The format of the FxEvaluateExpression method is:
1083 %
1084 %      double FxEvaluateChannelExpression(FxInfo *fx_info,
1085 %        const PixelChannel channel,const ssize_t x,const ssize_t y,
1086 %        double *alpha,Exceptioninfo *exception)
1087 %      double FxEvaluateExpression(FxInfo *fx_info,
1088 %        double *alpha,Exceptioninfo *exception)
1089 %
1090 %  A description of each parameter follows:
1091 %
1092 %    o fx_info: the fx info.
1093 %
1094 %    o channel: the channel.
1095 %
1096 %    o x,y: the pixel position.
1097 %
1098 %    o alpha: the result.
1099 %
1100 %    o exception: return any errors or warnings in this structure.
1101 %
1102 */
1103
1104 static double FxChannelStatistics(FxInfo *fx_info,Image *image,
1105   PixelChannel channel,const char *symbol,ExceptionInfo *exception)
1106 {
1107   ChannelType
1108     channel_mask;
1109
1110   char
1111     key[MagickPathExtent],
1112     statistic[MagickPathExtent];
1113
1114   const char
1115     *value;
1116
1117   register const char
1118     *p;
1119
1120   channel_mask=UndefinedChannel;
1121   for (p=symbol; (*p != '.') && (*p != '\0'); p++) ;
1122   if (*p == '.')
1123     {
1124       ssize_t
1125         option;
1126
1127       option=ParseCommandOption(MagickPixelChannelOptions,MagickTrue,p+1);
1128       if (option >= 0)
1129         {
1130           channel=(PixelChannel) option;
1131           channel_mask=SetPixelChannelMask(image,(ChannelType) (1 << channel));
1132         }
1133     }
1134   (void) FormatLocaleString(key,MagickPathExtent,"%p.%.20g.%s",(void *) image,
1135     (double) channel,symbol);
1136   value=(const char *) GetValueFromSplayTree(fx_info->symbols,key);
1137   if (value != (const char *) NULL)
1138     {
1139       if (channel_mask != UndefinedChannel)
1140         (void) SetPixelChannelMask(image,channel_mask);
1141       return(QuantumScale*StringToDouble(value,(char **) NULL));
1142     }
1143   (void) DeleteNodeFromSplayTree(fx_info->symbols,key);
1144   if (LocaleNCompare(symbol,"depth",5) == 0)
1145     {
1146       size_t
1147         depth;
1148
1149       depth=GetImageDepth(image,exception);
1150       (void) FormatLocaleString(statistic,MagickPathExtent,"%.20g",(double)
1151         depth);
1152     }
1153   if (LocaleNCompare(symbol,"kurtosis",8) == 0)
1154     {
1155       double
1156         kurtosis,
1157         skewness;
1158
1159       (void) GetImageKurtosis(image,&kurtosis,&skewness,exception);
1160       (void) FormatLocaleString(statistic,MagickPathExtent,"%g",kurtosis);
1161     }
1162   if (LocaleNCompare(symbol,"maxima",6) == 0)
1163     {
1164       double
1165         maxima,
1166         minima;
1167
1168       (void) GetImageRange(image,&minima,&maxima,exception);
1169       (void) FormatLocaleString(statistic,MagickPathExtent,"%g",maxima);
1170     }
1171   if (LocaleNCompare(symbol,"mean",4) == 0)
1172     {
1173       double
1174         mean,
1175         standard_deviation;
1176
1177       (void) GetImageMean(image,&mean,&standard_deviation,exception);
1178       (void) FormatLocaleString(statistic,MagickPathExtent,"%g",mean);
1179     }
1180   if (LocaleNCompare(symbol,"minima",6) == 0)
1181     {
1182       double
1183         maxima,
1184         minima;
1185
1186       (void) GetImageRange(image,&minima,&maxima,exception);
1187       (void) FormatLocaleString(statistic,MagickPathExtent,"%g",minima);
1188     }
1189   if (LocaleNCompare(symbol,"skewness",8) == 0)
1190     {
1191       double
1192         kurtosis,
1193         skewness;
1194
1195       (void) GetImageKurtosis(image,&kurtosis,&skewness,exception);
1196       (void) FormatLocaleString(statistic,MagickPathExtent,"%g",skewness);
1197     }
1198   if (LocaleNCompare(symbol,"standard_deviation",18) == 0)
1199     {
1200       double
1201         mean,
1202         standard_deviation;
1203
1204       (void) GetImageMean(image,&mean,&standard_deviation,exception);
1205       (void) FormatLocaleString(statistic,MagickPathExtent,"%g",
1206         standard_deviation);
1207     }
1208   if (channel_mask != UndefinedChannel)
1209     (void) SetPixelChannelMask(image,channel_mask);
1210   (void) AddValueToSplayTree(fx_info->symbols,ConstantString(key),
1211     ConstantString(statistic));
1212   return(QuantumScale*StringToDouble(statistic,(char **) NULL));
1213 }
1214
1215 static double
1216   FxEvaluateSubexpression(FxInfo *,const PixelChannel,const ssize_t,
1217     const ssize_t,const char *,size_t *,double *,ExceptionInfo *);
1218
1219 static MagickOffsetType FxGCD(MagickOffsetType alpha,MagickOffsetType beta)
1220 {
1221   if (beta != 0)
1222     return(FxGCD(beta,alpha % beta));
1223   return(alpha);
1224 }
1225
1226 static inline const char *FxSubexpression(const char *expression,
1227   ExceptionInfo *exception)
1228 {
1229   const char
1230     *subexpression;
1231
1232   register ssize_t
1233     level;
1234
1235   level=0;
1236   subexpression=expression;
1237   while ((*subexpression != '\0') &&
1238          ((level != 1) || (strchr(")",(int) *subexpression) == (char *) NULL)))
1239   {
1240     if (strchr("(",(int) *subexpression) != (char *) NULL)
1241       level++;
1242     else
1243       if (strchr(")",(int) *subexpression) != (char *) NULL)
1244         level--;
1245     subexpression++;
1246   }
1247   if (*subexpression == '\0')
1248     (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1249       "UnbalancedParenthesis","`%s'",expression);
1250   return(subexpression);
1251 }
1252
1253 static double FxGetSymbol(FxInfo *fx_info,const PixelChannel channel,
1254   const ssize_t x,const ssize_t y,const char *expression,
1255   ExceptionInfo *exception)
1256 {
1257   char
1258     *q,
1259     subexpression[MagickPathExtent],
1260     symbol[MagickPathExtent];
1261
1262   const char
1263     *p,
1264     *value;
1265
1266   Image
1267     *image;
1268
1269   PixelInfo
1270     pixel;
1271
1272   double
1273     alpha,
1274     beta;
1275
1276   PointInfo
1277     point;
1278
1279   register ssize_t
1280     i;
1281
1282   size_t
1283     depth,
1284     length,
1285     level;
1286
1287   p=expression;
1288   i=GetImageIndexInList(fx_info->images);
1289   depth=0;
1290   level=0;
1291   point.x=(double) x;
1292   point.y=(double) y;
1293   if (isalpha((int) ((unsigned char) *(p+1))) == 0)
1294     {
1295       if (strchr("suv",(int) *p) != (char *) NULL)
1296         {
1297           switch (*p)
1298           {
1299             case 's':
1300             default:
1301             {
1302               i=GetImageIndexInList(fx_info->images);
1303               break;
1304             }
1305             case 'u': i=0; break;
1306             case 'v': i=1; break;
1307           }
1308           p++;
1309           if (*p == '[')
1310             {
1311               level++;
1312               q=subexpression;
1313               for (p++; *p != '\0'; )
1314               {
1315                 if (*p == '[')
1316                   level++;
1317                 else
1318                   if (*p == ']')
1319                     {
1320                       level--;
1321                       if (level == 0)
1322                         break;
1323                     }
1324                 *q++=(*p++);
1325               }
1326               *q='\0';
1327               alpha=FxEvaluateSubexpression(fx_info,channel,x,y,subexpression,
1328                 &depth,&beta,exception);
1329               i=(ssize_t) alpha;
1330               p++;
1331             }
1332           if (*p == '.')
1333             p++;
1334         }
1335       if ((*p == 'p') && (isalpha((int) ((unsigned char) *(p+1))) == 0))
1336         {
1337           p++;
1338           if (*p == '{')
1339             {
1340               level++;
1341               q=subexpression;
1342               for (p++; *p != '\0'; )
1343               {
1344                 if (*p == '{')
1345                   level++;
1346                 else
1347                   if (*p == '}')
1348                     {
1349                       level--;
1350                       if (level == 0)
1351                         break;
1352                     }
1353                 *q++=(*p++);
1354               }
1355               *q='\0';
1356               alpha=FxEvaluateSubexpression(fx_info,channel,x,y,subexpression,
1357                 &depth,&beta,exception);
1358               point.x=alpha;
1359               point.y=beta;
1360               p++;
1361             }
1362           else
1363             if (*p == '[')
1364               {
1365                 level++;
1366                 q=subexpression;
1367                 for (p++; *p != '\0'; )
1368                 {
1369                   if (*p == '[')
1370                     level++;
1371                   else
1372                     if (*p == ']')
1373                       {
1374                         level--;
1375                         if (level == 0)
1376                           break;
1377                       }
1378                   *q++=(*p++);
1379                 }
1380                 *q='\0';
1381                 alpha=FxEvaluateSubexpression(fx_info,channel,x,y,subexpression,
1382                   &depth,&beta,exception);
1383                 point.x+=alpha;
1384                 point.y+=beta;
1385                 p++;
1386               }
1387           if (*p == '.')
1388             p++;
1389         }
1390     }
1391   length=GetImageListLength(fx_info->images);
1392   while (i < 0)
1393     i+=(ssize_t) length;
1394   if (length != 0)
1395     i%=length;
1396   image=GetImageFromList(fx_info->images,i);
1397   if (image == (Image *) NULL)
1398     {
1399       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1400         "NoSuchImage","`%s'",expression);
1401       return(0.0);
1402     }
1403   GetPixelInfo(image,&pixel);
1404   (void) InterpolatePixelInfo(image,fx_info->view[i],image->interpolate,
1405     point.x,point.y,&pixel,exception);
1406   if ((strlen(p) > 2) && (LocaleCompare(p,"intensity") != 0) &&
1407       (LocaleCompare(p,"luma") != 0) && (LocaleCompare(p,"luminance") != 0) &&
1408       (LocaleCompare(p,"hue") != 0) && (LocaleCompare(p,"saturation") != 0) &&
1409       (LocaleCompare(p,"lightness") != 0))
1410     {
1411       char
1412         name[MagickPathExtent];
1413
1414       (void) CopyMagickString(name,p,MagickPathExtent);
1415       for (q=name+(strlen(name)-1); q > name; q--)
1416       {
1417         if (*q == ')')
1418           break;
1419         if (*q == '.')
1420           {
1421             *q='\0';
1422             break;
1423           }
1424       }
1425       if ((strlen(name) > 2) &&
1426           (GetValueFromSplayTree(fx_info->symbols,name) == (const char *) NULL))
1427         {
1428           PixelInfo
1429             *color;
1430
1431           color=(PixelInfo *) GetValueFromSplayTree(fx_info->colors,name);
1432           if (color != (PixelInfo *) NULL)
1433             {
1434               pixel=(*color);
1435               p+=strlen(name);
1436             }
1437           else
1438             {
1439               MagickBooleanType
1440                 status;
1441
1442               status=QueryColorCompliance(name,AllCompliance,&pixel,
1443                 fx_info->exception);
1444               if (status != MagickFalse)
1445                 {
1446                   (void) AddValueToSplayTree(fx_info->colors,ConstantString(
1447                     name),ClonePixelInfo(&pixel));
1448                   p+=strlen(name);
1449                 }
1450             }
1451         }
1452     }
1453   (void) CopyMagickString(symbol,p,MagickPathExtent);
1454   StripString(symbol);
1455   if (*symbol == '\0')
1456     {
1457       switch (channel)
1458       {
1459         case RedPixelChannel: return(QuantumScale*pixel.red);
1460         case GreenPixelChannel: return(QuantumScale*pixel.green);
1461         case BluePixelChannel: return(QuantumScale*pixel.blue);
1462         case BlackPixelChannel:
1463         {
1464           if (image->colorspace != CMYKColorspace)
1465             {
1466               (void) ThrowMagickException(exception,GetMagickModule(),
1467                 ImageError,"ColorSeparatedImageRequired","`%s'",
1468                 image->filename);
1469               return(0.0);
1470             }
1471           return(QuantumScale*pixel.black);
1472         }
1473         case AlphaPixelChannel:
1474         {
1475           if (pixel.alpha_trait == UndefinedPixelTrait)
1476             return(1.0);
1477           alpha=(double) (QuantumScale*pixel.alpha);
1478           return(alpha);
1479         }
1480         case IndexPixelChannel:
1481           return(0.0);
1482         case IntensityPixelChannel:
1483         {
1484           Quantum
1485             quantum_pixel[MaxPixelChannels];
1486
1487           SetPixelViaPixelInfo(image,&pixel,quantum_pixel);
1488           return(QuantumScale*GetPixelIntensity(image,quantum_pixel));
1489         }
1490         default:
1491           break;
1492       }
1493       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1494         "UnableToParseExpression","`%s'",p);
1495       return(0.0);
1496     }
1497   switch (*symbol)
1498   {
1499     case 'A':
1500     case 'a':
1501     {
1502       if (LocaleCompare(symbol,"a") == 0)
1503         return((QuantumScale*pixel.alpha));
1504       break;
1505     }
1506     case 'B':
1507     case 'b':
1508     {
1509       if (LocaleCompare(symbol,"b") == 0)
1510         return(QuantumScale*pixel.blue);
1511       break;
1512     }
1513     case 'C':
1514     case 'c':
1515     {
1516       if (LocaleNCompare(symbol,"channel",7) == 0)
1517         {
1518           GeometryInfo
1519             channel_info;
1520
1521           MagickStatusType
1522             flags;
1523
1524           flags=ParseGeometry(symbol+7,&channel_info);
1525           if (image->colorspace == CMYKColorspace)
1526             switch (channel)
1527             {
1528               case CyanPixelChannel:
1529               {
1530                 if ((flags & RhoValue) == 0)
1531                   return(0.0);
1532                 return(channel_info.rho);
1533               }
1534               case MagentaPixelChannel:
1535               {
1536                 if ((flags & SigmaValue) == 0)
1537                   return(0.0);
1538                 return(channel_info.sigma);
1539               }
1540               case YellowPixelChannel:
1541               {
1542                 if ((flags & XiValue) == 0)
1543                   return(0.0);
1544                 return(channel_info.xi);
1545               }
1546               case BlackPixelChannel:
1547               {
1548                 if ((flags & PsiValue) == 0)
1549                   return(0.0);
1550                 return(channel_info.psi);
1551               }
1552               case AlphaPixelChannel:
1553               {
1554                 if ((flags & ChiValue) == 0)
1555                   return(0.0);
1556                 return(channel_info.chi);
1557               }
1558               default:
1559                 return(0.0);
1560             }
1561           switch (channel)
1562           {
1563             case RedPixelChannel:
1564             {
1565               if ((flags & RhoValue) == 0)
1566                 return(0.0);
1567               return(channel_info.rho);
1568             }
1569             case GreenPixelChannel:
1570             {
1571               if ((flags & SigmaValue) == 0)
1572                 return(0.0);
1573               return(channel_info.sigma);
1574             }
1575             case BluePixelChannel:
1576             {
1577               if ((flags & XiValue) == 0)
1578                 return(0.0);
1579               return(channel_info.xi);
1580             }
1581             case BlackPixelChannel:
1582             {
1583               if ((flags & ChiValue) == 0)
1584                 return(0.0);
1585               return(channel_info.chi);
1586             }
1587             case AlphaPixelChannel:
1588             {
1589               if ((flags & PsiValue) == 0)
1590                 return(0.0);
1591               return(channel_info.psi);
1592             }
1593             default:
1594               return(0.0);
1595           }
1596         }
1597       if (LocaleCompare(symbol,"c") == 0)
1598         return(QuantumScale*pixel.red);
1599       break;
1600     }
1601     case 'D':
1602     case 'd':
1603     {
1604       if (LocaleNCompare(symbol,"depth",5) == 0)
1605         return(FxChannelStatistics(fx_info,image,channel,symbol,exception));
1606       break;
1607     }
1608     case 'G':
1609     case 'g':
1610     {
1611       if (LocaleCompare(symbol,"g") == 0)
1612         return(QuantumScale*pixel.green);
1613       break;
1614     }
1615     case 'K':
1616     case 'k':
1617     {
1618       if (LocaleNCompare(symbol,"kurtosis",8) == 0)
1619         return(FxChannelStatistics(fx_info,image,channel,symbol,exception));
1620       if (LocaleCompare(symbol,"k") == 0)
1621         {
1622           if (image->colorspace != CMYKColorspace)
1623             {
1624               (void) ThrowMagickException(exception,GetMagickModule(),
1625                 OptionError,"ColorSeparatedImageRequired","`%s'",
1626                 image->filename);
1627               return(0.0);
1628             }
1629           return(QuantumScale*pixel.black);
1630         }
1631       break;
1632     }
1633     case 'H':
1634     case 'h':
1635     {
1636       if (LocaleCompare(symbol,"h") == 0)
1637         return((double) image->rows);
1638       if (LocaleCompare(symbol,"hue") == 0)
1639         {
1640           double
1641             hue,
1642             lightness,
1643             saturation;
1644
1645           ConvertRGBToHSL(pixel.red,pixel.green,pixel.blue,&hue,&saturation,
1646             &lightness);
1647           return(hue);
1648         }
1649       break;
1650     }
1651     case 'I':
1652     case 'i':
1653     {
1654       if ((LocaleCompare(symbol,"image.depth") == 0) ||
1655           (LocaleCompare(symbol,"image.minima") == 0) ||
1656           (LocaleCompare(symbol,"image.maxima") == 0) ||
1657           (LocaleCompare(symbol,"image.mean") == 0) ||
1658           (LocaleCompare(symbol,"image.kurtosis") == 0) ||
1659           (LocaleCompare(symbol,"image.skewness") == 0) ||
1660           (LocaleCompare(symbol,"image.standard_deviation") == 0))
1661         return(FxChannelStatistics(fx_info,image,channel,symbol+6,exception));
1662       if (LocaleCompare(symbol,"image.resolution.x") == 0)
1663         return(image->resolution.x);
1664       if (LocaleCompare(symbol,"image.resolution.y") == 0)
1665         return(image->resolution.y);
1666       if (LocaleCompare(symbol,"intensity") == 0)
1667         {
1668           Quantum
1669             quantum_pixel[MaxPixelChannels];
1670
1671           SetPixelViaPixelInfo(image,&pixel,quantum_pixel);
1672           return(QuantumScale*GetPixelIntensity(image,quantum_pixel));
1673         }
1674       if (LocaleCompare(symbol,"i") == 0)
1675         return((double) x);
1676       break;
1677     }
1678     case 'J':
1679     case 'j':
1680     {
1681       if (LocaleCompare(symbol,"j") == 0)
1682         return((double) y);
1683       break;
1684     }
1685     case 'L':
1686     case 'l':
1687     {
1688       if (LocaleCompare(symbol,"lightness") == 0)
1689         {
1690           double
1691             hue,
1692             lightness,
1693             saturation;
1694
1695           ConvertRGBToHSL(pixel.red,pixel.green,pixel.blue,&hue,&saturation,
1696             &lightness);
1697           return(lightness);
1698         }
1699       if (LocaleCompare(symbol,"luma") == 0)
1700         {
1701           double
1702             luma;
1703
1704           luma=0.212656*pixel.red+0.715158*pixel.green+0.072186*pixel.blue;
1705           return(QuantumScale*luma);
1706         }
1707       if (LocaleCompare(symbol,"luminance") == 0)
1708         {
1709           double
1710             luminence;
1711
1712           luminence=0.212656*pixel.red+0.715158*pixel.green+0.072186*pixel.blue;
1713           return(QuantumScale*luminence);
1714         }
1715       break;
1716     }
1717     case 'M':
1718     case 'm':
1719     {
1720       if (LocaleNCompare(symbol,"maxima",6) == 0)
1721         return(FxChannelStatistics(fx_info,image,channel,symbol,exception));
1722       if (LocaleNCompare(symbol,"mean",4) == 0)
1723         return(FxChannelStatistics(fx_info,image,channel,symbol,exception));
1724       if (LocaleNCompare(symbol,"minima",6) == 0)
1725         return(FxChannelStatistics(fx_info,image,channel,symbol,exception));
1726       if (LocaleCompare(symbol,"m") == 0)
1727         return(QuantumScale*pixel.green);
1728       break;
1729     }
1730     case 'N':
1731     case 'n':
1732     {
1733       if (LocaleCompare(symbol,"n") == 0)
1734         return((double) GetImageListLength(fx_info->images));
1735       break;
1736     }
1737     case 'O':
1738     case 'o':
1739     {
1740       if (LocaleCompare(symbol,"o") == 0)
1741         return(QuantumScale*pixel.alpha);
1742       break;
1743     }
1744     case 'P':
1745     case 'p':
1746     {
1747       if (LocaleCompare(symbol,"page.height") == 0)
1748         return((double) image->page.height);
1749       if (LocaleCompare(symbol,"page.width") == 0)
1750         return((double) image->page.width);
1751       if (LocaleCompare(symbol,"page.x") == 0)
1752         return((double) image->page.x);
1753       if (LocaleCompare(symbol,"page.y") == 0)
1754         return((double) image->page.y);
1755       break;
1756     }
1757     case 'Q':
1758     case 'q':
1759     {
1760       if (LocaleCompare(symbol,"quality") == 0)
1761         return((double) image->quality);
1762       break;
1763     }
1764     case 'R':
1765     case 'r':
1766     {
1767       if (LocaleCompare(symbol,"resolution.x") == 0)
1768         return(image->resolution.x);
1769       if (LocaleCompare(symbol,"resolution.y") == 0)
1770         return(image->resolution.y);
1771       if (LocaleCompare(symbol,"r") == 0)
1772         return(QuantumScale*pixel.red);
1773       break;
1774     }
1775     case 'S':
1776     case 's':
1777     {
1778       if (LocaleCompare(symbol,"saturation") == 0)
1779         {
1780           double
1781             hue,
1782             lightness,
1783             saturation;
1784
1785           ConvertRGBToHSL(pixel.red,pixel.green,pixel.blue,&hue,&saturation,
1786             &lightness);
1787           return(saturation);
1788         }
1789       if (LocaleNCompare(symbol,"skewness",8) == 0)
1790         return(FxChannelStatistics(fx_info,image,channel,symbol,exception));
1791       if (LocaleNCompare(symbol,"standard_deviation",18) == 0)
1792         return(FxChannelStatistics(fx_info,image,channel,symbol,exception));
1793       break;
1794     }
1795     case 'T':
1796     case 't':
1797     {
1798       if (LocaleCompare(symbol,"t") == 0)
1799         return((double) GetImageIndexInList(fx_info->images));
1800       break;
1801     }
1802     case 'W':
1803     case 'w':
1804     {
1805       if (LocaleCompare(symbol,"w") == 0)
1806         return((double) image->columns);
1807       break;
1808     }
1809     case 'Y':
1810     case 'y':
1811     {
1812       if (LocaleCompare(symbol,"y") == 0)
1813         return(QuantumScale*pixel.blue);
1814       break;
1815     }
1816     case 'Z':
1817     case 'z':
1818     {
1819       if (LocaleCompare(symbol,"z") == 0)
1820         return((double)GetImageDepth(image, fx_info->exception));
1821       break;
1822     }
1823     default:
1824       break;
1825   }
1826   value=(const char *) GetValueFromSplayTree(fx_info->symbols,symbol);
1827   if (value != (const char *) NULL)
1828     return(StringToDouble(value,(char **) NULL));
1829   (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1830     "UnableToParseExpression","`%s'",symbol);
1831   return(0.0);
1832 }
1833
1834 static const char *FxOperatorPrecedence(const char *expression,
1835   ExceptionInfo *exception)
1836 {
1837   typedef enum
1838   {
1839     UndefinedPrecedence,
1840     NullPrecedence,
1841     BitwiseComplementPrecedence,
1842     ExponentPrecedence,
1843     ExponentialNotationPrecedence,
1844     MultiplyPrecedence,
1845     AdditionPrecedence,
1846     ShiftPrecedence,
1847     RelationalPrecedence,
1848     EquivalencyPrecedence,
1849     BitwiseAndPrecedence,
1850     BitwiseOrPrecedence,
1851     LogicalAndPrecedence,
1852     LogicalOrPrecedence,
1853     TernaryPrecedence,
1854     AssignmentPrecedence,
1855     CommaPrecedence,
1856     SeparatorPrecedence
1857   } FxPrecedence;
1858
1859   FxPrecedence
1860     precedence,
1861     target;
1862
1863   register const char
1864     *subexpression;
1865
1866   register int
1867     c;
1868
1869   size_t
1870     level;
1871
1872   c=0;
1873   level=0;
1874   subexpression=(const char *) NULL;
1875   target=NullPrecedence;
1876   while (*expression != '\0')
1877   {
1878     precedence=UndefinedPrecedence;
1879     if ((isspace((int) ((unsigned char) *expression)) != 0) || (c == (int) '@'))
1880       {
1881         expression++;
1882         continue;
1883       }
1884     switch (*expression)
1885     {
1886       case 'A':
1887       case 'a':
1888       {
1889 #if defined(MAGICKCORE_HAVE_ACOSH)
1890         if (LocaleNCompare(expression,"acosh",5) == 0)
1891           {
1892             expression+=5;
1893             break;
1894           }
1895 #endif
1896 #if defined(MAGICKCORE_HAVE_ASINH)
1897         if (LocaleNCompare(expression,"asinh",5) == 0)
1898           {
1899             expression+=5;
1900             break;
1901           }
1902 #endif
1903 #if defined(MAGICKCORE_HAVE_ATANH)
1904         if (LocaleNCompare(expression,"atanh",5) == 0)
1905           {
1906             expression+=5;
1907             break;
1908           }
1909 #endif
1910         if (LocaleNCompare(expression,"atan2",5) == 0)
1911           {
1912             expression+=5;
1913             break;
1914           }
1915         break;
1916       }
1917       case 'E':
1918       case 'e':
1919       {
1920         if ((isdigit((int) ((unsigned char) c)) != 0) &&
1921             ((LocaleNCompare(expression,"E+",2) == 0) ||
1922              (LocaleNCompare(expression,"E-",2) == 0)))
1923           {
1924             expression+=2;  /* scientific notation */
1925             break;
1926           }
1927       }
1928       case 'J':
1929       case 'j':
1930       {
1931         if ((LocaleNCompare(expression,"j0",2) == 0) ||
1932             (LocaleNCompare(expression,"j1",2) == 0))
1933           {
1934             expression+=2;
1935             break;
1936           }
1937         break;
1938       }
1939       case '#':
1940       {
1941         while (isxdigit((int) ((unsigned char) *(expression+1))) != 0)
1942           expression++;
1943         break;
1944       }
1945       default:
1946         break;
1947     }
1948     if ((c == (int) '{') || (c == (int) '['))
1949       level++;
1950     else
1951       if ((c == (int) '}') || (c == (int) ']'))
1952         level--;
1953     if (level == 0)
1954       switch ((unsigned char) *expression)
1955       {
1956         case '~':
1957         case '!':
1958         {
1959           precedence=BitwiseComplementPrecedence;
1960           break;
1961         }
1962         case '^':
1963         case '@':
1964         {
1965           precedence=ExponentPrecedence;
1966           break;
1967         }
1968         default:
1969         {
1970           if (((c != 0) && ((isdigit((int) ((unsigned char) c)) != 0) ||
1971                (strchr(")",(int) ((unsigned char) c)) != (char *) NULL))) &&
1972               (((islower((int) ((unsigned char) *expression)) != 0) ||
1973                (strchr("(",(int) ((unsigned char) *expression)) != (char *) NULL)) ||
1974                ((isdigit((int) ((unsigned char) c)) == 0) &&
1975                 (isdigit((int) ((unsigned char) *expression)) != 0))) &&
1976               (strchr("xy",(int) ((unsigned char) *expression)) == (char *) NULL))
1977             precedence=MultiplyPrecedence;
1978           break;
1979         }
1980         case '*':
1981         case '/':
1982         case '%':
1983         {
1984           precedence=MultiplyPrecedence;
1985           break;
1986         }
1987         case '+':
1988         case '-':
1989         {
1990           if ((strchr("(+-/*%:&^|<>~,",c) == (char *) NULL) ||
1991               (isalpha(c) != 0))
1992             precedence=AdditionPrecedence;
1993           break;
1994         }
1995         case LeftShiftOperator:
1996         case RightShiftOperator:
1997         {
1998           precedence=ShiftPrecedence;
1999           break;
2000         }
2001         case '<':
2002         case LessThanEqualOperator:
2003         case GreaterThanEqualOperator:
2004         case '>':
2005         {
2006           precedence=RelationalPrecedence;
2007           break;
2008         }
2009         case EqualOperator:
2010         case NotEqualOperator:
2011         {
2012           precedence=EquivalencyPrecedence;
2013           break;
2014         }
2015         case '&':
2016         {
2017           precedence=BitwiseAndPrecedence;
2018           break;
2019         }
2020         case '|':
2021         {
2022           precedence=BitwiseOrPrecedence;
2023           break;
2024         }
2025         case LogicalAndOperator:
2026         {
2027           precedence=LogicalAndPrecedence;
2028           break;
2029         }
2030         case LogicalOrOperator:
2031         {
2032           precedence=LogicalOrPrecedence;
2033           break;
2034         }
2035         case ExponentialNotation:
2036         {
2037           precedence=ExponentialNotationPrecedence;
2038           break;
2039         }
2040         case ':':
2041         case '?':
2042         {
2043           precedence=TernaryPrecedence;
2044           break;
2045         }
2046         case '=':
2047         {
2048           precedence=AssignmentPrecedence;
2049           break;
2050         }
2051         case ',':
2052         {
2053           precedence=CommaPrecedence;
2054           break;
2055         }
2056         case ';':
2057         {
2058           precedence=SeparatorPrecedence;
2059           break;
2060         }
2061       }
2062     if ((precedence == BitwiseComplementPrecedence) ||
2063         (precedence == TernaryPrecedence) ||
2064         (precedence == AssignmentPrecedence))
2065       {
2066         if (precedence > target)
2067           {
2068             /*
2069               Right-to-left associativity.
2070             */
2071             target=precedence;
2072             subexpression=expression;
2073           }
2074       }
2075     else
2076       if (precedence >= target)
2077         {
2078           /*
2079             Left-to-right associativity.
2080           */
2081           target=precedence;
2082           subexpression=expression;
2083         }
2084     if (strchr("(",(int) *expression) != (char *) NULL)
2085       expression=FxSubexpression(expression,exception);
2086     c=(int) (*expression++);
2087   }
2088   return(subexpression);
2089 }
2090
2091 static double FxEvaluateSubexpression(FxInfo *fx_info,
2092   const PixelChannel channel,const ssize_t x,const ssize_t y,
2093   const char *expression,size_t *depth,double *beta,ExceptionInfo *exception)
2094 {
2095 #define FxMaxParenthesisDepth  58
2096
2097   char
2098     *q,
2099     subexpression[MagickPathExtent];
2100
2101   double
2102     alpha,
2103     gamma;
2104
2105   register const char
2106     *p;
2107
2108   *beta=0.0;
2109   if (exception->severity >= ErrorException)
2110     return(0.0);
2111   while (isspace((int) ((unsigned char) *expression)) != 0)
2112     expression++;
2113   if (*expression == '\0')
2114     return(0.0);
2115   *subexpression='\0';
2116   p=FxOperatorPrecedence(expression,exception);
2117   if (p != (const char *) NULL)
2118     {
2119       (void) CopyMagickString(subexpression,expression,(size_t)
2120         (p-expression+1));
2121       alpha=FxEvaluateSubexpression(fx_info,channel,x,y,subexpression,depth,
2122         beta,exception);
2123       switch ((unsigned char) *p)
2124       {
2125         case '~':
2126         {
2127           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2128             exception);
2129           *beta=(double) (~(size_t) *beta);
2130           return(*beta);
2131         }
2132         case '!':
2133         {
2134           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2135             exception);
2136           return(*beta == 0.0 ? 1.0 : 0.0);
2137         }
2138         case '^':
2139         {
2140           *beta=pow(alpha,FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,
2141             beta,exception));
2142           return(*beta);
2143         }
2144         case '*':
2145         case ExponentialNotation:
2146         {
2147           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2148             exception);
2149           return(alpha*(*beta));
2150         }
2151         case '/':
2152         {
2153           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2154             exception);
2155           if (*beta == 0.0)
2156             {
2157               (void) ThrowMagickException(exception,GetMagickModule(),
2158                 OptionError,"DivideByZero","`%s'",expression);
2159               return(0.0);
2160             }
2161           return(alpha/(*beta));
2162         }
2163         case '%':
2164         {
2165           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2166             exception);
2167           *beta=fabs(floor((*beta)+0.5));
2168           if (*beta == 0.0)
2169             {
2170               (void) ThrowMagickException(exception,GetMagickModule(),
2171                 OptionError,"DivideByZero","`%s'",expression);
2172               return(0.0);
2173             }
2174           return(fmod(alpha,*beta));
2175         }
2176         case '+':
2177         {
2178           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2179             exception);
2180           return(alpha+(*beta));
2181         }
2182         case '-':
2183         {
2184           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2185             exception);
2186           return(alpha-(*beta));
2187         }
2188         case LeftShiftOperator:
2189         {
2190           gamma=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2191             exception);
2192           *beta=(double) ((size_t) (alpha+0.5) << (size_t) (gamma+0.5));
2193           return(*beta);
2194         }
2195         case RightShiftOperator:
2196         {
2197           gamma=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2198             exception);
2199           *beta=(double) ((size_t) (alpha+0.5) >> (size_t) (gamma+0.5));
2200           return(*beta);
2201         }
2202         case '<':
2203         {
2204           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2205             exception);
2206           return(alpha < *beta ? 1.0 : 0.0);
2207         }
2208         case LessThanEqualOperator:
2209         {
2210           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2211             exception);
2212           return(alpha <= *beta ? 1.0 : 0.0);
2213         }
2214         case '>':
2215         {
2216           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2217             exception);
2218           return(alpha > *beta ? 1.0 : 0.0);
2219         }
2220         case GreaterThanEqualOperator:
2221         {
2222           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2223             exception);
2224           return(alpha >= *beta ? 1.0 : 0.0);
2225         }
2226         case EqualOperator:
2227         {
2228           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2229             exception);
2230           return(fabs(alpha-(*beta)) < MagickEpsilon ? 1.0 : 0.0);
2231         }
2232         case NotEqualOperator:
2233         {
2234           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2235             exception);
2236           return(fabs(alpha-(*beta)) >= MagickEpsilon ? 1.0 : 0.0);
2237         }
2238         case '&':
2239         {
2240           gamma=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2241             exception);
2242           *beta=(double) ((size_t) (alpha+0.5) & (size_t) (gamma+0.5));
2243           return(*beta);
2244         }
2245         case '|':
2246         {
2247           gamma=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2248             exception);
2249           *beta=(double) ((size_t) (alpha+0.5) | (size_t) (gamma+0.5));
2250           return(*beta);
2251         }
2252         case LogicalAndOperator:
2253         {
2254           p++;
2255           if (alpha <= 0.0)
2256             {
2257               *beta=0.0;
2258               return(*beta);
2259             }
2260           gamma=FxEvaluateSubexpression(fx_info,channel,x,y,p,depth,beta,
2261             exception);
2262           *beta=(gamma > 0.0) ? 1.0 : 0.0;
2263           return(*beta);
2264         }
2265         case LogicalOrOperator:
2266         {
2267           p++;
2268           if (alpha > 0.0)
2269             {
2270              *beta=1.0;
2271              return(*beta);
2272             }
2273           gamma=FxEvaluateSubexpression(fx_info,channel,x,y,p,depth,beta,
2274             exception);
2275           *beta=(gamma > 0.0) ? 1.0 : 0.0;
2276           return(*beta);
2277         }
2278         case '?':
2279         {
2280           (void) CopyMagickString(subexpression,++p,MagickPathExtent);
2281           q=subexpression;
2282           p=StringToken(":",&q);
2283           if (q == (char *) NULL)
2284             {
2285               (void) ThrowMagickException(exception,GetMagickModule(),
2286                 OptionError,"UnableToParseExpression","`%s'",subexpression);
2287               return(0.0);
2288             }
2289           if (fabs(alpha) >= MagickEpsilon)
2290             gamma=FxEvaluateSubexpression(fx_info,channel,x,y,p,depth,beta,
2291               exception);
2292           else
2293             gamma=FxEvaluateSubexpression(fx_info,channel,x,y,q,depth,beta,
2294               exception);
2295           return(gamma);
2296         }
2297         case '=':
2298         {
2299           char
2300             numeric[MagickPathExtent];
2301
2302           q=subexpression;
2303           while (isalpha((int) ((unsigned char) *q)) != 0)
2304             q++;
2305           if (*q != '\0')
2306             {
2307               (void) ThrowMagickException(exception,GetMagickModule(),
2308                 OptionError,"UnableToParseExpression","`%s'",subexpression);
2309               return(0.0);
2310             }
2311           ClearMagickException(exception);
2312           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2313             exception);
2314           (void) FormatLocaleString(numeric,MagickPathExtent,"%g",*beta);
2315           (void) DeleteNodeFromSplayTree(fx_info->symbols,subexpression);
2316           (void) AddValueToSplayTree(fx_info->symbols,ConstantString(
2317             subexpression),ConstantString(numeric));
2318           return(*beta);
2319         }
2320         case ',':
2321         {
2322           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2323             exception);
2324           return(alpha);
2325         }
2326         case ';':
2327         {
2328           *beta=FxEvaluateSubexpression(fx_info,channel,x,y,++p,depth,beta,
2329             exception);
2330           return(*beta);
2331         }
2332         default:
2333         {
2334           gamma=alpha*FxEvaluateSubexpression(fx_info,channel,x,y,p,depth,beta,
2335             exception);
2336           return(gamma);
2337         }
2338       }
2339     }
2340   if (strchr("(",(int) *expression) != (char *) NULL)
2341     {
2342       (*depth)++;
2343       if (*depth >= FxMaxParenthesisDepth)
2344         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
2345           "ParenthesisNestedTooDeeply","`%s'",expression);
2346       (void) CopyMagickString(subexpression,expression+1,MagickPathExtent);
2347       subexpression[strlen(subexpression)-1]='\0';
2348       gamma=FxEvaluateSubexpression(fx_info,channel,x,y,subexpression,depth,
2349         beta,exception);
2350       (*depth)--;
2351       return(gamma);
2352     }
2353   switch (*expression)
2354   {
2355     case '+':
2356     {
2357       gamma=FxEvaluateSubexpression(fx_info,channel,x,y,expression+1,depth,beta,
2358         exception);
2359       return(1.0*gamma);
2360     }
2361     case '-':
2362     {
2363       gamma=FxEvaluateSubexpression(fx_info,channel,x,y,expression+1,depth,beta,
2364         exception);
2365       return(-1.0*gamma);
2366     }
2367     case '~':
2368     {
2369       gamma=FxEvaluateSubexpression(fx_info,channel,x,y,expression+1,depth,beta,
2370         exception);
2371       return((double) (~(size_t) (gamma+0.5)));
2372     }
2373     case 'A':
2374     case 'a':
2375     {
2376       if (LocaleNCompare(expression,"abs",3) == 0)
2377         {
2378           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2379             beta,exception);
2380           return(fabs(alpha));
2381         }
2382 #if defined(MAGICKCORE_HAVE_ACOSH)
2383       if (LocaleNCompare(expression,"acosh",5) == 0)
2384         {
2385           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2386             beta,exception);
2387           return(acosh(alpha));
2388         }
2389 #endif
2390       if (LocaleNCompare(expression,"acos",4) == 0)
2391         {
2392           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2393             beta,exception);
2394           return(acos(alpha));
2395         }
2396 #if defined(MAGICKCORE_HAVE_J1)
2397       if (LocaleNCompare(expression,"airy",4) == 0)
2398         {
2399           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2400             beta,exception);
2401           if (alpha == 0.0)
2402             return(1.0);
2403           gamma=2.0*j1((MagickPI*alpha))/(MagickPI*alpha);
2404           return(gamma*gamma);
2405         }
2406 #endif
2407 #if defined(MAGICKCORE_HAVE_ASINH)
2408       if (LocaleNCompare(expression,"asinh",5) == 0)
2409         {
2410           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2411             beta,exception);
2412           return(asinh(alpha));
2413         }
2414 #endif
2415       if (LocaleNCompare(expression,"asin",4) == 0)
2416         {
2417           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2418             beta,exception);
2419           return(asin(alpha));
2420         }
2421       if (LocaleNCompare(expression,"alt",3) == 0)
2422         {
2423           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2424             beta,exception);
2425           return(((ssize_t) alpha) & 0x01 ? -1.0 : 1.0);
2426         }
2427       if (LocaleNCompare(expression,"atan2",5) == 0)
2428         {
2429           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2430             beta,exception);
2431           return(atan2(alpha,*beta));
2432         }
2433 #if defined(MAGICKCORE_HAVE_ATANH)
2434       if (LocaleNCompare(expression,"atanh",5) == 0)
2435         {
2436           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2437             beta,exception);
2438           return(atanh(alpha));
2439         }
2440 #endif
2441       if (LocaleNCompare(expression,"atan",4) == 0)
2442         {
2443           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2444             beta,exception);
2445           return(atan(alpha));
2446         }
2447       if (LocaleCompare(expression,"a") == 0)
2448         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2449       break;
2450     }
2451     case 'B':
2452     case 'b':
2453     {
2454       if (LocaleCompare(expression,"b") == 0)
2455         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2456       break;
2457     }
2458     case 'C':
2459     case 'c':
2460     {
2461       if (LocaleNCompare(expression,"ceil",4) == 0)
2462         {
2463           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2464             beta,exception);
2465           return(ceil(alpha));
2466         }
2467       if (LocaleNCompare(expression,"clamp",5) == 0)
2468         {
2469           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2470             beta,exception);
2471           if (alpha < 0.0)
2472             return(0.0);
2473           if (alpha > 1.0)
2474             return(1.0);
2475           return(alpha);
2476         }
2477       if (LocaleNCompare(expression,"cosh",4) == 0)
2478         {
2479           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2480             beta,exception);
2481           return(cosh(alpha));
2482         }
2483       if (LocaleNCompare(expression,"cos",3) == 0)
2484         {
2485           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2486             beta,exception);
2487           return(cos(alpha));
2488         }
2489       if (LocaleCompare(expression,"c") == 0)
2490         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2491       break;
2492     }
2493     case 'D':
2494     case 'd':
2495     {
2496       if (LocaleNCompare(expression,"debug",5) == 0)
2497         {
2498           const char
2499             *type;
2500
2501           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2502             beta,exception);
2503           if (fx_info->images->colorspace == CMYKColorspace)
2504             switch (channel)
2505             {
2506               case CyanPixelChannel: type="cyan"; break;
2507               case MagentaPixelChannel: type="magenta"; break;
2508               case YellowPixelChannel: type="yellow"; break;
2509               case AlphaPixelChannel: type="opacity"; break;
2510               case BlackPixelChannel: type="black"; break;
2511               default: type="unknown"; break;
2512             }
2513           else
2514             switch (channel)
2515             {
2516               case RedPixelChannel: type="red"; break;
2517               case GreenPixelChannel: type="green"; break;
2518               case BluePixelChannel: type="blue"; break;
2519               case AlphaPixelChannel: type="opacity"; break;
2520               default: type="unknown"; break;
2521             }
2522           (void) CopyMagickString(subexpression,expression+6,MagickPathExtent);
2523           if (strlen(subexpression) > 1)
2524             subexpression[strlen(subexpression)-1]='\0';
2525           if (fx_info->file != (FILE *) NULL)
2526             (void) FormatLocaleFile(fx_info->file,"%s[%.20g,%.20g].%s: "
2527                "%s=%.*g\n",fx_info->images->filename,(double) x,(double) y,type,
2528                subexpression,GetMagickPrecision(),alpha);
2529           return(0.0);
2530         }
2531       if (LocaleNCompare(expression,"drc",3) == 0)
2532         {
2533           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2534             beta,exception);
2535           return((alpha/(*beta*(alpha-1.0)+1.0)));
2536         }
2537       break;
2538     }
2539     case 'E':
2540     case 'e':
2541     {
2542       if (LocaleCompare(expression,"epsilon") == 0)
2543         return(MagickEpsilon);
2544 #if defined(MAGICKCORE_HAVE_ERF)
2545       if (LocaleNCompare(expression,"erf",3) == 0)
2546         {
2547           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2548             beta,exception);
2549           return(erf(alpha));
2550         }
2551 #endif
2552       if (LocaleNCompare(expression,"exp",3) == 0)
2553         {
2554           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2555             beta,exception);
2556           return(exp(alpha));
2557         }
2558       if (LocaleCompare(expression,"e") == 0)
2559         return(2.7182818284590452354);
2560       break;
2561     }
2562     case 'F':
2563     case 'f':
2564     {
2565       if (LocaleNCompare(expression,"floor",5) == 0)
2566         {
2567           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2568             beta,exception);
2569           return(floor(alpha));
2570         }
2571       break;
2572     }
2573     case 'G':
2574     case 'g':
2575     {
2576       if (LocaleNCompare(expression,"gauss",5) == 0)
2577         {
2578           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2579             beta,exception);
2580           gamma=exp((-alpha*alpha/2.0))/sqrt(2.0*MagickPI);
2581           return(gamma);
2582         }
2583       if (LocaleNCompare(expression,"gcd",3) == 0)
2584         {
2585           MagickOffsetType
2586             gcd;
2587
2588           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2589             beta,exception);
2590           gcd=FxGCD((MagickOffsetType) (alpha+0.5),(MagickOffsetType) (*beta+
2591             0.5));
2592           return((double) gcd);
2593         }
2594       if (LocaleCompare(expression,"g") == 0)
2595         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2596       break;
2597     }
2598     case 'H':
2599     case 'h':
2600     {
2601       if (LocaleCompare(expression,"h") == 0)
2602         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2603       if (LocaleCompare(expression,"hue") == 0)
2604         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2605       if (LocaleNCompare(expression,"hypot",5) == 0)
2606         {
2607           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2608             beta,exception);
2609           return(hypot(alpha,*beta));
2610         }
2611       break;
2612     }
2613     case 'K':
2614     case 'k':
2615     {
2616       if (LocaleCompare(expression,"k") == 0)
2617         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2618       break;
2619     }
2620     case 'I':
2621     case 'i':
2622     {
2623       if (LocaleCompare(expression,"intensity") == 0)
2624         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2625       if (LocaleNCompare(expression,"int",3) == 0)
2626         {
2627           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2628             beta,exception);
2629           return(floor(alpha));
2630         }
2631       if (LocaleNCompare(expression,"isnan",5) == 0)
2632         {
2633           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2634             beta,exception);
2635           return((double) !!IsNaN(alpha));
2636         }
2637       if (LocaleCompare(expression,"i") == 0)
2638         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2639       break;
2640     }
2641     case 'J':
2642     case 'j':
2643     {
2644       if (LocaleCompare(expression,"j") == 0)
2645         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2646 #if defined(MAGICKCORE_HAVE_J0)
2647       if (LocaleNCompare(expression,"j0",2) == 0)
2648         {
2649           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+2,depth,
2650             beta,exception);
2651           return(j0(alpha));
2652         }
2653 #endif
2654 #if defined(MAGICKCORE_HAVE_J1)
2655       if (LocaleNCompare(expression,"j1",2) == 0)
2656         {
2657           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+2,depth,
2658             beta,exception);
2659           return(j1(alpha));
2660         }
2661 #endif
2662 #if defined(MAGICKCORE_HAVE_J1)
2663       if (LocaleNCompare(expression,"jinc",4) == 0)
2664         {
2665           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2666             beta,exception);
2667           if (alpha == 0.0)
2668             return(1.0);
2669           gamma=(2.0*j1((MagickPI*alpha))/(MagickPI*alpha));
2670           return(gamma);
2671         }
2672 #endif
2673       break;
2674     }
2675     case 'L':
2676     case 'l':
2677     {
2678       if (LocaleNCompare(expression,"ln",2) == 0)
2679         {
2680           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+2,depth,
2681             beta,exception);
2682           return(log(alpha));
2683         }
2684       if (LocaleNCompare(expression,"logtwo",6) == 0)
2685         {
2686           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+6,depth,
2687             beta,exception);
2688           return(log10(alpha))/log10(2.0);
2689         }
2690       if (LocaleNCompare(expression,"log",3) == 0)
2691         {
2692           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2693             beta,exception);
2694           return(log10(alpha));
2695         }
2696       if (LocaleCompare(expression,"lightness") == 0)
2697         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2698       break;
2699     }
2700     case 'M':
2701     case 'm':
2702     {
2703       if (LocaleCompare(expression,"MaxRGB") == 0)
2704         return(QuantumRange);
2705       if (LocaleNCompare(expression,"maxima",6) == 0)
2706         break;
2707       if (LocaleNCompare(expression,"max",3) == 0)
2708         {
2709           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2710             beta,exception);
2711           return(alpha > *beta ? alpha : *beta);
2712         }
2713       if (LocaleNCompare(expression,"minima",6) == 0)
2714         break;
2715       if (LocaleNCompare(expression,"min",3) == 0)
2716         {
2717           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2718             beta,exception);
2719           return(alpha < *beta ? alpha : *beta);
2720         }
2721       if (LocaleNCompare(expression,"mod",3) == 0)
2722         {
2723           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2724             beta,exception);
2725           gamma=alpha-floor((alpha/(*beta)))*(*beta);
2726           return(gamma);
2727         }
2728       if (LocaleCompare(expression,"m") == 0)
2729         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2730       break;
2731     }
2732     case 'N':
2733     case 'n':
2734     {
2735       if (LocaleNCompare(expression,"not",3) == 0)
2736         {
2737           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2738             beta,exception);
2739           return((double) (alpha < MagickEpsilon));
2740         }
2741       if (LocaleCompare(expression,"n") == 0)
2742         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2743       break;
2744     }
2745     case 'O':
2746     case 'o':
2747     {
2748       if (LocaleCompare(expression,"Opaque") == 0)
2749         return(1.0);
2750       if (LocaleCompare(expression,"o") == 0)
2751         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2752       break;
2753     }
2754     case 'P':
2755     case 'p':
2756     {
2757       if (LocaleCompare(expression,"phi") == 0)
2758         return(MagickPHI);
2759       if (LocaleCompare(expression,"pi") == 0)
2760         return(MagickPI);
2761       if (LocaleNCompare(expression,"pow",3) == 0)
2762         {
2763           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2764             beta,exception);
2765           return(pow(alpha,*beta));
2766         }
2767       if (LocaleCompare(expression,"p") == 0)
2768         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2769       break;
2770     }
2771     case 'Q':
2772     case 'q':
2773     {
2774       if (LocaleCompare(expression,"QuantumRange") == 0)
2775         return(QuantumRange);
2776       if (LocaleCompare(expression,"QuantumScale") == 0)
2777         return(QuantumScale);
2778       break;
2779     }
2780     case 'R':
2781     case 'r':
2782     {
2783       if (LocaleNCompare(expression,"rand",4) == 0)
2784         {
2785 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2786         #pragma omp critical (MagickCore_FxEvaluateSubexpression)
2787 #endif
2788           alpha=GetPseudoRandomValue(fx_info->random_info);
2789           return(alpha);
2790         }
2791       if (LocaleNCompare(expression,"round",5) == 0)
2792         {
2793           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2794             beta,exception);
2795           return(floor(alpha+0.5));
2796         }
2797       if (LocaleCompare(expression,"r") == 0)
2798         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2799       break;
2800     }
2801     case 'S':
2802     case 's':
2803     {
2804       if (LocaleCompare(expression,"saturation") == 0)
2805         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2806       if (LocaleNCompare(expression,"sign",4) == 0)
2807         {
2808           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2809             beta,exception);
2810           return(alpha < 0.0 ? -1.0 : 1.0);
2811         }
2812       if (LocaleNCompare(expression,"sinc",4) == 0)
2813         {
2814           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2815             beta,exception);
2816           if (alpha == 0)
2817             return(1.0);
2818           gamma=sin((MagickPI*alpha))/(MagickPI*alpha);
2819           return(gamma);
2820         }
2821       if (LocaleNCompare(expression,"sinh",4) == 0)
2822         {
2823           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2824             beta,exception);
2825           return(sinh(alpha));
2826         }
2827       if (LocaleNCompare(expression,"sin",3) == 0)
2828         {
2829           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2830             beta,exception);
2831           return(sin(alpha));
2832         }
2833       if (LocaleNCompare(expression,"sqrt",4) == 0)
2834         {
2835           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2836             beta,exception);
2837           return(sqrt(alpha));
2838         }
2839       if (LocaleNCompare(expression,"squish",6) == 0)
2840         {
2841           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+6,depth,
2842             beta,exception);
2843           return((1.0/(1.0+exp(-alpha))));
2844         }
2845       if (LocaleCompare(expression,"s") == 0)
2846         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2847       break;
2848     }
2849     case 'T':
2850     case 't':
2851     {
2852       if (LocaleNCompare(expression,"tanh",4) == 0)
2853         {
2854           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,depth,
2855             beta,exception);
2856           return(tanh(alpha));
2857         }
2858       if (LocaleNCompare(expression,"tan",3) == 0)
2859         {
2860           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,depth,
2861             beta,exception);
2862           return(tan(alpha));
2863         }
2864       if (LocaleCompare(expression,"Transparent") == 0)
2865         return(0.0);
2866       if (LocaleNCompare(expression,"trunc",5) == 0)
2867         {
2868           alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,depth,
2869             beta,exception);
2870           if (alpha >= 0.0)
2871             return(floor(alpha));
2872           return(ceil(alpha));
2873         }
2874       if (LocaleCompare(expression,"t") == 0)
2875         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2876       break;
2877     }
2878     case 'U':
2879     case 'u':
2880     {
2881       if (LocaleCompare(expression,"u") == 0)
2882         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2883       break;
2884     }
2885     case 'V':
2886     case 'v':
2887     {
2888       if (LocaleCompare(expression,"v") == 0)
2889         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2890       break;
2891     }
2892     case 'W':
2893     case 'w':
2894     {
2895       if (LocaleNCompare(expression,"while",5) == 0)
2896         {
2897           do
2898           {
2899             alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+5,
2900               depth,beta,exception);
2901           } while (fabs(alpha) >= MagickEpsilon);
2902           return(*beta);
2903         }
2904       if (LocaleCompare(expression,"w") == 0)
2905         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2906       break;
2907     }
2908     case 'Y':
2909     case 'y':
2910     {
2911       if (LocaleCompare(expression,"y") == 0)
2912         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2913       break;
2914     }
2915     case 'Z':
2916     case 'z':
2917     {
2918       if (LocaleCompare(expression,"z") == 0)
2919         return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2920       break;
2921     }
2922     default:
2923       break;
2924   }
2925   q=(char *) expression;
2926   alpha=InterpretSiPrefixValue(expression,&q);
2927   if (q == expression)
2928     return(FxGetSymbol(fx_info,channel,x,y,expression,exception));
2929   return(alpha);
2930 }
2931
2932 MagickPrivate MagickBooleanType FxEvaluateExpression(FxInfo *fx_info,
2933   double *alpha,ExceptionInfo *exception)
2934 {
2935   MagickBooleanType
2936     status;
2937
2938   status=FxEvaluateChannelExpression(fx_info,GrayPixelChannel,0,0,alpha,
2939     exception);
2940   return(status);
2941 }
2942
2943 MagickExport MagickBooleanType FxPreprocessExpression(FxInfo *fx_info,
2944   double *alpha,ExceptionInfo *exception)
2945 {
2946   FILE
2947     *file;
2948
2949   MagickBooleanType
2950     status;
2951
2952   file=fx_info->file;
2953   fx_info->file=(FILE *) NULL;
2954   status=FxEvaluateChannelExpression(fx_info,GrayPixelChannel,0,0,alpha,
2955     exception);
2956   fx_info->file=file;
2957   return(status);
2958 }
2959
2960 MagickPrivate MagickBooleanType FxEvaluateChannelExpression(FxInfo *fx_info,
2961   const PixelChannel channel,const ssize_t x,const ssize_t y,
2962   double *alpha,ExceptionInfo *exception)
2963 {
2964   double
2965     beta;
2966
2967   size_t
2968     depth;
2969
2970   depth=0;
2971   beta=0.0;
2972   *alpha=FxEvaluateSubexpression(fx_info,channel,x,y,fx_info->expression,&depth,
2973     &beta,exception);
2974   return(exception->severity == OptionError ? MagickFalse : MagickTrue);
2975 }
2976 \f
2977 /*
2978 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2979 %                                                                             %
2980 %                                                                             %
2981 %                                                                             %
2982 %     F x I m a g e                                                           %
2983 %                                                                             %
2984 %                                                                             %
2985 %                                                                             %
2986 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2987 %
2988 %  FxImage() applies a mathematical expression to the specified image.
2989 %
2990 %  The format of the FxImage method is:
2991 %
2992 %      Image *FxImage(const Image *image,const char *expression,
2993 %        ExceptionInfo *exception)
2994 %
2995 %  A description of each parameter follows:
2996 %
2997 %    o image: the image.
2998 %
2999 %    o expression: A mathematical expression.
3000 %
3001 %    o exception: return any errors or warnings in this structure.
3002 %
3003 */
3004
3005 static FxInfo **DestroyFxThreadSet(FxInfo **fx_info)
3006 {
3007   register ssize_t
3008     i;
3009
3010   assert(fx_info != (FxInfo **) NULL);
3011   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
3012     if (fx_info[i] != (FxInfo *) NULL)
3013       fx_info[i]=DestroyFxInfo(fx_info[i]);
3014   fx_info=(FxInfo **) RelinquishMagickMemory(fx_info);
3015   return(fx_info);
3016 }
3017
3018 static FxInfo **AcquireFxThreadSet(const Image *image,const char *expression,
3019   ExceptionInfo *exception)
3020 {
3021   char
3022     *fx_expression;
3023
3024   FxInfo
3025     **fx_info;
3026
3027   double
3028     alpha;
3029
3030   register ssize_t
3031     i;
3032
3033   size_t
3034     number_threads;
3035
3036   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
3037   fx_info=(FxInfo **) AcquireQuantumMemory(number_threads,sizeof(*fx_info));
3038   if (fx_info == (FxInfo **) NULL)
3039     {
3040       (void) ThrowMagickException(exception,GetMagickModule(),
3041         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
3042       return((FxInfo **) NULL);
3043     }
3044   (void) ResetMagickMemory(fx_info,0,number_threads*sizeof(*fx_info));
3045   if (*expression != '@')
3046     fx_expression=ConstantString(expression);
3047   else
3048     fx_expression=FileToString(expression+1,~0UL,exception);
3049   for (i=0; i < (ssize_t) number_threads; i++)
3050   {
3051     MagickBooleanType
3052       status;
3053
3054     fx_info[i]=AcquireFxInfo(image,fx_expression,exception);
3055     if (fx_info[i] == (FxInfo *) NULL)
3056       break;
3057     status=FxPreprocessExpression(fx_info[i],&alpha,exception);
3058     if (status == MagickFalse)
3059       break;
3060   }
3061   fx_expression=DestroyString(fx_expression);
3062   if (i < (ssize_t) number_threads)
3063     fx_info=DestroyFxThreadSet(fx_info);
3064   return(fx_info);
3065 }
3066
3067 MagickExport Image *FxImage(const Image *image,const char *expression,
3068   ExceptionInfo *exception)
3069 {
3070 #define FxImageTag  "Fx/Image"
3071
3072   CacheView
3073     *fx_view,
3074     *image_view;
3075
3076   FxInfo
3077     **magick_restrict fx_info;
3078
3079   Image
3080     *fx_image;
3081
3082   MagickBooleanType
3083     status;
3084
3085   MagickOffsetType
3086     progress;
3087
3088   ssize_t
3089     y;
3090
3091   assert(image != (Image *) NULL);
3092   assert(image->signature == MagickCoreSignature);
3093   if (image->debug != MagickFalse)
3094     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3095   fx_info=AcquireFxThreadSet(image,expression,exception);
3096   if (fx_info == (FxInfo **) NULL)
3097     return((Image *) NULL);
3098   fx_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
3099   if (fx_image == (Image *) NULL)
3100     {
3101       fx_info=DestroyFxThreadSet(fx_info);
3102       return((Image *) NULL);
3103     }
3104   if (SetImageStorageClass(fx_image,DirectClass,exception) == MagickFalse)
3105     {
3106       fx_info=DestroyFxThreadSet(fx_info);
3107       fx_image=DestroyImage(fx_image);
3108       return((Image *) NULL);
3109     }
3110   /*
3111     Fx image.
3112   */
3113   status=MagickTrue;
3114   progress=0;
3115   image_view=AcquireVirtualCacheView(image,exception);
3116   fx_view=AcquireAuthenticCacheView(fx_image,exception);
3117 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3118   #pragma omp parallel for schedule(static,4) shared(progress,status) \
3119     magick_threads(image,fx_image,fx_image->rows,1)
3120 #endif
3121   for (y=0; y < (ssize_t) fx_image->rows; y++)
3122   {
3123     const int
3124       id = GetOpenMPThreadId();
3125
3126     register const Quantum
3127       *magick_restrict p;
3128
3129     register Quantum
3130       *magick_restrict q;
3131
3132     register ssize_t
3133       x;
3134
3135     if (status == MagickFalse)
3136       continue;
3137     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3138     q=QueueCacheViewAuthenticPixels(fx_view,0,y,fx_image->columns,1,exception);
3139     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3140       {
3141         status=MagickFalse;
3142         continue;
3143       }
3144     for (x=0; x < (ssize_t) fx_image->columns; x++)
3145     {
3146       register ssize_t
3147         i;
3148
3149       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3150       {
3151         double
3152           alpha;
3153
3154         PixelChannel channel=GetPixelChannelChannel(image,i);
3155         PixelTrait traits=GetPixelChannelTraits(image,channel);
3156         PixelTrait fx_traits=GetPixelChannelTraits(fx_image,channel);
3157         if ((traits == UndefinedPixelTrait) ||
3158             (fx_traits == UndefinedPixelTrait))
3159           continue;
3160         if (((fx_traits & CopyPixelTrait) != 0) ||
3161             (GetPixelWriteMask(image,p) == 0))
3162           {
3163             SetPixelChannel(fx_image,channel,p[i],q);
3164             continue;
3165           }
3166         alpha=0.0;
3167         (void) FxEvaluateChannelExpression(fx_info[id],channel,x,y,&alpha,
3168           exception);
3169         q[i]=ClampToQuantum(QuantumRange*alpha);
3170       }
3171       p+=GetPixelChannels(image);
3172       q+=GetPixelChannels(fx_image);
3173     }
3174     if (SyncCacheViewAuthenticPixels(fx_view,exception) == MagickFalse)
3175       status=MagickFalse;
3176     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3177       {
3178         MagickBooleanType
3179           proceed;
3180
3181 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3182         #pragma omp critical (MagickCore_FxImage)
3183 #endif
3184         proceed=SetImageProgress(image,FxImageTag,progress++,image->rows);
3185         if (proceed == MagickFalse)
3186           status=MagickFalse;
3187       }
3188   }
3189   fx_view=DestroyCacheView(fx_view);
3190   image_view=DestroyCacheView(image_view);
3191   fx_info=DestroyFxThreadSet(fx_info);
3192   if (status == MagickFalse)
3193     fx_image=DestroyImage(fx_image);
3194   return(fx_image);
3195 }
3196 \f
3197 /*
3198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3199 %                                                                             %
3200 %                                                                             %
3201 %                                                                             %
3202 %     I m p l o d e I m a g e                                                 %
3203 %                                                                             %
3204 %                                                                             %
3205 %                                                                             %
3206 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3207 %
3208 %  ImplodeImage() creates a new image that is a copy of an existing
3209 %  one with the image pixels "implode" by the specified percentage.  It
3210 %  allocates the memory necessary for the new Image structure and returns a
3211 %  pointer to the new image.
3212 %
3213 %  The format of the ImplodeImage method is:
3214 %
3215 %      Image *ImplodeImage(const Image *image,const double amount,
3216 %        const PixelInterpolateMethod method,ExceptionInfo *exception)
3217 %
3218 %  A description of each parameter follows:
3219 %
3220 %    o implode_image: Method ImplodeImage returns a pointer to the image
3221 %      after it is implode.  A null image is returned if there is a memory
3222 %      shortage.
3223 %
3224 %    o image: the image.
3225 %
3226 %    o amount:  Define the extent of the implosion.
3227 %
3228 %    o method: the pixel interpolation method.
3229 %
3230 %    o exception: return any errors or warnings in this structure.
3231 %
3232 */
3233 MagickExport Image *ImplodeImage(const Image *image,const double amount,
3234   const PixelInterpolateMethod method,ExceptionInfo *exception)
3235 {
3236 #define ImplodeImageTag  "Implode/Image"
3237
3238   CacheView
3239     *canvas_view,
3240     *implode_view,
3241     *interpolate_view;
3242
3243   Image
3244     *canvas,
3245     *implode_image;
3246
3247   MagickBooleanType
3248     status;
3249
3250   MagickOffsetType
3251     progress;
3252
3253   double
3254     radius;
3255
3256   PointInfo
3257     center,
3258     scale;
3259
3260   ssize_t
3261     y;
3262
3263   /*
3264     Initialize implode image attributes.
3265   */
3266   assert(image != (Image *) NULL);
3267   assert(image->signature == MagickCoreSignature);
3268   if (image->debug != MagickFalse)
3269     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3270   assert(exception != (ExceptionInfo *) NULL);
3271   assert(exception->signature == MagickCoreSignature);
3272   canvas=CloneImage(image,0,0,MagickTrue,exception);
3273   if (canvas == (Image *) NULL)
3274     return((Image *) NULL);
3275   if ((canvas->alpha_trait == UndefinedPixelTrait) &&
3276       (canvas->background_color.alpha != OpaqueAlpha))
3277     (void) SetImageAlphaChannel(canvas,OpaqueAlphaChannel,exception);
3278   implode_image=CloneImage(canvas,canvas->columns,canvas->rows,MagickTrue,
3279     exception);
3280   if (implode_image == (Image *) NULL)
3281     {
3282       canvas=DestroyImage(canvas);
3283       return((Image *) NULL);
3284     }
3285   if (SetImageStorageClass(implode_image,DirectClass,exception) == MagickFalse)
3286     {
3287       canvas=DestroyImage(canvas);
3288       implode_image=DestroyImage(implode_image);
3289       return((Image *) NULL);
3290     }
3291   /*
3292     Compute scaling factor.
3293   */
3294   scale.x=1.0;
3295   scale.y=1.0;
3296   center.x=0.5*canvas->columns;
3297   center.y=0.5*canvas->rows;
3298   radius=center.x;
3299   if (canvas->columns > canvas->rows)
3300     scale.y=(double) canvas->columns/(double) canvas->rows;
3301   else
3302     if (canvas->columns < canvas->rows)
3303       {
3304         scale.x=(double) canvas->rows/(double) canvas->columns;
3305         radius=center.y;
3306       }
3307   /*
3308     Implode image.
3309   */
3310   status=MagickTrue;
3311   progress=0;
3312   canvas_view=AcquireVirtualCacheView(canvas,exception);
3313   interpolate_view=AcquireVirtualCacheView(canvas,exception);
3314   implode_view=AcquireAuthenticCacheView(implode_image,exception);
3315 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3316   #pragma omp parallel for schedule(static,4) shared(progress,status) \
3317     magick_threads(canvas,implode_image,canvas->rows,1)
3318 #endif
3319   for (y=0; y < (ssize_t) canvas->rows; y++)
3320   {
3321     double
3322       distance;
3323
3324     PointInfo
3325       delta;
3326
3327     register const Quantum
3328       *magick_restrict p;
3329
3330     register ssize_t
3331       x;
3332
3333     register Quantum
3334       *magick_restrict q;
3335
3336     if (status == MagickFalse)
3337       continue;
3338     p=GetCacheViewVirtualPixels(canvas_view,0,y,canvas->columns,1,exception);
3339     q=QueueCacheViewAuthenticPixels(implode_view,0,y,implode_image->columns,1,
3340       exception);
3341     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3342       {
3343         status=MagickFalse;
3344         continue;
3345       }
3346     delta.y=scale.y*(double) (y-center.y);
3347     for (x=0; x < (ssize_t) canvas->columns; x++)
3348     {
3349       register ssize_t
3350         i;
3351
3352       /*
3353         Determine if the pixel is within an ellipse.
3354       */
3355       if (GetPixelWriteMask(canvas,p) == 0)
3356         {
3357           SetPixelBackgoundColor(implode_image,q);
3358           p+=GetPixelChannels(canvas);
3359           q+=GetPixelChannels(implode_image);
3360           continue;
3361         }
3362       delta.x=scale.x*(double) (x-center.x);
3363       distance=delta.x*delta.x+delta.y*delta.y;
3364       if (distance >= (radius*radius))
3365         for (i=0; i < (ssize_t) GetPixelChannels(canvas); i++)
3366         {
3367           PixelChannel channel = GetPixelChannelChannel(canvas,i);
3368           PixelTrait traits = GetPixelChannelTraits(canvas,channel);
3369           PixelTrait implode_traits = GetPixelChannelTraits(implode_image,
3370             channel);
3371           if ((traits == UndefinedPixelTrait) ||
3372               (implode_traits == UndefinedPixelTrait))
3373             continue;
3374           SetPixelChannel(implode_image,channel,p[i],q);
3375         }
3376       else
3377         {
3378           double
3379             factor;
3380
3381           /*
3382             Implode the pixel.
3383           */
3384           factor=1.0;
3385           if (distance > 0.0)
3386             factor=pow(sin(MagickPI*sqrt((double) distance)/radius/2),-amount);
3387           status=InterpolatePixelChannels(canvas,interpolate_view,implode_image,
3388             method,(double) (factor*delta.x/scale.x+center.x),(double) (factor*
3389             delta.y/scale.y+center.y),q,exception);
3390         }
3391       p+=GetPixelChannels(canvas);
3392       q+=GetPixelChannels(implode_image);
3393     }
3394     if (SyncCacheViewAuthenticPixels(implode_view,exception) == MagickFalse)
3395       status=MagickFalse;
3396     if (canvas->progress_monitor != (MagickProgressMonitor) NULL)
3397       {
3398         MagickBooleanType
3399           proceed;
3400
3401 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3402         #pragma omp critical (MagickCore_ImplodeImage)
3403 #endif
3404         proceed=SetImageProgress(canvas,ImplodeImageTag,progress++,
3405           canvas->rows);
3406         if (proceed == MagickFalse)
3407           status=MagickFalse;
3408       }
3409   }
3410   implode_view=DestroyCacheView(implode_view);
3411   interpolate_view=DestroyCacheView(interpolate_view);
3412   canvas_view=DestroyCacheView(canvas_view);
3413   canvas=DestroyImage(canvas);
3414   if (status == MagickFalse)
3415     implode_image=DestroyImage(implode_image);
3416   return(implode_image);
3417 }
3418 \f
3419 /*
3420 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3421 %                                                                             %
3422 %                                                                             %
3423 %                                                                             %
3424 %     M o r p h I m a g e s                                                   %
3425 %                                                                             %
3426 %                                                                             %
3427 %                                                                             %
3428 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3429 %
3430 %  The MorphImages() method requires a minimum of two images.  The first
3431 %  image is transformed into the second by a number of intervening images
3432 %  as specified by frames.
3433 %
3434 %  The format of the MorphImage method is:
3435 %
3436 %      Image *MorphImages(const Image *image,const size_t number_frames,
3437 %        ExceptionInfo *exception)
3438 %
3439 %  A description of each parameter follows:
3440 %
3441 %    o image: the image.
3442 %
3443 %    o number_frames:  Define the number of in-between image to generate.
3444 %      The more in-between frames, the smoother the morph.
3445 %
3446 %    o exception: return any errors or warnings in this structure.
3447 %
3448 */
3449 MagickExport Image *MorphImages(const Image *image,const size_t number_frames,
3450   ExceptionInfo *exception)
3451 {
3452 #define MorphImageTag  "Morph/Image"
3453
3454   double
3455     alpha,
3456     beta;
3457
3458   Image
3459     *morph_image,
3460     *morph_images;
3461
3462   MagickBooleanType
3463     status;
3464
3465   MagickOffsetType
3466     scene;
3467
3468   register const Image
3469     *next;
3470
3471   register ssize_t
3472     n;
3473
3474   ssize_t
3475     y;
3476
3477   /*
3478     Clone first frame in sequence.
3479   */
3480   assert(image != (Image *) NULL);
3481   assert(image->signature == MagickCoreSignature);
3482   if (image->debug != MagickFalse)
3483     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3484   assert(exception != (ExceptionInfo *) NULL);
3485   assert(exception->signature == MagickCoreSignature);
3486   morph_images=CloneImage(image,0,0,MagickTrue,exception);
3487   if (morph_images == (Image *) NULL)
3488     return((Image *) NULL);
3489   if (GetNextImageInList(image) == (Image *) NULL)
3490     {
3491       /*
3492         Morph single image.
3493       */
3494       for (n=1; n < (ssize_t) number_frames; n++)
3495       {
3496         morph_image=CloneImage(image,0,0,MagickTrue,exception);
3497         if (morph_image == (Image *) NULL)
3498           {
3499             morph_images=DestroyImageList(morph_images);
3500             return((Image *) NULL);
3501           }
3502         AppendImageToList(&morph_images,morph_image);
3503         if (image->progress_monitor != (MagickProgressMonitor) NULL)
3504           {
3505             MagickBooleanType
3506               proceed;
3507
3508             proceed=SetImageProgress(image,MorphImageTag,(MagickOffsetType) n,
3509               number_frames);
3510             if (proceed == MagickFalse)
3511               status=MagickFalse;
3512           }
3513       }
3514       return(GetFirstImageInList(morph_images));
3515     }
3516   /*
3517     Morph image sequence.
3518   */
3519   status=MagickTrue;
3520   scene=0;
3521   next=image;
3522   for ( ; GetNextImageInList(next) != (Image *) NULL; next=GetNextImageInList(next))
3523   {
3524     for (n=0; n < (ssize_t) number_frames; n++)
3525     {
3526       CacheView
3527         *image_view,
3528         *morph_view;
3529
3530       beta=(double) (n+1.0)/(double) (number_frames+1.0);
3531       alpha=1.0-beta;
3532       morph_image=ResizeImage(next,(size_t) (alpha*next->columns+beta*
3533         GetNextImageInList(next)->columns+0.5),(size_t) (alpha*next->rows+beta*
3534         GetNextImageInList(next)->rows+0.5),next->filter,exception);
3535       if (morph_image == (Image *) NULL)
3536         {
3537           morph_images=DestroyImageList(morph_images);
3538           return((Image *) NULL);
3539         }
3540       status=SetImageStorageClass(morph_image,DirectClass,exception);
3541       if (status == MagickFalse)
3542         {
3543           morph_image=DestroyImage(morph_image);
3544           return((Image *) NULL);
3545         }
3546       AppendImageToList(&morph_images,morph_image);
3547       morph_images=GetLastImageInList(morph_images);
3548       morph_image=ResizeImage(GetNextImageInList(next),morph_images->columns,
3549         morph_images->rows,GetNextImageInList(next)->filter,exception);
3550       if (morph_image == (Image *) NULL)
3551         {
3552           morph_images=DestroyImageList(morph_images);
3553           return((Image *) NULL);
3554         }
3555       image_view=AcquireVirtualCacheView(morph_image,exception);
3556       morph_view=AcquireAuthenticCacheView(morph_images,exception);
3557 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3558       #pragma omp parallel for schedule(static,4) shared(status) \
3559         magick_threads(morph_image,morph_image,morph_image->rows,1)
3560 #endif
3561       for (y=0; y < (ssize_t) morph_images->rows; y++)
3562       {
3563         MagickBooleanType
3564           sync;
3565
3566         register const Quantum
3567           *magick_restrict p;
3568
3569         register ssize_t
3570           x;
3571
3572         register Quantum
3573           *magick_restrict q;
3574
3575         if (status == MagickFalse)
3576           continue;
3577         p=GetCacheViewVirtualPixels(image_view,0,y,morph_image->columns,1,
3578           exception);
3579         q=GetCacheViewAuthenticPixels(morph_view,0,y,morph_images->columns,1,
3580           exception);
3581         if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3582           {
3583             status=MagickFalse;
3584             continue;
3585           }
3586         for (x=0; x < (ssize_t) morph_images->columns; x++)
3587         {
3588           register ssize_t
3589             i;
3590
3591           for (i=0; i < (ssize_t) GetPixelChannels(morph_image); i++)
3592           {
3593             PixelChannel channel=GetPixelChannelChannel(morph_image,i);
3594             PixelTrait traits=GetPixelChannelTraits(morph_image,channel);
3595             PixelTrait morph_traits=GetPixelChannelTraits(morph_images,channel);
3596             if ((traits == UndefinedPixelTrait) ||
3597                 (morph_traits == UndefinedPixelTrait))
3598               continue;
3599             if (((morph_traits & CopyPixelTrait) != 0) ||
3600                 (GetPixelWriteMask(morph_images,p) == 0))
3601               {
3602                 SetPixelChannel(morph_image,channel,p[i],q);
3603                 continue;
3604               }
3605             SetPixelChannel(morph_image,channel,ClampToQuantum(alpha*
3606               GetPixelChannel(morph_images,channel,q)+beta*p[i]),q);
3607           }
3608           p+=GetPixelChannels(morph_image);
3609           q+=GetPixelChannels(morph_images);
3610         }
3611         sync=SyncCacheViewAuthenticPixels(morph_view,exception);
3612         if (sync == MagickFalse)
3613           status=MagickFalse;
3614       }
3615       morph_view=DestroyCacheView(morph_view);
3616       image_view=DestroyCacheView(image_view);
3617       morph_image=DestroyImage(morph_image);
3618     }
3619     if (n < (ssize_t) number_frames)
3620       break;
3621     /*
3622       Clone last frame in sequence.
3623     */
3624     morph_image=CloneImage(GetNextImageInList(next),0,0,MagickTrue,exception);
3625     if (morph_image == (Image *) NULL)
3626       {
3627         morph_images=DestroyImageList(morph_images);
3628         return((Image *) NULL);
3629       }
3630     AppendImageToList(&morph_images,morph_image);
3631     morph_images=GetLastImageInList(morph_images);
3632     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3633       {
3634         MagickBooleanType
3635           proceed;
3636
3637 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3638         #pragma omp critical (MagickCore_MorphImages)
3639 #endif
3640         proceed=SetImageProgress(image,MorphImageTag,scene,
3641           GetImageListLength(image));
3642         if (proceed == MagickFalse)
3643           status=MagickFalse;
3644       }
3645     scene++;
3646   }
3647   if (GetNextImageInList(next) != (Image *) NULL)
3648     {
3649       morph_images=DestroyImageList(morph_images);
3650       return((Image *) NULL);
3651     }
3652   return(GetFirstImageInList(morph_images));
3653 }
3654 \f
3655 /*
3656 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3657 %                                                                             %
3658 %                                                                             %
3659 %                                                                             %
3660 %     P l a s m a I m a g e                                                   %
3661 %                                                                             %
3662 %                                                                             %
3663 %                                                                             %
3664 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3665 %
3666 %  PlasmaImage() initializes an image with plasma fractal values.  The image
3667 %  must be initialized with a base color and the random number generator
3668 %  seeded before this method is called.
3669 %
3670 %  The format of the PlasmaImage method is:
3671 %
3672 %      MagickBooleanType PlasmaImage(Image *image,const SegmentInfo *segment,
3673 %        size_t attenuate,size_t depth,ExceptionInfo *exception)
3674 %
3675 %  A description of each parameter follows:
3676 %
3677 %    o image: the image.
3678 %
3679 %    o segment:   Define the region to apply plasma fractals values.
3680 %
3681 %    o attenuate: Define the plasma attenuation factor.
3682 %
3683 %    o depth: Limit the plasma recursion depth.
3684 %
3685 %    o exception: return any errors or warnings in this structure.
3686 %
3687 */
3688
3689 static inline Quantum PlasmaPixel(RandomInfo *random_info,
3690   const double pixel,const double noise)
3691 {
3692   Quantum
3693     plasma;
3694
3695   plasma=ClampToQuantum(pixel+noise*GetPseudoRandomValue(random_info)-
3696     noise/2.0);
3697   if (plasma <= 0)
3698     return((Quantum) 0);
3699   if (plasma >= QuantumRange)
3700     return(QuantumRange);
3701   return(plasma);
3702 }
3703
3704 static MagickBooleanType PlasmaImageProxy(Image *image,CacheView *image_view,
3705   CacheView *u_view,CacheView *v_view,RandomInfo *random_info,
3706   const SegmentInfo *segment,size_t attenuate,size_t depth,
3707   ExceptionInfo *exception)
3708 {
3709   double
3710     plasma;
3711
3712   register const Quantum
3713     *magick_restrict u,
3714     *magick_restrict v;
3715
3716   register Quantum
3717     *magick_restrict q;
3718
3719   register ssize_t
3720     i;
3721
3722   ssize_t
3723     x,
3724     x_mid,
3725     y,
3726     y_mid;
3727
3728   if ((fabs(segment->x2-segment->x1) <= MagickEpsilon) &&
3729       (fabs(segment->y2-segment->y1) <= MagickEpsilon))
3730     return(MagickTrue);
3731   if (depth != 0)
3732     {
3733       MagickBooleanType
3734         status;
3735
3736       SegmentInfo
3737         local_info;
3738
3739       /*
3740         Divide the area into quadrants and recurse.
3741       */
3742       depth--;
3743       attenuate++;
3744       x_mid=(ssize_t) ceil((segment->x1+segment->x2)/2-0.5);
3745       y_mid=(ssize_t) ceil((segment->y1+segment->y2)/2-0.5);
3746       local_info=(*segment);
3747       local_info.x2=(double) x_mid;
3748       local_info.y2=(double) y_mid;
3749       (void) PlasmaImageProxy(image,image_view,u_view,v_view,random_info,
3750         &local_info,attenuate,depth,exception);
3751       local_info=(*segment);
3752       local_info.y1=(double) y_mid;
3753       local_info.x2=(double) x_mid;
3754       (void) PlasmaImageProxy(image,image_view,u_view,v_view,random_info,
3755         &local_info,attenuate,depth,exception);
3756       local_info=(*segment);
3757       local_info.x1=(double) x_mid;
3758       local_info.y2=(double) y_mid;
3759       (void) PlasmaImageProxy(image,image_view,u_view,v_view,random_info,
3760         &local_info,attenuate,depth,exception);
3761       local_info=(*segment);
3762       local_info.x1=(double) x_mid;
3763       local_info.y1=(double) y_mid;
3764       status=PlasmaImageProxy(image,image_view,u_view,v_view,random_info,
3765         &local_info,attenuate,depth,exception);
3766       return(status);
3767     }
3768   x_mid=(ssize_t) ceil((segment->x1+segment->x2)/2-0.5);
3769   y_mid=(ssize_t) ceil((segment->y1+segment->y2)/2-0.5);
3770   if ((fabs(segment->x1-x_mid) < MagickEpsilon) &&
3771       (fabs(segment->x2-x_mid) < MagickEpsilon) &&
3772       (fabs(segment->y1-y_mid) < MagickEpsilon) &&
3773       (fabs(segment->y2-y_mid) < MagickEpsilon))
3774     return(MagickFalse);
3775   /*
3776     Average pixels and apply plasma.
3777   */
3778   plasma=(double) QuantumRange/(2.0*attenuate);
3779   if ((fabs(segment->x1-x_mid) > MagickEpsilon) ||
3780       (fabs(segment->x2-x_mid) > MagickEpsilon))
3781     {
3782       /*
3783         Left pixel.
3784       */
3785       x=(ssize_t) ceil(segment->x1-0.5);
3786       u=GetCacheViewVirtualPixels(u_view,x,(ssize_t) ceil(segment->y1-0.5),1,1,
3787         exception);
3788       v=GetCacheViewVirtualPixels(v_view,x,(ssize_t) ceil(segment->y2-0.5),1,1,
3789         exception);
3790       q=QueueCacheViewAuthenticPixels(image_view,x,y_mid,1,1,exception);
3791       if ((u == (const Quantum *) NULL) || (v == (const Quantum *) NULL) ||
3792           (q == (Quantum *) NULL))
3793         return(MagickTrue);
3794       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3795       {
3796         PixelChannel channel=GetPixelChannelChannel(image,i);
3797         PixelTrait traits=GetPixelChannelTraits(image,channel);
3798         if (traits == UndefinedPixelTrait)
3799           continue;
3800         q[i]=PlasmaPixel(random_info,(u[i]+v[i])/2.0,plasma);
3801       }
3802       (void) SyncCacheViewAuthenticPixels(image_view,exception);
3803       if (fabs(segment->x1-segment->x2) > MagickEpsilon)
3804         {
3805           /*
3806             Right pixel.
3807           */
3808           x=(ssize_t) ceil(segment->x2-0.5);
3809           u=GetCacheViewVirtualPixels(u_view,x,(ssize_t) ceil(segment->y1-0.5),
3810             1,1,exception);
3811           v=GetCacheViewVirtualPixels(v_view,x,(ssize_t) ceil(segment->y2-0.5),
3812             1,1,exception);
3813           q=QueueCacheViewAuthenticPixels(image_view,x,y_mid,1,1,exception);
3814           if ((u == (const Quantum *) NULL) || (v == (const Quantum *) NULL) ||
3815               (q == (Quantum *) NULL))
3816             return(MagickTrue);
3817           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3818           {
3819             PixelChannel channel=GetPixelChannelChannel(image,i);
3820             PixelTrait traits=GetPixelChannelTraits(image,channel);
3821             if (traits == UndefinedPixelTrait)
3822               continue;
3823             q[i]=PlasmaPixel(random_info,(u[i]+v[i])/2.0,plasma);
3824           }
3825           (void) SyncCacheViewAuthenticPixels(image_view,exception);
3826         }
3827     }
3828   if ((fabs(segment->y1-y_mid) > MagickEpsilon) ||
3829       (fabs(segment->y2-y_mid) > MagickEpsilon))
3830     {
3831       if ((fabs(segment->x1-x_mid) > MagickEpsilon) ||
3832           (fabs(segment->y2-y_mid) > MagickEpsilon))
3833         {
3834           /*
3835             Bottom pixel.
3836           */
3837           y=(ssize_t) ceil(segment->y2-0.5);
3838           u=GetCacheViewVirtualPixels(u_view,(ssize_t) ceil(segment->x1-0.5),y,
3839             1,1,exception);
3840           v=GetCacheViewVirtualPixels(v_view,(ssize_t) ceil(segment->x2-0.5),y,
3841             1,1,exception);
3842           q=QueueCacheViewAuthenticPixels(image_view,x_mid,y,1,1,exception);
3843           if ((u == (const Quantum *) NULL) || (v == (const Quantum *) NULL) ||
3844               (q == (Quantum *) NULL))
3845             return(MagickTrue);
3846           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3847           {
3848             PixelChannel channel=GetPixelChannelChannel(image,i);
3849             PixelTrait traits=GetPixelChannelTraits(image,channel);
3850             if (traits == UndefinedPixelTrait)
3851               continue;
3852             q[i]=PlasmaPixel(random_info,(u[i]+v[i])/2.0,plasma);
3853           }
3854           (void) SyncCacheViewAuthenticPixels(image_view,exception);
3855         }
3856       if (fabs(segment->y1-segment->y2) > MagickEpsilon)
3857         {
3858           /*
3859             Top pixel.
3860           */
3861           y=(ssize_t) ceil(segment->y1-0.5);
3862           u=GetCacheViewVirtualPixels(u_view,(ssize_t) ceil(segment->x1-0.5),y,
3863             1,1,exception);
3864           v=GetCacheViewVirtualPixels(v_view,(ssize_t) ceil(segment->x2-0.5),y,
3865             1,1,exception);
3866           q=QueueCacheViewAuthenticPixels(image_view,x_mid,y,1,1,exception);
3867           if ((u == (const Quantum *) NULL) || (v == (const Quantum *) NULL) ||
3868               (q == (Quantum *) NULL))
3869             return(MagickTrue);
3870           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3871           {
3872             PixelChannel channel=GetPixelChannelChannel(image,i);
3873             PixelTrait traits=GetPixelChannelTraits(image,channel);
3874             if (traits == UndefinedPixelTrait)
3875               continue;
3876             q[i]=PlasmaPixel(random_info,(u[i]+v[i])/2.0,plasma);
3877           }
3878           (void) SyncCacheViewAuthenticPixels(image_view,exception);
3879         }
3880     }
3881   if ((fabs(segment->x1-segment->x2) > MagickEpsilon) ||
3882       (fabs(segment->y1-segment->y2) > MagickEpsilon))
3883     {
3884       /*
3885         Middle pixel.
3886       */
3887       x=(ssize_t) ceil(segment->x1-0.5);
3888       y=(ssize_t) ceil(segment->y1-0.5);
3889       u=GetCacheViewVirtualPixels(u_view,x,y,1,1,exception);
3890       x=(ssize_t) ceil(segment->x2-0.5);
3891       y=(ssize_t) ceil(segment->y2-0.5);
3892       v=GetCacheViewVirtualPixels(v_view,x,y,1,1,exception);
3893       q=QueueCacheViewAuthenticPixels(image_view,x_mid,y_mid,1,1,exception);
3894       if ((u == (const Quantum *) NULL) || (v == (const Quantum *) NULL) ||
3895           (q == (Quantum *) NULL))
3896         return(MagickTrue);
3897       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3898       {
3899         PixelChannel channel=GetPixelChannelChannel(image,i);
3900         PixelTrait traits=GetPixelChannelTraits(image,channel);
3901         if (traits == UndefinedPixelTrait)
3902           continue;
3903         q[i]=PlasmaPixel(random_info,(u[i]+v[i])/2.0,plasma);
3904       }
3905       (void) SyncCacheViewAuthenticPixels(image_view,exception);
3906     }
3907   if ((fabs(segment->x2-segment->x1) < 3.0) &&
3908       (fabs(segment->y2-segment->y1) < 3.0))
3909     return(MagickTrue);
3910   return(MagickFalse);
3911 }
3912
3913 MagickExport MagickBooleanType PlasmaImage(Image *image,
3914   const SegmentInfo *segment,size_t attenuate,size_t depth,
3915   ExceptionInfo *exception)
3916 {
3917   CacheView
3918     *image_view,
3919     *u_view,
3920     *v_view;
3921
3922   MagickBooleanType
3923     status;
3924
3925   RandomInfo
3926     *random_info;
3927
3928   if (image->debug != MagickFalse)
3929     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
3930   assert(image != (Image *) NULL);
3931   assert(image->signature == MagickCoreSignature);
3932   if (image->debug != MagickFalse)
3933     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
3934   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
3935     return(MagickFalse);
3936   image_view=AcquireAuthenticCacheView(image,exception);
3937   u_view=AcquireVirtualCacheView(image,exception);
3938   v_view=AcquireVirtualCacheView(image,exception);
3939   random_info=AcquireRandomInfo();
3940   status=PlasmaImageProxy(image,image_view,u_view,v_view,random_info,segment,
3941     attenuate,depth,exception);
3942   random_info=DestroyRandomInfo(random_info);
3943   v_view=DestroyCacheView(v_view);
3944   u_view=DestroyCacheView(u_view);
3945   image_view=DestroyCacheView(image_view);
3946   return(status);
3947 }
3948 \f
3949 /*
3950 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3951 %                                                                             %
3952 %                                                                             %
3953 %                                                                             %
3954 %   P o l a r o i d I m a g e                                                 %
3955 %                                                                             %
3956 %                                                                             %
3957 %                                                                             %
3958 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3959 %
3960 %  PolaroidImage() simulates a Polaroid picture.
3961 %
3962 %  The format of the PolaroidImage method is:
3963 %
3964 %      Image *PolaroidImage(const Image *image,const DrawInfo *draw_info,
3965 %        const char *caption,const double angle,
3966 %        const PixelInterpolateMethod method,ExceptionInfo exception)
3967 %
3968 %  A description of each parameter follows:
3969 %
3970 %    o image: the image.
3971 %
3972 %    o draw_info: the draw info.
3973 %
3974 %    o caption: the Polaroid caption.
3975 %
3976 %    o angle: Apply the effect along this angle.
3977 %
3978 %    o method: the pixel interpolation method.
3979 %
3980 %    o exception: return any errors or warnings in this structure.
3981 %
3982 */
3983 MagickExport Image *PolaroidImage(const Image *image,const DrawInfo *draw_info,
3984   const char *caption,const double angle,const PixelInterpolateMethod method,
3985   ExceptionInfo *exception)
3986 {
3987   Image
3988     *bend_image,
3989     *caption_image,
3990     *flop_image,
3991     *picture_image,
3992     *polaroid_image,
3993     *rotate_image,
3994     *trim_image;
3995
3996   size_t
3997     height;
3998
3999   ssize_t
4000     quantum;
4001
4002   /*
4003     Simulate a Polaroid picture.
4004   */
4005   assert(image != (Image *) NULL);
4006   assert(image->signature == MagickCoreSignature);
4007   if (image->debug != MagickFalse)
4008     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4009   assert(exception != (ExceptionInfo *) NULL);
4010   assert(exception->signature == MagickCoreSignature);
4011   quantum=(ssize_t) MagickMax(MagickMax((double) image->columns,(double)
4012     image->rows)/25.0,10.0);
4013   height=image->rows+2*quantum;
4014   caption_image=(Image *) NULL;
4015   if (caption != (const char *) NULL)
4016     {
4017       char
4018         geometry[MagickPathExtent],
4019         *text;
4020
4021       DrawInfo
4022         *annotate_info;
4023
4024       ImageInfo
4025         *image_info;
4026
4027       MagickBooleanType
4028         status;
4029
4030       ssize_t
4031         count;
4032
4033       TypeMetric
4034         metrics;
4035
4036       /*
4037         Generate caption image.
4038       */
4039       caption_image=CloneImage(image,image->columns,1,MagickTrue,exception);
4040       if (caption_image == (Image *) NULL)
4041         return((Image *) NULL);
4042       image_info=AcquireImageInfo();
4043       annotate_info=CloneDrawInfo((const ImageInfo *) NULL,draw_info);
4044       text=InterpretImageProperties(image_info,(Image *) image,caption,
4045         exception);
4046       image_info=DestroyImageInfo(image_info);
4047       (void) CloneString(&annotate_info->text,text);
4048       count=FormatMagickCaption(caption_image,annotate_info,MagickTrue,&metrics,
4049         &text,exception);
4050       status=SetImageExtent(caption_image,image->columns,(size_t) ((count+1)*
4051         (metrics.ascent-metrics.descent)+0.5),exception);
4052       if (status == MagickFalse)
4053         caption_image=DestroyImage(caption_image);
4054       else
4055         {
4056           caption_image->background_color=image->border_color;
4057           (void) SetImageBackgroundColor(caption_image,exception);
4058           (void) CloneString(&annotate_info->text,text);
4059           (void) FormatLocaleString(geometry,MagickPathExtent,"+0+%g",
4060             metrics.ascent);
4061           if (annotate_info->gravity == UndefinedGravity)
4062             (void) CloneString(&annotate_info->geometry,AcquireString(
4063               geometry));
4064           (void) AnnotateImage(caption_image,annotate_info,exception);
4065           height+=caption_image->rows;
4066         }
4067       annotate_info=DestroyDrawInfo(annotate_info);
4068       text=DestroyString(text);
4069     }
4070   picture_image=CloneImage(image,image->columns+2*quantum,height,MagickTrue,
4071     exception);
4072   if (picture_image == (Image *) NULL)
4073     {
4074       if (caption_image != (Image *) NULL)
4075         caption_image=DestroyImage(caption_image);
4076       return((Image *) NULL);
4077     }
4078   picture_image->background_color=image->border_color;
4079   (void) SetImageBackgroundColor(picture_image,exception);
4080   (void) CompositeImage(picture_image,image,OverCompositeOp,MagickTrue,quantum,
4081     quantum,exception);
4082   if (caption_image != (Image *) NULL)
4083     {
4084       (void) CompositeImage(picture_image,caption_image,OverCompositeOp,
4085         MagickTrue,quantum,(ssize_t) (image->rows+3*quantum/2),exception);
4086       caption_image=DestroyImage(caption_image);
4087     }
4088   (void) QueryColorCompliance("none",AllCompliance,
4089     &picture_image->background_color,exception);
4090   (void) SetImageAlphaChannel(picture_image,OpaqueAlphaChannel,exception);
4091   rotate_image=RotateImage(picture_image,90.0,exception);
4092   picture_image=DestroyImage(picture_image);
4093   if (rotate_image == (Image *) NULL)
4094     return((Image *) NULL);
4095   picture_image=rotate_image;
4096   bend_image=WaveImage(picture_image,0.01*picture_image->rows,2.0*
4097     picture_image->columns,method,exception);
4098   picture_image=DestroyImage(picture_image);
4099   if (bend_image == (Image *) NULL)
4100     return((Image *) NULL);
4101   picture_image=bend_image;
4102   rotate_image=RotateImage(picture_image,-90.0,exception);
4103   picture_image=DestroyImage(picture_image);
4104   if (rotate_image == (Image *) NULL)
4105     return((Image *) NULL);
4106   picture_image=rotate_image;
4107   picture_image->background_color=image->background_color;
4108   polaroid_image=ShadowImage(picture_image,80.0,2.0,quantum/3,quantum/3,
4109     exception);
4110   if (polaroid_image == (Image *) NULL)
4111     {
4112       picture_image=DestroyImage(picture_image);
4113       return(picture_image);
4114     }
4115   flop_image=FlopImage(polaroid_image,exception);
4116   polaroid_image=DestroyImage(polaroid_image);
4117   if (flop_image == (Image *) NULL)
4118     {
4119       picture_image=DestroyImage(picture_image);
4120       return(picture_image);
4121     }
4122   polaroid_image=flop_image;
4123   (void) CompositeImage(polaroid_image,picture_image,OverCompositeOp,
4124     MagickTrue,(ssize_t) (-0.01*picture_image->columns/2.0),0L,exception);
4125   picture_image=DestroyImage(picture_image);
4126   (void) QueryColorCompliance("none",AllCompliance,
4127     &polaroid_image->background_color,exception);
4128   rotate_image=RotateImage(polaroid_image,angle,exception);
4129   polaroid_image=DestroyImage(polaroid_image);
4130   if (rotate_image == (Image *) NULL)
4131     return((Image *) NULL);
4132   polaroid_image=rotate_image;
4133   trim_image=TrimImage(polaroid_image,exception);
4134   polaroid_image=DestroyImage(polaroid_image);
4135   if (trim_image == (Image *) NULL)
4136     return((Image *) NULL);
4137   polaroid_image=trim_image;
4138   return(polaroid_image);
4139 }
4140 \f
4141 /*
4142 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4143 %                                                                             %
4144 %                                                                             %
4145 %                                                                             %
4146 %     S e p i a T o n e I m a g e                                             %
4147 %                                                                             %
4148 %                                                                             %
4149 %                                                                             %
4150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4151 %
4152 %  MagickSepiaToneImage() applies a special effect to the image, similar to the
4153 %  effect achieved in a photo darkroom by sepia toning.  Threshold ranges from
4154 %  0 to QuantumRange and is a measure of the extent of the sepia toning.  A
4155 %  threshold of 80% is a good starting point for a reasonable tone.
4156 %
4157 %  The format of the SepiaToneImage method is:
4158 %
4159 %      Image *SepiaToneImage(const Image *image,const double threshold,
4160 %        ExceptionInfo *exception)
4161 %
4162 %  A description of each parameter follows:
4163 %
4164 %    o image: the image.
4165 %
4166 %    o threshold: the tone threshold.
4167 %
4168 %    o exception: return any errors or warnings in this structure.
4169 %
4170 */
4171 MagickExport Image *SepiaToneImage(const Image *image,const double threshold,
4172   ExceptionInfo *exception)
4173 {
4174 #define SepiaToneImageTag  "SepiaTone/Image"
4175
4176   CacheView
4177     *image_view,
4178     *sepia_view;
4179
4180   Image
4181     *sepia_image;
4182
4183   MagickBooleanType
4184     status;
4185
4186   MagickOffsetType
4187     progress;
4188
4189   ssize_t
4190     y;
4191
4192   /*
4193     Initialize sepia-toned image attributes.
4194   */
4195   assert(image != (const Image *) NULL);
4196   assert(image->signature == MagickCoreSignature);
4197   if (image->debug != MagickFalse)
4198     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4199   assert(exception != (ExceptionInfo *) NULL);
4200   assert(exception->signature == MagickCoreSignature);
4201   sepia_image=CloneImage(image,0,0,MagickTrue,exception);
4202   if (sepia_image == (Image *) NULL)
4203     return((Image *) NULL);
4204   if (SetImageStorageClass(sepia_image,DirectClass,exception) == MagickFalse)
4205     {
4206       sepia_image=DestroyImage(sepia_image);
4207       return((Image *) NULL);
4208     }
4209   /*
4210     Tone each row of the image.
4211   */
4212   status=MagickTrue;
4213   progress=0;
4214   image_view=AcquireVirtualCacheView(image,exception);
4215   sepia_view=AcquireAuthenticCacheView(sepia_image,exception);
4216 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4217   #pragma omp parallel for schedule(static,4) shared(progress,status) \
4218     magick_threads(image,sepia_image,image->rows,1)
4219 #endif
4220   for (y=0; y < (ssize_t) image->rows; y++)
4221   {
4222     register const Quantum
4223       *magick_restrict p;
4224
4225     register ssize_t
4226       x;
4227
4228     register Quantum
4229       *magick_restrict q;
4230
4231     if (status == MagickFalse)
4232       continue;
4233     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
4234     q=GetCacheViewAuthenticPixels(sepia_view,0,y,sepia_image->columns,1,
4235       exception);
4236     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
4237       {
4238         status=MagickFalse;
4239         continue;
4240       }
4241     for (x=0; x < (ssize_t) image->columns; x++)
4242     {
4243       double
4244         intensity,
4245         tone;
4246
4247       intensity=GetPixelIntensity(image,p);
4248       tone=intensity > threshold ? (double) QuantumRange : intensity+
4249         (double) QuantumRange-threshold;
4250       SetPixelRed(sepia_image,ClampToQuantum(tone),q);
4251       tone=intensity > (7.0*threshold/6.0) ? (double) QuantumRange :
4252         intensity+(double) QuantumRange-7.0*threshold/6.0;
4253       SetPixelGreen(sepia_image,ClampToQuantum(tone),q);
4254       tone=intensity < (threshold/6.0) ? 0 : intensity-threshold/6.0;
4255       SetPixelBlue(sepia_image,ClampToQuantum(tone),q);
4256       tone=threshold/7.0;
4257       if ((double) GetPixelGreen(image,q) < tone)
4258         SetPixelGreen(sepia_image,ClampToQuantum(tone),q);
4259       if ((double) GetPixelBlue(image,q) < tone)
4260         SetPixelBlue(sepia_image,ClampToQuantum(tone),q);
4261       SetPixelAlpha(sepia_image,GetPixelAlpha(image,p),q);
4262       p+=GetPixelChannels(image);
4263       q+=GetPixelChannels(sepia_image);
4264     }
4265     if (SyncCacheViewAuthenticPixels(sepia_view,exception) == MagickFalse)
4266       status=MagickFalse;
4267     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4268       {
4269         MagickBooleanType
4270           proceed;
4271
4272 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4273         #pragma omp critical (MagickCore_SepiaToneImage)
4274 #endif
4275         proceed=SetImageProgress(image,SepiaToneImageTag,progress++,
4276           image->rows);
4277         if (proceed == MagickFalse)
4278           status=MagickFalse;
4279       }
4280   }
4281   sepia_view=DestroyCacheView(sepia_view);
4282   image_view=DestroyCacheView(image_view);
4283   (void) NormalizeImage(sepia_image,exception);
4284   (void) ContrastImage(sepia_image,MagickTrue,exception);
4285   if (status == MagickFalse)
4286     sepia_image=DestroyImage(sepia_image);
4287   return(sepia_image);
4288 }
4289 \f
4290 /*
4291 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4292 %                                                                             %
4293 %                                                                             %
4294 %                                                                             %
4295 %     S h a d o w I m a g e                                                   %
4296 %                                                                             %
4297 %                                                                             %
4298 %                                                                             %
4299 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4300 %
4301 %  ShadowImage() simulates a shadow from the specified image and returns it.
4302 %
4303 %  The format of the ShadowImage method is:
4304 %
4305 %      Image *ShadowImage(const Image *image,const double alpha,
4306 %        const double sigma,const ssize_t x_offset,const ssize_t y_offset,
4307 %        ExceptionInfo *exception)
4308 %
4309 %  A description of each parameter follows:
4310 %
4311 %    o image: the image.
4312 %
4313 %    o alpha: percentage transparency.
4314 %
4315 %    o sigma: the standard deviation of the Gaussian, in pixels.
4316 %
4317 %    o x_offset: the shadow x-offset.
4318 %
4319 %    o y_offset: the shadow y-offset.
4320 %
4321 %    o exception: return any errors or warnings in this structure.
4322 %
4323 */
4324 MagickExport Image *ShadowImage(const Image *image,const double alpha,
4325   const double sigma,const ssize_t x_offset,const ssize_t y_offset,
4326   ExceptionInfo *exception)
4327 {
4328 #define ShadowImageTag  "Shadow/Image"
4329
4330   CacheView
4331     *image_view;
4332
4333   ChannelType
4334     channel_mask;
4335
4336   Image
4337     *border_image,
4338     *clone_image,
4339     *shadow_image;
4340
4341   MagickBooleanType
4342     status;
4343
4344   PixelInfo
4345     background_color;
4346
4347   RectangleInfo
4348     border_info;
4349
4350   ssize_t
4351     y;
4352
4353   assert(image != (Image *) NULL);
4354   assert(image->signature == MagickCoreSignature);
4355   if (image->debug != MagickFalse)
4356     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4357   assert(exception != (ExceptionInfo *) NULL);
4358   assert(exception->signature == MagickCoreSignature);
4359   clone_image=CloneImage(image,0,0,MagickTrue,exception);
4360   if (clone_image == (Image *) NULL)
4361     return((Image *) NULL);
4362   if (IsGrayColorspace(image->colorspace) != MagickFalse)
4363     (void) SetImageColorspace(clone_image,sRGBColorspace,exception);
4364   (void) SetImageVirtualPixelMethod(clone_image,EdgeVirtualPixelMethod,
4365     exception);
4366   border_info.width=(size_t) floor(2.0*sigma+0.5);
4367   border_info.height=(size_t) floor(2.0*sigma+0.5);
4368   border_info.x=0;
4369   border_info.y=0;
4370   (void) QueryColorCompliance("none",AllCompliance,&clone_image->border_color,
4371     exception);
4372   clone_image->alpha_trait=BlendPixelTrait;
4373   border_image=BorderImage(clone_image,&border_info,OverCompositeOp,exception);
4374   clone_image=DestroyImage(clone_image);
4375   if (border_image == (Image *) NULL)
4376     return((Image *) NULL);
4377   if (border_image->alpha_trait == UndefinedPixelTrait)
4378     (void) SetImageAlphaChannel(border_image,OpaqueAlphaChannel,exception);
4379   /*
4380     Shadow image.
4381   */
4382   status=MagickTrue;
4383   background_color=border_image->background_color;
4384   background_color.alpha_trait=BlendPixelTrait;
4385   image_view=AcquireAuthenticCacheView(border_image,exception);
4386   for (y=0; y < (ssize_t) border_image->rows; y++)
4387   {
4388     register Quantum
4389       *magick_restrict q;
4390
4391     register ssize_t
4392       x;
4393
4394     if (status == MagickFalse)
4395       continue;
4396     q=QueueCacheViewAuthenticPixels(image_view,0,y,border_image->columns,1,
4397       exception);
4398     if (q == (Quantum *) NULL)
4399       {
4400         status=MagickFalse;
4401         continue;
4402       }
4403     for (x=0; x < (ssize_t) border_image->columns; x++)
4404     {
4405       if (border_image->alpha_trait != UndefinedPixelTrait)
4406         background_color.alpha=GetPixelAlpha(border_image,q)*alpha/100.0;
4407       SetPixelViaPixelInfo(border_image,&background_color,q);
4408       q+=GetPixelChannels(border_image);
4409     }
4410     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
4411       status=MagickFalse;
4412   }
4413   image_view=DestroyCacheView(image_view);
4414   if (status == MagickFalse)
4415     {
4416       border_image=DestroyImage(border_image);
4417       return((Image *) NULL);
4418     }
4419   channel_mask=SetImageChannelMask(border_image,AlphaChannel);
4420   shadow_image=BlurImage(border_image,0.0,sigma,exception);
4421   border_image=DestroyImage(border_image);
4422   if (shadow_image == (Image *) NULL)
4423     return((Image *) NULL);
4424   (void) SetPixelChannelMask(shadow_image,channel_mask);
4425   if (shadow_image->page.width == 0)
4426     shadow_image->page.width=shadow_image->columns;
4427   if (shadow_image->page.height == 0)
4428     shadow_image->page.height=shadow_image->rows;
4429   shadow_image->page.width+=x_offset-(ssize_t) border_info.width;
4430   shadow_image->page.height+=y_offset-(ssize_t) border_info.height;
4431   shadow_image->page.x+=x_offset-(ssize_t) border_info.width;
4432   shadow_image->page.y+=y_offset-(ssize_t) border_info.height;
4433   return(shadow_image);
4434 }
4435 \f
4436 /*
4437 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4438 %                                                                             %
4439 %                                                                             %
4440 %                                                                             %
4441 %     S k e t c h I m a g e                                                   %
4442 %                                                                             %
4443 %                                                                             %
4444 %                                                                             %
4445 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4446 %
4447 %  SketchImage() simulates a pencil sketch.  We convolve the image with a
4448 %  Gaussian operator of the given radius and standard deviation (sigma).  For
4449 %  reasonable results, radius should be larger than sigma.  Use a radius of 0
4450 %  and SketchImage() selects a suitable radius for you.  Angle gives the angle
4451 %  of the sketch.
4452 %
4453 %  The format of the SketchImage method is:
4454 %
4455 %    Image *SketchImage(const Image *image,const double radius,
4456 %      const double sigma,const double angle,ExceptionInfo *exception)
4457 %
4458 %  A description of each parameter follows:
4459 %
4460 %    o image: the image.
4461 %
4462 %    o radius: the radius of the Gaussian, in pixels, not counting the
4463 %      center pixel.
4464 %
4465 %    o sigma: the standard deviation of the Gaussian, in pixels.
4466 %
4467 %    o angle: apply the effect along this angle.
4468 %
4469 %    o exception: return any errors or warnings in this structure.
4470 %
4471 */
4472 MagickExport Image *SketchImage(const Image *image,const double radius,
4473   const double sigma,const double angle,ExceptionInfo *exception)
4474 {
4475   CacheView
4476     *random_view;
4477
4478   Image
4479     *blend_image,
4480     *blur_image,
4481     *dodge_image,
4482     *random_image,
4483     *sketch_image;
4484
4485   MagickBooleanType
4486     status;
4487
4488   RandomInfo
4489     **magick_restrict random_info;
4490
4491   ssize_t
4492     y;
4493
4494 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4495   unsigned long
4496     key;
4497 #endif
4498
4499   /*
4500     Sketch image.
4501   */
4502   random_image=CloneImage(image,image->columns << 1,image->rows << 1,
4503     MagickTrue,exception);
4504   if (random_image == (Image *) NULL)
4505     return((Image *) NULL);
4506   status=MagickTrue;
4507   random_info=AcquireRandomInfoThreadSet();
4508   random_view=AcquireAuthenticCacheView(random_image,exception);
4509 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4510   key=GetRandomSecretKey(random_info[0]);
4511   #pragma omp parallel for schedule(static,4) shared(status) \
4512     magick_threads(random_image,random_image,random_image->rows,key == ~0UL)
4513 #endif
4514   for (y=0; y < (ssize_t) random_image->rows; y++)
4515   {
4516     const int
4517       id = GetOpenMPThreadId();
4518
4519     register Quantum
4520       *magick_restrict q;
4521
4522     register ssize_t
4523       x;
4524
4525     if (status == MagickFalse)
4526       continue;
4527     q=QueueCacheViewAuthenticPixels(random_view,0,y,random_image->columns,1,
4528       exception);
4529     if (q == (Quantum *) NULL)
4530       {
4531         status=MagickFalse;
4532         continue;
4533       }
4534     for (x=0; x < (ssize_t) random_image->columns; x++)
4535     {
4536       double
4537         value;
4538
4539       register ssize_t
4540         i;
4541
4542       if (GetPixelWriteMask(random_image,q) == 0)
4543         {
4544           q+=GetPixelChannels(random_image);
4545           continue;
4546         }
4547       value=GetPseudoRandomValue(random_info[id]);
4548       for (i=0; i < (ssize_t) GetPixelChannels(random_image); i++)
4549       {
4550         PixelChannel channel=GetPixelChannelChannel(image,i);
4551         PixelTrait traits=GetPixelChannelTraits(image,channel);
4552         if (traits == UndefinedPixelTrait)
4553           continue;
4554         q[i]=ClampToQuantum(QuantumRange*value);
4555       }
4556       q+=GetPixelChannels(random_image);
4557     }
4558     if (SyncCacheViewAuthenticPixels(random_view,exception) == MagickFalse)
4559       status=MagickFalse;
4560   }
4561   random_view=DestroyCacheView(random_view);
4562   random_info=DestroyRandomInfoThreadSet(random_info);
4563   if (status == MagickFalse)
4564     {
4565       random_image=DestroyImage(random_image);
4566       return(random_image);
4567     }
4568   blur_image=MotionBlurImage(random_image,radius,sigma,angle,exception);
4569   random_image=DestroyImage(random_image);
4570   if (blur_image == (Image *) NULL)
4571     return((Image *) NULL);
4572   dodge_image=EdgeImage(blur_image,radius,exception);
4573   blur_image=DestroyImage(blur_image);
4574   if (dodge_image == (Image *) NULL)
4575     return((Image *) NULL);
4576   (void) NormalizeImage(dodge_image,exception);
4577   (void) NegateImage(dodge_image,MagickFalse,exception);
4578   (void) TransformImage(&dodge_image,(char *) NULL,"50%",exception);
4579   sketch_image=CloneImage(image,0,0,MagickTrue,exception);
4580   if (sketch_image == (Image *) NULL)
4581     {
4582       dodge_image=DestroyImage(dodge_image);
4583       return((Image *) NULL);
4584     }
4585   (void) CompositeImage(sketch_image,dodge_image,ColorDodgeCompositeOp,
4586     MagickTrue,0,0,exception);
4587   dodge_image=DestroyImage(dodge_image);
4588   blend_image=CloneImage(image,0,0,MagickTrue,exception);
4589   if (blend_image == (Image *) NULL)
4590     {
4591       sketch_image=DestroyImage(sketch_image);
4592       return((Image *) NULL);
4593     }
4594   if (blend_image->alpha_trait != BlendPixelTrait)
4595     (void) SetImageAlpha(blend_image,TransparentAlpha,exception);
4596   (void) SetImageArtifact(blend_image,"compose:args","20x80");
4597   (void) CompositeImage(sketch_image,blend_image,BlendCompositeOp,MagickTrue,
4598     0,0,exception);
4599   blend_image=DestroyImage(blend_image);
4600   return(sketch_image);
4601 }
4602 \f
4603 /*
4604 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4605 %                                                                             %
4606 %                                                                             %
4607 %                                                                             %
4608 %     S o l a r i z e I m a g e                                               %
4609 %                                                                             %
4610 %                                                                             %
4611 %                                                                             %
4612 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4613 %
4614 %  SolarizeImage() applies a special effect to the image, similar to the effect
4615 %  achieved in a photo darkroom by selectively exposing areas of photo
4616 %  sensitive paper to light.  Threshold ranges from 0 to QuantumRange and is a
4617 %  measure of the extent of the solarization.
4618 %
4619 %  The format of the SolarizeImage method is:
4620 %
4621 %      MagickBooleanType SolarizeImage(Image *image,const double threshold,
4622 %        ExceptionInfo *exception)
4623 %
4624 %  A description of each parameter follows:
4625 %
4626 %    o image: the image.
4627 %
4628 %    o threshold:  Define the extent of the solarization.
4629 %
4630 %    o exception: return any errors or warnings in this structure.
4631 %
4632 */
4633 MagickExport MagickBooleanType SolarizeImage(Image *image,
4634   const double threshold,ExceptionInfo *exception)
4635 {
4636 #define SolarizeImageTag  "Solarize/Image"
4637
4638   CacheView
4639     *image_view;
4640
4641   MagickBooleanType
4642     status;
4643
4644   MagickOffsetType
4645     progress;
4646
4647   ssize_t
4648     y;
4649
4650   assert(image != (Image *) NULL);
4651   assert(image->signature == MagickCoreSignature);
4652   if (image->debug != MagickFalse)
4653     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4654   if (IsGrayColorspace(image->colorspace) != MagickFalse)
4655     (void) SetImageColorspace(image,sRGBColorspace,exception);
4656   if (image->storage_class == PseudoClass)
4657     {
4658       register ssize_t
4659         i;
4660
4661       /*
4662         Solarize colormap.
4663       */
4664       for (i=0; i < (ssize_t) image->colors; i++)
4665       {
4666         if ((double) image->colormap[i].red > threshold)
4667           image->colormap[i].red=QuantumRange-image->colormap[i].red;
4668         if ((double) image->colormap[i].green > threshold)
4669           image->colormap[i].green=QuantumRange-image->colormap[i].green;
4670         if ((double) image->colormap[i].blue > threshold)
4671           image->colormap[i].blue=QuantumRange-image->colormap[i].blue;
4672       }
4673     }
4674   /*
4675     Solarize image.
4676   */
4677   status=MagickTrue;
4678   progress=0;
4679   image_view=AcquireAuthenticCacheView(image,exception);
4680 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4681   #pragma omp parallel for schedule(static,4) shared(progress,status) \
4682     magick_threads(image,image,image->rows,1)
4683 #endif
4684   for (y=0; y < (ssize_t) image->rows; y++)
4685   {
4686     register ssize_t
4687       x;
4688
4689     register Quantum
4690       *magick_restrict q;
4691
4692     if (status == MagickFalse)
4693       continue;
4694     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
4695     if (q == (Quantum *) NULL)
4696       {
4697         status=MagickFalse;
4698         continue;
4699       }
4700     for (x=0; x < (ssize_t) image->columns; x++)
4701     {
4702       register ssize_t
4703         i;
4704
4705       if (GetPixelWriteMask(image,q) == 0)
4706         {
4707           q+=GetPixelChannels(image);
4708           continue;
4709         }
4710       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
4711       {
4712         PixelChannel channel=GetPixelChannelChannel(image,i);
4713         PixelTrait traits=GetPixelChannelTraits(image,channel);
4714         if ((traits & UpdatePixelTrait) == 0)
4715           continue;
4716         if ((double) q[i] > threshold)
4717           q[i]=QuantumRange-q[i];
4718       }
4719       q+=GetPixelChannels(image);
4720     }
4721     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
4722       status=MagickFalse;
4723     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4724       {
4725         MagickBooleanType
4726           proceed;
4727
4728 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4729         #pragma omp critical (MagickCore_SolarizeImage)
4730 #endif
4731         proceed=SetImageProgress(image,SolarizeImageTag,progress++,image->rows);
4732         if (proceed == MagickFalse)
4733           status=MagickFalse;
4734       }
4735   }
4736   image_view=DestroyCacheView(image_view);
4737   return(status);
4738 }
4739 \f
4740 /*
4741 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4742 %                                                                             %
4743 %                                                                             %
4744 %                                                                             %
4745 %   S t e g a n o I m a g e                                                   %
4746 %                                                                             %
4747 %                                                                             %
4748 %                                                                             %
4749 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4750 %
4751 %  SteganoImage() hides a digital watermark within the image.  Recover
4752 %  the hidden watermark later to prove that the authenticity of an image.
4753 %  Offset defines the start position within the image to hide the watermark.
4754 %
4755 %  The format of the SteganoImage method is:
4756 %
4757 %      Image *SteganoImage(const Image *image,Image *watermark,
4758 %        ExceptionInfo *exception)
4759 %
4760 %  A description of each parameter follows:
4761 %
4762 %    o image: the image.
4763 %
4764 %    o watermark: the watermark image.
4765 %
4766 %    o exception: return any errors or warnings in this structure.
4767 %
4768 */
4769 MagickExport Image *SteganoImage(const Image *image,const Image *watermark,
4770   ExceptionInfo *exception)
4771 {
4772 #define GetBit(alpha,i) ((((size_t) (alpha) >> (size_t) (i)) & 0x01) != 0)
4773 #define SetBit(alpha,i,set) (Quantum) ((set) != 0 ? (size_t) (alpha) \
4774   | (one << (size_t) (i)) : (size_t) (alpha) & ~(one << (size_t) (i)))
4775 #define SteganoImageTag  "Stegano/Image"
4776
4777   CacheView
4778     *stegano_view,
4779     *watermark_view;
4780
4781   Image
4782     *stegano_image;
4783
4784   int
4785     c;
4786
4787   MagickBooleanType
4788     status;
4789
4790   PixelInfo
4791     pixel;
4792
4793   register Quantum
4794     *q;
4795
4796   register ssize_t
4797     x;
4798
4799   size_t
4800     depth,
4801     one;
4802
4803   ssize_t
4804     i,
4805     j,
4806     k,
4807     y;
4808
4809   /*
4810     Initialize steganographic image attributes.
4811   */
4812   assert(image != (const Image *) NULL);
4813   assert(image->signature == MagickCoreSignature);
4814   if (image->debug != MagickFalse)
4815     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4816   assert(watermark != (const Image *) NULL);
4817   assert(watermark->signature == MagickCoreSignature);
4818   assert(exception != (ExceptionInfo *) NULL);
4819   assert(exception->signature == MagickCoreSignature);
4820   one=1UL;
4821   stegano_image=CloneImage(image,0,0,MagickTrue,exception);
4822   if (stegano_image == (Image *) NULL)
4823     return((Image *) NULL);
4824   stegano_image->depth=MAGICKCORE_QUANTUM_DEPTH;
4825   if (SetImageStorageClass(stegano_image,DirectClass,exception) == MagickFalse)
4826     {
4827       stegano_image=DestroyImage(stegano_image);
4828       return((Image *) NULL);
4829     }
4830   /*
4831     Hide watermark in low-order bits of image.
4832   */
4833   c=0;
4834   i=0;
4835   j=0;
4836   depth=stegano_image->depth;
4837   k=stegano_image->offset;
4838   status=MagickTrue;
4839   watermark_view=AcquireVirtualCacheView(watermark,exception);
4840   stegano_view=AcquireAuthenticCacheView(stegano_image,exception);
4841   for (i=(ssize_t) depth-1; (i >= 0) && (j < (ssize_t) depth); i--)
4842   {
4843     for (y=0; (y < (ssize_t) watermark->rows) && (j < (ssize_t) depth); y++)
4844     {
4845       for (x=0; (x < (ssize_t) watermark->columns) && (j < (ssize_t) depth); x++)
4846       {
4847         ssize_t
4848           offset;
4849
4850         (void) GetOneCacheViewVirtualPixelInfo(watermark_view,x,y,&pixel,
4851           exception);
4852         offset=k/(ssize_t) stegano_image->columns;
4853         if (offset >= (ssize_t) stegano_image->rows)
4854           break;
4855         q=GetCacheViewAuthenticPixels(stegano_view,k % (ssize_t)
4856           stegano_image->columns,k/(ssize_t) stegano_image->columns,1,1,
4857           exception);
4858         if (q == (Quantum *) NULL)
4859           break;
4860         switch (c)
4861         {
4862           case 0:
4863           {
4864             SetPixelRed(stegano_image,SetBit(GetPixelRed(stegano_image,q),j,
4865               GetBit(GetPixelInfoIntensity(stegano_image,&pixel),i)),q);
4866             break;
4867           }
4868           case 1:
4869           {
4870             SetPixelGreen(stegano_image,SetBit(GetPixelGreen(stegano_image,q),j,
4871               GetBit(GetPixelInfoIntensity(stegano_image,&pixel),i)),q);
4872             break;
4873           }
4874           case 2:
4875           {
4876             SetPixelBlue(stegano_image,SetBit(GetPixelBlue(stegano_image,q),j,
4877               GetBit(GetPixelInfoIntensity(stegano_image,&pixel),i)),q);
4878             break;
4879           }
4880         }
4881         if (SyncCacheViewAuthenticPixels(stegano_view,exception) == MagickFalse)
4882           break;
4883         c++;
4884         if (c == 3)
4885           c=0;
4886         k++;
4887         if (k == (ssize_t) (stegano_image->columns*stegano_image->columns))
4888           k=0;
4889         if (k == stegano_image->offset)
4890           j++;
4891       }
4892     }
4893     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4894       {
4895         MagickBooleanType
4896           proceed;
4897
4898         proceed=SetImageProgress(image,SteganoImageTag,(MagickOffsetType)
4899           (depth-i),depth);
4900         if (proceed == MagickFalse)
4901           status=MagickFalse;
4902       }
4903   }
4904   stegano_view=DestroyCacheView(stegano_view);
4905   watermark_view=DestroyCacheView(watermark_view);
4906   if (status == MagickFalse)
4907     stegano_image=DestroyImage(stegano_image);
4908   return(stegano_image);
4909 }
4910 \f
4911 /*
4912 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4913 %                                                                             %
4914 %                                                                             %
4915 %                                                                             %
4916 %   S t e r e o A n a g l y p h I m a g e                                     %
4917 %                                                                             %
4918 %                                                                             %
4919 %                                                                             %
4920 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4921 %
4922 %  StereoAnaglyphImage() combines two images and produces a single image that
4923 %  is the composite of a left and right image of a stereo pair.  Special
4924 %  red-green stereo glasses are required to view this effect.
4925 %
4926 %  The format of the StereoAnaglyphImage method is:
4927 %
4928 %      Image *StereoImage(const Image *left_image,const Image *right_image,
4929 %        ExceptionInfo *exception)
4930 %      Image *StereoAnaglyphImage(const Image *left_image,
4931 %        const Image *right_image,const ssize_t x_offset,const ssize_t y_offset,
4932 %        ExceptionInfo *exception)
4933 %
4934 %  A description of each parameter follows:
4935 %
4936 %    o left_image: the left image.
4937 %
4938 %    o right_image: the right image.
4939 %
4940 %    o exception: return any errors or warnings in this structure.
4941 %
4942 %    o x_offset: amount, in pixels, by which the left image is offset to the
4943 %      right of the right image.
4944 %
4945 %    o y_offset: amount, in pixels, by which the left image is offset to the
4946 %      bottom of the right image.
4947 %
4948 %
4949 */
4950 MagickExport Image *StereoImage(const Image *left_image,
4951   const Image *right_image,ExceptionInfo *exception)
4952 {
4953   return(StereoAnaglyphImage(left_image,right_image,0,0,exception));
4954 }
4955
4956 MagickExport Image *StereoAnaglyphImage(const Image *left_image,
4957   const Image *right_image,const ssize_t x_offset,const ssize_t y_offset,
4958   ExceptionInfo *exception)
4959 {
4960 #define StereoImageTag  "Stereo/Image"
4961
4962   const Image
4963     *image;
4964
4965   Image
4966     *stereo_image;
4967
4968   MagickBooleanType
4969     status;
4970
4971   ssize_t
4972     y;
4973
4974   assert(left_image != (const Image *) NULL);
4975   assert(left_image->signature == MagickCoreSignature);
4976   if (left_image->debug != MagickFalse)
4977     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
4978       left_image->filename);
4979   assert(right_image != (const Image *) NULL);
4980   assert(right_image->signature == MagickCoreSignature);
4981   assert(exception != (ExceptionInfo *) NULL);
4982   assert(exception->signature == MagickCoreSignature);
4983   assert(right_image != (const Image *) NULL);
4984   image=left_image;
4985   if ((left_image->columns != right_image->columns) ||
4986       (left_image->rows != right_image->rows))
4987     ThrowImageException(ImageError,"LeftAndRightImageSizesDiffer");
4988   /*
4989     Initialize stereo image attributes.
4990   */
4991   stereo_image=CloneImage(left_image,left_image->columns,left_image->rows,
4992     MagickTrue,exception);
4993   if (stereo_image == (Image *) NULL)
4994     return((Image *) NULL);
4995   if (SetImageStorageClass(stereo_image,DirectClass,exception) == MagickFalse)
4996     {
4997       stereo_image=DestroyImage(stereo_image);
4998       return((Image *) NULL);
4999     }
5000   (void) SetImageColorspace(stereo_image,sRGBColorspace,exception);
5001   /*
5002     Copy left image to red channel and right image to blue channel.
5003   */
5004   status=MagickTrue;
5005   for (y=0; y < (ssize_t) stereo_image->rows; y++)
5006   {
5007     register const Quantum
5008       *magick_restrict p,
5009       *magick_restrict q;
5010
5011     register ssize_t
5012       x;
5013
5014     register Quantum
5015       *magick_restrict r;
5016
5017     p=GetVirtualPixels(left_image,-x_offset,y-y_offset,image->columns,1,
5018       exception);
5019     q=GetVirtualPixels(right_image,0,y,right_image->columns,1,exception);
5020     r=QueueAuthenticPixels(stereo_image,0,y,stereo_image->columns,1,exception);
5021     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL) ||
5022         (r == (Quantum *) NULL))
5023       break;
5024     for (x=0; x < (ssize_t) stereo_image->columns; x++)
5025     {
5026       SetPixelRed(image,GetPixelRed(left_image,p),r);
5027       SetPixelGreen(image,GetPixelGreen(right_image,q),r);
5028       SetPixelBlue(image,GetPixelBlue(right_image,q),r);
5029       if ((GetPixelAlphaTraits(stereo_image) & CopyPixelTrait) != 0)
5030         SetPixelAlpha(image,(GetPixelAlpha(left_image,p)+
5031           GetPixelAlpha(right_image,q))/2,r);
5032       p+=GetPixelChannels(left_image);
5033       q+=GetPixelChannels(right_image);
5034       r+=GetPixelChannels(stereo_image);
5035     }
5036     if (SyncAuthenticPixels(stereo_image,exception) == MagickFalse)
5037       break;
5038     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5039       {
5040         MagickBooleanType
5041           proceed;
5042
5043         proceed=SetImageProgress(image,StereoImageTag,(MagickOffsetType) y,
5044           stereo_image->rows);
5045         if (proceed == MagickFalse)
5046           status=MagickFalse;
5047       }
5048   }
5049   if (status == MagickFalse)
5050     stereo_image=DestroyImage(stereo_image);
5051   return(stereo_image);
5052 }
5053 \f
5054 /*
5055 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5056 %                                                                             %
5057 %                                                                             %
5058 %                                                                             %
5059 %     S w i r l I m a g e                                                     %
5060 %                                                                             %
5061 %                                                                             %
5062 %                                                                             %
5063 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5064 %
5065 %  SwirlImage() swirls the pixels about the center of the image, where
5066 %  degrees indicates the sweep of the arc through which each pixel is moved.
5067 %  You get a more dramatic effect as the degrees move from 1 to 360.
5068 %
5069 %  The format of the SwirlImage method is:
5070 %
5071 %      Image *SwirlImage(const Image *image,double degrees,
5072 %        const PixelInterpolateMethod method,ExceptionInfo *exception)
5073 %
5074 %  A description of each parameter follows:
5075 %
5076 %    o image: the image.
5077 %
5078 %    o degrees: Define the tightness of the swirling effect.
5079 %
5080 %    o method: the pixel interpolation method.
5081 %
5082 %    o exception: return any errors or warnings in this structure.
5083 %
5084 */
5085 MagickExport Image *SwirlImage(const Image *image,double degrees,
5086   const PixelInterpolateMethod method,ExceptionInfo *exception)
5087 {
5088 #define SwirlImageTag  "Swirl/Image"
5089
5090   CacheView
5091     *canvas_view,
5092     *interpolate_view,
5093     *swirl_view;
5094
5095   Image
5096     *canvas,
5097     *swirl_image;
5098
5099   MagickBooleanType
5100     status;
5101
5102   MagickOffsetType
5103     progress;
5104
5105   double
5106     radius;
5107
5108   PointInfo
5109     center,
5110     scale;
5111
5112   ssize_t
5113     y;
5114
5115   /*
5116     Initialize swirl image attributes.
5117   */
5118   assert(image != (const Image *) NULL);
5119   assert(image->signature == MagickCoreSignature);
5120   if (image->debug != MagickFalse)
5121     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5122   assert(exception != (ExceptionInfo *) NULL);
5123   assert(exception->signature == MagickCoreSignature);
5124   canvas=CloneImage(image,0,0,MagickTrue,exception);
5125   if (canvas == (Image *) NULL)
5126     return((Image *) NULL);
5127   if ((canvas->alpha_trait == UndefinedPixelTrait) &&
5128       (canvas->background_color.alpha != OpaqueAlpha))
5129     (void) SetImageAlphaChannel(canvas,OpaqueAlphaChannel,exception);
5130   swirl_image=CloneImage(canvas,canvas->columns,canvas->rows,MagickTrue,
5131     exception);
5132   if (swirl_image == (Image *) NULL)
5133     {
5134       canvas=DestroyImage(canvas);
5135       return((Image *) NULL);
5136     }
5137   if (SetImageStorageClass(swirl_image,DirectClass,exception) == MagickFalse)
5138     {
5139       canvas=DestroyImage(canvas);
5140       swirl_image=DestroyImage(swirl_image);
5141       return((Image *) NULL);
5142     }
5143   /*
5144     Compute scaling factor.
5145   */
5146   center.x=(double) canvas->columns/2.0;
5147   center.y=(double) canvas->rows/2.0;
5148   radius=MagickMax(center.x,center.y);
5149   scale.x=1.0;
5150   scale.y=1.0;
5151   if (canvas->columns > canvas->rows)
5152     scale.y=(double) canvas->columns/(double) canvas->rows;
5153   else
5154     if (canvas->columns < canvas->rows)
5155       scale.x=(double) canvas->rows/(double) canvas->columns;
5156   degrees=(double) DegreesToRadians(degrees);
5157   /*
5158     Swirl image.
5159   */
5160   status=MagickTrue;
5161   progress=0;
5162   canvas_view=AcquireVirtualCacheView(canvas,exception);
5163   interpolate_view=AcquireVirtualCacheView(image,exception);
5164   swirl_view=AcquireAuthenticCacheView(swirl_image,exception);
5165 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5166   #pragma omp parallel for schedule(static,4) shared(progress,status) \
5167     magick_threads(canvas,swirl_image,canvas->rows,1)
5168 #endif
5169   for (y=0; y < (ssize_t) canvas->rows; y++)
5170   {
5171     double
5172       distance;
5173
5174     PointInfo
5175       delta;
5176
5177     register const Quantum
5178       *magick_restrict p;
5179
5180     register ssize_t
5181       x;
5182
5183     register Quantum
5184       *magick_restrict q;
5185
5186     if (status == MagickFalse)
5187       continue;
5188     p=GetCacheViewVirtualPixels(canvas_view,0,y,canvas->columns,1,exception);
5189     q=QueueCacheViewAuthenticPixels(swirl_view,0,y,swirl_image->columns,1,
5190       exception);
5191     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
5192       {
5193         status=MagickFalse;
5194         continue;
5195       }
5196     delta.y=scale.y*(double) (y-center.y);
5197     for (x=0; x < (ssize_t) canvas->columns; x++)
5198     {
5199       /*
5200         Determine if the pixel is within an ellipse.
5201       */
5202       if (GetPixelWriteMask(canvas,p) == 0)
5203         {
5204           SetPixelBackgoundColor(swirl_image,q);
5205           p+=GetPixelChannels(canvas);
5206           q+=GetPixelChannels(swirl_image);
5207           continue;
5208         }
5209       delta.x=scale.x*(double) (x-center.x);
5210       distance=delta.x*delta.x+delta.y*delta.y;
5211       if (distance >= (radius*radius))
5212         {
5213           register ssize_t
5214             i;
5215
5216           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
5217           {
5218             PixelChannel channel = GetPixelChannelChannel(canvas,i);
5219             PixelTrait traits = GetPixelChannelTraits(canvas,channel);
5220             PixelTrait swirl_traits = GetPixelChannelTraits(swirl_image,
5221               channel);
5222             if ((traits == UndefinedPixelTrait) ||
5223                 (swirl_traits == UndefinedPixelTrait))
5224               continue;
5225             SetPixelChannel(swirl_image,channel,p[i],q);
5226           }
5227         }
5228       else
5229         {
5230           double
5231             cosine,
5232             factor,
5233             sine;
5234
5235           /*
5236             Swirl the pixel.
5237           */
5238           factor=1.0-sqrt((double) distance)/radius;
5239           sine=sin((double) (degrees*factor*factor));
5240           cosine=cos((double) (degrees*factor*factor));
5241           status=InterpolatePixelChannels(canvas,interpolate_view,swirl_image,
5242             method,((cosine*delta.x-sine*delta.y)/scale.x+center.x),(double)
5243             ((sine*delta.x+cosine*delta.y)/scale.y+center.y),q,exception);
5244         }
5245       p+=GetPixelChannels(canvas);
5246       q+=GetPixelChannels(swirl_image);
5247     }
5248     if (SyncCacheViewAuthenticPixels(swirl_view,exception) == MagickFalse)
5249       status=MagickFalse;
5250     if (canvas->progress_monitor != (MagickProgressMonitor) NULL)
5251       {
5252         MagickBooleanType
5253           proceed;
5254
5255 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5256         #pragma omp critical (MagickCore_SwirlImage)
5257 #endif
5258         proceed=SetImageProgress(canvas,SwirlImageTag,progress++,canvas->rows);
5259         if (proceed == MagickFalse)
5260           status=MagickFalse;
5261       }
5262   }
5263   swirl_view=DestroyCacheView(swirl_view);
5264   interpolate_view=DestroyCacheView(interpolate_view);
5265   canvas_view=DestroyCacheView(canvas_view);
5266   canvas=DestroyImage(canvas);
5267   if (status == MagickFalse)
5268     swirl_image=DestroyImage(swirl_image);
5269   return(swirl_image);
5270 }
5271 \f
5272 /*
5273 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5274 %                                                                             %
5275 %                                                                             %
5276 %                                                                             %
5277 %     T i n t I m a g e                                                       %
5278 %                                                                             %
5279 %                                                                             %
5280 %                                                                             %
5281 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5282 %
5283 %  TintImage() applies a color vector to each pixel in the image.  The length
5284 %  of the vector is 0 for black and white and at its maximum for the midtones.
5285 %  The vector weighting function is f(x)=(1-(4.0*((x-0.5)*(x-0.5))))
5286 %
5287 %  The format of the TintImage method is:
5288 %
5289 %      Image *TintImage(const Image *image,const char *blend,
5290 %        const PixelInfo *tint,ExceptionInfo *exception)
5291 %
5292 %  A description of each parameter follows:
5293 %
5294 %    o image: the image.
5295 %
5296 %    o blend: A color value used for tinting.
5297 %
5298 %    o tint: A color value used for tinting.
5299 %
5300 %    o exception: return any errors or warnings in this structure.
5301 %
5302 */
5303 MagickExport Image *TintImage(const Image *image,const char *blend,
5304   const PixelInfo *tint,ExceptionInfo *exception)
5305 {
5306 #define TintImageTag  "Tint/Image"
5307
5308   CacheView
5309     *image_view,
5310     *tint_view;
5311
5312   double
5313     intensity;
5314
5315   GeometryInfo
5316     geometry_info;
5317
5318   Image
5319     *tint_image;
5320
5321   MagickBooleanType
5322     status;
5323
5324   MagickOffsetType
5325     progress;
5326
5327   PixelInfo
5328     color_vector;
5329
5330   MagickStatusType
5331     flags;
5332
5333   ssize_t
5334     y;
5335
5336   /*
5337     Allocate tint image.
5338   */
5339   assert(image != (const Image *) NULL);
5340   assert(image->signature == MagickCoreSignature);
5341   if (image->debug != MagickFalse)
5342     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5343   assert(exception != (ExceptionInfo *) NULL);
5344   assert(exception->signature == MagickCoreSignature);
5345   tint_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
5346   if (tint_image == (Image *) NULL)
5347     return((Image *) NULL);
5348   if (SetImageStorageClass(tint_image,DirectClass,exception) == MagickFalse)
5349     {
5350       tint_image=DestroyImage(tint_image);
5351       return((Image *) NULL);
5352     }
5353   if ((IsGrayColorspace(image->colorspace) != MagickFalse) &&
5354       (IsPixelInfoGray(tint) == MagickFalse))
5355     (void) SetImageColorspace(tint_image,sRGBColorspace,exception);
5356   if (blend == (const char *) NULL)
5357     return(tint_image);
5358   /*
5359     Determine RGB values of the color.
5360   */
5361   GetPixelInfo(image,&color_vector);
5362   flags=ParseGeometry(blend,&geometry_info);
5363   color_vector.red=geometry_info.rho;
5364   color_vector.green=geometry_info.rho;
5365   color_vector.blue=geometry_info.rho;
5366   color_vector.alpha=(MagickRealType) OpaqueAlpha;
5367   if ((flags & SigmaValue) != 0)
5368     color_vector.green=geometry_info.sigma;
5369   if ((flags & XiValue) != 0)
5370     color_vector.blue=geometry_info.xi;
5371   if ((flags & PsiValue) != 0)
5372     color_vector.alpha=geometry_info.psi;
5373   if (image->colorspace == CMYKColorspace)
5374     {
5375       color_vector.black=geometry_info.rho;
5376       if ((flags & PsiValue) != 0)
5377         color_vector.black=geometry_info.psi;
5378       if ((flags & ChiValue) != 0)
5379         color_vector.alpha=geometry_info.chi;
5380     }
5381   intensity=(double) GetPixelInfoIntensity((const Image *) NULL,tint);
5382   color_vector.red=(double) (color_vector.red*tint->red/100.0-intensity);
5383   color_vector.green=(double) (color_vector.green*tint->green/100.0-intensity);
5384   color_vector.blue=(double) (color_vector.blue*tint->blue/100.0-intensity);
5385   color_vector.black=(double) (color_vector.black*tint->black/100.0-intensity);
5386   color_vector.alpha=(double) (color_vector.alpha*tint->alpha/100.0-intensity);
5387   /*
5388     Tint image.
5389   */
5390   status=MagickTrue;
5391   progress=0;
5392   image_view=AcquireVirtualCacheView(image,exception);
5393   tint_view=AcquireAuthenticCacheView(tint_image,exception);
5394 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5395   #pragma omp parallel for schedule(static,4) shared(progress,status) \
5396     magick_threads(image,tint_image,image->rows,1)
5397 #endif
5398   for (y=0; y < (ssize_t) image->rows; y++)
5399   {
5400     register const Quantum
5401       *magick_restrict p;
5402
5403     register Quantum
5404       *magick_restrict q;
5405
5406     register ssize_t
5407       x;
5408
5409     if (status == MagickFalse)
5410       continue;
5411     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
5412     q=QueueCacheViewAuthenticPixels(tint_view,0,y,tint_image->columns,1,
5413       exception);
5414     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
5415       {
5416         status=MagickFalse;
5417         continue;
5418       }
5419     for (x=0; x < (ssize_t) image->columns; x++)
5420     {
5421       PixelInfo
5422         pixel;
5423
5424       double
5425         weight;
5426
5427       register ssize_t
5428         i;
5429
5430       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
5431       {
5432         PixelChannel channel=GetPixelChannelChannel(image,i);
5433         PixelTrait traits=GetPixelChannelTraits(image,channel);
5434         PixelTrait tint_traits=GetPixelChannelTraits(tint_image,channel);
5435         if ((traits == UndefinedPixelTrait) ||
5436             (tint_traits == UndefinedPixelTrait))
5437           continue;
5438         if (((tint_traits & CopyPixelTrait) != 0) ||
5439             (GetPixelWriteMask(image,p) == 0))
5440           {
5441             SetPixelChannel(tint_image,channel,p[i],q);
5442             continue;
5443           }
5444       }
5445       GetPixelInfo(image,&pixel);
5446       weight=QuantumScale*GetPixelRed(image,p)-0.5;
5447       pixel.red=(double) GetPixelRed(image,p)+color_vector.red*(1.0-(4.0*
5448         (weight*weight)));
5449       weight=QuantumScale*GetPixelGreen(image,p)-0.5;
5450       pixel.green=(double) GetPixelGreen(image,p)+color_vector.green*(1.0-(4.0*
5451         (weight*weight)));
5452       weight=QuantumScale*GetPixelBlue(image,p)-0.5;
5453       pixel.blue=(double) GetPixelBlue(image,p)+color_vector.blue*(1.0-(4.0*
5454         (weight*weight)));
5455       weight=QuantumScale*GetPixelBlack(image,p)-0.5;
5456       pixel.black=(double) GetPixelBlack(image,p)+color_vector.black*(1.0-(4.0*
5457         (weight*weight)));
5458       SetPixelViaPixelInfo(tint_image,&pixel,q);
5459       p+=GetPixelChannels(image);
5460       q+=GetPixelChannels(tint_image);
5461     }
5462     if (SyncCacheViewAuthenticPixels(tint_view,exception) == MagickFalse)
5463       status=MagickFalse;
5464     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5465       {
5466         MagickBooleanType
5467           proceed;
5468
5469 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5470         #pragma omp critical (MagickCore_TintImage)
5471 #endif
5472         proceed=SetImageProgress(image,TintImageTag,progress++,image->rows);
5473         if (proceed == MagickFalse)
5474           status=MagickFalse;
5475       }
5476   }
5477   tint_view=DestroyCacheView(tint_view);
5478   image_view=DestroyCacheView(image_view);
5479   if (status == MagickFalse)
5480     tint_image=DestroyImage(tint_image);
5481   return(tint_image);
5482 }
5483 \f
5484 /*
5485 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5486 %                                                                             %
5487 %                                                                             %
5488 %                                                                             %
5489 %     V i g n e t t e I m a g e                                               %
5490 %                                                                             %
5491 %                                                                             %
5492 %                                                                             %
5493 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5494 %
5495 %  VignetteImage() softens the edges of the image in vignette style.
5496 %
5497 %  The format of the VignetteImage method is:
5498 %
5499 %      Image *VignetteImage(const Image *image,const double radius,
5500 %        const double sigma,const ssize_t x,const ssize_t y,
5501 %        ExceptionInfo *exception)
5502 %
5503 %  A description of each parameter follows:
5504 %
5505 %    o image: the image.
5506 %
5507 %    o radius: the radius of the pixel neighborhood.
5508 %
5509 %    o sigma: the standard deviation of the Gaussian, in pixels.
5510 %
5511 %    o x, y:  Define the x and y ellipse offset.
5512 %
5513 %    o exception: return any errors or warnings in this structure.
5514 %
5515 */
5516 MagickExport Image *VignetteImage(const Image *image,const double radius,
5517   const double sigma,const ssize_t x,const ssize_t y,ExceptionInfo *exception)
5518 {
5519   char
5520     ellipse[MagickPathExtent];
5521
5522   DrawInfo
5523     *draw_info;
5524
5525   Image
5526     *canvas,
5527     *blur_image,
5528     *oval_image,
5529     *vignette_image;
5530
5531   assert(image != (Image *) NULL);
5532   assert(image->signature == MagickCoreSignature);
5533   if (image->debug != MagickFalse)
5534     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5535   assert(exception != (ExceptionInfo *) NULL);
5536   assert(exception->signature == MagickCoreSignature);
5537   canvas=CloneImage(image,0,0,MagickTrue,exception);
5538   if (canvas == (Image *) NULL)
5539     return((Image *) NULL);
5540   if ((canvas->alpha_trait == UndefinedPixelTrait) &&
5541       (canvas->background_color.alpha != OpaqueAlpha))
5542     (void) SetImageAlphaChannel(canvas,OpaqueAlphaChannel,exception);
5543   if (SetImageStorageClass(canvas,DirectClass,exception) == MagickFalse)
5544     {
5545       canvas=DestroyImage(canvas);
5546       return((Image *) NULL);
5547     }
5548   oval_image=CloneImage(canvas,canvas->columns,canvas->rows,MagickTrue,
5549     exception);
5550   if (oval_image == (Image *) NULL)
5551     {
5552       canvas=DestroyImage(canvas);
5553       return((Image *) NULL);
5554     }
5555   (void) QueryColorCompliance("#000000",AllCompliance,
5556     &oval_image->background_color,exception);
5557   (void) SetImageBackgroundColor(oval_image,exception);
5558   draw_info=CloneDrawInfo((const ImageInfo *) NULL,(const DrawInfo *) NULL);
5559   (void) QueryColorCompliance("#ffffff",AllCompliance,&draw_info->fill,
5560     exception);
5561   (void) QueryColorCompliance("#ffffff",AllCompliance,&draw_info->stroke,
5562     exception);
5563   (void) FormatLocaleString(ellipse,MagickPathExtent,"ellipse %g,%g,%g,%g,"
5564     "0.0,360.0",image->columns/2.0,image->rows/2.0,image->columns/2.0-x,
5565     image->rows/2.0-y);
5566   draw_info->primitive=AcquireString(ellipse);
5567   (void) DrawImage(oval_image,draw_info,exception);
5568   draw_info=DestroyDrawInfo(draw_info);
5569   blur_image=BlurImage(oval_image,radius,sigma,exception);
5570   oval_image=DestroyImage(oval_image);
5571   if (blur_image == (Image *) NULL)
5572     {
5573       canvas=DestroyImage(canvas);
5574       return((Image *) NULL);
5575     }
5576   blur_image->alpha_trait=UndefinedPixelTrait;
5577   (void) CompositeImage(canvas,blur_image,IntensityCompositeOp,MagickTrue,
5578     0,0,exception);
5579   blur_image=DestroyImage(blur_image);
5580   vignette_image=MergeImageLayers(canvas,FlattenLayer,exception);
5581   canvas=DestroyImage(canvas);
5582   if (vignette_image != (Image *) NULL)
5583     (void) TransformImageColorspace(vignette_image,image->colorspace,exception);
5584   return(vignette_image);
5585 }
5586 \f
5587 /*
5588 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5589 %                                                                             %
5590 %                                                                             %
5591 %                                                                             %
5592 %     W a v e I m a g e                                                       %
5593 %                                                                             %
5594 %                                                                             %
5595 %                                                                             %
5596 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5597 %
5598 %  WaveImage() creates a "ripple" effect in the image by shifting the pixels
5599 %  vertically along a sine wave whose amplitude and wavelength is specified
5600 %  by the given parameters.
5601 %
5602 %  The format of the WaveImage method is:
5603 %
5604 %      Image *WaveImage(const Image *image,const double amplitude,
5605 %        const double wave_length,const PixelInterpolateMethod method,
5606 %        ExceptionInfo *exception)
5607 %
5608 %  A description of each parameter follows:
5609 %
5610 %    o image: the image.
5611 %
5612 %    o amplitude, wave_length:  Define the amplitude and wave length of the
5613 %      sine wave.
5614 %
5615 %    o interpolate: the pixel interpolation method.
5616 %
5617 %    o exception: return any errors or warnings in this structure.
5618 %
5619 */
5620 MagickExport Image *WaveImage(const Image *image,const double amplitude,
5621   const double wave_length,const PixelInterpolateMethod method,
5622   ExceptionInfo *exception)
5623 {
5624 #define WaveImageTag  "Wave/Image"
5625
5626   CacheView
5627     *canvas_view,
5628     *wave_view;
5629
5630   Image
5631     *canvas,
5632     *wave_image;
5633
5634   MagickBooleanType
5635     status;
5636
5637   MagickOffsetType
5638     progress;
5639
5640   double
5641     *sine_map;
5642
5643   register ssize_t
5644     i;
5645
5646   ssize_t
5647     y;
5648
5649   /*
5650     Initialize wave image attributes.
5651   */
5652   assert(image != (Image *) NULL);
5653   assert(image->signature == MagickCoreSignature);
5654   if (image->debug != MagickFalse)
5655     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5656   assert(exception != (ExceptionInfo *) NULL);
5657   assert(exception->signature == MagickCoreSignature);
5658   canvas=CloneImage(image,0,0,MagickTrue,exception);
5659   if (canvas == (Image *) NULL)
5660     return((Image *) NULL);
5661   if ((canvas->alpha_trait == UndefinedPixelTrait) &&
5662       (canvas->background_color.alpha != OpaqueAlpha))
5663     (void) SetImageAlpha(canvas,OpaqueAlpha,exception);
5664   wave_image=CloneImage(canvas,canvas->columns,(size_t) (canvas->rows+2.0*
5665     fabs(amplitude)),MagickTrue,exception);
5666   if (wave_image == (Image *) NULL)
5667     {
5668       canvas=DestroyImage(canvas);
5669       return((Image *) NULL);
5670     }
5671   if (SetImageStorageClass(wave_image,DirectClass,exception) == MagickFalse)
5672     {
5673       canvas=DestroyImage(canvas);
5674       wave_image=DestroyImage(wave_image);
5675       return((Image *) NULL);
5676     }
5677   /*
5678     Allocate sine map.
5679   */
5680   sine_map=(double *) AcquireQuantumMemory((size_t) wave_image->columns,
5681     sizeof(*sine_map));
5682   if (sine_map == (double *) NULL)
5683     {
5684       canvas=DestroyImage(canvas);
5685       wave_image=DestroyImage(wave_image);
5686       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
5687     }
5688   for (i=0; i < (ssize_t) wave_image->columns; i++)
5689     sine_map[i]=fabs(amplitude)+amplitude*sin((double) ((2.0*MagickPI*i)/
5690       wave_length));
5691   /*
5692     Wave image.
5693   */
5694   status=MagickTrue;
5695   progress=0;
5696   canvas_view=AcquireVirtualCacheView(canvas,exception);
5697   wave_view=AcquireAuthenticCacheView(wave_image,exception);
5698   (void) SetCacheViewVirtualPixelMethod(canvas_view,
5699     BackgroundVirtualPixelMethod);
5700 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5701   #pragma omp parallel for schedule(static,4) shared(progress,status) \
5702     magick_threads(canvas,wave_image,wave_image->rows,1)
5703 #endif
5704   for (y=0; y < (ssize_t) wave_image->rows; y++)
5705   {
5706     register Quantum
5707       *magick_restrict q;
5708
5709     register ssize_t
5710       x;
5711
5712     if (status == MagickFalse)
5713       continue;
5714     q=QueueCacheViewAuthenticPixels(wave_view,0,y,wave_image->columns,1,
5715       exception);
5716     if (q == (Quantum *) NULL)
5717       {
5718         status=MagickFalse;
5719         continue;
5720       }
5721     for (x=0; x < (ssize_t) wave_image->columns; x++)
5722     {
5723       status=InterpolatePixelChannels(canvas,canvas_view,wave_image,method,
5724         (double) x,(double) (y-sine_map[x]),q,exception);
5725       q+=GetPixelChannels(wave_image);
5726     }
5727     if (SyncCacheViewAuthenticPixels(wave_view,exception) == MagickFalse)
5728       status=MagickFalse;
5729     if (image->progress_monitor != (MagickProgressMonitor) NULL)
5730       {
5731         MagickBooleanType
5732           proceed;
5733
5734 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5735         #pragma omp critical (MagickCore_WaveImage)
5736 #endif
5737         proceed=SetImageProgress(canvas,WaveImageTag,progress++,canvas->rows);
5738         if (proceed == MagickFalse)
5739           status=MagickFalse;
5740       }
5741   }
5742   wave_view=DestroyCacheView(wave_view);
5743   canvas_view=DestroyCacheView(canvas_view);
5744   canvas=DestroyImage(canvas);
5745   sine_map=(double *) RelinquishMagickMemory(sine_map);
5746   if (status == MagickFalse)
5747     wave_image=DestroyImage(wave_image);
5748   return(wave_image);
5749 }
5750 \f
5751 /*
5752 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5753 %                                                                             %
5754 %                                                                             %
5755 %                                                                             %
5756 %     W a v e l e t D e n o i s e I m a g e                                   %
5757 %                                                                             %
5758 %                                                                             %
5759 %                                                                             %
5760 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5761 %
5762 %  WaveletDenoiseImage() removes noise from the image using a wavelet
5763 %  transform.  The wavelet transform is a fast hierarchical scheme for
5764 %  processing an image using a set of consecutive lowpass and high_pass filters,
5765 %  followed by a decimation.  This results in a decomposition into different
5766 %  scales which can be regarded as different “frequency bands”, determined by
5767 %  the mother wavelet.  Adapted from dcraw.c by David Coffin.
5768 %
5769 %  The format of the WaveletDenoiseImage method is:
5770 %
5771 %      Image *WaveletDenoiseImage(const Image *image,const double threshold,
5772 %        const double softness,ExceptionInfo *exception)
5773 %
5774 %  A description of each parameter follows:
5775 %
5776 %    o image: the image.
5777 %
5778 %    o threshold: set the threshold for smoothing.
5779 %
5780 %    o softness: attenuate the smoothing threshold.
5781 %
5782 %    o exception: return any errors or warnings in this structure.
5783 %
5784 */
5785
5786 static inline void HatTransform(const float *magick_restrict pixels,
5787   const size_t stride,const size_t extent,const size_t scale,float *kernel)
5788 {
5789   const float
5790     *magick_restrict p,
5791     *magick_restrict q,
5792     *magick_restrict r;
5793
5794   register ssize_t
5795     i;
5796
5797   p=pixels;
5798   q=pixels+scale*stride;
5799   r=pixels+scale*stride;
5800   for (i=0; i < (ssize_t) scale; i++)
5801   {
5802     kernel[i]=0.25f*(*p+(*p)+(*q)+(*r));
5803     p+=stride;
5804     q-=stride;
5805     r+=stride;
5806   }
5807   for ( ; i < (ssize_t) (extent-scale); i++)
5808   {
5809     kernel[i]=0.25f*(2.0f*(*p)+*(p-scale*stride)+*(p+scale*stride));
5810     p+=stride;
5811   }
5812   q=p-scale*stride;
5813   r=pixels+stride*(extent-2);
5814   for ( ; i < (ssize_t) extent; i++)
5815   {
5816     kernel[i]=0.25f*(*p+(*p)+(*q)+(*r));
5817     p+=stride;
5818     q+=stride;
5819     r-=stride;
5820   }
5821 }
5822
5823 MagickExport Image *WaveletDenoiseImage(const Image *image,
5824   const double threshold,const double softness,ExceptionInfo *exception)
5825 {
5826   CacheView
5827     *image_view,
5828     *noise_view;
5829
5830   float
5831     *kernel,
5832     *pixels;
5833
5834   Image
5835     *noise_image;
5836
5837   MagickBooleanType
5838     status;
5839
5840   MagickSizeType
5841     number_pixels;
5842
5843   MemoryInfo
5844     *pixels_info;
5845
5846   ssize_t
5847     channel;
5848
5849   static const float
5850     noise_levels[] = { 0.8002f, 0.2735f, 0.1202f, 0.0585f, 0.0291f, 0.0152f,
5851       0.0080f, 0.0044f };
5852
5853   /*
5854     Initialize noise image attributes.
5855   */
5856   assert(image != (const Image *) NULL);
5857   assert(image->signature == MagickCoreSignature);
5858   if (image->debug != MagickFalse)
5859     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
5860   assert(exception != (ExceptionInfo *) NULL);
5861   assert(exception->signature == MagickCoreSignature);
5862 #if defined(MAGICKCORE_OPENCL_SUPPORT)
5863   noise_image=AccelerateWaveletDenoiseImage(image,threshold,exception);
5864   if (noise_image != (Image *) NULL)
5865     return(noise_image);
5866 #endif
5867   noise_image=CloneImage(image,0,0,MagickTrue,exception);
5868   if (noise_image == (Image *) NULL)
5869     return((Image *) NULL);
5870   if (SetImageStorageClass(noise_image,DirectClass,exception) == MagickFalse)
5871     {
5872       noise_image=DestroyImage(noise_image);
5873       return((Image *) NULL);
5874     }
5875   if (AcquireMagickResource(WidthResource,4*image->columns) == MagickFalse)
5876     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
5877   pixels_info=AcquireVirtualMemory(3*image->columns,image->rows*
5878     sizeof(*pixels));
5879   kernel=(float *) AcquireQuantumMemory(MagickMax(image->rows,image->columns),
5880     GetOpenMPMaximumThreads()*sizeof(*kernel));
5881   if ((pixels_info == (MemoryInfo *) NULL) || (kernel == (float *) NULL))
5882     {
5883       if (kernel != (float *) NULL)
5884         kernel=(float *) RelinquishMagickMemory(kernel);
5885       if (pixels_info != (MemoryInfo *) NULL)
5886         pixels_info=RelinquishVirtualMemory(pixels_info);
5887       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
5888     }
5889   pixels=(float *) GetVirtualMemoryBlob(pixels_info);
5890   status=MagickTrue;
5891   number_pixels=(MagickSizeType) image->columns*image->rows;
5892   image_view=AcquireAuthenticCacheView(image,exception);
5893   noise_view=AcquireAuthenticCacheView(noise_image,exception);
5894   for (channel=0; channel < (ssize_t) GetPixelChannels(image); channel++)
5895   {
5896     register ssize_t
5897       i;
5898
5899     size_t
5900       high_pass,
5901       low_pass;
5902
5903     ssize_t
5904       level,
5905       y;
5906
5907     PixelChannel
5908       pixel_channel;
5909
5910     PixelTrait
5911       traits;
5912
5913     if (status == MagickFalse)
5914       continue;
5915     traits=GetPixelChannelTraits(image,(PixelChannel) channel);
5916     if (traits == UndefinedPixelTrait)
5917       continue;
5918     pixel_channel=GetPixelChannelChannel(image,channel);
5919     if ((pixel_channel != RedPixelChannel) &&
5920         (pixel_channel != GreenPixelChannel) &&
5921         (pixel_channel != BluePixelChannel))
5922       continue;
5923     /*
5924       Copy channel from image to wavelet pixel array.
5925     */
5926     i=0;
5927     for (y=0; y < (ssize_t) image->rows; y++)
5928     {
5929       register const Quantum
5930         *magick_restrict p;
5931
5932       ssize_t
5933         x;
5934
5935       p=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
5936       if (p == (const Quantum *) NULL)
5937         {
5938           status=MagickFalse;
5939           break;
5940         }
5941       for (x=0; x < (ssize_t) image->columns; x++)
5942       {
5943         pixels[i++]=(float) p[channel];
5944         p+=GetPixelChannels(image);
5945       }
5946     }
5947     /*
5948       Low pass filter outputs are called approximation kernel & high pass
5949       filters are referred to as detail kernel. The detail kernel
5950       have high values in the noisy parts of the signal.
5951     */
5952     high_pass=0;
5953     for (level=0; level < 5; level++)
5954     {
5955       double
5956         magnitude;
5957
5958       ssize_t
5959         x,
5960         y;
5961
5962       low_pass=(size_t) (number_pixels*((level & 0x01)+1));
5963 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5964       #pragma omp parallel for schedule(static,1) \
5965         magick_threads(image,image,image->rows,1)
5966 #endif
5967       for (y=0; y < (ssize_t) image->rows; y++)
5968       {
5969         const int
5970           id = GetOpenMPThreadId();
5971
5972         register float
5973           *magick_restrict p,
5974           *magick_restrict q;
5975
5976         register ssize_t
5977           x;
5978
5979         p=kernel+id*image->columns;
5980         q=pixels+y*image->columns;
5981         HatTransform(q+high_pass,1,image->columns,(size_t) (1 << level),p);
5982         q+=low_pass;
5983         for (x=0; x < (ssize_t) image->columns; x++)
5984           *q++=(*p++);
5985       }
5986 #if defined(MAGICKCORE_OPENMP_SUPPORT)
5987       #pragma omp parallel for schedule(static,1) \
5988         magick_threads(image,image,image->columns,1)
5989 #endif
5990       for (x=0; x < (ssize_t) image->columns; x++)
5991       {
5992         const int
5993           id = GetOpenMPThreadId();
5994
5995         register float
5996           *magick_restrict p,
5997           *magick_restrict q;
5998
5999         register ssize_t
6000           y;
6001
6002         p=kernel+id*image->rows;
6003         q=pixels+x+low_pass;
6004         HatTransform(q,image->columns,image->rows,(size_t) (1 << level),p);
6005         for (y=0; y < (ssize_t) image->rows; y++)
6006         {
6007           *q=(*p++);
6008           q+=image->columns;
6009         }
6010       }
6011       /*
6012         To threshold, each coefficient is compared to a threshold value and
6013         attenuated / shrunk by some factor.
6014       */
6015       magnitude=threshold*noise_levels[level];
6016       for (i=0; i < (ssize_t) number_pixels; ++i)
6017       {
6018         pixels[high_pass+i]-=pixels[low_pass+i];
6019         if (pixels[high_pass+i] < -magnitude)
6020           pixels[high_pass+i]+=magnitude-softness*magnitude;
6021         else
6022           if (pixels[high_pass+i] > magnitude)
6023             pixels[high_pass+i]-=magnitude-softness*magnitude;
6024           else
6025             pixels[high_pass+i]*=softness;
6026         if (high_pass != 0)
6027           pixels[i]+=pixels[high_pass+i];
6028       }
6029       high_pass=low_pass;
6030     }
6031     /*
6032       Reconstruct image from the thresholded wavelet kernel.
6033     */
6034     i=0;
6035     for (y=0; y < (ssize_t) image->rows; y++)
6036     {
6037       MagickBooleanType
6038         sync;
6039
6040       register Quantum
6041         *magick_restrict q;
6042
6043       register ssize_t
6044         x;
6045
6046       ssize_t
6047         offset;
6048
6049       q=GetCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
6050         exception);
6051       if (q == (Quantum *) NULL)
6052         {
6053           status=MagickFalse;
6054           break;
6055         }
6056       offset=GetPixelChannelOffset(noise_image,pixel_channel);
6057       for (x=0; x < (ssize_t) image->columns; x++)
6058       {
6059         MagickRealType
6060           pixel;
6061
6062         pixel=(MagickRealType) pixels[i]+pixels[low_pass+i];
6063         q[offset]=ClampToQuantum(pixel);
6064         i++;
6065         q+=GetPixelChannels(noise_image);
6066       }
6067       sync=SyncCacheViewAuthenticPixels(noise_view,exception);
6068       if (sync == MagickFalse)
6069         status=MagickFalse;
6070     }
6071     if (image->progress_monitor != (MagickProgressMonitor) NULL)
6072       {
6073         MagickBooleanType
6074           proceed;
6075
6076         proceed=SetImageProgress(image,AddNoiseImageTag,(MagickOffsetType)
6077           channel,GetPixelChannels(image));
6078         if (proceed == MagickFalse)
6079           status=MagickFalse;
6080       }
6081   }
6082   noise_view=DestroyCacheView(noise_view);
6083   image_view=DestroyCacheView(image_view);
6084   kernel=(float *) RelinquishMagickMemory(kernel);
6085   pixels_info=RelinquishVirtualMemory(pixels_info);
6086   if (status == MagickFalse)
6087     noise_image=DestroyImage(noise_image);
6088   return(noise_image);
6089 }