Shaders

// Compulsory fragment shader: initial chunk
uniform float width; 
uniform float height; 
uniform vec3 ambient_color;
uniform vec3 int_diff_color;
uniform vec3 ext_diff_color;
uniform vec3 specular_color;
uniform vec3 sing_color;
uniform float shininess;
//uniform float opacity;
uniform float transparency;
uniform float reflection;
uniform float curv_scale;
uniform int   princ_curv;
uniform vec3 light1_direction;
uniform vec3 light1_color;
uniform vec3 light2_direction;
uniform vec3 light2_color;
uniform int   bg_mode;
uniform vec3  bg_color;
uniform float singularity;
uniform float accuracy;
uniform float iso;
uniform float a;
uniform float b;
uniform float c;      
uniform float scale;
uniform float center_x;
uniform float center_y;
uniform float center_z;
uniform float grid_size;
uniform float grid_thickness;
uniform float grid_attenuation;
uniform int   shape;
uniform int   color_model;
uniform float camera_distance;
uniform float view_radius;
uniform float time;
uniform samplerCube cubeMap;

// constants
const float exterior_index = 1.0;

// new uniforms
uniform int   nb_bissections; //< default to 8
uniform float interior_index;
uniform int   nb_reflections;
uniform int   nb_refractions;
uniform vec3  dir_color_on;
uniform vec3  dir_color_off;
uniform float dir_on_opacity;
uniform float dir_off_opacity;
    in vec3 vertPos;
    in vec3 projPos;
    in mat4 projViewInv;
    in mat4 myModelMatrix;
    in mat4 myModelViewMatrix;
    in vec3 surfNormal;
out vec4 outColor;
// // 3D function in cube
// float vwidth  = view_radius; // cube width  
// float vwidth2 = vwidth*vwidth; // cube width

// Estimates the Hessian of f at position p
mat3 hessf( vec3 p ) {
    mat3 H;
    float eps  = 0.001 * scale; // 0.1 / scale;
    float eps2 = eps * eps;
    // numerical centered derivatives
    vec3  dx   = vec3( eps, 0.0, 0.0 );
    vec3  dy   = vec3( 0.0, eps, 0.0 );
    vec3  dz   = vec3( 0.0, 0.0, eps );
    float fp   = f( p );
    float fppdx= f( p+dx );
    float fpmdx= f( p-dx );        
    float fppdy= f( p+dy );
    float fpmdy= f( p-dy );        
    float fppdz= f( p+dz );
    float fpmdz= f( p-dz );        
    H[ 0 ][ 0 ] = ( fppdx - 2.*fp + fpmdx ) / eps2;
    H[ 1 ][ 1 ] = ( fppdy - 2.*fp + fpmdy ) / eps2;
    H[ 2 ][ 2 ] = ( fppdz - 2.*fp + fpmdz ) / eps2;
    H[ 0 ][ 1 ] = ( ( f(p+dx+dy) - f(p-dx+dy) ) - ( f(p+dx-dy) - f(p-dx-dy) ) ) / (4.*eps2);
    H[ 0 ][ 2 ] = ( ( f(p+dx+dz) - f(p-dx+dz) ) - ( f(p+dx-dz) - f(p-dx-dz) ) ) / (4.*eps2);
    H[ 1 ][ 2 ] = ( ( f(p+dy+dz) - f(p-dy+dz) ) - ( f(p+dy-dz) - f(p-dy-dz) ) ) / (4.*eps2);
    H[ 1 ][ 0 ] = H[ 0 ][ 1 ];
    H[ 2 ][ 0 ] = H[ 0 ][ 2 ];
    H[ 2 ][ 1 ] = H[ 1 ][ 2 ];        
    return H;
}

