Bingo, Computer Graphics & Game Developer

CentOS 7

修改SSH默认端口22

SSH远程登录部分可参见SSH原理与运用。具体终端命令可参见怎样修改 CentOS 7 SSH 端口, 本质为修改开放端口+防火墙修改

如若和我一样为阿里云服务器, 则需要找到控制台中的安全组规则,自定义的开放想要SSH连接的端口,否则在Xshell等工具中仍然无法连接。

QQ图片20180321174851

关闭账号密码登录

vi /etc/ssh/sshd_config

找到PasswordAuthentication yes,更改为no即可。此时再用一未信任的设备连接,可用终端命令登陆

ssh -p 22 hostName@address

提示以下错误信息则证明已无法直接通过账号密码直接登陆

Permission denied(publickey,gssapi-keyex,gssapi-with-mic).

Nginx搭建开发

Nginx的环境搭建可见Nginx Documentation,在配置SSH默认端口过程中可能将Nginx使用的默认80端口禁用,可参考Centos 7 firewall 命令在防火墙中打开对80的访问。

以免在服务器上预览开发效率过低,因此在本地也同时配置了Nginx,同时搭配Browsersync可以保证只要Web文件有修改时浏览器自动刷新,提高开发效率。

Nginx默认在80端口开放,因此启用代理模式可以实现类似Hexo主题开发中的效果。

browser-sync start --proxy "localhost:80" --files "**"

引号中的内容可更换为**.css等glob语法的正则表达式.

Git的自动化部署

Git的环境搭建可参见搭建Git服务器。由于更换了端口号,因此所有Git命令的对象都变为了类似git clone ssh://host@ip:port/xxx/xxx

在配置完Nginx后,为了能git push后将整个repo拷贝到Nginx的指定网页文件夹中(默认为Nginx/html), 可参考VPS服务器搭建Hexo博客教程中的安装配置git部分实现Git Hook的功能。

如若为Windows平台,可能遇到push后post-receive文件无法自动执行的问题,因为行尾Unix与Windows格式不一,可在Sublime或其他文本编辑器中修改行尾为UNIX即可解决(其他无法执行原因可直接执行该脚本查看)

具体可以参考bad interpreter: No such file or directory



理解来自CSDN中的递归删除算法,个人补充细节描述

题为「设计递归算法,删除无头结点的单链表L中值为x的结点」

void deleteX(LinkList &L, ElemType x)  
{  
    LNode *p;  
    if(L==NULL) return;  
    if(L->data==x)  
    {  
        p=L;  
        L=L->next; // (1)单链表保持连续 
        free(p);  
        deleteX(L,x); // (2)递归调用
    }
    else
        deleteX(L->next,x); //(3)向后遍历
}  

理解难度在于似乎仅仅只是释放了当前结点并使得L向前推进,而导致了断链。

例如要在单链表[1,3,5,7][1, 3, 5, 7]中删除元素xx(例如x=3x=3),这里递归第二次调用函数讲会遇到data==3data==3

核心理解在于这里的LinkList L是引用, 观察堆栈会非常清晰的说明这个问题。

deleteX

当代码中(2)(2)处递归调用deleteX时, 参数作为L->next的引用传入(引用为别名,相当于本体操作),即递归中参数L为上一函数的L->next本体,此时的代码相当于

deleteX(L->next, x)  
{  
    // ...
    p = L->next;  
    L->next = L->next->next; // 单链表保持连续 
    free(p);  
    deleteX(L->next, x);
    // ...
}  

如若将引用更换为指针,相较于引用的实现较为繁琐,但更能观察到上述代码工作的本质。

// 仅有传入指针的指针才能改变上一函数中指针的值
void deleteX(LinkList **L, ElemType x)  
{  
    LNode *p;  
    if((*L)==NULL) return;  
    if((*L)->data==x)  // L->next指向的结点为x
    {  
        p=*L;  // p指向L->next指向的结点
        *L=(*L)->next; // L->next=L->next->next
        free(p);  
        deleteX(L,x); // 传入L->next的地址
    }
    else
        deleteX(&((*L)->next),x); // 传入L->next指向的结点的next指针的地址
}  



