demotool/test/lib/raymarching.glsl

//! type library
//! include utils.glsl

vec2 RM_Map( in vec3 pos );

// -----------------------------------------------------------------------------

vec3 RM_GetNormal(
	in vec3 pos )
{
	const vec2 e = vec2( 1 , -1 ) * .0001;
	const vec4 f = vec4(
		RM_Map( pos + e.xyy ).x ,
		RM_Map( pos + e.yyx ).x ,
		RM_Map( pos + e.yxy ).x ,
		RM_Map( pos + e.xxx ).x );
	return normalize( f.x * e.xyy
			+ f.y * e.yyx
			+ f.z * e.yxy
			+ f.w * e.xxx );
}


// -----------------------------------------------------------------------------

vec3 RM_Basic(
	in vec3 origin ,
	in vec3 direction ,
	in int steps ,
	in float factor ,
	in float epsilon ,
	in float dMin ,
	in float dMax )
{
	int i = 0;
	float dist = dMin , mat = -1;

	for ( ; i < steps ; ++ i ) {
		vec2 res = RM_Map( origin + direction * dist );
		if ( abs( res.x ) < epsilon || dist > dMax ) {
			break;
		}
		dist += res.x * factor;
		mat = res.y;
	}

	return vec3( dist , dist >= dMax ? -1 : mat , i );
}

vec3 RM_Advanced(
	in vec3 origin ,
	in vec3 direction ,
	in int steps ,
	in float factor ,
	in float epsilon ,
	in float dMin ,
	in float dMax )
{
	int i = 0;
	float dist = dMin ,
	      omega = factor ,
	      cError = 1 / 0. ,
	      cDist = dMin ,
	      pRad = 0 ,
	      sLen = 0;

	for ( ; i < steps ; ++ i ) {
		vec2 res = RM_Map( origin + direction * dist );
		float rad = abs( res.x );

		bool sorFail = omega > 1 && ( rad + pRad ) < sLen;
		if ( sorFail ) {
			sLen -= omega * sLen;
			omega = 1;
		} else {
			sLen = res.x * omega;
		}
		pRad = rad;

		float error = rad / dist;
		if ( !sorFail && error < cError ) {
			cError = error;
			cDist = dist;
		}

		if ( !sorFail && error < epsilon || dist > dMax ) {
			break;
		}

		dist += sLen;
	}
	if ( dist <= dMax && cError <= epsilon ) {
		return vec3( cDist , RM_Map( origin + direction * cDist ).y , i );
	}
	return vec3( cDist , -1 , steps );
}

float RM_TanPhi(
	in vec2 uv ,
	in float pixelRadius )
{
	// 1 / tan( FOV * PI / 360 )
	// u_NearPlane = 2 * tan( PI * ( 180 - FOV ) / 360 );
	// atan( u_NearPlane / 2 ) = PI * ( 180 - FOV ) / 360
	//			= PI / 2 - PI * FOV / 360
	// atan( u_NearPlane / 2 ) - PI / 2 = - PI * FOV / 360
	// PI / 2 - atan( u_NearPlane / 2 ) = PI * FOV / 360
	// ...
	// ...
	// or we could just pass it as a fucking uniform
	const float hfov = ( PI * .5 - atan( u_NearPlane * .5 ) ) * .5;

	vec2 e = ( abs( uv.x ) > abs( uv.y ) )
		? vec2( uv.x < 0 ? -1 : 1 , 0 )
		: vec2( 0 , uv.y < 0 ? -1 : 1  );

	vec3	d = normalize( vec3( uv , -1 / tan( hfov ) ) ) ,
		a = vec3( ( e + e.yx ) * pixelRadius + uv ,
			-1 / tan( hfov ) ) ,
		b = vec3( ( e - e.yx ) * pixelRadius + uv ,
			-1 / tan( hfov ) ) ,
		p2ca = a / dot( a , d ) ,
		p2cb = a / dot( b , d ) ,
		diff = p2cb - p2ca ,
		de = mix( p2ca , p2cb , M_Saturate(
			dot( d - p2ca , diff ) / M_LengthSqr( diff ) ) );

	float sinPhiSquared = M_LengthSqr( cross( d , normalize( de ) ) );
	return sqrt( sinPhiSquared / ( 1 - sinPhiSquared ) );
}

vec3 RM_ReduceDiscontinuity(
	in vec3 origin ,
	in vec3 rayDir ,
	in float hitDist ,
	in float tanPhi ,
	in int iterations )
{
	vec3 hitPos = origin + hitDist * rayDir;
	for ( int i = 0 ; i < iterations ; i ++ ) {
		hitPos = hitPos + rayDir * (
			RM_Map( hitPos ).x - tanPhi * length( origin - hitPos )
		);
	}
	return hitPos;
}