// Estimates the curvatures of f at position p
vec2 mean_gaussian_curvatures( vec3 p ) {
    vec3  G = gradf( p );
    float g = length( G );
    mat3  H = hessf( p );
    mat4  C = mat4( H );
//    for ( int i = 0; i < 3; i ++ )
//	for ( int j = 0; j < 3; j ++ )
//	    C[ i ][ j ] = H[ i ][ j ];
    C[0][3] = G[ 0 ];
    C[1][3] = G[ 1 ];
    C[2][3] = G[ 2 ];
    C[3][0] = G[ 0 ];
    C[3][1] = G[ 1 ];
    C[3][2] = G[ 2 ];
    C[3][3] = 0.0;
    vec2 curv_hg;
    curv_hg.y = - determinant( C ) / (g*g*g*g);
    curv_hg.x = - 0.5*( dot( H*G, G ) - g * g * (H[0][0]+H[1][1]+H[2][2]) ) / (g*g*g);
    return curv_hg;
}
// Estimates the principal curvatures of f at position p
// k1 : first principal curvature (highest value k1 >= k2)
// k2 : second principal curvature (lowest value k2 <= k1)
void principal_curvatures( in vec3 p, out float k1, out float k2 ) {
    vec2 curv_hg  = mean_gaussian_curvatures( p );
    float dk = sqrt( max( 0.0, curv_hg.x*curv_hg.x - curv_hg.y ) );
    k2       = curv_hg.x - dk;
    k1       = curv_hg.x + dk;
}

// Estimates the principal directions of curvatures of f at position p
/*
 * Copyright 2015 
 * Hélène Perrier <helene.perrier@liris.cnrs.fr>
 * Jérémy Levallois <jeremy.levallois@liris.cnrs.fr>
 * David Coeurjolly <david.coeurjolly@liris.cnrs.fr>
 * Jacques-Olivier Lachaud <jacques-olivier.lachaud@univ-savoie.fr>
 * Jean-Philippe Farrugia <jean-philippe.farrugia@liris.cnrs.fr>
 * Jean-Claude Iehl <jean-claude.iehl@liris.cnrs.fr>
 * 
 * This file is part of ICTV.
 * 
 * ICTV is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * ICTV is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with ICTV.  If not, see <http://www.gnu.org/licenses/>
 */
