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

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %                              GGGG  EEEEE  M   M                             %
   7 %                             G      E      MM MM                             %
   8 %                             G GG   EEE    M M M                             %
   9 %                             G   G  E      M   M                             %
  10 %                              GGGG  EEEEE  M   M                             %
  11 %                                                                             %
  12 %                                                                             %
  13 %                    Graphic Gems - Graphic Support Methods                   %
  14 %                                                                             %
  15 %                               Software Design                               %
  16 %                                 John Cristy                                 %
  17 %                                 August 1996                                 %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "MagickCore/studio.h"
  44 #include "MagickCore/color-private.h"
  45 #include "MagickCore/draw.h"
  46 #include "MagickCore/gem.h"
  47 #include "MagickCore/gem-private.h"
  48 #include "MagickCore/image.h"
  49 #include "MagickCore/image-private.h"
  50 #include "MagickCore/log.h"
  51 #include "MagickCore/memory_.h"
  52 #include "MagickCore/pixel-accessor.h"
  53 #include "MagickCore/pixel-private.h"
  54 #include "MagickCore/quantum.h"
  55 #include "MagickCore/quantum-private.h"
  56 #include "MagickCore/random_.h"
  57 #include "MagickCore/resize.h"
  58 #include "MagickCore/transform.h"
  59 #include "MagickCore/signature-private.h"
  60 \f
  61 /*
  62 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  63 %                                                                             %
  64 %                                                                             %
  65 %                                                                             %
  66 %   C o n v e r t H S B T o R G B                                             %
  67 %                                                                             %
  68 %                                                                             %
  69 %                                                                             %
  70 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  71 %
  72 %  ConvertHSBToRGB() transforms a (hue, saturation, brightness) to a (red,
  73 %  green, blue) triple.
  74 %
  75 %  The format of the ConvertHSBToRGBImage method is:
  76 %
  77 %      void ConvertHSBToRGB(const double hue,const double saturation,
  78 %        const double brightness,double *red,double *green,double *blue)
  79 %
  80 %  A description of each parameter follows:
  81 %
  82 %    o hue, saturation, brightness: A double value representing a
  83 %      component of the HSB color space.
  84 %
  85 %    o red, green, blue: A pointer to a pixel component of type Quantum.
  86 %
  87 */
  88 MagickPrivate void ConvertHSBToRGB(const double hue,const double saturation,
  89   const double brightness,double *red,double *green,double *blue)
  90 {
  91   double
  92     f,
  93     h,
  94     p,
  95     q,
  96     t;
  97
  98   /*
  99     Convert HSB to RGB colorspace.
 100   */
 101   assert(red != (double *) NULL);
 102   assert(green != (double *) NULL);
 103   assert(blue != (double *) NULL);
 104   if (saturation == 0.0)
 105     {
 106       *red=QuantumRange*brightness;
 107       *green=(*red);
 108       *blue=(*red);
 109       return;
 110     }
 111   h=6.0*(hue-floor(hue));
 112   f=h-floor((double) h);
 113   p=brightness*(1.0-saturation);
 114   q=brightness*(1.0-saturation*f);
 115   t=brightness*(1.0-(saturation*(1.0-f)));
 116   switch ((int) h)
 117   {
 118     case 0:
 119     default:
 120     {
 121       *red=QuantumRange*brightness;
 122       *green=QuantumRange*t;
 123       *blue=QuantumRange*p;
 124       break;
 125     }
 126     case 1:
 127     {
 128       *red=QuantumRange*q;
 129       *green=QuantumRange*brightness;
 130       *blue=QuantumRange*p;
 131       break;
 132     }
 133     case 2:
 134     {
 135       *red=QuantumRange*p;
 136       *green=QuantumRange*brightness;
 137       *blue=QuantumRange*t;
 138       break;
 139     }
 140     case 3:
 141     {
 142       *red=QuantumRange*p;
 143       *green=QuantumRange*q;
 144       *blue=QuantumRange*brightness;
 145       break;
 146     }
 147     case 4:
 148     {
 149       *red=QuantumRange*t;
 150       *green=QuantumRange*p;
 151       *blue=QuantumRange*brightness;
 152       break;
 153     }
 154     case 5:
 155     {
 156       *red=QuantumRange*brightness;
 157       *green=QuantumRange*p;
 158       *blue=QuantumRange*q;
 159       break;
 160     }
 161   }
 162 }
 163 \f
 164 /*
 165 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 166 %                                                                             %
 167 %                                                                             %
 168 %                                                                             %
 169 %   C o n v e r t H S L T o R G B                                             %
 170 %                                                                             %
 171 %                                                                             %
 172 %                                                                             %
 173 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 174 %
 175 %  ConvertHSLToRGB() transforms a (hue, saturation, lightness) to a (red,
 176 %  green, blue) triple.
 177 %
 178 %  The format of the ConvertHSLToRGBImage method is:
 179 %
 180 %      void ConvertHSLToRGB(const double hue,const double saturation,
 181 %        const double lightness,double *red,double *green,double *blue)
 182 %
 183 %  A description of each parameter follows:
 184 %
 185 %    o hue, saturation, lightness: A double value representing a
 186 %      component of the HSL color space.
 187 %
 188 %    o red, green, blue: A pointer to a pixel component of type Quantum.
 189 %
 190 */
 191
 192 static inline double ConvertHueToRGB(double m1,double m2,double hue)
 193 {
 194   if (hue < 0.0)
 195     hue+=1.0;
 196   if (hue > 1.0)
 197     hue-=1.0;
 198   if ((6.0*hue) < 1.0)
 199     return(m1+6.0*(m2-m1)*hue);
 200   if ((2.0*hue) < 1.0)
 201     return(m2);
 202   if ((3.0*hue) < 2.0)
 203     return(m1+6.0*(m2-m1)*(2.0/3.0-hue));
 204   return(m1);
 205 }
 206
 207 MagickExport void ConvertHSLToRGB(const double hue,const double saturation,
 208   const double lightness,double *red,double *green,double *blue)
 209 {
 210   double
 211     b,
 212     g,
 213     r,
 214     m1,
 215     m2;
 216
 217   /*
 218     Convert HSL to RGB colorspace.
 219   */
 220   assert(red != (double *) NULL);
 221   assert(green != (double *) NULL);
 222   assert(blue != (double *) NULL);
 223   if (saturation == 0)
 224     {
 225       *red=QuantumRange*lightness;
 226       *green=(*red);
 227       *blue=(*red);
 228       return;
 229     }
 230   if (lightness < 0.5)
 231     m2=lightness*(saturation+1.0);
 232   else
 233     m2=(lightness+saturation)-(lightness*saturation);
 234   m1=2.0*lightness-m2;
 235   r=ConvertHueToRGB(m1,m2,hue+1.0/3.0);
 236   g=ConvertHueToRGB(m1,m2,hue);
 237   b=ConvertHueToRGB(m1,m2,hue-1.0/3.0);
 238   *red=QuantumRange*r;
 239   *green=QuantumRange*g;
 240   *blue=QuantumRange*b;
 241 }
 242 \f
 243 /*
 244 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 245 %                                                                             %
 246 %                                                                             %
 247 %                                                                             %
 248 %   C o n v e r t H W B T o R G B                                             %
 249 %                                                                             %
 250 %                                                                             %
 251 %                                                                             %
 252 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 253 %
 254 %  ConvertHWBToRGB() transforms a (hue, whiteness, blackness) to a (red, green,
 255 %  blue) triple.
 256 %
 257 %  The format of the ConvertHWBToRGBImage method is:
 258 %
 259 %      void ConvertHWBToRGB(const double hue,const double whiteness,
 260 %        const double blackness,double *red,double *green,double *blue)
 261 %
 262 %  A description of each parameter follows:
 263 %
 264 %    o hue, whiteness, blackness: A double value representing a
 265 %      component of the HWB color space.
 266 %
 267 %    o red, green, blue: A pointer to a pixel component of type Quantum.
 268 %
 269 */
 270 MagickPrivate void ConvertHWBToRGB(const double hue,const double whiteness,
 271   const double blackness,double *red,double *green,double *blue)
 272 {
 273   double
 274     b,
 275     f,
 276     g,
 277     n,
 278     r,
 279     v;
 280
 281   register ssize_t
 282     i;
 283
 284   /*
 285     Convert HWB to RGB colorspace.
 286   */
 287   assert(red != (double *) NULL);
 288   assert(green != (double *) NULL);
 289   assert(blue != (double *) NULL);
 290   v=1.0-blackness;
 291   if (hue == -1.0)
 292     {
 293       *red=QuantumRange*v;
 294       *green=QuantumRange*v;
 295       *blue=QuantumRange*v;
 296       return;
 297     }
 298   i=(ssize_t) floor(6.0*hue);
 299   f=6.0*hue-i;
 300   if ((i & 0x01) != 0)
 301     f=1.0-f;
 302   n=whiteness+f*(v-whiteness);  /* linear interpolation */
 303   switch (i)
 304   {
 305     default:
 306     case 6:
 307     case 0: r=v; g=n; b=whiteness; break;
 308     case 1: r=n; g=v; b=whiteness; break;
 309     case 2: r=whiteness; g=v; b=n; break;
 310     case 3: r=whiteness; g=n; b=v; break;
 311     case 4: r=n; g=whiteness; b=v; break;
 312     case 5: r=v; g=whiteness; b=n; break;
 313   }
 314   *red=QuantumRange*r;
 315   *green=QuantumRange*g;
 316   *blue=QuantumRange*b;
 317 }
 318 \f
 319 /*
 320 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 321 %                                                                             %
 322 %                                                                             %
 323 %                                                                             %
 324 %   C o n v e r t R G B T o H S B                                             %
 325 %                                                                             %
 326 %                                                                             %
 327 %                                                                             %
 328 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 329 %
 330 %  ConvertRGBToHSB() transforms a (red, green, blue) to a (hue, saturation,
 331 %  brightness) triple.
 332 %
 333 %  The format of the ConvertRGBToHSB method is:
 334 %
 335 %      void ConvertRGBToHSB(const double red,const double green,
 336 %        const double blue,double *hue,double *saturation,double *brightness)
 337 %
 338 %  A description of each parameter follows:
 339 %
 340 %    o red, green, blue: A Quantum value representing the red, green, and
 341 %      blue component of a pixel..
 342 %
 343 %    o hue, saturation, brightness: A pointer to a double value representing a
 344 %      component of the HSB color space.
 345 %
 346 */
 347 MagickPrivate void ConvertRGBToHSB(const double red,const double green,
 348   const double blue,double *hue,double *saturation,double *brightness)
 349 {
 350   double
 351     b,
 352     delta,
 353     g,
 354     max,
 355     min,
 356     r;
 357
 358   /*
 359     Convert RGB to HSB colorspace.
 360   */
 361   assert(hue != (double *) NULL);
 362   assert(saturation != (double *) NULL);
 363   assert(brightness != (double *) NULL);
 364   *hue=0.0;
 365   *saturation=0.0;
 366   *brightness=0.0;
 367   r=red;
 368   g=green;
 369   b=blue;
 370   min=r < g ? r : g;
 371   if (b < min)
 372     min=b;
 373   max=r > g ? r : g;
 374   if (b > max)
 375     max=b;
 376   if (max == 0.0)
 377     return;
 378   delta=max-min;
 379   *saturation=delta/max;
 380   *brightness=QuantumScale*max;
 381   if (delta == 0.0)
 382     return;
 383   if (r == max)
 384     *hue=(g-b)/delta;
 385   else
 386     if (g == max)
 387       *hue=2.0+(b-r)/delta;
 388     else
 389       *hue=4.0+(r-g)/delta;
 390   *hue/=6.0;
 391   if (*hue < 0.0)
 392     *hue+=1.0;
 393 }
 394 \f
 395 /*
 396 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 397 %                                                                             %
 398 %                                                                             %
 399 %                                                                             %
 400 %   C o n v e r t R G B T o H S L                                             %
 401 %                                                                             %
 402 %                                                                             %
 403 %                                                                             %
 404 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 405 %
 406 %  ConvertRGBToHSL() transforms a (red, green, blue) to a (hue, saturation,
 407 %  lightness) triple.
 408 %
 409 %  The format of the ConvertRGBToHSL method is:
 410 %
 411 %      void ConvertRGBToHSL(const double red,const double green,
 412 %        const double blue,double *hue,double *saturation,double *lightness)
 413 %
 414 %  A description of each parameter follows:
 415 %
 416 %    o red, green, blue: A Quantum value representing the red, green, and
 417 %      blue component of a pixel..
 418 %
 419 %    o hue, saturation, lightness: A pointer to a double value representing a
 420 %      component of the HSL color space.
 421 %
 422 */
 423
 424 static inline double MagickMax(const double x,const double y)
 425 {
 426   if (x > y)
 427     return(x);
 428   return(y);
 429 }
 430
 431 static inline double MagickMin(const double x,const double y)
 432 {
 433   if (x < y)
 434     return(x);
 435   return(y);
 436 }
 437
 438 MagickExport void ConvertRGBToHSL(const double red,const double green,
 439   const double blue,double *hue,double *saturation,double *lightness)
 440 {
 441   double
 442     b,
 443     delta,
 444     g,
 445     max,
 446     min,
 447     r;
 448
 449   /*
 450     Convert RGB to HSL colorspace.
 451   */
 452   assert(hue != (double *) NULL);
 453   assert(saturation != (double *) NULL);
 454   assert(lightness != (double *) NULL);
 455   r=QuantumScale*red;
 456   g=QuantumScale*green;
 457   b=QuantumScale*blue;
 458   max=MagickMax(r,MagickMax(g,b));
 459   min=MagickMin(r,MagickMin(g,b));
 460   *lightness=(double) ((min+max)/2.0);
 461   delta=max-min;
 462   if (delta == 0.0)
 463     {
 464       *hue=0.0;
 465       *saturation=0.0;
 466       return;
 467     }
 468   if (*lightness < 0.5)
 469     *saturation=(double) (delta/(min+max));
 470   else
 471     *saturation=(double) (delta/(2.0-max-min));
 472   if (r == max)
 473     *hue=((((max-b)/6.0)+(delta/2.0))-(((max-g)/6.0)+(delta/2.0)))/delta;
 474   else
 475     if (g == max)
 476       *hue=(1.0/3.0)+((((max-r)/6.0)+(delta/2.0))-(((max-b)/6.0)+(delta/2.0)))/
 477         delta;
 478     else
 479       if (b == max)
 480         *hue=(2.0/3.0)+((((max-g)/6.0)+(delta/2.0))-(((max-r)/6.0)+
 481           (delta/2.0)))/delta;
 482   if (*hue < 0.0)
 483     *hue+=1.0;
 484   if (*hue > 1.0)
 485     *hue-=1.0;
 486 }
 487 \f
 488 /*
 489 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 490 %                                                                             %
 491 %                                                                             %
 492 %                                                                             %
 493 %   C o n v e r t R G B T o H W B                                             %
 494 %                                                                             %
 495 %                                                                             %
 496 %                                                                             %
 497 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 498 %
 499 %  ConvertRGBToHWB() transforms a (red, green, blue) to a (hue, whiteness,
 500 %  blackness) triple.
 501 %
 502 %  The format of the ConvertRGBToHWB method is:
 503 %
 504 %      void ConvertRGBToHWB(const double red,const double green,
 505 %        const double blue,double *hue,double *whiteness,double *blackness)
 506 %
 507 %  A description of each parameter follows:
 508 %
 509 %    o red, green, blue: A Quantum value representing the red, green, and
 510 %      blue component of a pixel.
 511 %
 512 %    o hue, whiteness, blackness: A pointer to a double value representing a
 513 %      component of the HWB color space.
 514 %
 515 */
 516 MagickPrivate void ConvertRGBToHWB(const double red,const double green,
 517   const double blue,double *hue,double *whiteness,double *blackness)
 518 {
 519   double
 520     b,
 521     f,
 522     g,
 523     p,
 524     r,
 525     v,
 526     w;
 527
 528   /*
 529     Convert RGB to HWB colorspace.
 530   */
 531   assert(hue != (double *) NULL);
 532   assert(whiteness != (double *) NULL);
 533   assert(blackness != (double *) NULL);
 534   r=red;
 535   g=green;
 536   b=blue;
 537   w=MagickMin(r,MagickMin(g,b));
 538   v=MagickMax(r,MagickMax(g,b));
 539   *blackness=1.0-QuantumScale*v;
 540   *whiteness=QuantumScale*w;
 541   if (v == w)
 542     {
 543       *hue=(-1.0);
 544       return;
 545     }
 546   f=(r == w) ? g-b : ((g == w) ? b-r : r-g);
 547   p=(r == w) ? 3.0 : ((g == w) ? 5.0 : 1.0);
 548   *hue=(p-f/(v-1.0*w))/6.0;
 549 }
 550 \f
 551 /*
 552 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 553 %                                                                             %
 554 %                                                                             %
 555 %                                                                             %
 556 %   E x p a n d A f f i n e                                                   %
 557 %                                                                             %
 558 %                                                                             %
 559 %                                                                             %
 560 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 561 %
 562 %  ExpandAffine() computes the affine's expansion factor, i.e. the square root
 563 %  of the factor by which the affine transform affects area. In an affine
 564 %  transform composed of scaling, rotation, shearing, and translation, returns
 565 %  the amount of scaling.
 566 %
 567 %  The format of the ExpandAffine method is:
 568 %
 569 %      double ExpandAffine(const AffineMatrix *affine)
 570 %
 571 %  A description of each parameter follows:
 572 %
 573 %    o expansion: Method ExpandAffine returns the affine's expansion factor.
 574 %
 575 %    o affine: A pointer the affine transform of type AffineMatrix.
 576 %
 577 */
 578 MagickExport double ExpandAffine(const AffineMatrix *affine)
 579 {
 580   assert(affine != (const AffineMatrix *) NULL);
 581   return(sqrt(fabs(affine->sx*affine->sy-affine->rx*affine->ry)));
 582 }
 583 \f
 584 /*
 585 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 586 %                                                                             %
 587 %                                                                             %
 588 %                                                                             %
 589 %   G e n e r a t e D i f f e r e n t i a l N o i s e                         %
 590 %                                                                             %
 591 %                                                                             %
 592 %                                                                             %
 593 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 594 %
 595 %  GenerateDifferentialNoise() generates differentual noise.
 596 %
 597 %  The format of the GenerateDifferentialNoise method is:
 598 %
 599 %      double GenerateDifferentialNoise(RandomInfo *random_info,
 600 %        const Quantum pixel,const NoiseType noise_type,const double attenuate)
 601 %
 602 %  A description of each parameter follows:
 603 %
 604 %    o random_info: the random info.
 605 %
 606 %    o pixel: noise is relative to this pixel value.
 607 %
 608 %    o noise_type: the type of noise.
 609 %
 610 %    o attenuate:  attenuate the noise.
 611 %
 612 */
 613 MagickPrivate double GenerateDifferentialNoise(RandomInfo *random_info,
 614   const Quantum pixel,const NoiseType noise_type,const double attenuate)
 615 {
 616 #define SigmaUniform  (attenuate*0.015625)
 617 #define SigmaGaussian  (attenuate*0.015625)
 618 #define SigmaImpulse  (attenuate*0.1)
 619 #define SigmaLaplacian (attenuate*0.0390625)
 620 #define SigmaMultiplicativeGaussian  (attenuate*0.5)
 621 #define SigmaPoisson  (attenuate*12.5)
 622 #define SigmaRandom  (attenuate)
 623 #define TauGaussian  (attenuate*0.078125)
 624
 625   double
 626     alpha,
 627     beta,
 628     noise,
 629     sigma;
 630
 631   alpha=GetPseudoRandomValue(random_info);
 632   switch (noise_type)
 633   {
 634     case UniformNoise:
 635     default:
 636     {
 637       noise=(double) (pixel+QuantumRange*SigmaUniform*(alpha-0.5));
 638       break;
 639     }
 640     case GaussianNoise:
 641     {
 642       double
 643         gamma,
 644         tau;
 645
 646       if (alpha == 0.0)
 647         alpha=1.0;
 648       beta=GetPseudoRandomValue(random_info);
 649       gamma=sqrt(-2.0*log(alpha));
 650       sigma=gamma*cos((double) (2.0*MagickPI*beta));
 651       tau=gamma*sin((double) (2.0*MagickPI*beta));
 652       noise=(double) (pixel+sqrt((double) pixel)*SigmaGaussian*sigma+
 653         QuantumRange*TauGaussian*tau);
 654       break;
 655     }
 656     case ImpulseNoise:
 657     {
 658       if (alpha < (SigmaImpulse/2.0))
 659         noise=0.0;
 660       else
 661         if (alpha >= (1.0-(SigmaImpulse/2.0)))
 662           noise=(double) QuantumRange;
 663         else
 664           noise=(double) pixel;
 665       break;
 666     }
 667     case LaplacianNoise:
 668     {
 669       if (alpha <= 0.5)
 670         {
 671           if (alpha <= MagickEpsilon)
 672             noise=(double) (pixel-QuantumRange);
 673           else
 674             noise=(double) (pixel+QuantumRange*SigmaLaplacian*log(2.0*alpha)+
 675               0.5);
 676           break;
 677         }
 678       beta=1.0-alpha;
 679       if (beta <= (0.5*MagickEpsilon))
 680         noise=(double) (pixel+QuantumRange);
 681       else
 682         noise=(double) (pixel-QuantumRange*SigmaLaplacian*log(2.0*beta)+0.5);
 683       break;
 684     }
 685     case MultiplicativeGaussianNoise:
 686     {
 687       sigma=1.0;
 688       if (alpha > MagickEpsilon)
 689         sigma=sqrt(-2.0*log(alpha));
 690       beta=GetPseudoRandomValue(random_info);
 691       noise=(double) (pixel+pixel*SigmaMultiplicativeGaussian*sigma*
 692         cos((double) (2.0*MagickPI*beta))/2.0);
 693       break;
 694     }
 695     case PoissonNoise:
 696     {
 697       double
 698         poisson;
 699
 700       register ssize_t
 701         i;
 702
 703       poisson=exp(-SigmaPoisson*QuantumScale*pixel);
 704       for (i=0; alpha > poisson; i++)
 705       {
 706         beta=GetPseudoRandomValue(random_info);
 707         alpha*=beta;
 708       }
 709       noise=(double) (QuantumRange*i/SigmaPoisson);
 710       break;
 711     }
 712     case RandomNoise:
 713     {
 714       noise=(double) (QuantumRange*SigmaRandom*alpha);
 715       break;
 716     }
 717   }
 718   return(noise);
 719 }
 720 \f
 721 /*
 722 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 723 %                                                                             %
 724 %                                                                             %
 725 %                                                                             %
 726 %   G e t O p t i m a l K e r n e l W i d t h                                 %
 727 %                                                                             %
 728 %                                                                             %
 729 %                                                                             %
 730 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 731 %
 732 %  GetOptimalKernelWidth() computes the optimal kernel radius for a convolution
 733 %  filter.  Start with the minimum value of 3 pixels and walk out until we drop
 734 %  below the threshold of one pixel numerical accuracy.
 735 %
 736 %  The format of the GetOptimalKernelWidth method is:
 737 %
 738 %      size_t GetOptimalKernelWidth(const double radius,
 739 %        const double sigma)
 740 %
 741 %  A description of each parameter follows:
 742 %
 743 %    o width: Method GetOptimalKernelWidth returns the optimal width of
 744 %      a convolution kernel.
 745 %
 746 %    o radius: the radius of the Gaussian, in pixels, not counting the center
 747 %      pixel.
 748 %
 749 %    o sigma: the standard deviation of the Gaussian, in pixels.
 750 %
 751 */
 752 MagickPrivate size_t GetOptimalKernelWidth1D(const double radius,
 753   const double sigma)
 754 {
 755   double
 756     alpha,
 757     beta,
 758     gamma,
 759     normalize,
 760     value;
 761
 762   register ssize_t
 763     i;
 764
 765   size_t
 766     width;
 767
 768   ssize_t
 769     j;
 770
 771   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
 772   if (radius > MagickEpsilon)
 773     return((size_t) (2.0*ceil(radius)+1.0));
 774   gamma=fabs(sigma);
 775   if (gamma <= MagickEpsilon)
 776     return(3UL);
 777   alpha=MagickEpsilonReciprocal(2.0*gamma*gamma);
 778   beta=(double) MagickEpsilonReciprocal((MagickRealType) MagickSQ2PI*gamma);
 779   for (width=5; ; )
 780   {
 781     normalize=0.0;
 782     j=(ssize_t) width/2;
 783     for (i=(-j); i <= j; i++)
 784       normalize+=exp(-((double) (i*i))*alpha)*beta;
 785     value=exp(-((double) (j*j))*alpha)*beta/normalize;
 786     if ((value < QuantumScale) || (value < MagickEpsilon))
 787       break;
 788     width+=2;
 789   }
 790   return((size_t) (width-2));
 791 }
 792
 793 MagickPrivate size_t GetOptimalKernelWidth2D(const double radius,
 794   const double sigma)
 795 {
 796   double
 797     alpha,
 798     beta,
 799     gamma,
 800     normalize,
 801     value;
 802
 803   size_t
 804     width;
 805
 806   ssize_t
 807     j,
 808     u,
 809     v;
 810
 811   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
 812   if (radius > MagickEpsilon)
 813     return((size_t) (2.0*ceil(radius)+1.0));
 814   gamma=fabs(sigma);
 815   if (gamma <= MagickEpsilon)
 816     return(3UL);
 817   alpha=MagickEpsilonReciprocal(2.0*gamma*gamma);
 818   beta=(double) MagickEpsilonReciprocal((MagickRealType) Magick2PI*gamma*gamma);
 819   for (width=5; ; )
 820   {
 821     normalize=0.0;
 822     j=(ssize_t) width/2;
 823     for (v=(-j); v <= j; v++)
 824       for (u=(-j); u <= j; u++)
 825         normalize+=exp(-((double) (u*u+v*v))*alpha)*beta;
 826     value=exp(-((double) (j*j))*alpha)*beta/normalize;
 827     if ((value < QuantumScale) || (value < MagickEpsilon))
 828       break;
 829     width+=2;
 830   }
 831   return((size_t) (width-2));
 832 }
 833
 834 MagickPrivate size_t  GetOptimalKernelWidth(const double radius,
 835   const double sigma)
 836 {
 837   return(GetOptimalKernelWidth1D(radius,sigma));
 838 }