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
TrigonometricMomentProblemShadow.fx
Go to the documentation of this file.
1 
6 
8 static const float PI=3.1415926535897932384626433832795f;
9 
18 float CircleToParameter(float2 UnitCirclePoint){
19  return (UnitCirclePoint.y>=0.0f)?(-UnitCirclePoint.x):(2.0f+UnitCirclePoint.x);
20 }
21 float2 CircleToParameter(float2x2 UnitCirclePoint){
22  return (UnitCirclePoint._m01_m11>=0.0f)?(-UnitCirclePoint._m00_m10):(2.0f+UnitCirclePoint._m00_m10);
23 }
24 float3 CircleToParameter(float3x2 UnitCirclePoint){
25  return (UnitCirclePoint._m01_m11_m21>=0.0f)?(-UnitCirclePoint._m00_m10_m20):(2.0f+UnitCirclePoint._m00_m10_m20);
26 }
27 float4 CircleToParameter(float4x2 UnitCirclePoint){
28  return (UnitCirclePoint._m01_m11_m21_m31>=0.0f)?(-UnitCirclePoint._m00_m10_m20_m30):(2.0f+UnitCirclePoint._m00_m10_m20_m30);
29 }
30 void CircleToParameter(out float pOutParameter[5],float2 pUnitCirclePoint[5]){
31  [unroll] for(int i=0;i!=5;++i){
32  pOutParameter[i]=((pUnitCirclePoint[i].y>=0.0f)?(-pUnitCirclePoint[i].x):(2.0f+pUnitCirclePoint[i].x));
33  }
34 }
36 
37 
42 float GetWeightExclusive(float3 pWeight,float3x2 pRoot){
43  float3 pParameter=CircleToParameter(pRoot);
44  return dot((pParameter.yz<pParameter[0])?pWeight.yz:float2(0.0f,0.0f),float2(1.0f,1.0f));
45 }
46 
47 
50 float3x2 MomentGeneratingFunction2(float Point){
51  float3x2 MomentVector;
52  MomentVector[0]=float2(1.0f,0.0f);
53  sincos(2.0f*PI*Point,MomentVector[1].y,MomentVector[1].x);
54  MomentVector[2]=Square(MomentVector[1]);
55  return MomentVector;
56 }
57 
58 
61 float3x2 MomentGeneratingFunction2(float2 CirclePoint){
62  float3x2 MomentVector;
63  MomentVector[0]=float2(1.0f,0.0f);
64  MomentVector[1]=CirclePoint;
65  MomentVector[2]=Square(MomentVector[1]);
66  return MomentVector;
67 }
68 
69 
75 void Get2TMSMCanonicalDistribution(out float3 pOutWeight,out float3x2 pOutRoot,float2 pTrigonometricMoment[2],float PrescribedRoot){
76  // Gather inputs in convenient ways
77  float2 b1=pTrigonometricMoment[0];
78  float2 b2=pTrigonometricMoment[1];
79  // Compute the right-hand side of the first linear system
80  float3x2 bz0=MomentGeneratingFunction2(PrescribedRoot);
81  // Construct the quadratic polynomial that has points of support as roots
82  float3x2 c=Conjugate(SolveToeplitzSystem(b1,b2,bz0[1],bz0[2]));
83  // Solve the polynomial equation to get all roots of the distribution
84  float2x2 pRoot=GetRoots(c);
85  pOutRoot[0]=bz0[1];
86  pOutRoot[1]=pRoot[0];
87  pOutRoot[2]=pRoot[1];
88  // Compute the weights
89  float InvDenominator=1.0f/dot(pOutRoot._m21_m01_m11-pOutRoot._m11_m21_m01,pOutRoot._m00_m10_m20);
90  pOutWeight[0]= dot(float3(pOutRoot[1].y-b1.y,b1.y-pOutRoot[2].y,pOutRoot[2].y-pOutRoot[1].y),float3(pOutRoot[2].x,pOutRoot[1].x,b1.x))*InvDenominator;
91  pOutWeight[1]=-dot(float3(pOutRoot[0].y-b1.y,b1.y-pOutRoot[2].y,pOutRoot[2].y-pOutRoot[0].y),float3(pOutRoot[2].x,pOutRoot[0].x,b1.x))*InvDenominator;
92  pOutWeight[2]=1.0f-pOutWeight[0]-pOutWeight[1];
93 }
94 
95 
99 float Get2TMSMCanonicalWeightExclusive(float2 pTrigonometricMoment[2],float PrescribedRoot){
100  // Compute roots of the canonical representation
101  float3x2 bz0=MomentGeneratingFunction2(PrescribedRoot);
102  float2x2 pAdditionalRoot=GetRoots(Conjugate(SolveToeplitzSystem(pTrigonometricMoment[0],pTrigonometricMoment[1],bz0[1],bz0[2])));
103  // Sort the roots
104  float2 pAdditionalParameter=CircleToParameter(pAdditionalRoot);
105  bool SwapRoots=(pAdditionalParameter.x>pAdditionalParameter.y);
106  float3x2 pRoot;
107  pRoot[0]=bz0[1];
108  pRoot[1]=(SwapRoots?pAdditionalRoot[1]:pAdditionalRoot[0]);
109  pRoot[2]=(SwapRoots?pAdditionalRoot[0]:pAdditionalRoot[1]);
110  pAdditionalParameter=SwapRoots?pAdditionalParameter.yx:pAdditionalParameter.xy;
111  // Figure out which combination of weights is needed
112  float PrescribedParameter=CircleToParameter(pRoot[0]);
113  float4 Switch=
114  (pAdditionalParameter.y<PrescribedParameter)?float4(pRoot[1],1.0f,1.0f):(
115  (pAdditionalParameter.x<PrescribedParameter)?float4(pRoot[0],0.0f,1.0f):
116  float4(0.0f,0.0f,0.0f,0.0f));
117  // Evaluate the combination of weights
118  float2 b1=pTrigonometricMoment[0];
119  float Quotient=dot(float3(Switch.y-b1.y,b1.y-pRoot[2].y,pRoot[2].y-Switch.y),float3(pRoot[2].x,Switch.x,b1.x))
120  /dot(pRoot._m21_m01_m11-pRoot._m11_m21_m01,pRoot._m00_m10_m20);
121  return Switch[2]-Switch[3]*Quotient;
122 }
123 
124 
127 void GetExtremalWeightHelpers2(out float3x2 OutMomentEnd,out float3x2 OutInvertedMomentEnd,out float3x2 OutInvertedMomentOne,out float2 OutEndDotOne,out float OutOneDotOne,float2 pTrigonometricMoment[2],float IntervalEnd){
128  // Evaluate the moment generating function at IntervalEnd and compute its
129  // inverse with respect to the Toeplitz matrix
130  OutMomentEnd=MomentGeneratingFunction2(IntervalEnd);
131  OutInvertedMomentEnd=SolveToeplitzSystem(pTrigonometricMoment[0],pTrigonometricMoment[1],OutMomentEnd[1],OutMomentEnd[2]);
132  // Do the same for the moment generating function at one
133  OutInvertedMomentOne=SolveToeplitzSystem(pTrigonometricMoment[0],pTrigonometricMoment[1],float2(1.0f,0.0f),float2(1.0f,0.0f));
134  // Compute various dot products, which are needed for both remaining cases
135  OutEndDotOne=Conjugate(OutInvertedMomentEnd[0]+OutInvertedMomentEnd[1]+OutInvertedMomentEnd[2]);
136  OutOneDotOne=dot(OutInvertedMomentOne._m00_m10_m20,float3(1.0f,1.0f,1.0f));
137 }
138 
139 
146 float4x2 GetWeightMapCriticalPoints(float3x2 MomentEnd,float3x2 InvertedMomentEnd,float2 EndDotOne,float OneDotOne,float2 pTrigonometricMoment[2],float IntervalEnd){
147  // Construct the polynomial
148  float2 pCoefficient[5]={float2(0.0f,0.0f),float2(0.0f,0.0f),float2(0.0f,0.0f),float2(0.0f,0.0f),float2(0.0f,0.0f)};
149  // Compute the inverse Toeplitz matrix
150  float3x3 ToeplitzInverseReal,ToeplitzInverseImaginary;
151  GetToeplitzInverse(ToeplitzInverseReal,ToeplitzInverseImaginary,pTrigonometricMoment);
152  // Start with the terms, which do not include the interval ends, but need
153  // to be multiplied by EndDotOne afterwards
154  [unroll] for(int i=0;i!=3;++i){
155  [unroll] for(int j=0;j!=3;++j){
156  float2 InverseToeplitzEntry=float2(ToeplitzInverseReal[i][j],ToeplitzInverseImaginary[i][j]);
157  pCoefficient[2+j-i]+=(j-i)*InverseToeplitzEntry;
158  }
159  }
160  // Multiply the first part by EndDotOne
161  for(int i=0;i!=5;++i){
162  pCoefficient[i]=Multiply(pCoefficient[i],EndDotOne);
163  }
164  // Invert the moment vector for the other interval end
165  float3x2 InvertedMomentOne;
166  InvertedMomentOne._m00_m10_m20=mul(ToeplitzInverseReal,float3(1.0f,1.0f,1.0f));
167  InvertedMomentOne._m01_m11_m21=mul(ToeplitzInverseImaginary,float3(1.0f,1.0f,1.0f));
168  // Construct the parts of the polynomial, which are formed by outer products
169  // of vectors depending on the variable, multiplied by inverted interval end
170  // vectors
171  [unroll] for(int i=0;i!=3;++i){
172  [unroll] for(int j=0;j!=3;++j){
173  float2 EndTimesOne=Multiply(InvertedMomentOne[i],Conjugate(InvertedMomentEnd[j]));
174  pCoefficient[2+j-i]-=(j-i)*EndTimesOne;
175  }
176  }
177  // Compute the roots of the polynomial to obtain critical points of the
178  // weight map
179  float4x2 pCriticalPoint=GetRoots(pCoefficient);
180  // The critical points are known to lie on the unit circle, but due to numerical
181  // inaccurracies this may not hold in practice. Renormalize.
182  [unroll] for(int i=0;i!=4;++i){
183  pCriticalPoint[i]=normalize(pCriticalPoint[i]);
184  }
185  return pCriticalPoint;
186 }
187 
188 
192 float GetMinimizingWeightAtOne(float2 pTrigonometricMoment[2],float IntervalEnd){
193  // Compute the canonical distribution at one
194  float3x2 pRootOne;
195  float3 pWeightOne;
196  Get2TMSMCanonicalDistribution(pWeightOne,pRootOne,pTrigonometricMoment,1.0f);
197  float3 pRootOneParameter=CircleToParameter(pRootOne);
198  // Sort the two free roots of this canonical distribution
199  bool SwapRootsOne=(pRootOneParameter[2]<pRootOneParameter[1]);
200  pRootOneParameter.yz=SwapRootsOne?pRootOneParameter.zy:pRootOneParameter.yz;
201  pWeightOne.yz=SwapRootsOne?pWeightOne.zy:pWeightOne.yz;
202  // Determine the present interval
203  float2 SegmentEnd;
204  sincos(2.0f*PI*IntervalEnd,SegmentEnd.y,SegmentEnd.x);
205  float IntervalEndParameter=CircleToParameter(SegmentEnd);
206  // Any canonical distribution will not have a root between 0.0 and its
207  // prescribed root IntervalEnd, so any weight at one is fine. We use 0.0 for
208  // reasons of numerical stability.
209  if(IntervalEndParameter<pRootOneParameter[1]){
210  return 0.0f;
211  }
212 
213  // Compute some shared, intermediate values
214  float3x2 MomentEnd,InvertedMomentEnd,InvertedMomentOne;
215  float2 EndDotOne;
216  float OneDotOne;
217  GetExtremalWeightHelpers2(MomentEnd,InvertedMomentEnd,InvertedMomentOne,EndDotOne,OneDotOne,pTrigonometricMoment,IntervalEnd);
218 
219  // We need to maximize the weight at the only root in the relevant interval.
220  // This can be done by solving a quartic equation for potentially optimal
221  // locations of this root.
222  if(IntervalEndParameter<pRootOneParameter[2]){
223  // Get critical points of the relevant weight map
224  float4x2 pCriticalPoint=GetWeightMapCriticalPoints(MomentEnd,InvertedMomentEnd,EndDotOne,OneDotOne,pTrigonometricMoment,IntervalEnd);
225  // Compute the canonical distribution at IntervalEnd to obtain the interval
226  // containing the relevant critical point
227  float3x2 pRootEnd;
228  float3 pWeightEnd;
229  Get2TMSMCanonicalDistribution(pWeightEnd,pRootEnd,pTrigonometricMoment,IntervalEnd);
230  float3 pRootEndParameter=CircleToParameter(pRootEnd);
231  bool SwapRootsEnd=(pRootEndParameter[2]<pRootEndParameter[1]);
232  float MinFirstRootParameter=SwapRootsEnd?pRootEndParameter[2]:pRootEndParameter[1];
233  float MinFirstWeight=SwapRootsEnd?pWeightEnd[2]:pWeightEnd[1];
234  // Find the critical point in this interval
235  float2 FirstRoot;
236  bool RootFound=false;
237  float4 pCriticalPointParameter=CircleToParameter(pCriticalPoint);
238  [unroll] for(int i=0;i!=4;++i){
239  if(MinFirstRootParameter<=pCriticalPointParameter[i] && pCriticalPointParameter[i]<=pRootOneParameter[1]){
240  FirstRoot=pCriticalPoint[i];
241  RootFound=true;
242  }
243  }
244  if(!RootFound){
245  // There is no extremum in the relevant interval, so the correct
246  // solution has to be to use one of the two already computed canonical
247  // distributions. Find out which one.
248  return (MinFirstWeight>pWeightOne[1])?pWeightOne[0]:0.0f;
249  }
250  // Compute the corresponding weight at one
251  float3x2 MomentFirst;
252  MomentFirst[0]=float2(1.0f,0.0f);
253  MomentFirst[1]=FirstRoot;
254  MomentFirst[2]=Square(FirstRoot);
255  float2 FirstDotEnd=Dot(MomentFirst,InvertedMomentEnd);
256  float2 FirstDotOne=Dot(MomentFirst,InvertedMomentOne);
257  float Weight=Divide(FirstDotEnd,OneDotOne*FirstDotEnd-Multiply(FirstDotOne,Conjugate(EndDotOne))).x;
258  // Clamp and return it
259  return max(0.0f,min(pWeightOne[0],Weight));
260  }
261  // We need to maximize the sum of the weight at one and at IntervalEnd. This can
262  // be done by means of a rather explicit closed-form.
263  else{
264  float EndDotEnd=Dot(MomentEnd,InvertedMomentEnd).x;
265  float Weight=(Magnitude(EndDotOne)-EndDotEnd)/(dot(EndDotOne,EndDotOne)-OneDotOne*EndDotEnd);
266  // Clamp and return it
267  return max(0.0f,min(pWeightOne[0],Weight));
268  }
269 }
270 
271 
275 float Get2TMSMMinimalWeight(float2 pTrigonometricMoment[2],float IntervalEnd){
276  // Compute the weight at one
277  float WeightAtOne=GetMinimizingWeightAtOne(pTrigonometricMoment,IntervalEnd);
278  // Update the moment vector and renormalize it
279  float2 pModifiedMoment[2]={
280  (pTrigonometricMoment[0]-float2(WeightAtOne,0.0f))/(1.0f-WeightAtOne),
281  (pTrigonometricMoment[1]-float2(WeightAtOne,0.0f))/(1.0f-WeightAtOne)};
282  // Compute the weight of the corresponding canonical distribution
283  float CanonicalWeight=Get2TMSMCanonicalWeightExclusive(pModifiedMoment,IntervalEnd);
284  // Scale it down to account for the weight at one
285  return saturate(CanonicalWeight*(1.0f-WeightAtOne));
286 }
float2 Conjugate(float2 Z)
void GetToeplitzInverse(out float3x3 OutToeplitzInverseReal, out float3x3 OutToeplitzInverseImaginary, float2 pTrigonometricMoment[2])
float3 GetRoots(float4 Coefficient)
float GetMinimizingWeightAtOne(float2 pTrigonometricMoment[2], float IntervalEnd)
float2 Divide(float2 Numerator, float2 Denominator)
void Get2TMSMCanonicalDistribution(out float3 pOutWeight, out float3x2 pOutRoot, float2 pTrigonometricMoment[2], float PrescribedRoot)
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)
float GetWeightExclusive(float3 pWeight, float3x2 pRoot)
float3x2 MomentGeneratingFunction2(float Point)
float CircleToParameter(float2 UnitCirclePoint)
float Magnitude(float2 Z)
float Get2TMSMMinimalWeight(float2 pTrigonometricMoment[2], float IntervalEnd)
float4x2 GetWeightMapCriticalPoints(float3x2 MomentEnd, float3x2 InvertedMomentEnd, float2 EndDotOne, float OneDotOne, float2 pTrigonometricMoment[2], float IntervalEnd)
void GetExtremalWeightHelpers2(out float3x2 OutMomentEnd, out float3x2 OutInvertedMomentEnd, out float3x2 OutInvertedMomentOne, out float2 OutEndDotOne, out float OutOneDotOne, float2 pTrigonometricMoment[2], float IntervalEnd)
float Get2TMSMCanonicalWeightExclusive(float2 pTrigonometricMoment[2], float PrescribedRoot)