void getEigenValuesVectors( in mat3 mat_data,
			    out mat3 vectors,
			    out vec3 values )
{ 
    vec3 e = vec3(0);
    int dimension = 3;
    int dimensionMinusOne = 2;
    for( int j = 0; j < dimension; ++j )
	values[ j ] =  mat_data[dimensionMinusOne][ j ];
    // Householder reduction to tridiagonal form.
    for( int i = dimensionMinusOne; i > 0 && i <= dimensionMinusOne; --i ) {
	// Scale to avoid under/overflow.
	float scale = 0.0;
	float h =  0.0;
	for( int k = 0; k < i; ++k )
            scale += abs( values[ k ] );
	if ( scale ==  0.0 ) {
            e[ i ] = values[ i - 1 ];
            for( int j = 0; j < i; ++j ) {
		values[ j ] = mat_data[ ( i - 1 ) ] [  j  ];
		mat_data[ i ][ j ] = 0.0;
		mat_data[ j ][ i ] = 0.0;
            }
	} else {
            // Generate Householder vector.
            for ( int k = 0; k < i; ++k ) {
		values[ k ] /= scale;
		h += values[ k ] * values[ k ];
            }
            float f = values[ i - 1 ];
            float g = sqrt( h );
            if ( f >  0.0 ) g = -g;
            e[ i ] = scale * g;
            h -= f * g;
            values[ i - 1 ] = f - g;
            for ( int j = 0; j < i; ++j) e[ j ] =  0.0;
            // Apply similarity transformation to remaining columns.
            for ( int j = 0; j < i; ++j ) {
		f = values[ j ];
		mat_data[ j ][ i ] = f;
		g = e[ j ] +  mat_data[ j ][ j ] * f;
		for ( int k = j + 1; k <= i - 1; ++k ) {
		    g +=  mat_data[ k ][ j ] * values[ k ];
		    e[ k ] +=  mat_data[ k ][ j ] * f;
		}
		e[ j ] = g;
            }
            f = 0.0;
            for ( int j = 0; j < i; ++j ) {
		e[ j ] /= h;
		f += e[ j ] * values[ j ];
            }
            float hh = f / ( h + h );
            for ( int j = 0; j < i; ++j )
		e[ j ] -= hh * values[ j ];
            for ( int j = 0; j < i; ++j ) {
		f = values[ j ];
		g = e[ j ];
		for ( int k = j; k <= i - 1; ++k )
                    mat_data[ k ][ j ] = mat_data[ k ][ j ] - (f * e[ k ] + g * values[ k ]);
		values[ j ] =  mat_data[ i - 1][ j ];
		mat_data[ i ][ j ] = 0.0;
            }
	}
	values[ i ] = h;
    }
    // Accumulate transformations.
    for ( int i = 0; i < dimensionMinusOne; ++i ) {
	mat_data[dimensionMinusOne][ i ] =  mat_data[ i ][ i ];
	mat_data[ i ][ i ] = 1.0;
	float h = values[ i + 1 ];
	if ( h != 0.0 ) {
            for ( int k = 0; k <= i; ++k )
		values[ k ] =  mat_data[ k ][ i + 1 ] / h;
            for ( int j = 0; j <= i; ++j ) {
		float g = 0.0;
		for ( int k = 0; k <= i; ++k )
		    g +=  mat_data[ k ][ i + 1 ] *  mat_data[ k ][ j ];
		for ( int k = 0; k <= i; ++k )
                    mat_data[ k ][ j ] =   mat_data[k][ j ] - ( g * values[ k ] );
            }
	}
	for ( int k = 0; k <= i; ++k )
            mat_data[ k ][ i + 1 ] =  0.0;
    }
    for ( int j = 0; j < dimension; ++j ) {
        values[ j ] =  mat_data[ dimensionMinusOne ][ j ];
	mat_data[ dimensionMinusOne ][ j ] = 0.0;
    }
    mat_data[ dimensionMinusOne ][ dimensionMinusOne ] =  1.0;
    e[ 0 ] =  0.0;
    for ( int i = 1; i < dimension; ++i )
	e[ i - 1 ] = e[ i ];
    e[ dimensionMinusOne ] = 0.0;
    float f = float( 0.0 );
    float tst1 = float( 0.0 );
    float eps = float( pow( 2.0, -52.0 ));
    //float eps = float( pow( 2.0, -52.0 ));
    for( int l = 0; l < dimension; ++l ) {
	// Find small subdiagonal element
	tst1 = float( max( tst1, abs ( values[ l ] ) + abs( e[ l ] )));
	int m = l;
	while ( m < dimension ) {
          if ( abs ( e[ m ] ) <= eps * tst1 ) break;
            ++m;
        }
	// If m == l, d[l] is an eigenvalue,
	// otherwise, iterate.
	if( m > l && l<2 ) {
            int iter = 0;
            do {
		++iter;  // (Could check iteration count here.)
		// Compute implicit shift
		float g = values[ l ];
		float p = ( values[ l + 1 ] - g ) / ( float( 2.0 ) * e[ l ] );
		float r = float( sqrt ( p * p + float( 1.0 ) * float( 1.0 )));
		if ( p < 0.0 ) r = -r;
		values[ l ] = e[ l ] / ( p + r );
		values[ l + 1 ] = e[ l ] * ( p + r );
		float dl1 = values[ l + 1 ];
		float h = g - values[ l ];
		for( int i = l + 2; i < dimension; ++i )
                    values[ i ] -= h;
		f = f + h;
		// Implicit QL transformation.
		p = values[ m ];
		float c = float( 1.0 );
		float c2 = c;
		float c3 = c;
		float el1 = e[ l + 1 ];
		float s = float( 0.0 );
		float s2 = float( 0.0 );
		for ( int i = m - 1; i >= l && i <= m - 1; --i ) {
                    c3 = c2;
                    c2 = c;
                    s2 = s;
                    g = c * e[ i ];
                    h = c * p;
                    r = float( sqrt ( p * p + e[ i ] * e[ i ] ));
                    e[ i + 1 ] = s * r;
                    s = e[ i ] / r;
                    c = p / r;
                    p = c * values[ i ] - s * g;
                    values[ i + 1 ] = h + s * ( c * g + s * values[ i ] );
                    // Accumulate transformation.
                    for( int k = 0; k < dimension; ++k ) {
			h =  mat_data[ k ][ i + 1 ];
			mat_data[ k ][ i + 1 ] =  ( s *  mat_data[ k ][ i ] + c * h );
			mat_data[ k ][ i ] = ( c *  mat_data[ k ][ i ] - s * h );
                    }
                }
		p = - s * s2 * c3 * el1 * e[ l ] / dl1;
		e[ l ] = s * p;
		values[ l ] = c * p;
		// Check for convergence.
            } while ( abs ( e[ l ] ) > eps * tst1 && iter < 30);
        }
	values[ l ] = values[ l ] + f;
	e[ l ] = float( 0.0 );
    }
    // Sort eigenvalues and corresponding vectors.
    for ( int i = 0; i < dimensionMinusOne; ++i ) {
	int k = i;
	float p = values[ i ];
	for ( int j = i + 1; j < dimension; ++j ) {
            if ( values[ j ] < p ) {
		k = j;
		p = values[ j ];
            }
        }
	if ( k != i ) {
            values[ k ] = values[ i ];
            values[ i ] = p;
            for ( int j = 0; j < dimension; ++j ) {
		p =  mat_data[ j ][ i ];
		mat_data[ j ][ i ] =  mat_data[ j ][ k ];
		mat_data[ j ][ k ] = p;
            }
        }
    }
    for(int i=0; i<3; i++)
	for(int j=0; j<3; j++)
	    vectors[i][j] = mat_data[i][j];
}

