Shaders for Shadows
Shader code for rendering shadows of translucent occluders, soft shadows and single scattering with moment shadow maps.
 All Classes Files Functions Variables Pages
TrigonometricMomentProblemShadowAlgebra.fx
Go to the documentation of this file.
1 
6 #include "ComplexAlgebra.fx"
7 
10 float2 Dot(float3x2 LHS,float3x2 RHS){
11  return float2(
12  dot(transpose(LHS)[0],transpose(RHS)[0])+dot(transpose(LHS)[1],transpose(RHS)[1]),
13  dot(transpose(LHS)[0],transpose(RHS)[1])-dot(transpose(LHS)[1],transpose(RHS)[0]));
14 }
15 float2 Dot(float4x2 LHS,float4x2 RHS){
16  return float2(
17  dot(transpose(LHS)[0],transpose(RHS)[0])+dot(transpose(LHS)[1],transpose(RHS)[1]),
18  dot(transpose(LHS)[0],transpose(RHS)[1])-dot(transpose(LHS)[1],transpose(RHS)[0]));
19 }
20 float2 Dot(float2 LHS[5],float2 RHS[5]){
21  float2 Result=float2(0.0f,0.0f);
22  [unroll] for(int i=0;i!=5;++i){
23  Result+=Multiply(LHS[i],RHS[i]);
24  }
25  return Result;
26 }
28 
29 
36 float3x2 Multiply(float3x3 MatrixRealPart,float3x3 MatrixImaginaryPart,float3x2 RHSVector){
37  float3x2 Result;
38  Result._m00_m10_m20=mul(MatrixRealPart,RHSVector._m00_m10_m20)-mul(MatrixImaginaryPart,RHSVector._m01_m11_m21);
39  Result._m01_m11_m21=mul(MatrixRealPart,RHSVector._m01_m11_m21)+mul(MatrixImaginaryPart,RHSVector._m00_m10_m20);
40  return Result;
41 }
42 float4x2 Multiply(float4x4 MatrixRealPart,float4x4 MatrixImaginaryPart,float4x2 RHSVector){
43  float4x2 Result;
44  Result._m00_m10_m20_m30=mul(MatrixRealPart,RHSVector._m00_m10_m20_m30)-mul(MatrixImaginaryPart,RHSVector._m01_m11_m21_m31);
45  Result._m01_m11_m21_m31=mul(MatrixRealPart,RHSVector._m01_m11_m21_m31)+mul(MatrixImaginaryPart,RHSVector._m00_m10_m20_m30);
46  return Result;
47 }
48 void Multiply(out float2 OutProduct[5],float MatrixRealPart[5][5],float MatrixImaginaryPart[5][5],float2 RHSVector[5]){
49  [unroll] for(int i=0;i!=5;++i){
50  OutProduct[i]=float2(0.0f,0.0f);
51  [unroll] for(int j=0;j!=5;++j){
52  OutProduct[i]+=Multiply(float2(MatrixRealPart[i][j],MatrixImaginaryPart[i][j]),RHSVector[j]);
53  }
54  }
55 }
57 
58 
59 
64 float2x2 GetRoots(float3x2 pCoefficient){
65  float2 aInv=Reciprocal(pCoefficient[2]);
66  float2 p=Multiply(pCoefficient[1],aInv);
67  float2 q=Multiply(pCoefficient[0],aInv);
68  float2 Discriminant=0.25f*Square(p)-q;
69  float2 Root=SquareRoot(Discriminant);
70  float2 pHalf=-0.5f*p;
71  float2x2 pRoot;
72  pRoot[0]=pHalf-Root;
73  pRoot[1]=pHalf+Root;
74  return pRoot;
75 }
76 
78 float3x2 GetRoots(float4x2 pCoefficient){
79  // This implementation is based upon the algorithm and notions on Wikipedia:
80  // http://en.wikipedia.org/wiki/Cubic_polynomial#General_formula_for_roots
81  float2 bb=Square(pCoefficient[2]);
82  float2 ac=Multiply(pCoefficient[1],pCoefficient[3]);
83  float2 Delta0=bb-3.0f*ac;
84  float2 bbb=Multiply(bb,pCoefficient[2]);
85  float2 abc=Multiply(ac,pCoefficient[2]);
86  float2 aad=Multiply(Multiply(pCoefficient[3],pCoefficient[3]),pCoefficient[0]);
87  float2 Delta1=2.0f*bbb-9.0f*abc+27.0f*aad;
88  float2 C=CubicRoot((Delta1+SquareRoot(Multiply(Delta1,Delta1)-4.0f*Cube(Delta0)))*0.5f);
89  float2 aInv=Reciprocal(pCoefficient[3]);
90  float2 Delta0OverC=Divide(Delta0,C);
91  float3x2 pRoot;
92  pRoot[0]=-(1.0f/3.0f)*Multiply(aInv,pCoefficient[2]+C+Delta0OverC);
93  pRoot[1]=-(1.0f/3.0f)*Multiply(aInv,pCoefficient[2]+Multiply(float2(-1.0f, sqrt(3.0f))/2.0f,C)+Multiply(float2(-1.0f,-sqrt(3.0f))/2.0f,Delta0OverC));
94  pRoot[2]=-(1.0f/3.0f)*Multiply(aInv,pCoefficient[2]+Multiply(float2(-1.0f,-sqrt(3.0f))/2.0f,C)+Multiply(float2(-1.0f, sqrt(3.0f))/2.0f,Delta0OverC));
95  return pRoot;
96 }
97 
99 float4x2 GetRoots(float2 pCoefficient[5]){
100  // This implementation is based upon the algorithm and notions on Wikipedia:
101  // http://en.wikipedia.org/wiki/Quartic_equation#General_formula_for_roots
102  float2 cc=Square(pCoefficient[2]);
103  float2 bd=Multiply(pCoefficient[1],pCoefficient[3]);
104  float2 ae=Multiply(pCoefficient[0],pCoefficient[4]);
105  float2 Delta0=cc-3.0f*bd+12.0f*ae;
106 
107  float2 ccc=Multiply(cc,pCoefficient[2]);
108  float2 bcd=Multiply(bd,pCoefficient[2]);
109  float2 bb=Square(pCoefficient[3]);
110  float2 bbe=Multiply(bb,pCoefficient[0]);
111  float2 add=Multiply(Square(pCoefficient[1]),pCoefficient[4]);
112  float2 ace=Multiply(ae,pCoefficient[2]);
113  float2 Delta1=2.0f*ccc-9.0f*bcd+27.0f*bbe+27.0f*add-72.0f*ace;
114 
115  float2 Q=CubicRoot((Delta1+SquareRoot(Square(Delta1)-4.0f*Cube(Delta0)))*0.5f);
116 
117  float2 ac=Multiply(pCoefficient[2],pCoefficient[4]);
118  float2 aa=Square(pCoefficient[4]);
119  float2 p=Divide(8.0f*ac-3.0f*bb,8.0f*aa);
120  float2 aaa=Multiply(aa,pCoefficient[4]);
121  float2 aad=Multiply(aa,pCoefficient[1]);
122  float2 abc=Multiply(pCoefficient[2],Multiply(pCoefficient[3],pCoefficient[4]));
123  float2 bbb=Multiply(bb,pCoefficient[3]);
124  float2 q=Divide(bbb-4.0f*abc+8.0f*aad,8.0f*aaa);
125 
126  float2 aInv=Reciprocal(pCoefficient[4]);
127  float2 S=0.5f*SquareRoot((-2.0f/3.0f)*p+(1.0f/3.0f)*Multiply(aInv,Q+Divide(Delta0,Q)));
128 
129  float2 bOver4a=0.25f*Multiply(pCoefficient[3],aInv);
130  float2 SS4Plus2p=4.0f*Square(S)+2.0f*p;
131  float2 qOverS=Divide(q,S);
132  float2 FirstSquareRoot=SquareRoot(-SS4Plus2p+qOverS);
133  float2 SecondSquareRoot=SquareRoot(-SS4Plus2p-qOverS);
134 
135  float4x2 pRoot;
136  pRoot[0]=-bOver4a-S+0.5f*FirstSquareRoot;
137  pRoot[1]=-bOver4a-S-0.5f*FirstSquareRoot;
138  pRoot[2]=-bOver4a+S+0.5f*SecondSquareRoot;
139  pRoot[3]=-bOver4a+S-0.5f*SecondSquareRoot;
140  return pRoot;
141 }
142 
143 
146 float3 GetRoots(float p,float q){
147  // Use Vieta's substitution
148  float2 w=float2(-0.5f*q,sqrt(-(0.25f*q*q+p*p*p/27.0f)));
149  // Compute the third roots
150  float2 pThirdRoot[3];
151  pThirdRoot[0]=CubicRoot(w);
152  pThirdRoot[1]=Multiply(pThirdRoot[0],float2(-0.5f,0.5f*sqrt(3.0f)));
153  pThirdRoot[2]=Multiply(pThirdRoot[0],float2(-0.5f,-0.5f*sqrt(3.0f)));
154  // Revert the substitution and extract the real part
155  float3 pRoot;
156  [unroll] for(int i=0;i!=3;++i){
157  pRoot[i]=pThirdRoot[i].x-p*(1.0f/3.0f)*Reciprocal(pThirdRoot[i]).x;
158  }
159  return pRoot;
160 }
161 
162 
169 void GetToeplitzInverse(out float3x3 OutToeplitzInverseReal,out float3x3 OutToeplitzInverseImaginary,float2 pTrigonometricMoment[2]){
170  float2 b1=pTrigonometricMoment[0];
171  float2 b2=pTrigonometricMoment[1];
172  float m00=1.0f-dot(b1,b1);
173  float2 m10=-b1+Multiply(Conjugate(b1),b2);
174  float m11=1.0f-dot(b2,b2);
175  float2 m20=Multiply(b1,b1)-b2;
176  OutToeplitzInverseReal=float3x3(
177  m00.x,m10.x,m20.x,
178  m10.x,m11.x,m10.x,
179  m20.x,m10.x,m00.x
180  );
181  OutToeplitzInverseImaginary=float3x3(
182  0.0f,-m10.y,-m20.y,
183  m10.y, 0.0f,-m10.y,
184  m20.y, m10.y, 0.0f
185  );
186  // Compute the real part of the (non-Hermitian) dot product of the first row of
187  // the scaled inverse with the first column of the Toeplitz matrix
188  float Scaling=dot(OutToeplitzInverseReal[0],float3(1.0f,pTrigonometricMoment[0].x,pTrigonometricMoment[1].x))-dot(OutToeplitzInverseImaginary[0],float3(0.0f,pTrigonometricMoment[0].y,pTrigonometricMoment[1].y));
189  // It should be one, so we need to divide by that
190  OutToeplitzInverseReal/=Scaling;
191  OutToeplitzInverseImaginary/=Scaling;
192 }
193 
194 
200 float3x2 SolveToeplitzSystem(float2 TrigonometricMoment1,float2 TrigonometricMoment2,float2 RHS1,float2 RHS2){
201  float2 b1=TrigonometricMoment1;
202  float2 b2=TrigonometricMoment2;
203  // Construct a Cholesky decomposition of the Toeplitz matrix
204  float D11=1.0f-dot(b1,b1);
205  float InvD11=1.0f/D11;
206  float2 L21=(b1-Multiply(Conjugate(b1),b2))*InvD11;
207  float D22=1.0f-dot(b2,b2)-dot(L21,L21)*D11;
208  // Perform forward and backward substitution to solve the linear system
209  float3x2 c;
210  c[0]=float2(1.0f,0.0f);
211  c[1]=RHS1;
212  c[2]=RHS2;
213  // Forward substitution to solve L*c1=bz
214  c[1]-=b1;
215  c[2]-=b2+Multiply(L21,c[1]);
216  // Scaling to solve D*c2=c1
217  c[1]*=InvD11;
218  c[2]*=1.0f/D22;
219  // Backward substitution to solve L^H*c3=c2
220  c[1]-=Multiply(Conjugate(L21),c[2]);
221  c[0]-=Multiply(Conjugate(b1),c[1])+Multiply(Conjugate(b2),c[2]);
222  return c;
223 }
224 
225 
228 float GetDeterminant(float2 pTrigonometricMoment[2]){
229  return 1.0f-(2.0f*dot(pTrigonometricMoment[0],pTrigonometricMoment[0])+dot(pTrigonometricMoment[1],pTrigonometricMoment[1]))+2.0f*Multiply(Square(pTrigonometricMoment[0]),Conjugate(pTrigonometricMoment[1])).x;
230 }
float2 Cube(float2 Z)
float2 Conjugate(float2 Z)
void GetToeplitzInverse(out float3x3 OutToeplitzInverseReal, out float3x3 OutToeplitzInverseImaginary, float2 pTrigonometricMoment[2])
float3 GetRoots(float4 Coefficient)
float2 Divide(float2 Numerator, float2 Denominator)
float2 Square(float2 Z)
float2 Dot(float3x2 LHS, float3x2 RHS)
float2 Multiply(float2 LHS, float2 RHS)
float3x2 SolveToeplitzSystem(float2 TrigonometricMoment1, float2 TrigonometricMoment2, float2 RHS1, float2 RHS2)
float2 Reciprocal(float2 Z)
float2 CubicRoot(float2 Z)
float2 SquareRoot(float2 Z)
float GetDeterminant(float2 pTrigonometricMoment[2])