坐标变换

通常摄像机动画,模型动画中会使用旋转平移变换来完成坐标变换。这里使用基变换来完成通用的坐标系转换。

现要求得世界坐标PworldP_{world}对应的代求坐标系Ω\Omega下坐标PΩP_{\Omega},已知Ω\Omega坐标轴向量分别为X,Y,Z\overset{\rightharpoonup}{X}, \overset{\rightharpoonup}{Y}, \overset{\rightharpoonup}{Z},坐标原点为OO

worldToObject[X.xY.xZ.xO.xX.xY.yZ.yO.yX.zY.zZ.zO.z0001]worldToObject
\left[
\begin{matrix}
\overset{\rightharpoonup}{X}.x & \overset{\rightharpoonup}{Y}.x & \overset{\rightharpoonup}{Z}.x & O.x \
\overset{\rightharpoonup}{X}.x & \overset{\rightharpoonup}{Y}.y & \overset{\rightharpoonup}{Z}.y & O.y \\
\overset{\rightharpoonup}{X}.z & \overset{\rightharpoonup}{Y}.z & \overset{\rightharpoonup}{Z}.z & O.z \\
0 & 0 & 0 & 1
\end{matrix}
\right]

即有PΩ=worldToObjectPworldP_{\Omega}=worldToObject * P_{world}.(该矩阵含义可展开其含义自明)

PBRT以及Atmos中关于相机各类坐标系的转换

PBRT中的NDC坐标系在Atmos中未使用,因为其本身在PBRT中并不直接体现

具体的转换管线流程如下,实现公式可由下图各自推导出来,不做赘述。

rasterPosition(imageX, imageY);

// raster to screen
screenPosition(rasterPosition.x - image->width / 2.0f, -rasterPosition.y + image->height / 2.0f);

// screen to camera
// 这里pbrt中使用投影矩阵的逆变换得到cameraPos,不过由于透视的特殊性,其方向向量是一致的。
cameraPosition.x = 2 * canvasSize.x * screenPosition.x / image->width;
cameraPosition.y = 2 * canvasSize.y * screenPosition.y / image->height;
cameraPosition.z = canvasDistance;

// camera to world
ray->direction = (cameraToWorld * cameraPosition).getNormalized();
ray->origin = cameraToWorld * zeroVector3;

之后可以用worldToCamera矩阵进行转换。

worldToCamera[right.xup.xdirection.xorigin.xright.xup.ydirection.yorigin.yright.zup.zdirection.zorigin.z0001]worldToCamera
\left[
\begin{matrix}
right.x & up.x & direction.x & origin.x \
right.x & up.y & direction.y & origin.y \\
right.z & up.z & direction.z & origin.z \\
0 & 0 & 0 & 1
\end{matrix}
\right]



本文运用上文中多维随机变量分布转换规律,来推导出PBRT中一系列重要结论


极坐标与笛卡尔坐标

已知p(x,y)=p(r,θ)Jt(r,θ)p(x, y)=\frac{p(r, \theta)}{\vert J_t(r, \theta) \vert},其中
T1(r,θ)=rcosθ=xT2(r,θ)=rsinθ=y
\begin{array}{l}
T_1(r,\theta)=rcos\theta=x \\
T_2(r,\theta)=rsin\theta=y
\end{array}

Jt(r,θ)=[cosθrsinθsinθrcosθ]=rJ_t(r, \theta)=\left[
\begin{matrix}
cos\theta & -rsin\theta \\
sin\theta & rcos\theta \\
\end{matrix}
\right] = r

则有p(r,θ)=rp(x,y)p(r, \theta) = rp(x, y)


球坐标与笛卡尔坐标

