/* **************************************************************** * D.c * **************************************************************** * MODULE PURPOSE: * This module contains microfacet distribution routines and * various support routines. One of support functions includes * computing the coefficients corresponding to beta for each * function (as described in Section 4.3.1-Improved Specular * Reflection). Beta is the angle between H and N where the * function is equal to half the value at N = H. Beta is in * radians. * * MODULE CONTENTS: * D_phong_init - compute Ns given beta * D_phong - Phong cosine power function * D_blinn_init - compute Ns given beta * D_blinn - Blinn cosine power function * D_gaussian_init - compute C1 given beta * D_gaussian - Gaussian distribution * D_reitz_init - compute c2 given beta * D_reitz - Trowbridge-Reitz * D_cook_init - compute m given beta * D_cook - Cook model * D_g - compute apparent roughness * D_tau - compute correlation distance * D_sigma - compute rms roughness * D_m - compute slope * D_coherent - coherent roughness attenuation * (per Beckmann) * D_incoherent - incoherent microfacet attenuation * (per Beckmann and Bahar) * D_Vxy - Vxy used in D_incoherent() * * ASSUMPTIONS: * > The following are defined in math.h * HUGE largest floating pont value * M_PI pi * M_PI_2 pi/2 * M_PI_4 pi/4 * M_SQRT2 root of 2 * * > The direction vectors N, V, and L are assumed to be oriented * to the same side of the surface. */ #include #include #include "geo.h" #include "D.h" /* **************************************************************** * double D_phong_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns Ns given beta for the Phong (1975) specular function, * per footnote to Eq.(\*(rH) */ double D_phong_init (beta) double beta; { if (beta <= 0.0) return HUGE; if (beta >= M_PI_4) return 0.0; return -(log(2.0) / log(cos(2.0 * beta))); } /* **************************************************************** * double D_phong (N, L, V, Ns) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double Ns (in) - specular exponent * * Returns the value of the Phong (1975) specular function given * surface normal, incident light direction, view direction, and * the specular exponent. Refer to Eqs.(\*(rU) and (\*(rI). */ double D_phong (N, L, V, Ns) DIR_VECT *N, *L, *V; double Ns; { double Rv_dot_L; if ((Rv_dot_L = geo_dot(geo_rfl(V,N),L)) < 0.0) return 0.0; return pow (Rv_dot_L, Ns); } /* **************************************************************** * double D_blinn_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns Ns given beta for the cosine power distribution * function presented by Blinn (1977). Refer to Eq.(\*(rHa) */ double D_blinn_init (beta) double beta; { if (beta <= 0.0) return HUGE; if (beta >= M_PI_2) return 0.0; return -(log(2.0) / log(cos(beta))); } /* **************************************************************** * double D_blinn (N, L, V, Ns) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double Ns (in) - specular exponent * * Returns the cosine power distribution function presented by * Blinn (1977) given surface normal, incident light direction, * view direction, and the specular exponent. Refer to Eq.(\*(e5) */ double D_blinn (N, L, V, Ns) DIR_VECT *N, *L, *V; double Ns; { double N_dot_H; if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0; return pow (N_dot_H, Ns); } /* **************************************************************** * double D_gaussian_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns C1 given beta for the Gaussian distribution presented * by Blinn (1977). Refer to Eq.(\*(rHb) */ double D_gaussian_init (beta) double beta; { if (beta <= 0.0) return HUGE; return sqrt(log(2.0)) / beta; } /* **************************************************************** * double D_gaussian (N, L, V, C1) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double C1 (in) - shape coefficient * * Returns the Gaussian distribution function presented by * Blinn (1977) given surface normal, incident light direction, * view direction, and the specular exponent. Refer to Eq.(\*(eA) */ double D_gaussian (N, L, V, C1) DIR_VECT *N, *L, *V; double C1; { double tmp; double N_dot_H; if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0; tmp = acos(N_dot_H) * C1; return exp (-(tmp * tmp)); } /* **************************************************************** * double D_reitz_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns C2 squared given beta for the Trowbridge-Reitz * distribution function presented by Blinn (1977). Refer to * Eq.(\*(rHc). C2 is squared for computational efficiency later. */ double D_reitz_init (beta) double beta; { double cos_beta; if (beta <= 0.0) return 0.0; cos_beta = cos(beta); return ((cos_beta * cos_beta) - 1.0) / ((cos_beta * cos_beta) - sqrt(2.0)); } /* **************************************************************** * double D_reitz (N, L, V, C2_2) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double C2_2 (in) - C2 squared * * Returns the Trowbridge-Reitz distribution function presented * by Blinn (1977) given surface normal, incident light * direction, view direction, and the specular exponent. Refer * to Eq.(\*(eB) */ double D_reitz (N, L, V, C2_2) DIR_VECT *N, *L, *V; double C2_2; { double tmp; double N_dot_H; if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0; tmp = C2_2 / ((N_dot_H * N_dot_H * (C2_2 - 1.0)) + 1.0); return tmp * tmp; } /* **************************************************************** * double D_cook_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns m squared corresponding to beta, refer to Eq.(\*(rT). * m is squared for computational efficiency in later use. */ double D_cook_init (beta) double beta; { double tan_beta; if (beta <= 0.0) return 0.0; if (beta >= M_PI_2) return HUGE; tan_beta = tan(beta); return -(tan_beta * tan_beta) / log(pow(cos(beta),4.0)/2.0); } /* **************************************************************** * double D_cook (N, L, V, m_2) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double m_2 (in) - m squared * * Returns microfacet distribution probability (Cook 1982) derived * from (Beckmann 1963). Refer to Eq.(\*(rX). */ double D_cook (N, L, V, m_2) DIR_VECT *N, *L, *V; double m_2; { double tmp; double N_dot_H; if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0; tmp = -(1.0 - (N_dot_H * N_dot_H)) / (N_dot_H * N_dot_H * m_2); return exp(tmp)/(4.0 * M_PI * m_2 * pow(N_dot_H,4.0)); } /* **************************************************************** * double D_g (N, L, V, sigma, lambda) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double sigma (in) RMS roughness, nm * double lambda (in) wavelength, nm * * Returns the apparent roughness, g, as given by Beckmann (1963). * Refer to Eq.(\*(rY). */ double D_g (N, L, V, sigma, lambda) DIR_VECT *N, *L, *V; double sigma, lambda; { double tmp; tmp = geo_dot(N,L) + geo_dot(N,V); return (4.0 * M_PI * M_PI * sigma * sigma * tmp * tmp) / (lambda * lambda); } /* **************************************************************** * double D_tau (sigma, m) * double sigma (in) rms roughness (nm) * double m (in) rms slope * * Returns tau (correlation distance) in nm. The relationship of * roughness, slope, and correlation distance given by * Beckmann (1963) is assumed, refer to Eq.(\*(rD). */ double D_tau (sigma, m) double sigma, m; { return (2.0 * sigma) / m; } /* **************************************************************** * double D_sigma (m, tau) * double m (in) rms slope * double tau (in) correlation distance (nm) * * Returns sigma (rms roughness). The relationship of roughness, * slope, and correlation distance given by Beckmann (1963) is * assumed, refer to Eq.(\*(rD). */ double D_sigma (m, tau) double m, tau; { return (m * tau)/ 2.0; } /* **************************************************************** * double D_m (sigma, tau) * double sigma (in) rms roughness (nm) * double tau (in) correlation distance (nm) * * Returns m (rms slope). The relationship of roughness, slope, * and correlation distance given by Beckmann (1963) is assumed, * refer to Eq.(\*(rD). */ double D_m (sigma, tau) double sigma, tau; { return (2.0 * sigma) / tau; } /* **************************************************************** * double D_coherent (g) * double g (in) apparent roughness * * Returns the coherent roughness attenuation given * by Beckmann (1963), refer to Eq.(\*(rZ). */ double D_coherent (g) double g; { return exp(-g); } /* **************************************************************** * double D_incoherent (N, L, V, m, g, lambda, tau) * DIR_VECT *N, *L, *V (in) N, L, V vectors * double m (in) m (slope) * double g (in) apparent roughness * double lambda (in) wavelength (nm) * double tau (in) correlation distance (nm) * * Returns the microfacet distribution function given by * Beckmann (1963) as interpreted by Koestner (1986), refer to * Eq.(\*(rC). */ double D_incoherent (N, L, V, m, g, lambda, tau) DIR_VECT *N, *L, *V; double m, g, lambda, tau; { DIR_VECT *H; double D1, D2; double denom; double N_dot_H, N_dot_H_2; if (g <= 0.0) return 0.0; if ((H = geo_H(V,L)) == NULL) return 0.0; if ((N_dot_H = geo_dot(N,H)) <= 0.0) return 0.0; N_dot_H_2 = N_dot_H * N_dot_H; denom = 4.0 * M_PI * m * m * N_dot_H_2 * N_dot_H_2; /* if g < 8.0, evaluate the series expansion, Eq.(\*(rCb) * terminating when a term falls below 0.00001. For small g * the evaluation terminates quickly. For large g the series * converges slowly. If g is near 8.0 the first terms of the * series will be less that the termination criteria. This * requires an additional termination criteria that the values * of the terms are decreasing at the time the termination * criteria is reached. */ if (g < 8.0) { double inc, ct, ct_factorial, g_pow; double last_inc, Vxz_t_over_4; Vxz_t_over_4 = (tau * tau * D_Vxz(N,L,V,lambda)) / 4.0; D1 = 0.0; ct = 1.0; ct_factorial = 1.0; g_pow = g * g; inc = 0.0; do { last_inc = inc; inc = (g_pow / (ct * ct_factorial)) * exp(-(g + Vxz_t_over_4/ct)); D1 += inc; ct += 1.0; ct_factorial *= ct; g_pow *= g; } while (((inc/denom) > 0.00001) || (inc > last_inc)); D1 /= denom; /* if g < 5.0, only the series expansion is used. If * 5.0 < g < 8.0, then we are in a transition range * for interpolation between the series expansion and * the convergence expression. */ if (g < 5.0) return D1; } /* if g > 5.0, evaluate the convergence expression, Eq.(\*(rCc) * (this routine would have returned earlier if G < 5.0). */ { double tmp; tmp = (N_dot_H_2 - 1.0) / (N_dot_H_2 * m * m); D2 = exp(tmp) / denom; /* if g > 8.0, only the convergence expression is used, * otherwise we are in a transition zone between the * convergent expression and series expansion. */ if (g > 8.0) return D2; } /* linear interpolation between D1 and D2 for * 5.0 < g < 8.0 */ return (((8.0 - g) * D1) + ((g - 5.0) * D2)) / 3.0; } /* **************************************************************** * double D_Vxz (N, L, V, lambda) * DIR_VECT *N, *L, *V (in) N, L, V vectors * double lambda (in) wavelength (nm) * * Returns the value of Vxz per Eq.(\*(rCd). */ double D_Vxz (N, L, V, lambda) DIR_VECT *N, *L, *V; double lambda; { double N_dot_L, N_dot_V, cos_phi; DIR_VECT sin_i, sin_r; double len_2_i, len_2_r; /* compute the sine squared of the incident angle */ N_dot_L = geo_dot(N, L); len_2_i = 1.0 - (N_dot_L * N_dot_L); /* compute the sine squared of the reflected angle */ N_dot_V = geo_dot(N, V); len_2_r = 1.0 - (N_dot_V * N_dot_V); /* if the sine of either the incident or the * reflected angle is zero, then the middle * term is zero. */ if ((len_2_i > 0.0) && (len_2_r > 0.0)) { /* the two sine vectors are of length equal to the sine of * the angles. The dot product is the cosine times the * lengths of the vectors (sines of the angles). */ sin_i.i = L->i - (N_dot_L * N->i); sin_i.j = L->j - (N_dot_L * N->j); sin_i.k = L->k - (N_dot_L * N->k); sin_r.i = V->i - (N_dot_V * N->i); sin_r.j = V->j - (N_dot_V * N->j); sin_r.k = V->k - (N_dot_V * N->k); cos_phi = geo_dot (&sin_i, &sin_r); return (4.0 * M_PI * M_PI * (len_2_i + (2.0 * cos_phi) + len_2_r)) / (lambda * lambda); } else { return (4.0 * M_PI * M_PI * (len_2_i + len_2_r)) / (lambda * lambda); } } /* ************************************************************* */