// Compulsory fragment shader

// Ray structure
struct Ray {
    vec3 pos;
    vec3 dir;
};

// @return the step of the ray tracing
float ray_trace_step()
{
    return 0.1 * accuracy / scale;
}
// @return the maximal number of steps of the ray tracing
int   ray_trace_n()
{
    return int( ceil( 20.0 * view_radius / accuracy ) );
}

// @return 'true' iff the ray is coming to the surface from the exterior
bool ray_is_from_exterior( in Ray ray, in vec3 n )
{
    return dot( ray.dir, n ) <= 0.0;
}

// Approximates the solution f(x)=0 by bissection around two positions.
// Invariant: f(p) < iso <= f(q)
vec3 bissection( in vec3 p, in vec3 q )
{   // 10: precision is 0.001
    // 8: precision is 0.004
    for ( int i = 0; i < nb_bissections; i++ )
    {
	vec3 m = 0.5*(p + q);
	float fm = f( m );
	if ( fm < iso ) p = m;
	else            q = m;
    }
    return 0.5*(p + q);
}

// Roughly approximates a position p, such that f(p)=iso.
// In practice, finds the first change of sign.
vec3 ray_trace( in Ray ray, out bool found )
{
    vec3 d = ray_trace_step() * ray.dir;
    vec3 p = ray.pos;
    vec3 q = p + d;
    int  n = ray_trace_n();
    found  = true;
    if ( f( p ) < iso ) {
	for ( int i = 0; i < n; i++ ) {
	    if ( f( q ) >= iso ) return bissection( q - d, q );
	    q  += d;
	}
    } else {
	for ( int i = 0; i < n; i++ ) {
	    if ( f( q ) <  iso ) return bissection( q, q - d );
	    q  += d;
	}
    }
    found  = false;
    // Did not find intersection, return a point outside the view shape.
    return q + normalize( d ) * 2. * view_radius;
}

// @return the reflected ray at given position according to the given normal
Ray ray_reflection( in Ray ray, in vec3 n )
{
    float eps = 0.01*ray_trace_step();
    Ray rray  = Ray( ray.pos, reflect( ray.dir, n ) );
    rray.pos += rray.dir * eps;
    return rray;
}