T1(r,θ,ϕ)=rsinθcosϕ=xT2(r,θ,ϕ)=rsinθsinϕ=yT3(r,θ,ϕ)=rcosθ=z
\begin{array}{l}
T_1(r,\theta, \phi)=rsin\theta cos\phi = x \\
T_2(r,\theta, \phi)=rsin\theta sin\phi = y \\
T_3(r,\theta, \phi)=rcos\theta = z
\end{array}

p(x,y,z)=p(r,θ,ϕ)Jt(r,θ,ϕ)p(x, y, z)=\frac{p(r, \theta, \phi)}{\vert J_t(r, \theta, \phi) \vert},其中

Jt(r,θ,ϕ)=[sinθcosϕrcosθcosϕrsinθsinϕsinθsinϕrcosθsinϕrsinθcosϕcosθrsinθ0]=r2sinθJ_t(r, \theta, \phi)=\left[
\begin{matrix}
sin\theta cos\phi & rcos\theta cos\phi & -rsin\theta sin\phi \\
sin\theta sin\phi & rcos\theta sin\phi & rsin\theta cos\phi \\
cos\theta & -rsin\theta & 0\\
\end{matrix}
\right] = r^2sin\theta

则有p(r,θ,ϕ)=r2sinθp(x,y,z)p(r, \theta, \phi) = r^2sin\theta p(x, y, z)


立体角与球坐标

立体角的定义为在单位圆上投影的面积,在球坐标中有dω=sinθdθdϕd\omega=sin\theta d\theta d\phi成立

因此基于物理Lambertion中的Ωf(ωi,ωo,p)cosθdω\int_{\Omega}{f(\omega_i, \omega_o, p)cos\theta d\omega}=1, 可得f(ωi,ωo,p)=102πdϕ0π2cosθsinθdθ=1πf(\omega_i, \omega_o, p) = \frac{1}{\int_{0}^{2\pi}d\phi \int_{0}^{\frac{\pi}{2}}{cos\theta sin\theta d\theta}} = \frac{1}{\pi}

若定义在某区域Ω\Omega上的概率为Pr{ωΩ}=Ωp(ω)dωPr\{ \omega \in \Omega \} = \int_{\Omega}{p(\omega)d\omega},也有Pr{(θ,ϕ)Ω}=Ωp(θ,ϕ)dθdϕPr\{ (\theta ,\phi) \in \Omega '\} = \int_{\Omega '}{p(\theta, \phi)d\theta d\phi}

因此
p(θ,ϕ)dθdϕ=p(ω)dωp(\theta, \phi)d\theta d\phi=p(\omega) d\omega
p(θ,ϕ)=sinθp(ω)p(\theta, \phi)=sin\theta p(\omega)


半球上立体角均匀采样

要对单位半球做关于立体角的均匀采样,则p(ω)=12πp(\omega)=\frac{1}{2\pi},根据上文中与球坐标关系,p(θ,ϕ)=sinθ2πp(\theta, \phi)=\frac{sin\theta}{2\pi}

则有两独立同均匀分布随机变量(ξ,ψ)[0,1](\xi,\psi)\in [0,1],要转换成满足半球上均匀分布的球坐标表示。

p(θ)=0π2p(θ,ϕ)dϕ=sinθ,p(ϕ)=02πp(θ,ϕ)dθ=12πp(\theta) = \int_{0}^{\frac{\pi}{2}}p(\theta, \phi)d\phi = sin\theta, \quad p(\phi) = \int_{0}^{2\pi}p(\theta, \phi)d\theta = \frac{1}{2\pi}

p(θϕ)=p(θ,ϕ)p(ϕ)=sinθ,p(ϕθ)=p(θ,ϕ)p(θ)=12πp(\theta \vert \phi) = \frac{p(\theta, \phi)}{p(\phi)}=sin\theta, \quad p(\phi \vert \theta) = \frac{p(\theta, \phi)}{p(\theta)}=\frac{1}{2\pi}

这里很明显可以看出来θ,ϕ\theta,\phi为两个独立随机变量

