00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifndef LUX_MCDISTRIBUTION_H
00024 #define LUX_MCDISTRIBUTION_H
00025
00026 #include "mc.h"
00027
00028 namespace lux {
00029
00033 class Function1D {
00034 public:
00043 Function1D(float *f, int n) {
00044 func = new float[n];
00045 count = n;
00046 memcpy(func, f, n*sizeof(float));
00047 }
00048 ~Function1D() {
00049 delete[] func;
00050 }
00051
00059 float Eval(float x) const {
00060 float pos = Clamp(x, 0.f, 1.f) * count + .5f;
00061 int off1 = (int)pos;
00062 int off2 = min(count-1, off1 + 1);
00063 float d = pos - off1;
00064 return func[off1] * (1.f - d) * func[off2] * d;
00065 }
00066
00067
00068
00069
00070
00071 float *func;
00072
00073
00074
00075 int count;
00076 };
00077
00081 class Distribution1D {
00082 public:
00091 Distribution1D(float *f, int n) {
00092 func = new float[n];
00093 cdf = new float[n+1];
00094 count = n;
00095 memcpy(func, f, n*sizeof(float));
00096 ComputeStep1dCDF(func, n, &funcInt, cdf);
00097 invFuncInt = 1.f / funcInt;
00098 invCount = 1.f / count;
00099 }
00100 ~Distribution1D() {
00101 delete[] func;
00102 delete[] cdf;
00103 }
00104
00114 float Sample(float u, float *pdf) const {
00115
00116 float *ptr = std::lower_bound(cdf, cdf+count+1, u);
00117
00118 int offset = max(0, (int)(ptr-cdf-1));
00119
00120 u = (u - cdf[offset]) / (cdf[offset+1] - cdf[offset]);
00121 *pdf = func[offset] * invFuncInt;
00122 return offset + u;
00123 }
00124
00125
00126
00127
00128
00129 float *func, *cdf;
00134 float funcInt, invFuncInt, invCount;
00135
00136
00137
00138 int count;
00139 };
00140
00144 class IrregularFunction1D {
00145 public:
00156 IrregularFunction1D(float *aX, float *aFx, int aN) {
00157 count = aN;
00158 xFunc = new float[count];
00159 yFunc = new float[count];
00160 memcpy(xFunc, aX, aN*sizeof(float));
00161 memcpy(yFunc, aFx, aN*sizeof(float));
00162 }
00163
00164 ~IrregularFunction1D() {
00165 delete[] xFunc;
00166 delete[] yFunc;
00167 }
00168
00176 float Eval(float x) const {
00177 if( x <= xFunc[0] )
00178 return yFunc[ 0 ];
00179 else if( x >= xFunc[ count - 1 ] )
00180 return yFunc[ count - 1 ];
00181
00182 float *ptr = std::upper_bound(xFunc, xFunc+count, x);
00183 int offset = max(0, (int)(ptr-xFunc-1));
00184
00185 float d = (x - xFunc[offset]) / (xFunc[offset+1] - xFunc[offset]);
00186
00187 return (1.f - d) * yFunc[offset] + d * yFunc[offset+1];
00188 }
00189
00198 int IndexOf(float x, float *d) const {
00199 if( x <= xFunc[0] ) {
00200 *d = 0.f;
00201 return 0;
00202 }
00203 else if( x >= xFunc[ count - 1 ] ) {
00204 *d = 0.f;
00205 return count - 1;
00206 }
00207
00208 float *ptr = std::upper_bound(xFunc, xFunc+count, x);
00209 int offset = (int) (ptr-xFunc-1);
00210
00211 *d = (x - xFunc[offset]) / (xFunc[offset+1] - xFunc[offset]);
00212 return offset;
00213 }
00214
00215 private:
00216
00217
00218
00219
00220 float *xFunc, *yFunc;
00221
00222
00223
00224 int count;
00225 };
00226
00230 class IrregularDistribution1D {
00231 public:
00243 IrregularDistribution1D(float aX0, float aX1, float *aX, float *aFx, int aN) {
00244 count = aN;
00245 x0 = aX0;
00246 x1 = aX1;
00247 xFunc = new float[count];
00248 yFunc = new float[count];
00249 xCdf = new float[count+1];
00250 yCdf = new float[count+1];
00251 memcpy(xFunc, aX, count*sizeof(float));
00252 memcpy(yFunc, aFx, count*sizeof(float));
00253
00254
00255 xCdf[0] = aX0;
00256 for (int i = 1; i < count; ++i)
00257 xCdf[i] = ( xFunc[i-1] + xFunc[i] ) * .5f;
00258 xCdf[count] = aX1;
00259 yCdf[0] = 0.f;
00260 for (int i = 1; i < count+1; ++i) {
00261 yCdf[i] = yCdf[i-1] + max( 1e-3f, yFunc[i-1] ) * ( xCdf[i] - xCdf[i-1] );
00262 }
00263 funcInt = yCdf[count];
00264
00265 for (int i = 1; i < count+1; ++i)
00266 yCdf[i] /= funcInt;
00267
00268 invFuncInt = 1.f / funcInt;
00269 invCount = 1.f / count;
00270 }
00271
00272 ~IrregularDistribution1D() {
00273 delete[] xFunc;
00274 delete[] yFunc;
00275 delete[] xCdf;
00276 delete[] yCdf;
00277 }
00278
00288 float Sample(float u, float *pdf) const {
00289
00290 float *ptr = std::upper_bound(yCdf, yCdf+count+1, u);
00291 int offset = Clamp((int) (ptr-yCdf-1), 0, count - 1);
00292
00293 float du = (u - yCdf[offset]) / (yCdf[offset+1] - yCdf[offset]);
00294 *pdf = xFunc[offset] * invFuncInt;
00295 return xCdf[offset] + du * (xCdf[offset+1] - xCdf[offset]);
00296 }
00297
00305 float Eval(float x) const {
00306 if( x <= xFunc[0] )
00307 return yFunc[ 0 ];
00308 else if( x >= xFunc[ count - 1 ] )
00309 return yFunc[ count - 1 ];
00310
00311 float *ptr = std::upper_bound(xFunc, xFunc+count, x);
00312 int offset = (int) (ptr-xFunc-1);
00313
00314 float d = (x - xFunc[offset]) / (xFunc[offset+1] - xFunc[offset]);
00315
00316 return (1.f - d) * yFunc[offset] + d * yFunc[offset+1];
00317 }
00318
00327 int IndexOf(float x, float *d) const {
00328 if( x <= xFunc[0] ) {
00329 *d = 0.f;
00330 return 0;
00331 }
00332 else if( x >= xFunc[ count - 1 ] ) {
00333 *d = 0.f;
00334 return count - 1;
00335 }
00336
00337 float *ptr = std::upper_bound(xFunc, xFunc+count, x);
00338 int offset = (int) (ptr-xFunc-1);
00339
00340 *d = (x - xFunc[offset]) / (xFunc[offset+1] - xFunc[offset]);
00341 return offset;
00342 }
00343
00344
00348 float x0, x1;
00349
00350
00351
00352 float *xFunc, *yFunc;
00353
00354
00355
00356 float *xCdf, *yCdf;
00361 float funcInt, invFuncInt, invCount;
00362
00363
00364
00365 int count;
00366 };
00367
00368 }
00369
00370 #endif //LUX_MCDISTRIBUTION_H