// @return the refracted ray at given position according to the given normal
// nb: depends on int/ext and uniform interior_index (exterior_index is 1.0)
Ray ray_refraction( in Ray ray, in vec3 n )
{
    float eps    = 0.01*ray_trace_step();
    float ct1    = dot( -ray.dir, n );
    bool  ext    = ct1 >= 0.0;  // is incident ray coming from exterior
    float ratio  = ext ? 1.0 / interior_index : interior_index; //< n1/n2
    // Ray   rray;
    // rray.dir     = refract( ray.dir, ext ? n : -n, ratio );
    // rray.pos     = ray.pos + rray.dir * eps;
    // return rray;
    float ct2_sq = 1.0 - ratio*ratio*( 1.0 - ct1 * ct1 );
    float ct2    = ct2_sq >= 0.0 ? sqrt( ct2_sq ) : - sqrt( -ct2_sq );
    Ray   rray;
    rray.dir     = normalize( ratio * ray.dir + ( ratio * ct1 + ( ext ? -ct2 : ct2 ) ) * n );
    rray.pos     = ray.pos + rray.dir * eps;
    return rray;
}

// @return the ray coming from the camera and starting from the surface of the view shape
// if the observer is inside the view shape, the starting point is moved backward.
// @pre use uniforms camera_distance, projViewInv, surfNormal, scale
Ray ray_from_camera() {
    // camera is behind
    vec3 focal  = vec3( 0.0, 0.0, -camera_distance );
    vec4 viewp4 = projViewInv * vec4( focal, 1.0 );
    vec3 viewp  = normalize( vertPos - viewp4.xyz / viewp4.w );
    vec3 tracePos = ( dot( surfNormal, focal ) > 0.0 )
        ? ( vertPos - 2.0*dot( vertPos, viewp ) * viewp ) / scale
        : vertPos / scale;
    return Ray( tracePos, viewp );
}

// Computes the Phong illumination model
// @param ray      the incoming ray from the observer
// @param color    the material color
// @param normal   the normal to the surface (unoriented unit vector)
// @param lightDir the direction toward the light
// @param colDir   the color of the light
// @return the color as given by the Phong illumination model
// @pre uses uniform specular_color
//
//     ray       ^ normal  _> lightDir
//       \__     |      __/
//          \__  |  ___/
//             \>|_/
// --------------------------
//phongModel( ray, c, n, vec3(0.0,0.0,-1.0), vec3(1.,1.,1.) );
vec3 phong_model( in Ray ray,
		  in vec3 color,    in vec3 normal,
		  in vec3 lightDir, in vec3 lightCol )
{
    float sirradiance = dot( lightDir, normal );
    float irradiance  = abs( sirradiance ); //max(dot(lightDir, normal), 0.0);
    vec3  reflectDir  = reflect(-lightDir, normal);
    float  specDot    = max(dot(reflectDir, -ray.dir ), 0.0);
    color            += pow(specDot, shininess) * specular_color;
    color            *= lightCol * irradiance;
    return color;
}

// Compute a background color from a z in [-1,1]
vec3 background_color( float z ) {
    //if ( bg_mode == 0 ) return bg_color;
    z = clamp( z, -1., 1. );
    vec3 c0 = vec3( 1.0,0.0,0.0);
    vec3 c1 = vec3( 1.0,1.0,0.0);
    vec3 c2 = vec3( 1.0,1.0,1.0);
    vec3 c3 = vec3( 0.0,1.0,1.0);
    vec3 c4 = vec3( 0.0,0.0,1.0);
    vec3 c5 = vec3( 0.0,0.0,0.0);
    if ( z < -0.5 )
	return (1. - 2.*(z+1.)) * c0 + 2.*(z+1.) * c1;
    else if ( z < 0. )
	return (1. - 2.*(z+0.5)) * c1 + 2.*(z+0.5) * c2;
    else if ( z < 0.5 )
	return (1. - 2.*z) * c2 + 2.*z * c3;
    else if ( z < 0.75 )
	return (1. - 4.*(z-0.5)) * c3 + 4.*(z-0.5) * c4;
    else 
	return (1. - 4.*(z-0.75)) * c4 + 4.*(z-0.75) * c5;
}