Pr{θθ}=0θp(θ)dθ=1cosθ,Pr{ϕϕ}=0ϕp(ϕ)dϕ=2πϕPr\{\theta' \leqslant \theta\} = \int_{0}^{\theta}p(\theta)d\theta=1-cos\theta , Pr\{\phi' \leqslant \phi\} = \int_{0}^{\phi}p(\phi)d\phi=2\pi \phi

应用反CDF变换法则,P1(θ)=1arccosθ,P1(ϕ)=2πϕP^{-1}(\theta)=1-arccos\theta, P^{-1}(\phi)=2\pi \phi,因此
θ=arccos(1ξ)=arccosξϕ=2πψ
\begin{array}{l}
\theta=arccos(1-\xi)=arccos\xi \\
\phi=2\pi \psi
\end{array}

由于球坐标不容易在计算机中表示,因此转换为笛卡尔表示

x=sinθcosϕ=1ξ2cos2πψy=sinθsinϕ=1ξ2sin2πψz=cosθ=ξ
\begin{array}{l}
x=sin\theta cos\phi=\sqrt{1-\xi^2}cos2\pi \psi \\
y = sin\theta sin\phi = \sqrt{1-\xi^2}sin2\pi \psi \\
z = cos\theta = \xi
\end{array}

同理可推证球上立体角均匀采样

x=2ξ(1ξ)cos2πψy=2ξ(1ξ)sin2πψz=12ξ
\begin{array}{l}
x=2\sqrt{\xi(1-\xi)}cos2\pi \psi \\
y = 2\sqrt{\xi(1-\xi)}sin2\pi \psi \\
z = 1-2\xi
\end{array}

下图为半球上均匀采样的结果与完整球上均匀采样的结果

hemisphere and sphere


单位圆上的均匀采样

要根据面积对单位圆上进行均匀采样,p(x,y)=1πp(x,y)=\frac{1}{\pi},则p(r,θ)=rp(x,y)p(r, \theta)=rp(x,y)

p(θ)=01rπdr=12π,p(rθ)=p(r,θ)p(θ)=2rp(\theta)=\int_{0}^{1}\frac{r}{\pi}dr=\frac{1}{2\pi}, \quad p(r\vert \theta)=\frac{p(r,\theta)}{p(\theta)}=2r

则他们对于的CDF为P(θ)=0θp(θ)dθ=θ2πP(\theta)=\int_{0}^{\theta}p(\theta)d\theta=\frac{\theta}{2\pi}P(rθ)=0rp(rθ)dr=r2P(r\vert \theta)=\int_{0}^{r}p(r\vert \theta)dr=r^2

求出其对应反函数P1(θ)=2πθP^{-1}(\theta)=2\pi\thetaP1(rθ)=rP^{-1}(r\vert \theta)=\sqrt{r}

若有两满足在[0,1][0, 1]上均匀分布的随机变量ξ,ψ\xi, \psiθ=2πξ,r=ψ\theta=2\pi \xi, r=\sqrt{\psi}

则转换回笛卡尔坐标结果为

x=ψcos2πξy=ψsin2πξ\begin{array}{l}
x = \sqrt{\psi}cos2\pi \xi\\
y = \sqrt{\psi}sin2\pi \xi
\end{array}

这里r,θr, \theta也是两个独立的随机变量

PBRT中论述的Concentric Disk Sampling可以参见前文Depth of Field中的实现与可视化


基于Cosine函数的半球采样

这里使用基本的转换思路,p(ω)cosθp(\omega)\propto cos\theta,因此p(ω)=cosθπp(\omega)=\frac{cos\theta}{\pi}

p(θ,ϕ)=sinθcosθπp(\theta, \phi)=\frac{sin\theta cos\theta}{\pi},则有

p(θ)=02πp(θ,ϕ)dϕ=2sinθcosθ,p(ϕθ)=p(θ,ϕ)p(θ)=12πp(\theta)=\int_{0}^{2\pi}p(\theta, \phi)d\phi=2sin\theta cos\theta, \quad p(\phi \vert \theta)=\frac{p(\theta, \phi)}{p(\theta)}=\frac{1}{2\pi}

