home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / Libraries / VideoToolbox 94.11.17 / VideoToolboxSources / Normal.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-10-08  |  7.2 KB  |  233 lines  |  [TEXT/KAHL]
  1. /*
  2. Normal.c
  3. statistical functions related to the normal distribution.
  4. Also see: Binomial.c,ChiSquare.c, Exponential.c, Uniform.c
  5.  
  6.     f=NormalPdf(x);
  7.     p=Normal(x);
  8.     x=InverseNormal(p);
  9.     x=NormalSample();
  10.     f=Normal2DPdf(double r);
  11.     p=Normal2D(r);
  12.     r=InverseNormal2D(p);
  13.     r=Normal2DSample();
  14. Normal2D is a Gaussian pdf over two dimensions, integrated over all orientations,
  15. 0 to 2π. r is the distance from the origin, [0,Inf]. 
  16.  
  17.     BoundedNormalIntegers(distribution,n,mean,sd,min,max);
  18. Fills the "distribution" array with n ordered integers so that random samples
  19. from the array,
  20.     i=distribution[nrand(n)];
  21. will have, as nearly as possible, the specified distribution, i.e. they will be
  22. samples (rounded to integer) drawn from the interval [min-.5,max+.5], where min
  23. and max are integers, of a normal distribution with the specified mean and
  24. variance. Once the distribution array has been filled, random samples from it
  25. are a fast way of generating bounded Gaussian noise to be added to each pixel of
  26. an image. Note that "mean" and "sd" apply only to the unclipped underlying
  27. normal distribution. Call this
  28.     mean=MeanW(distribution,n,&sd);
  29. to compute the mean and sd of your integer distribution. The runtime of 
  30. BoundedNormalIntegers will be proportional to max-min+1, and nearly independent
  31. of n. It takes about 0.3 s on a Mac IIci for max-min+1=256.
  32.  
  33. Copyright (c) 1989,1990,1991,1992 Denis G. Pelli
  34. HISTORY:
  35. 1989    dgp wrote it.
  36. 4/8/90    dgp    changed the names of the routines. 
  37.             Made sure that domain error produces NAN.
  38. 6/90    dgp    added NormalSample()
  39. 7/30/91    dgp    now use NAN defined in VideoToolbox.h
  40. 12/28/91 dgp sped up NormalPdf() by calculating the scale factor only once
  41. 12/29/91 dgp extracted code to create new routine UniformSample.c
  42. 1/11/92    dgp    rewrote Normal()'s polynomial evaluation to halve the number of multiplies
  43.             Renamed NormalPDF() to NormalPdf().
  44. 1/19/92    dgp    defined the constants LOG2 and LOGPI in VideoToolbox.h
  45.             Added more domain tests, returning NAN if outside. 
  46.             Added more checks to main().
  47.             Wrote Normal2DPdf(),Normal2D(),InverseNormal2D(),and Normal2DSample().
  48. 2/1/92    dgp    Redefined Normal2DPdf(r) to now, more sensibly, treat r as the random
  49.             variable, rather than the implicit x and y, where r=sqrt(x^2+y^2). 
  50.             Previously, Normal2D(R)=Integrate[2 Pi r Normal2DPdf(r),{r,0,R}]
  51.             now Normal2D(R)=Integrate[Normal2DPdf(r),{r,0,R}]. There is no change
  52.             in Normal2D(). I suspect that Normal2D is in fact the Raleigh distribution,
  53.             and I will rename it if this is in fact the case.
  54. 11/13/92 dgp InverseNormal(0) now returns -INF, and InverseNormal(1) returns INF.
  55. 12/15/93 dgp declared some arguments "register".
  56. 3/26/94    dgp    added BoundedNormalIntegers().
  57.             Asked compiler to use 68881 instructions for transcendentals.
  58. */
  59. #include "VideoToolbox.h"
  60. #include <assert.h>
  61. #include <math.h>
  62. #include "mc68881.h"
  63. void BoundedNormalIntegers(short *distribution,long n,double mean,double sd
  64.     ,short min,short max);
  65. #if THINK_C && mc68881
  66.     #define exp _exp    /* use fast 68881 instruction instead of SANE */
  67.     #define log _log    /* use fast 68881 instruction instead of SANE */
  68.     #define sqrt _sqrt    /* use fast 68881 instruction instead of SANE */
  69. #endif
  70.  
  71. double NormalPdf(register double x)
  72. /* Gaussian pdf. Zero mean and unit variance. */
  73. {
  74.     if(IsNan(x))return x;
  75.     return exp(-0.5*(x*x+(LOG2+LOGPI)));
  76. }
  77.  
  78. double Normal(register double x)
  79. /*
  80. Cumulative normal distribution. From Abramowitz and Stegun Eq. (26.2.17).
  81. Error |e|<7.5 10^-8
  82. */
  83. {
  84.     register double P,t;
  85.     
  86.     if(x<0.0) return 1.0-Normal(-x);
  87.     t=1.0/(1.0+0.2316419*x);
  88.     P=(0.319381530+(-0.356563782+(1.781477937+(-1.821255978+1.330274429*t)*t)*t)*t)*t;
  89.     return 1.0-NormalPdf(x)*P;
  90. }
  91.  
  92. double InverseNormal(register double p)
  93. /*
  94. Inverse of Normal(), based on Abramowitz and Stegun Eq. 26.2.23.
  95. Error |e|<4.5 10^-4.
  96. */
  97. {
  98.     register double t,x;
  99.     
  100.     if(IsNan(p))return p;
  101.     if(p<0.0)return NAN;
  102.     if(p==0.0)return -INF;
  103.     if(p>0.5) return -InverseNormal(1.0-p);
  104.     t=sqrt(-2.0*log(p));
  105.     x=t-(2.515517+(0.802853+0.010328*t)*t)/(1.0+(1.432788+(0.189269+0.001308*t)*t)*t);
  106.     return -x;
  107. }
  108.  
  109. double NormalSample(void)
  110. {
  111.     return InverseNormal(UniformSample());
  112. }
  113.  
  114. double Normal2DPdf(double r)
  115. /* Gaussian pdf over two dimensions, integrated over all orientations, 0 to 2π. */
  116. /* The argument is taken to be the distance from the origin, [0,Inf]. */
  117. /* The rms is 1 */
  118. {
  119.     if(IsNan(r))return r;
  120.     if(r<=0.0)return 0.0;
  121.     return 2*r*exp(-r*r);
  122. }
  123.  
  124. double Normal2D(double r)
  125. /* Integral of Normal2DPdf() from zero to r. */
  126. {
  127.     if(IsNan(r))return r;
  128.     if(r<=0.0)return 0.0;
  129.     return 1.0-exp(-r*r);
  130. }
  131.  
  132. double InverseNormal2D(double p)
  133. {
  134.     if(IsNan(p))return p;
  135.     if(p<0.0 || p>1.0)return NAN;
  136.     return sqrt(-log(1.0-p));
  137. }
  138.  
  139. double Normal2DSample(void)
  140. /* rms is 1 */
  141. {
  142.     return InverseNormal2D(UniformSample());
  143. }
  144.  
  145. /*
  146.     Integrate[exp(-.5*u^2),{u,a,a+1/sd}]
  147.     =exp(-.5*a^2)*Integrate[exp(-.5*((a+e)^2-a^2)),{e,0,1/sd}]
  148.     =exp(-.5*a^2)*Integrate[exp(-.5*e*e)*exp(-a*e),{e,0,1/sd}]
  149.     =exp(-.5*a^2)*Integrate[(1-.5*e*e)*exp(-a*e),{e,0,1/sd}]
  150.     =exp(-.5*a^2)*((exp(-a/sd) - 1)/(-a)-.5*Integrate[e*e*exp(-a*e),{e,0,1/sd}])
  151.     =exp(-.5*a^2)*(1-exp(-a/sd))/a
  152. */
  153. void BoundedNormalIntegers(register short *distribution,long n,double mean,double sd
  154.     ,short min,short max)
  155. {
  156.     register short i;
  157.     register long j,count;
  158.     double p,p0,p1,x;
  159.     short shortcut;
  160.     
  161.     shortcut=sd*sd>n;    // guarantees count will err by at most ±0.5
  162.     j=0;
  163.     p=p0=Normal((min-.5-mean)/sd);
  164.     p1=Normal((max+.5-mean)/sd)-p0;
  165.     for(i=min;i<max;i++){
  166.         x=(i+.5-mean)/sd;
  167.         if(shortcut)p+=NormalPdf(x)*(exp(x/sd)-1)/x;
  168.         else p=Normal(x);
  169.         count=0.5+n*(p-p0)/p1;
  170.         while(j<count)distribution[j++]=i;
  171.     }
  172.     while(j<n)distribution[j++]=max;
  173. }
  174.  
  175.  
  176. #if 0 /* A test program. */
  177.     void main()
  178.     {
  179.         double x,y,sum,dx,a,b,mean,sd;
  180.         static double z[1000];
  181.         int i;
  182.         
  183.         Require(0);
  184.         srand(clock());
  185.         printf("%4s%15s%15s%20s%15s\n","x","NormalPdf(x)","Normal(x)","InverseNormal","Error");
  186.         for(x=-4.0;x<=4.0;x+=2.0){
  187.             printf("%4.1f%15.8f%15.8f%20.4f%15.4f\n",
  188.             x,NormalPdf(x),Normal(x),InverseNormal(Normal(x)),InverseNormal(Normal(x))-x);
  189.         }
  190.         sum=0.0;
  191.         dx=0.001;
  192.         for(x=-1.;x<0.;x+=dx)sum+=NormalPdf(x);
  193.         sum*=dx;
  194.         sum-=Normal(0.0)-Normal(-1.0);
  195.         printf("Partial integral of NormalPdf error %.5f\n",sum);
  196.         for(i=0;i<1000;i++)z[i]=NormalSample();
  197.         mean=Mean(z,1000,&sd);
  198.         printf("1000 samples mean %.2f sd %.2f\n",mean,sd);
  199.         printf("\n");
  200.     
  201.         printf("%4s%15s%15s%20s%15s\n","x","Normal2DPdf(x)","Normal2D(x)","InverseNormal2D","Error");
  202.         for(x=-1.;x<=5.0;x+=1.0){
  203.             printf("%4.1f%15.8f%15.8f%20.4f%15.4f\n",
  204.             x,Normal2DPdf(x),Normal2D(x),InverseNormal2D(Normal2D(x)),InverseNormal2D(Normal2D(x))-x);
  205.         }
  206.         sum=0.0;
  207.         dx=0.0001;
  208.         for(x=0;x<1.;x+=dx)sum+=Normal2DPdf(x);
  209.         sum*=dx;
  210.         sum-=Normal2D(1.0);
  211.         printf("Partial integral of Normal2DPdf error %.5f\n",sum);
  212.         for(i=0;i<1000;i++)z[i]=Normal2DSample();
  213.         mean=Mean(z,1000,&sd);
  214.         printf("1000 samples rms %.2f\n",sqrt(mean*mean+sd*sd));
  215.         printf("\n");
  216.         for(i=0;i<1000;i++){
  217.             x=NormalSample();
  218.             y=NormalSample();
  219.             z[i]=sqrt((x*x+y*y)/2.);
  220.         }
  221.         mean=Mean(z,1000,&sd);
  222.         printf("1000 (x,y) normal samples with sd 2^-0.5 have rms hypotenuse of %.2f\n",sqrt(mean*mean+sd*sd));
  223.         printf("\n");
  224.     
  225.         a=4.0*atan(1.0);
  226.         if(a!=PI)printf("4*atan(1)-PI %.19f\n",a-PI);
  227.         a=log(a);
  228.         if(a!=LOGPI)printf("Error: log(PI) %.19f, error in LOGPI %.19f\n",a,LOGPI-a);
  229.         a=log(2.0);
  230.         if(a!=LOG2)printf("Error: log(2) %.19f, error in LOG2 %.19f\n",a,LOG2-a);
  231.     }
  232. #endif
  233.