vec3 background_color( vec3 dir ) {
    // static backgrounds
    vec3 vdir  = (myModelMatrix * vec4( dir.x, dir.y, dir.z, 0.0 )).xyz; 
    if ( bg_mode == 0 ) return bg_color;
    else if ( bg_mode == 1 ) {
	//vec3 c     = textureCube( cubeMap, vec3( -vdir.x, vdir.y, vdir.z) ).rgb;
	//	float z       = 2.*gl_FragCoord.y / height - 1.0;
	return background_color( vdir.y );
    } else if ( bg_mode == 2 ) {
	if ( vdir.y >= 0.0 ) return background_color( vdir.y );
	float x = -0.5 * vdir.x / vdir.y;
	float y = -0.5 * vdir.z / vdir.y;
	float d = sqrt( x*x + y*y );
	float t = min( d, 30.0 ) / 30.0;
	x -= floor( x );
	y -= floor( y );
	if ( ( ( x >= 0.5 ) && ( y >= 0.5 ) ) || ( ( x < 0.5 ) && ( y < 0.5 ) ) )
	    return (1.0 - t)*vec3( 0.2, 0.2, 0.2 ) + t * vec3( 1.0, 1.0, 1.0 );
	else
	    return (1.0 - t)*vec3( 0.4, 0.4, 0.4 ) + t * vec3( 1.0, 1.0, 1.0 );
    } else {
	// dynamic background
	// Take into account the possible rotation of the object.
	vec3 c     = textureCube( cubeMap, vec3( -vdir.x, vdir.y, vdir.z) ).rgb;
	// Gamma correction needed for background photo.
	const float gamma = 1./2.2; 
	return vec3( pow( c.r, gamma ), pow( c.g, gamma ), pow( c.b, gamma ) );
    }
}

// Display convex-concave zones
vec4 display_color( in Ray ray, in vec3 g, in vec3 n )
{
    vec3 pos = ray.pos;
    const vec3 white    = vec3( 1.0, 1.0, 1.0 );
    const vec3 red      = vec3( 1.0, 0.0, 0.0 );
    const vec3 green    = vec3( 0.0, 1.0, 0.0 );
    const vec3 blue     = vec3( 0.0, 0.0, 1.0 );
    vec2 curv_hg  = mean_gaussian_curvatures( pos );
    float dk = sqrt( max( 0.0, curv_hg.x*curv_hg.x - curv_hg.y ) );
    float k1 = curv_hg.x - dk;
    float k2 = curv_hg.x + dk;
    float force = min( max( abs(k1), abs(k2) ) / curv_scale, 1.0 );
    float theta = atan( k2, k1 );
    vec3 radiance = (1.0 - force) * white;
    const float pi_2  = 1.57079632679;
    const float pi    = 2.0*pi_2;
    const float pi_4  = 0.5*pi_2;
    const float pi3_4 = 3.0*pi_4;
    //const float band = 1.57079632679 / 100.0;
    if ( theta <  0.0 ) theta += 4.0*pi_2;
    if ( theta <= pi_2 ) {
	float l   = (theta - pi_4) / pi_4;
	radiance += force * vec3( 1.0, l, 0.0 );
    } else if ( theta <= pi3_4 ) {
	float l   = (theta - pi_2) / pi_4;
	radiance += force * vec3( 1.0-l, 1.0, 0.0 );
    } else if ( theta <= pi ) {
	float l   = (theta - pi3_4) / pi_4;
	radiance += force * vec3( 0.0, 1.0, l );
    } else {
	float l   = (theta - pi) / pi_4;
	radiance += force * vec3( 0.0, 1.0 - l, 1.0 );
    }
#if GRID_MODE == 3
    float  bound  = grid_size * grid_thickness * grid_thickness;
    if ( mod( force, grid_size ) < bound )
        radiance *= grid_attenuation;
#endif
    return vec4( radiance, 1.0 );
}