可得CDF为

P(θ)=0θp(θ)dθ=1cos2θ2,P(ϕθ)=0ϕp(ϕθ)dϕ=ϕ2πP(\theta)=\int_{0}^{\theta}p(\theta)d\theta=\frac{1-cos2\theta}{2}, \quad P(\phi \vert \theta)=\int_{0}^{\phi}p(\phi \vert \theta)d\phi=\frac{\phi}{2\pi}

应用反CDF得P1(θ)=12arccos(12θ),P1(ϕθ)=2πϕP^{-1}(\theta)=\frac{1}{2}arccos(1-2\theta), P^{-1}(\phi \vert \theta)=2\pi \phi

θ=12arccos(12ξ)ϕ=2πψ
\begin{array}{l}
\theta = \frac{1}{2}arccos(1-2\xi)\\
\phi = 2\pi \psi
\end{array}

转换为笛卡尔下的表示为
x=sinθcosϕ=sin(12arccos(12ξ)))cos2πψy=sinθsinϕ=sin(12arccos(12ξ)))sin2πψz=cosθ=cos(12arccos(12ξ))\begin{array}{l}
x=sin\theta cos\phi = sin(\frac{1}{2}arccos(1-2\xi)))cos2\pi\psi \\
y=sin\theta sin\phi = sin(\frac{1}{2}arccos(1-2\xi)))sin2\pi\psi\
z=cos\theta = cos(\frac{1}{2}arccos(1-2\xi))
\end{array}

PBRT中给出了另一种计算办法,将均匀分布在圆盘上的点投影到半球上,其结果就满足Cosine-Weighted。

以下给出证明

假定在圆盘上的极坐标分布为p(r,ϕ)=rπp(r, \phi)=\frac{r}{\pi}(使用ϕ\phi方便后续将rr对应为θ\theta相关),p(r,ϕ)=p(θ,ϕ)Jt(r,θ)p(r, \phi)=\frac{p(\theta, \phi)}{\vert J_t(r, \theta)\vert}

r=sinθ=T1(r)ϕ=ϕ=T2(ϕ)
\begin{array}{l}
r=sin\theta = T_1(r) \\
\phi=\phi = T_2(\phi)
\end{array}

Jt(r,θ)=[cosθ001]=cosθJ_t(r, \theta)=\left[
\begin{matrix}
cos\theta & 0 \\
0 & 1 \\
\end{matrix}
\right] = cos\theta

即有p(θ,ϕ)=cosθp(r,ϕ)=rcosθπ=sinθcosθπp(\theta, \phi)=cos\theta p(r, \phi)=\frac{r cos\theta}{\pi} = \frac{sin\theta cos\theta}{\pi},和上述中得到的结果一样。

下图为两种方法产生的Cosine-Weighted对比

cartisan and projected


三角形均匀采样

PBRT中采用了等腰切两边为1的三角形特殊情况,不过以下计算办法可以变换回任何三角形

p(u,v)=2p(u,v)=2(面积倒数),则有

p(u)=01up(u,v)dv=22u,p(vu)=p(u,v)p(u)=11uP(u)=0up(u)du=2uu2,P(vu)=0vp(vu)dv=v1u
\begin{array}{l}
p(u)=\int_{0}^{1-u}p(u,v)dv=2-2u, \quad p(v\vert u)=\frac{p(u,v)}{p(u)}=\frac{1}{1-u} \\
P(u)=\int_{0}^{u}p(u)du=2u-u^2, \quad P(v\vert u)=\int_{0}^{v}p(v\vert u)dv=\frac{v}{1-u}
\end{array}

这里反变换CDF开根号时要注意到u,v[0,1]u,v\in [0,1]

因此有u=1ξ,v=ξψu=1-\sqrt{\xi},v=\sqrt{\xi}\psi

