Discontinuity reduction by Mercury et al.

shaders/lib/raymarching.glsl

@@ -1,4 +1,5 @@
 //! type library
+//! include lib/utils.glsl
 
 vec2 RM_Map( in vec3 pos );
 
@@ -88,3 +89,60 @@ vec3 RM_Advanced(
 	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 fov = PI * .5 - atan( u_NearPlane * .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( fov ) ) ) ,
+		a = vec3( ( e + e.yx ) * pixelRadius + uv ,
+			-1 / tan( fov ) ) ,
+		b = vec3( ( e - e.yx ) * pixelRadius + uv ,
+			-1 / tan( fov ) ) ,
+		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 ) );
+}
+
+float RM_ErrFunc(
+	in vec3 origin ,
+	in vec3 pos )
+{
+	return 0; // FIXME
+}
+
+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 ++ ) {
+		const float err = length( origin - hitPos ) * tanPhi;
+		hitPos = hitPos + rayDir * (
+			RM_Map( hitPos ).x - tanPhi
+		);
+	}
+	return hitPos;
+}

shaders/lib/utils.glsl

@@ -28,3 +28,17 @@ vec3 M_NormalizeColor(
 	const float l = M_Luminosity( color );
 	return l > 0 ? ( color / l ) : vec3( 1 );
 }
+
+// -----------------------------------------------------------------------------
+
+float M_Saturate(
+	in float x )
+{
+	return clamp( x , 0 , 1 );
+}
+
+float M_LengthSqr(
+	in vec3 x )
+{
+	return dot( x , x );
+}

shaders/scene.f.glsl

@@ -94,24 +94,26 @@ void main( )
 			+ uv.x * side * u_Resolution.x / u_Resolution.y
 			+ uv.y * up );
 
+	float tanPhi = RM_TanPhi( uv , u_Render.z );
+
 	vec3 r = RM_Advanced( u_CamPos , rayDir ,
 			int( u_Render.x ) , u_Render.y ,
-			u_Render.z , .001 , u_Render.w );
+			/*u_Render.z */ tanPhi , .001 , u_Render.w );
-	vec3 hitPos = u_CamPos + rayDir * r.x;
 	const vec3 background = vec3( .005 );
 
 	vec3 bc = background;
 	if ( r.y >= 0. ) {
-		const int midx = int( r.y );
+		const vec3 hitPos = RM_ReduceDiscontinuity(
+				u_CamPos , rayDir , r.x ,
+				tanPhi , 10 );
 		const vec3 normal = RM_GetNormal( hitPos );
 
+		const int midx = int( r.y );
 		const vec3 sunlightDir = normalize( vec3( 0 , 1 , 1 ) );
 		const vec3 sunlightColor = vec3( 1 , 1 , .99 ) * .1;
-
 		const vec3 pl1Pos = vec3( 3 , 3 , 5 );
 		const vec3 pl1Color = vec3( 1 , .7 , .8 );
 		const float pl1Power = 20;
-
 		const vec3 pl2Pos = vec3( -3 , -3 , 5 );
 		const vec3 pl2Color = vec3( .6 , 1 , .8 );
 		const float pl2Power = 20;