home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / renaisnc / rcnstrct.lha / Reconstruct / Cdf.c next >
Encoding:
C/C++ Source or Header  |  1992-01-04  |  5.0 KB  |  210 lines
  1. /*
  2. ** Author:      James Painter
  3. ** Purpose:     Data structure and access routines for discrete integration
  4. **              (1D)
  5. **
  6. ** Copyright (c), 1988 GRAIL, University of Washington
  7. **
  8. ** $Revision: 1.3 $
  9. ** $Date: 89/01/25 20:17:02 $
  10. ** $Locker: jamie $
  11. ** $Log:    Cdf.c,v $
  12.  * Revision 1.3  89/01/25  20:17:02  jamie
  13.  * *** empty log message ***
  14.  * 
  15.  * Revision 1.3  89/01/11  21:27:02  jamie
  16.  * *** empty log message ***
  17.  * 
  18.  * Revision 1.2  88/10/20  20:39:55  jamie
  19.  * *** empty log message ***
  20.  * 
  21.  * Revision 1.1  88/10/18  21:38:34  jamie
  22.  * Initial revision
  23.  * 
  24. */
  25.  
  26. #include <stdio.h>
  27. #include "Math.h"
  28.  
  29. /* -------------------------------------------------------------------
  30.    Define the size of the table
  31.    ------------------------------------------------------------------- */
  32. #define SIZE 4096
  33.  
  34. /* ---------------------------------------------------------------------
  35.   Define the Cdf table itself.  It is allocated dynamically and accessed
  36.   via an indirection table.
  37.   ---------------------------------------------------------------------- */
  38. static double *Cdf;
  39.  
  40.  
  41. /* ----------------------------------------------------------------------
  42.    The number of radius of the filter
  43.    The filter covers -FilterSupport to +FilterSupport
  44.    ---------------------------------------------------------------------- */
  45. static double FilterSupport = 3.0;
  46.  
  47. /************************
  48.  a univariate Gaussian pdf 
  49.  ************************/
  50. double
  51. Gauss(x)
  52. double x;
  53. {
  54.     static double normalizing_factor =  0.398942280401433;
  55.                                     /* 1.0 / sqrt( 2.0*M_PI) */ ;
  56.  
  57.     return ( normalizing_factor*exp(-(x*x)/2.0) );
  58. }
  59.  
  60.  
  61. double sinc(x)
  62. double x;
  63. {
  64.   double z = M_PI*x;
  65.  
  66.   return (fabs(z) < 1.0E-20) ? 1.0 : sin(z) / z;
  67. }
  68.  
  69. /*
  70.  * univariate scaled Sinc function
  71.  */
  72. double
  73. ScaledSincFunction(x)
  74. double x;
  75. {
  76.   extern double CutOff;
  77.   double scale = 2.0*CutOff;
  78.  
  79.   return scale*sinc( scale*x );
  80. }
  81.  
  82. /*
  83.  * univariate windowed, truncated Sinc function
  84.  */
  85. double
  86. WindowedSincFunction(x)
  87. double x;
  88. {
  89.   extern int window;
  90.   double weight = (window) ? sinc(x/FilterSupport) : 1.0;
  91.  
  92.   return (fabs(x) > FilterSupport) ?  0 : (ScaledSincFunction(x) * weight);
  93. }
  94.  
  95.  
  96. void
  97. SetFilterSupport( support )
  98. double support;
  99. {
  100.   FilterSupport = support;
  101. }
  102.  
  103. /*
  104.    Ten point gaussian-legendre quadrature.
  105.    Reference "Numerical Recipes": p 122"
  106. */
  107. double
  108. qguass( Fcn, a, b )
  109. double (*Fcn)();
  110. double a, b;
  111. {
  112.   register int j;
  113.  
  114.   double xm, xr, dx, result = 0.0;
  115.   static double x[] = {        /* Abscissas */
  116.     .1488743389, .4333953941, .6794095682, .8650633666, .9739065285,
  117.   };
  118.   static double w[] = {        /* Weights */
  119.     .2955242247, .2692667193, .2190863625, .1494513491, .0666713443, 
  120.   };
  121.  
  122.   xm = 0.5*(b+a);
  123.   xr = 0.5*(b-a);
  124.   for(j=0; j<5; j++) {
  125.     dx = xr*x[j];
  126.     result += w[j]*( (*Fcn)(xm+dx) + (*Fcn)(xm-dx) );
  127.   }
  128.   return result*xr;
  129. }
  130.  
  131. /* -----------------------------------------------------------------------
  132.    Build the cumulative density function and store it as a table.
  133.    Integration is approximated by Gauss-Legendre integration
  134.    ------------------------------------------------------------------------ */
  135. void
  136. BuildCdf( Fcn)
  137. double (*Fcn)();
  138. {
  139.   double x,xDelta, previous;
  140.   int c;
  141.  
  142.   /* Allocate memory for the table */
  143.   Cdf = (double *) malloc( SIZE * sizeof(double) );
  144.  
  145.   if (Cdf == (double *) 0L )
  146.     {
  147.       fprintf( stderr, "BuildCdf:  No memory!\n" );
  148.       return;
  149.     }
  150.  
  151.   /* We get an extra entry in our table for free sine there is no point
  152.   ** in storing the zero entry of the Cdf.  
  153.   */
  154.   x = -FilterSupport;
  155.   xDelta = (2*FilterSupport) / (SIZE + 1.0); 
  156.   previous = 0.0;
  157.   for(c=0; c < SIZE; c++)
  158.     {
  159.       /* Integrate from x to x+delta and add to previous total */
  160.       previous +=  qguass( Fcn, x, x+xDelta );
  161.       x = -FilterSupport + c*xDelta;
  162.       Cdf[c] = previous;
  163.     }
  164. }
  165.  
  166.  
  167. /* Access the Cdf with bounds checking. */
  168. #define TableAccess(c)                                                 \
  169. (                                                                      \
  170.   (c < 0)    ?                                                         \
  171.        0.0           :                                                 \
  172.        Cdf [ (c >= SIZE) ? SIZE-1 : c ]                                \
  173. )
  174.  
  175. /* ------------------------------------------------------------------------
  176.    Lookup a value in the (1-d) Cdf.  linear interpolation is used to get
  177.    values in-between table values.  The range of the table is 
  178.    from -FilterSupport to FilterSupport in each direction.
  179.    Extrapolation is not done.  Instead the last entry is held.
  180.   ------------------------------------------------------------------------- */
  181. double
  182. CdfLookup( x )
  183. double x;
  184. {
  185.   double u, um;
  186.   int iu, iv;
  187.   double a, b;
  188.  
  189.   /* Map the input into table row and column indicies. */
  190.   u = SIZE * (x + FilterSupport) / (2.0*FilterSupport)  - 1;
  191.  
  192.   /* Take the floor */
  193.   iu = (int) u;
  194.  
  195.   /* Find the frational part */
  196.   u = u - iu;
  197.   um = 1.0 - u;
  198.  
  199.   /* Get the two values we will interpolate between */
  200.   a = TableAccess(iu) ;
  201.   b = TableAccess(iu+1);
  202.  
  203.   /* Do the interpolation */
  204.   return (um*a + u*b);
  205. }
  206.  
  207.  
  208.  
  209.  
  210.