下图中即为三角形上的均匀采样

triangle sample


根据二维分段函数分布采样

这里分段函数的采样即为离散随机变量,若有横向纵向各nu,nvn_u,n_v个元素,函数上指定位置一点值为f(u,v)f(u,v),那么就有

If=vuf(u,v)dudv=1nunvi=0nu1j=0nv1f(ui,vj)p(u,v)=f(u,v)If
\begin{array}{c}
I_f=\int_v \int_u f(u,v)dudv=\frac{1}{n_un_v}\sum_{i=0}^{n_u-1}\sum_{j=0}^{n_v-1}f(u_i,v_j) \\
p(u,v)=\frac{f(u,v)}{I_f}
\end{array}

而边缘概率密度可直接得到
p(u)=vp(u,v)dv=1nuif(u,vi)If,p(vu)=p(u,v)p(u)=f(u,v)1nuif(u,vi)p(u)=\int_v p(u,v)dv=\frac{\frac{1}{n_u}\sum_i f(u, v_i)}{I_f} , \quad p(v\vert u)=\frac{p(u,v)}{p(u)}=\frac{f(u,v)}{\frac{1}{n_u}\sum_i f(u, v_i)}

之后利用CDF反变换,即PBRT中distribution1d的实现办法可完成对于p(u),p(vu)p(u),p(v\vert u)上的采样

此处效果可见前文中的二维分段函数采样的结果


圆锥均匀采样

类似于球上的均匀采样, 易证随机变量θ,ϕ\theta, \phi相互独立

已知p(ω)=12π(1cosθmax)p(\omega)=\frac{1}{2\pi(1-cos\theta_{max})},则p(θ,ϕ)=sinθ2π(1cosθmax)p(\theta, \phi)=\frac{sin\theta}{2\pi(1-cos\theta_{max})}

p(θ)=02πp(θ,ϕ)dϕ=sinθ1cosθmax,p(ϕ)=02πsinθ2π(1cosθmax)dθ=12πp(\theta)=\int_{0}^{2\pi}p(\theta, \phi)d\phi=\frac{sin\theta}{1-cos\theta_{max}}, \quad p(\phi)=\int_{0}^{2\pi}\frac{sin\theta}{2\pi(1-cos\theta_{max})}d\theta=\frac{1}{2\pi}

P(θ)=0θsinθ1cosθmaxdθ=1cosθ1cosθmax,P(ϕ)=0ϕ12πdϕ=ϕ2πP(\theta)=\int_{0}^{\theta}\frac{sin\theta}{1-cos\theta_{max}}d\theta=\frac{1-cos\theta}{1-cos\theta_{max}}, \quad P(\phi)=\int_{0}^{\phi}\frac{1}{2\pi}d\phi=\frac{\phi}{2\pi}

若有两随机变量xi,ψxi, \psi,运用反变换可有(此处省略了θ,ϕ\theta, \phi到笛卡尔的转换)

cosθ=1ξ+ξcosθmax,ϕ=2πψcos\theta=1-\xi +\xi cos\theta_{max}, \quad \phi=2\pi\psi

下图中即为圆锥上的均匀采样

cone sampling



关于连续型随机变量的分布转换,除了使用反CDF函数将均匀分布函数转换为其他分布函数以外,存在一种通用办法实现不同分布之间的转换,不只局限于均匀分布。

假定存在随机变量XX,其累积分布函数为FX(x)F_X(x)。若存在单调函数gg,使得Y=g(X)Y=g(X)(因为是单调函数,因此XXYY的关系为单射),则为了计算YY的累积分布函数FY(y)F_Y(y)

g(x)g(x)此时为严格单调递增

FY(y)=Pr{Yy}=Pr{g(X)y}=Pr{Xg1(y)}=FX(g1(y))
\begin{array}{l}
F_Y(y) & = & Pr\{ Y\leqslant y \} \\
& = & Pr\{g(X)\leqslant y\} \\
& = & Pr\{ X\leqslant g^{-1}(y) \} \\
& = & F_X(g^{-1}(y))
\end{array}

