float G( vector V, M, N; float alpha ){
float a = 1 / ( alpha * tan( acos(V.N) ) );
float x = step( 0, (V.M) / (V.N) );
if( a < 1.6 ){
return x*(3.535*a+2.181*a*a)/(1+2.276*a+2.577*a*a);
} else {
return x*1;
}
}
float sqr( float x ){ return x*x; }
extern normal N;
extern vector I;
extern point P;
extern varying vector dPdu;
normal NN = normalize( N );
vector V = normalize( -I );
color C = 0;
color trs_reflected = 1;
uniform float i;
for( i=0; i<Samples; i+= 1 )
{
uniform float a = ceil(sqrt(Samples));
float x1 = (floor(i/a)+random())/a;
float x2 = (mod(i,a)+random())/a;
float theta = atan( sqrt( -sqr( Roughness )*log( 1-x1 ) ) );
float phi = 2*PI*x2;
matrix rot = rotate( matrix 1, theta, dPdu );
rot = rotate( rot, phi, NN );
vector MM = normalize( transform( rot, NN ) );
float Kr, Kt;
vector R, T;
fresnel( -V, MM, 1/IoR, Kr, Kt, R, T );
R = normalize(R);
float Gi = G( V, MM, NN, Roughness );
float Go = G( R, MM, NN, Roughness );
float weight = abs(V.MM) * Gi * Go * Kr / ( abs(V.NN) * abs(MM.NN) );
C += trace( P, R, "transmission", trs_reflected ) * weight;
}
C *= ReflectionStrength / Samples;
output_reflection = C;