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