此时对终式左右对x求导(累积分布函数导数为概率密度函数)

fY(y)dydx=fX(g1(y))fY(y)=(dydx)1fX(g1(y))
\begin{array}{l}
f_Y(y)\frac{dy}{dx}=f_X(g^{-1}(y)) \\
f_Y(y) = (\frac{dy}{dx})^{-1}f_X(g^{-1}(y)) \\
\end{array}

再设g(x)g(x)此时为严格单调递减

FY(y)=Pr{Yy}=Pr{g(X)y}=1Pr{Xg1(y)}=1FX(g1(y))
\begin{array}{l}
F_Y(y) & = & Pr\{ Y\leqslant y \} \\
& = & Pr\{g(X)\leqslant y\} \\
& = & 1 - Pr\{ X\leqslant g^{-1}(y) \} \\
& = & 1 - F_X(g^{-1}(y))
\end{array}

同样的有
fY(y)=(dydx)1fX(g1(y))f_Y(y) = -(\frac{dy}{dx})^{-1}f_X(g^{-1}(y))

则有fY(y)=(dydx)1fX(g1(y))f_Y(y) = \vert (\frac{dy}{dx})^{-1} \vert f_X(g^{-1}(y))

PBRT中讲到,当我们有一概率密度为px(x)p_x(x)的随机变量XX,想要转换为满足已知分布为py(y)p_y(y)的随机变量YY,那么即有下式成立(此过程与上式正好相反)

FX(x)=Pr{Xx}=Pr{g(X)g(x)}=Pr{Yg(x)}=FY(g(x))
\begin{array}{l}
F_X(x) & = & Pr\{ X\leqslant x \} \\
& = & Pr\{g(X)\leqslant g(x)\} \\
& = & Pr\{ Y\leqslant g(x) \} \\
& = & F_Y(g(x))
\end{array}

g(x)=FY1(FX(x))g(x) = F_Y^{-1}(F_X(x))或当g(X)g(X)为单调减函数时g(x)=FY1(1FX(x))g(x) = F_Y^{-1}(1 - F_X(x)),PBRT中未对此处进行详细说明。

因此当XX为满足[0,1][0, 1]区间均匀分布的随机变量时,FX(x)=xF_X(x)=x,那么上式中FX(x))F_X(x))或是1FX(x)1 - F_X(x)都满足[0,1][0, 1]区间均匀分布,则g(x)=FY1(x)g(x)=F_Y^{-1}(x),即Inverse Transform Sampling得证。

多维随机变量的转换T(x)T(x)为双射变换,则满足fY(y)=fY(T(x))=fX(x)JT(x)f_Y(y)=f_Y(T(x))=\frac{f_X(x)}{\vert J_T(x)\vert },且T(x)=(T1(x),Tn(x))T(x)=(T_1(x),\cdots T_n(x))

JT(x)=[T1x1T1xnTnx1Tnxn]J_T(x)=\left[
\begin{matrix}
\frac{\partial T_1}{\partial x_1} & \cdots & \frac{\partial T_1}{\partial x_n} \\
\vdots & \ddots & \vdots \\
\frac{\partial T_n}{\partial x_1} & \cdots & \frac{\partial T_n}{\partial x_n} \\
\end{matrix}
\right]



蒙特卡洛方法求积分以及部分所需的概率论术语可以在PBRT以及概率论与数理统计中找到

PBRT中给出了蒙特卡洛估计量,其中随机变量Xi[a,b]X_i \in [a,b]且独立同分布,分布满足概率密度函数 p(x)p(x)

FN=1Ni=1Nf(Xi)p(Xi)F_N = \frac{1}{N}\sum^{N}_{i=1}\frac{f(X_i)}{p(X_i)}

其期望为

