Originally Posted by

**iMalc**
This really isn't all that hard to figure out:

Code:

switch ((int)(_nextafter(x, -999)*10)) {
case 0: t = calculation; break;
case 1: t = calculation2; break;
}

The _nextafter bit can be replaced with just x if you aren't that interested in the boundary cases afterall.

Except that this is not possible the values which are being compared to are not fixed. I am actually generating the C code in matlab, which I later plug into a dll which is being fed to a ray tracer. This raytracer is used to determine for example light distributions of sources (measured using near field gonio photometers) from which the rays start and interact with a certain geometry (reflectors) of which I have measured the BSDF (or to be more specific BRDF as I am only dealing with reflections caused by surface scattering).

Anyhow for example a distribution of a BRDF might look like Q(theta_s,phi_s,theta_i,phi_i,lambda)=cos(theta_i) In this case we assume a isotropic surface (hence phi_s is not present and takes a random value between 0 and 2*Pi), theta_i goes from 0 to Pi/2. Actually its a surface which scatters in a lambertian way.

Anyhow, with this function it is quite easy to analytically compute the cumulative function, normalise this and inverse it. But as soon as the iniital function gets more complicated we can no longer analytically compute the inverse of the normalised cumulative function. This is why I am numerically integrating the initial function, normalising this and then flipping the x - y values, thus giving the ICDF (inverse cumulative distribution function). This ICDF takes a uniformly distributed random variable between 0 and 1 (because of the normalization) and returns a value between 0 and Pi/2 representing the theta_s value according to the initial distribution function.

Now to get to where I need the if structure. Suppose I numerically calculate the ICDF using 1000 samples, then I could either store all these samples in a table and use some kind of interpolation. Right now we are assuming that we have an isotropic function, when we are dealing with data from real measurements we will of course no longer have isotropic functions. Furthermore the BRDF value will also depend on the incoming angles phi_i and theta_i along with lambda. So when dealing with real data we will have to do more then one interpolation to get our result. This is why I am fitting a cubic spline to the ICDF in which I can plug my variables and get proper results for phi_s and theta_s.

So to get to the point where I can say why I can not use fixed intervals is because I want to sample the ICDF not on the x-axis which goes from 0 to 1 using equidistant samples, I want to have equidistant samples on the y-axis. If I would not do this then my final results in the raytrace would not be accurate (it would not follow the initial distribution function) as there would be places in the ICDF which are undersampled.

An example of how part of the structure looks like right now (beware it doesnt look too nice as it is being generated by matlab):

Code:

if(x>=0 && x<0.00034947)
t = ((double)4986107.6052 * pow(x-0.0000000000e+000,3)) + ((double)-1.9448905812e+004 * pow(x-0.0000000000e+000,2)) + ((double)2.8774584596e+001*(x-0.0000000000e+000)) +0.0000000000e+000;
else if(x>=0.00034947 && x<0.0013835)
t = ((double)4.9861076052e+006 * pow(x-3.4947388154e-004,3)) + ((double)-1.4221362676e+004 * pow(x-3.4947388154e-004,2)) + ((double)1.7007705175e+001*(x-3.4947388154e-004)) +7.8934488784e-003;
else if(x>=0.0013835 && x<0.0031001)
t = ((double)-3.8408721761e+005 * pow(x-1.3835180461e-003,3)) + ((double)1.2462037430e+003 * pow(x-1.3835180461e-003,2)) + ((double)3.5908177960e+000*(x-1.3835180461e-003)) +1.5786897757e-002;
else if(x>=0.0031001 && x<0.005496)
t = ((double)9.9996059239e+004 * pow(x-3.1001367606e-003,3)) + ((double)-7.3179017428e+002 * pow(x-3.1001367606e-003,2)) + ((double)4.4738697550e+000*(x-3.1001367606e-003)) +2.3680346635e-002;
else if(x>=0.005496 && x<0.0085666)
t = .....
//somewhere in the middle of the structure the if structures look like
else if(x>=0.69757 && x<0.71036)
t = ((double)1.3754739284e+000 * pow(x-6.9756733940e-001,3)) + ((double)5.3877102783e-001 * pow(x-6.9756733940e-001,2)) + ((double)6.0992720295e-001*(x-6.9756733940e-001)) +4.5782003495e-001;
else if(x>=0.71036 && x<0.72285)
t = ((double)1.5444259520e+000 * pow(x-7.1035969012e-001,3)) + ((double)5.9155766253e-001 * pow(x-7.1035969012e-001,2)) + ((double)6.2438676399e-001*(x-7.1035969012e-001)) +4.6571348382e-001;
else if(x>=0.72285 && x<0.73503)
t = ....

I am also attaching some images to claridy things visually. The initial distribution is being showed (cos(theta)*sin(theta) - as we are in spherical coordinates and need to take the Jacobian into account, the radius is 1 by the way).

You can also see the ICDF (which is already the cubic spline) which exists by taking 10 and 100 samples from the numerically computed cumulative normalised distribution function, inverse these values and fit 10 and 100 cubic splines. Also attached are the histograms of the values which are returned when feeding 10^6 uniformly distributed random variables between 0 and 1. It is clear that 10 samples are not enough, 100 seem to work, even though when used in a real raytracing simulation it seems that it is still not enough (250 samples do the trick though).