// light shader that explores a strategy for soft shadows
// Bertrand Bry-Marfaing
// www.3db-site.com
// Mai 2008

light
softShadowSpot(	float intensity = 1;
				color lightColor = 1;
				float coneAngle = 35; // in degrees
				float penumbraAngle = 5; // outside the coneAngle
				string shadowMap = "";
				float samples = 64;
				float shadowBias = 0.01;
				float minRadius = 0.001;
				float maxRadius = 0.1;
				float biasConeSlope = 2;
				float shadowDecay = 0; // shadow decreases with occluder distance
				float shadDecayDistMin = 10; 
				float shadDecayDistMax = 25;)
{
	
	uniform float illumAngle = radians( coneAngle/2 + penumbraAngle );
	uniform float cosoutside = cos( illumAngle );
	uniform float cosinside = cos( radians( coneAngle/2 )  );
	uniform vector lightDir = vector "shader" ( 0,0,1 );
	
	
	/*
	_______ generate a series of random coordinates ______
			(implementation of Hammersley point set)
	*/
	uniform float sampleCoordX[ samples ];
	uniform float sampleCoordY[ samples ];
	// note: init of variable size array (although uniform variable in this case) does not compile with prman.
	uniform float iter = 0;
	
	// find the num of decimals
	uniform float numOfDigits = ceil( log( samples, 2 ) );
	uniform float revBinValues[ numOfDigits ];
	
	for( iter = 0; iter < samples; iter += 1 )
	{
		
		// generate binary fractions and their reversed
		uniform float digit = 0;
		uniform float remainder = iter;
		for( digit = (numOfDigits - 1); digit >= 0; digit -= 1 )
		{
			uniform float power = pow( 2, digit );
			if ( remainder >= power  )
			{
				revBinValues[ digit ] = 1;
				remainder -= power;
			}
			else
				revBinValues[ digit ] = 0;
			
		}
		
		// convert back to floats
		uniform float x = 0;
		uniform float y = 0;
		uniform float index = 0;
		for( index = 1; index <=  numOfDigits; index += 1 )
		{
			x += ( pow( 2, -index ) * revBinValues[ numOfDigits-index ] );
			y += ( pow( 2, -index ) * revBinValues[ index-1 ] );
		}
		
		// scale x to compensate for the non-power-of-two sample numbers
		x *= pow( 2, numOfDigits ) / samples ;
		
		// store coordinates
		sampleCoordX[ iter ] = x;
		sampleCoordY[ iter ] = y;
	}
	// ___________________________________
	
	
	
	
	/* 
		Function that interprets the previously generated x,y coordinates
		as angle and radius, and retruns the new coordinates.
		
		Output samples may be denser around 0,0 (for averaging occluder depth)
		or evenly spread (for shadowing sampling )
	*/
	float randomNumber = float random();
	void HammersleyToRadial( float shadowSampling; output float sampleCoordX; output float sampleCoordY )
	{
		
		;
		// if these samples are for light shadowing sampling
		if ( shadowSampling == 1 ){
			// samples should be evenly distributed
			sampleCoordY = sqrt( sampleCoordY );
			
			// vary the samples positions per shaded point
			extern float randomNumber;
			sampleCoordX = mod( sampleCoordX + randomNumber, 1 );
		} 
		
		float angleOfSample = sampleCoordX * 2 * PI;
		sampleCoordX = cos( angleOfSample ) * sampleCoordY; // cos = adj / hypot
		sampleCoordY = sqrt( pow( sampleCoordY, 2 ) - pow( sampleCoordX, 2 ) ); // a2 + b2 = c2 
		if( angleOfSample > PI )
			sampleCoordY *= -1;
	}
	// ___________________________________
	
	
	
	
	/*
		Function that samples the shadow map,
		returns sampled depth and normalised raster distance of sample
	*/
	float sampleShadowMap( float pShadS; float pShadT; float sampleX; float sampleY; float radius; float uniformSampling; output float distFromPs )
	{
		extern string shadowMap;
		
		// sample coord
		float ss = sampleX;
		float tt = sampleY;
		HammersleyToRadial( uniformSampling, ss, tt );
		ss *= radius;
		tt *= radius;
		ss += pShadS;
		tt += pShadT;
		
		// find distance from Ps coord
		distFromPs = distance( point( pShadS, pShadT, 0), point( ss, tt, 0 ) );
		
		// sample map without filtering
		return texture( shadowMap, ss, tt, "width", 0.00001 , "samples", 1 ); // this will generate warnings at compilation time with 3Delight
	}
	// ___________________________________
	
	
	
	
	
	illuminate( point "shader" ( 0,0,0 ), lightDir, illumAngle )
	{
		Cl = 0;
		float cosangle = normalize( L ).lightDir;
		float atten = smoothstep( cosoutside, cosinside, cosangle );
		float shadowValue = 0;
		
		
		if ( shadowMap != "" )
		{
			/*
			_____ Compute soft shadow _____
				
			*/
				
				
			// get the s and t coord of Ps on shadowmap
			float pShadS, pShadT;
			uniform matrix shadProjSpace;
			textureinfo( shadowMap, "projectionmatrix", shadProjSpace );
			point shadProjP = transform( shadProjSpace, Ps );
			pShadS = ( 1 + xcomp( shadProjP ) ) * 0.5;
			pShadT = ( 1 - ycomp( shadProjP ) ) * 0.5; 
			
			
			
			// get Ps distance from shadow plane
			uniform matrix shadCamSpace;
			textureinfo(  shadowMap, "viewingmatrix", shadCamSpace );
			point shadCamP = transform( shadCamSpace, Ps );
			float PsDepth = zcomp( shadCamP );
			
			
			
			// get the average occluder depth in shadowmap around Ps
			float avrgedDepth = 0;
			float totalWeight = 0;
			float iterator = 0;
			for ( iterator = 0; iterator < samples; iterator += 1 )
			{
				float sampleDepth, distFromPs;
				sampleDepth = sampleShadowMap( pShadS, pShadT, sampleCoordX[ iterator ], sampleCoordY[ iterator ], maxRadius, 0, distFromPs );
				 
				// weight according to raster distance from Ps coord
				if ( (sampleDepth+shadowBias) < PsDepth )
				{
					float sampleWeight = 1 - 0.9 * ( distFromPs / maxRadius );
					totalWeight += sampleWeight;
					avrgedDepth += ( sampleDepth * sampleWeight ); 
				}
			}
			if ( totalWeight > 0 ) avrgedDepth /= totalWeight;
			// _________________________________________
			
			
			
			
			if( avrgedDepth > 0 ) // there is an occluder
			{
					
				// compute the radius of the filter
				float filterRadius = minRadius + (( maxRadius-minRadius ) * ( 1 - ( avrgedDepth / PsDepth )));
						
					
				// sample the shadow map with varying filter width
				avrgedDepth = 0;
				totalWeight = 0;
				
				// get the shadow occlusion of Ps
				for ( iterator = 0; iterator < samples; iterator += 1 )
				{
					
					float sampleDepth, distFromPs;
					sampleDepth = sampleShadowMap( pShadS, pShadT, sampleCoordX[ iterator ], sampleCoordY[ iterator ], filterRadius, 1, distFromPs );
							
					// add bias cone
					sampleDepth += shadowBias + ( distFromPs * biasConeSlope ); 
							
					if ( sampleDepth < PsDepth )
					{ 
						shadowValue += 1;
						
						// if shadow decay
						if( shadowDecay > 0 )
						{
							// store the occluder depth of the new sampling disc
							float sampleWeight = 1 - 0.9 * ( distFromPs / filterRadius );
							totalWeight += sampleWeight;
							avrgedDepth += ( sampleDepth * sampleWeight ); 
						}
					}
				}
				shadowValue /= samples;
				// __________________________________________
				
				
				
				
				
				// decrease light attenuation as Ps moves further from occluder
				if ( shadowDecay > 0 )
				{
					if ( totalWeight > 0 ) avrgedDepth /= totalWeight;
					float distFromOccluder = PsDepth - avrgedDepth;
					shadowValue *= 1 - ( smoothstep( shadDecayDistMin, shadDecayDistMax, distFromOccluder ) * shadowDecay );
				}
				// __________________________________________
				
			
			}	
			
			atten *= ( 1 - shadowValue );
		}
		
		Cl = intensity * lightColor * atten;
	}
}