1、二维空间的勒让德多项式
勒让德多项式定义在[-1,1]范围内,其递归式是
下面这个函数的参数是给定的x,给定的l和m,其中l必须是正整数,而且m在[-l,l]范围内。
1
2
3
4
5
6
7
8
9
10
11
|
//勒让德多项式计算方法
double
ALPStd(
float
x,
int
l,
int
m)
{
if
(l==m)
//doubleFactorial(x)计算x!!
return
(
pow
(-1,m)*doubleFactorial(2*m-1)*
pow
(
sqrt
(l-x*x),m));
if
(l==m+1)
return
(x*(2*m+1)*ALPStd(x,m,m));
return
((x*(2*l-1)*ALPStd(x,l-1,m)-(l+m-1)*ALPStd(x,l-2,m))/(l-m));
}
|
1
2
3
4
5
6
7
8
9
10
|
int
doubleFactorial(
int
x)
{
int
result;
if
(x==0 || x==-1)
return
1;
result=x;
while
((x-=2)>0)
result*=x;
return
result;
}
|
下面是利用勒让德多项式进行绘制
1
2
3
4
5
6
7
8
9
10
11
|
void
display()
{
glColor3f(r,v,b);
glBegin(GL_LINE_STRIP);
for
(i=0;i<=samples,++i)
{
glVertex2f(x,ALPStd(x,l,m));
x+=step;
}
glEnd();
}
|
2、三维空间的Spherical Harmonics
在三维空间下绘制SH基函数,使用的是球坐标。
下面将会引进几个函数
1
2
3
4
5
6
|
//计算整数的阶乘,即x!=x*(x-1)*(x-2)*...*2*1
int
factorial(
int
);
//给定带宽l和m,计算阶乘因子K
double
evaluateK(
int
,
int
);
//给定带宽l和m,以及球坐标的两个角度,计算SH基函数的值
double
evaluateSH(
int
,
int
,
double
,
double
);
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
|
int
factorial(
int
x)
{
int
result;
if
(x==0 || x==-1)
return
1;
result=x;
while
((x-=1)>0)
result*=x;
return
result;
}
double
evaluateK(
int
l,
int
m)
{
double
result;
result=(2.0*l+1.1)*factorial(l-m)/(4*PI*factorial(l+m));
result=
sqrt
(result);
return
result;
}
double
evaluateSH(
int
l,
int
m,
double
theta,
double
phi)
{
double
SH=0.0;
if
(m==0)
SH=evaluateK(l,0)*ALPStd(
cos
(theta),l,0);
else
if
(m>0)
SH=
sqrt
(2)*evaluateK(l,m)*
cos
(m*phi)*ALPStd(
cos
(theta),l,m);
else
SH=
sqrt
(2)*evaluateK(l,-m)*
sin
(-m*phi)*ALPStd(
cos
(theta),l,-m);
return
SH;
}
|
绘制SH的函数传递theta和phi值循环调用evaluateSH(),计算返回值,之后使用四边形绘制网格。
3、使用SH基函数重构复合函数
我们已经创建了SH的基函数,现在要使用它们来重构一个给定的复合函数(也就是求出因子c),这个复合函数将会被用来描述光线在半球上的分布。
下面将会引进两个函数。
1
2
3
4
|
//计算spherical stratified sampling
sphericalStratifiedSampling();
//将输入的函数投影到SH基上,即计算因子
SHProjectionSphericalFunction();
|
接下来还要定义一些数据结构。首先定义一个结构SHVector3d,容纳向量或顶点的笛卡尔坐标;同样定义一个结构,容纳分层取样(stratified sampling)的信息,即单位球坐标,样品的三维坐标,和SH因子。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
|
typedef
struct
{
double
x,y,z;
}SHVector3d;
typedef
struct
{
SHVector3d sph;
SHVector3d vec;
double
*coeff;
}SHSample;
void
sphericalStratifiedSampling(SHSample* samples,
int
sqrtNumSamples,
int
nBands)
{
//Indexes
int
a,b,index,l,m;
//Loop Index
int
i=0;
//Inverse of the number of samples
double
invNumSamples=1.0/sqrtNumSamples;
//Cartesian and spherical coordinates
double
x,y,theta,phi;
//Loop on width of grid
for
(a=0;a<sqrtNumSamples,++a)
{
//Loop on height of grid
for
(b=0;b<sqrtNumSamples;++b)
{
//jitter center of current cell
x=(a+rnd())*invNumSamples;
y=(b+rnd())*invNumSamples;
//Calculate corresponding spherical angles
theta=2.0*
acos
(
sqrt
(1.0-x));
phi=2.0*PI*y;
//Sample spherical coordinates
samples[i].sph.x=theta;
samples[i].shp.y=phi;
samples[i].shp.z=1.0;
//Sample normal
samples[i].vec.x=
sin
(theta)*
cos
(phi);
samples[i].vec.y=
sin
(theta)*
sin
(phi);
samples[i].vec.z=
cos
(theta);
//Calculate SH coefficients of current sample
for
(l=0;l<nBands;++l)
{
for
(m=-l;m<=l;++m)
{
index=l*(l+1)+m;
samples[i].coeff[index]=evaluateSH(l,m,theta,phi);
}
}
++i;
}
}
}
//复合函数的定义
typedef
double
(*SphericalFunction)(
double
theta,
double
phi);
void
SHProjectSphericalFunction(SphericalFunction mySPFunc,
SHSample* samples,
double
*result,
int
numSamples,
int
numCoeffs)
{
//Weighting factor
double
dWeight=4.0*PI;
double
factor=dWeight/numSamples;
//Spherical angles
double
theta,phi;
//Loop indexes
int
i,n;
//Loop through all samples
for
(i=0;i<numSamples;++i)
{
//Get current sample spherical coordinates
theta=samples[i].sph.x;
phi=samples[i].sph.y;
//Update SH coefficients
for
(n=0;n<numCoeffs;++n)
{
result[n]+=mySPFunc(theta,phi)*samples[i].coeff[n];
}
}
for
(n=0;n<numCoeffs;++n)
result[n]*=factor;
}
|