E[FN]=E[1Ni=1Nf(Xi)p(Xi)]=1Ni=1Nabf(x)p(x)p(x)dx=1Ni=1Nabf(x)dx=abf(x)dx\begin{array}{l}
E[F_N] & = E[\frac{1}{N}\sum_{i=1}^{N}\frac{f(X_i)}{p(X_i)}] \\
& = \frac{1}{N}\sum_{i=1}^{N}\int_{a}^{b}\frac{f(x)}{p(x)}p(x)dx \\
& = \frac{1}{N}\sum_{i=1}^{N}\int_{a}^{b}f(x)dx \\
& = \int_{a}^{b}f(x)dx \\
\end{array}

方差为

σ2=1N(f(x)p(x)I)2p(x)dx\sigma^2=\frac{1}{N}\int(\frac{f(x)}{p(x)}-I)^2p(x)dx

结果的误差与标准差成正比,因此随着样本数增加误差缩小速度仅与N\sqrt{N}相关(即常说的增加4倍采样数目才能缩小一半的误差)

PBRT以及许多其他文章书籍等直接给出了这个估计量的期望计算(其中第二步到第三步并未直接说明来由,尽管PBRT在对期望的简介中已给出式子),这一步只能证明估计量本身无偏*

以下将给出证明

引用一个知识点: Law of the unconscious statistician,简称LOTUS。用法是已知随机变量X的分布fXf_X,但并不知道函数g(X)的分布fg(X)f_{g(X)}。那么此时函数g(X)g(X)的期望为

E[g(X)]=xg(x)fX(x)E[g(X)]=\sum_x g(x)f_X(x)(X为离散型随机变量)

E[g(X)]=g(x)fX(x)dxE[g(X)]=\int_{-\infty}^{\infty} g(x)f_X(x)dx(X为连续型随机变量)

此处细节介绍也可以看Wyman的技术博客,其中就提到了许多地方都没有涉及到的LOTUS。

E[f(Xi)pdf(Xi)]=E[f(x)pdf(x)]=abf(x)pdf(x)pdf(x)dx=abf(x)dx\begin{array}{l}
E[ \frac {f(X_{i})}{pdf(X_{i})} ] & = E[ \frac {f(x)}{pdf(x)} ] \\
& = \int _{a}^{b}\frac {f(x)}{pdf(x)}pdf(x)dx \\
& = \int _{a}^{b}f(x)dx
\end{array}

这里积分区间变为[,][-\infty, \infty]也是如此。以下将会使用到大数定律的知识背景。

引用自维基百科:设 a1,a2,...,an,...a_{1},a_{2},…,a_{n},…为相互独立的随机变量,其数学期望为: E(ai)=μ(i=1,2,...)E(a_{i})=\mu (i=1,2,…),方差为: Var(ai)=σ2(i=1,2,...)Var(a_{i})=\sigma ^{2}(i=1,2,…)
则序列a=1ni=1n\overline{a} = \frac{1}{n} \sum_{i=1}^{n}aia_{i} 依概率收敛于 μ\mu(即收敛于此数列的数学期望 E(ai)E(a_{i}))。

取一组独立同分布的随机变量{ξ}\{\xi\},且ξ\xi[a,b][a,b]内满足分布律p(x)p(x),则令p(x)=f(x)p(x)p^(x)=\frac{f(x)}{p(x)},**则{p(ξi)}\{ p^(\xi_i) \}也是一组独立同分布的随机变量。**,那么计算其得到的期望其实就是积分本身(见下)。

由大数定理

Pr(limN1Ni=1Np(ξi)=I)=1Pr(\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}p^*(\xi_i)=I)=1

那么这里已经有随机变量{p(ξi)}\{ p^(\xi_i) \}存在,且期望等于积分值,**那么根据大数定理可知,这里的{p(ξi)}\{ p^(\xi_i) \}就是会依概率收敛于积分值。**证毕。



Bingo

@BentleyJobs

Graduated from JNU, interested in cg & game developing, once worked in Aurogon to develope GuJianOL. See about for more info.