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