From 6bb35ae40e5198f6bbd7fd0a193c2fb7b226dd9c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Emmanuel=20Beno=C3=AEt?= <tseeker@nocternity.net>
Date: Thu, 5 Oct 2017 17:03:31 +0200
Subject: [PATCH] Discontinuity reduction by Mercury et al.

---
 shaders/lib/raymarching.glsl | 58 ++++++++++++++++++++++++++++++++++++
 shaders/lib/utils.glsl       | 18 +++++++++--
 shaders/scene.f.glsl         | 12 ++++----
 3 files changed, 81 insertions(+), 7 deletions(-)

diff --git a/shaders/lib/raymarching.glsl b/shaders/lib/raymarching.glsl
index 645b294..21a574d 100644
--- a/shaders/lib/raymarching.glsl
+++ b/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;
+}
diff --git a/shaders/lib/utils.glsl b/shaders/lib/utils.glsl
index 3ef367d..9d64c09 100644
--- a/shaders/lib/utils.glsl
+++ b/shaders/lib/utils.glsl
@@ -17,14 +17,28 @@ float M_Hash(
 // -----------------------------------------------------------------------------
 
 float M_Luminosity(
-		in vec3 color )
+	in vec3 color )
 {
 	return dot( color , vec3( .299 , .587 , .114 ) );
 }
 
 vec3 M_NormalizeColor(
-		in vec3 color )
+	in vec3 color )
 {
 	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 );
+}
diff --git a/shaders/scene.f.glsl b/shaders/scene.f.glsl
index fd92d75..1609ed8 100644
--- a/shaders/scene.f.glsl
+++ b/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 );
-	vec3 hitPos = u_CamPos + rayDir * r.x;
+			/*u_Render.z */ tanPhi , .001 , u_Render.w );
 	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;