1 files changed, 118 insertions, 0 deletions
diff --git a/demo/d3d/shaders/ellipsoidDepthVS.hlsl b/demo/d3d/shaders/ellipsoidDepthVS.hlsl
new file mode 100644
index 0000000..f26317e
--- /dev/null
+++ b/demo/d3d/shaders/ellipsoidDepthVS.hlsl
@@ -0,0 +1,118 @@
+#include "shaderCommon.h"
+
+cbuffer constBuf : register(b0)
+{
+	FluidShaderConst gParams;
+};
+
+// returns 1.0 for x==0.0 (unlike glsl)
+float Sign(float x) { return x < 0.0 ? -1.0 : 1.0; }
+
+bool solveQuadratic(float a, float b, float c, out float minT, out float maxT)
+{
+#if 0
+	// for debugging
+	minT = -0.5;
+	maxT = 0.5;
+	return true;
+#else
+	//minT = 0.0f;
+	//maxT = 0.0f;
+#endif
+
+	if (a == 0.0 && b == 0.0)
+	{
+		minT = maxT = 0.0;
+		return false;
+	}
+
+	float discriminant = b*b - 4.0*a*c;
+
+	if (discriminant < 0.0)
+	{
+		return false;
+	}
+
+	float t = -0.5*(b + Sign(b)*sqrt(discriminant));
+	minT = t / a;
+	maxT = c / t;
+
+	if (minT > maxT)
+	{
+		float tmp = minT;
+		minT = maxT;
+		maxT = tmp;
+	}
+
+	return true;
+}
+
+float DotInvW(float4 a, float4 b) { return a.x*b.x + a.y*b.y + a.z*b.z - a.w*b.w; }
+
+FluidVertexOut ellipsoidDepthVS(FluidVertexIn input, uint instance : SV_VertexID)
+{
+	const float4 q1 = input.q1;
+	const float4 q2 = input.q2;
+	const float4 q3 = input.q3;
+
+	const float4x4 modelViewProjectionMatrix = gParams.modelViewProjection;
+	const float4x4 modelViewMatrixInverse = gParams.inverseModelView;
+
+	float3 worldPos = input.position.xyz;
+	
+	// construct quadric matrix
+	float4x4 q;
+	q._m00_m10_m20_m30 = float4(q1.xyz*q1.w, 0.0);
+	q._m01_m11_m21_m31 = float4(q2.xyz*q2.w, 0.0);
+	q._m02_m12_m22_m32 = float4(q3.xyz*q3.w, 0.0);
+	q._m03_m13_m23_m33 = float4(worldPos, 1.0);
+
+	// transforms a normal to parameter space (inverse transpose of (q*modelview)^-T)
+	float4x4 invClip = mul(modelViewProjectionMatrix, q);
+
+	// solve for the right hand bounds in homogenous clip space
+	float a1 = DotInvW(invClip[3], invClip[3]);
+	float b1 = -2.0f*DotInvW(invClip[0], invClip[3]);
+	float c1 = DotInvW(invClip[0], invClip[0]);
+
+	float xmin;
+	float xmax;
+	solveQuadratic(a1, b1, c1, xmin, xmax);
+
+	// solve for the right hand bounds in homogenous clip space
+	float a2 = DotInvW(invClip[3], invClip[3]);
+	float b2 = -2.0f*DotInvW(invClip[1], invClip[3]);
+	float c2 = DotInvW(invClip[1], invClip[1]);
+
+	float ymin;
+	float ymax;
+	solveQuadratic(a2, b2, c2, ymin, ymax);
+
+	FluidVertexOut output;
+	output.position = float4(worldPos.xyz, 1.0);
+	output.bounds = float4(xmin, xmax, ymin, ymax);
+
+	// construct inverse quadric matrix (used for ray-casting in parameter space)
+	float4x4 invq;
+	invq._m00_m10_m20_m30 = float4(q1.xyz / q1.w, 0.0);
+	invq._m01_m11_m21_m31 = float4(q2.xyz / q2.w, 0.0);
+	invq._m02_m12_m22_m32 = float4(q3.xyz / q3.w, 0.0);
+	invq._m03_m13_m23_m33 = float4(0.0, 0.0, 0.0, 1.0);
+
+	invq = transpose(invq);
+	invq._m03_m13_m23_m33 = -(mul(invq, output.position));
+
+	// transform a point from view space to parameter space
+	invq = mul(invq, modelViewMatrixInverse);
+
+	// pass down
+	output.invQ0 = invq._m00_m10_m20_m30;
+	output.invQ1 = invq._m01_m11_m21_m31;
+	output.invQ2 = invq._m02_m12_m22_m32;
+	output.invQ3 = invq._m03_m13_m23_m33;
+	
+	// compute ndc pos for frustrum culling in GS
+	float4 projPos = mul(modelViewProjectionMatrix, float4(worldPos.xyz, 1.0));
+	output.ndcPos = projPos / projPos.w;
+	return output;
+}