// Display principal directions
vec4 display_color( in Ray ray, in vec3 g, in vec3 n )
{
    const float pi_2  = 1.57079632679;
    const float pi    = 2.0*pi_2;
    vec3  d1; // first principal direction (highest curvature)
    float k1; // highest curvature
    vec3  d2; // second principal direction (lowest curvature)
    float k2; // lowest curvature
    vec3  d3;
    vec3  pos = ray.pos;
    principal_directions( pos, d1, d2, d3 );
    principal_curvatures( pos, k1, k2 );
    float curv1 = min( abs( 10.0 * k1 / curv_scale ), 1.0 );
    float curv2 = min( abs( 10.0 * k2 / curv_scale ), 1.0 );
    float freq  = 10.0 / grid_size;
    float print1 = grid_thickness+(1.0+6.*grid_thickness)*abs(sin( freq * dot( d1, pos ) ) );
    float print2 = grid_thickness+(1.0+6.*grid_thickness)*abs(sin( freq * dot( d2, pos ) ) );
    float lc     = max( curv1,curv2);
    vec4  col_in = vec4( dir_color_on,  dir_on_opacity );
    vec4  col_out= vec4( dir_color_off, dir_off_opacity );
    float lambda = 3.0 * ( 0.5 + curv1 - print1 );
    if ( princ_curv == 1 )      lambda = 10.0 * ( 0.5 + curv2 - print2 );
    else if ( princ_curv == 2 ) lambda = max( lambda, 10.0 * ( 0.5 + curv2 - print2 ) );
    vec4 lambda_c = lambda * col_in + (1.-lambda) * col_out;
    // float lambda = 0.8;
    vec4 radiance = clamp( lc * lambda_c + (1.-lc)*lambda_c,
			   vec4( min( col_in.rgb, col_out.rgb ), 0.0),
			   vec4( max( col_in.rgb, col_out.rgb ), 1.0) );
    return radiance;
}

// main common to all color models

// Compute color at given ray position on the surface
vec4 compute_local_color( in Ray ray, out vec3 n ) {
    vec3  g       = gradf( ray.pos );
    n             = normalize( g );
    vec4  c       = display_color( ray, g, n );
    vec3 radiance = illumination( ray, c.rgb, n );
    vec3 sing_rad = apply_singularity( radiance, g );
    return vec4( apply_fog( ray, sing_rad ), c.a );
}

void main() {
    // we must compute the direction of the tracing ray within the cube
    outColor     = vec4( 0.0, 0.0, 0.0, 1.0 );
    Ray  ray     = ray_from_camera();
    float coef   = 1.0;
    bool found   = false;
    vec3 pos     = ray_trace( ray, found );
    vec3 normal;
    float m      = 0.0001 + transparency + reflection;
    float k_t    = transparency * transparency / m;
    float k_r    = reflection   * reflection   / m;
    float k_d    = 1.0 - k_t - k_r;
    if ( found ) {
	// First color
	ray.pos       = pos;
	vec4 lcolor   = compute_local_color( ray, normal );
	outColor.rgb += (k_d * lcolor.a) * lcolor.rgb;
	Ray first     = ray;
	// Reflections
	for ( int i = 0; i < nb_reflections; i++ ) {
	    ray           = ray_reflection( ray, normal );
	    coef         *= dot( ray.dir, ray.dir ) < 0.1 ? 0.0 : k_r;
	    if ( coef <= 0.02 ) break;
	    pos           = ray_trace( ray, found );
	    if ( ! found ) break;
	    ray.pos       = pos;
	    lcolor        = compute_local_color( ray, normal );
	    outColor.rgb += (coef * k_d * lcolor.a) * lcolor.rgb;
	}
	outColor.rgb += coef * background_color( ray.dir );
	// Refractions
	coef = 1.0;
	ray = first;
	for ( int i = 0; i < nb_refractions; i++ ) {
	    ray           = ray_refraction( ray, normal );
	    coef         *= dot( ray.dir, ray.dir ) < 0.1 ? 0.0 : k_t + (1.0-lcolor.a);
	    if ( coef <= 0.02 ) break;
	    pos           = ray_trace( ray, found );
	    if ( ! found ) break;
	    ray.pos       = pos;
	    lcolor        = compute_local_color( ray, normal );
	    outColor.rgb += (coef * k_d * lcolor.a ) * lcolor.rgb;
	    //k_t          *= lcolor.a;
	}
	outColor.rgb += coef * background_color( ray.dir );
    } else {
	outColor.rgb = background_color( ray.dir );
    }