光子传输的蒙特卡罗方法并行实现.doc

光子传输的蒙特卡罗方法在 GPU 上的实现L‘ aszl 'o Szirmay-Kalos 等翻译：吴涛这篇文章提出一个快速的蒙特卡罗方法来解决激光和 gamma 射线光子在不均匀 媒介中辐射传输方程问题激光传输在计算机图表算法中起着重大作用，比如高 能 gamma 射线就在医疗与物理中扮演重要角色实时图像的应用在速度上受到 限制，因为我们通常想得到每秒至少 20 张图来获取运动的动态幻像快速模拟 在人机交互系统中也十分重要，比如像放射疗法或者物理实验设计，这些地方用 户可以放置一个源，并且希望得到辐射分配的迅速反馈多重散射模拟的速度在 迭代断层摄像术重建方面也很重要，此时散射辐射要从实际的重建数据中，且要 移除测量部分来估计，而重建仍然继续针对那些假设只代表不散射的部分剩余成 分17.1 光子传输的物理特性计算多重散射以及提出一种在不均匀媒介的实际方法更是困难重重最精确的逼 近方法是依据蒙特卡罗方法的积分运算以及通过跟踪随机光子的物理模拟现象， 来计算他们的总贡献光子与媒介粒子在该点碰撞的条件概率密度由消光系数C （v,X）定义，如果t粒子密度不均匀，它则由光子频率V以及位置X来确定。

消光系数通常可以用材 料密度c（ X）和一个仅与频率有关的因子的乘积来表达光子在碰撞之外还会被反射与吸收反射的概率称为反照率，用a（v ）表示随机反射方向用两个球坐标轴描述：散射角0和方位角p （图17.1）源点与新的方向的散射角0的概率密度可能与频率有关故用物理模型来描述（我们用瑞利散射【12】模型来模拟激光光子，用Klein-Nishina公式【11】【13】 来计算gamma光子）方位角p的概率密度鉴于轴向对称分布故而是均衡的光子在媒介中的散射能量为 E 的入射光子与材料的粒子碰撞因为碰撞，光子的方向和 0能量将改变新的方向由散射角0和方位角9指定新的光子能量E由入射光子能量和散 射角确定光子能量E = hh正比于它的频率，比值是普朗克常量h能量（也就是稍高 能量的 gamma 光子的频率）会因为散射而改变由康普顿定律，散射角与光子 能量的相对变化有独特的对应关系：EE 01 + 1(1 一 cos0 )m c 2其中，E是散射光子能量，E是入射光子能量,0m 是电子静止质量， e光速检查康普顿定律公式，我们得出当E关于电子能量（mc2）不可忽略时, 0e能量在散 射时 的变 化很大 。

但 是对于可见 光而言，就不再适用 了， 此时E =h □ mc ；因而，它们的散射能量将与入射光相似结果，我们能够按不0e同频率分解低能量（激光）辐射模拟，但在高能量（gamma射线），频率将成为 每次模拟光子的特征属性在计算机模拟中，光子根据目前提到的基本定律来跟踪模拟过程的输入数 据是媒介的密度，其由三维标量场密度C（X）定义，还有仅与频率有关的反照率 函数，以及散射角的概率密度场密度假定是离散数据像三维组织立体像素一样 仅与频率有关的反照率函数，以及散射角的概率密度函数即可以像表格一样存放 在组织中也可以用代数公式表示为了模拟光子传输过程，碰撞，吸收和散射方向需要根据物理模型定义的概 率密度来生成一种方法是根据给定的概率密度生产随机变量，称为反演法反 演法首先计算累计概率分配（CDF）,作为概率密度的积分，然后通过反演单位 间隔均匀分布的变量CDF来生成离散抽样在形式上,如果一个随机变量的CDF 是P（x），则随机抽样x可以通过解方程来产生：P（ x） r这里r是一个单位间隔内均匀分布的随机变量，通常可以调用随机数产生器 来获得17.2光子传输在GPU上的实现光子在空间的传输是相互独立的，因而可以用并行模拟，因此显示出其固有 的物理模型的并行性（自然是大规模并行的机器）。

模拟单个粒子路径的算法是 嵌套循环外循环的主体是负责寻找下一个散射位置，于是抽样吸收，如果还有 残存就抽样新的散射方向这层循环在光子被吸收或者离开模拟区域时结束不 同粒子的散射事件发生数目将会在一个大的时间间隔内波动，因而外层循环也会 以同样的次数执行内层循环一直访问体元直到下一个散射点确定；这个过程取 决于拉伸的随机数同样与沿着光子路径的媒介密度有关随机数与密度变化，因 此对于不同的抽样，访问的体元数目可能也是不同的图 17.2）图 17.2在光子传输模拟中控制路径的原因光子路径有随机数目散射点而且在两个散射点间又有不同数目的体元在并行线程中执行不同数目的迭代的问题在 GPU 中变得起决定性作用如 果我们将一个光子分配到GPU线程的技术单元，那么在warp中的最长的光子路 径与所有路径中最长的射线将会限制计算的速度为了保持 warp 的连续性，线 程执行的长度必须相似，所有的线程优先执行相同的指令序列为了克服这个困 难，我们根据 GPU 执行效率的需求调整传输过程的模型；也就是说，我们发展 了能在最少条件下执行的算法为了解决不同数目的散射事件，我们重复光子步骤，也就是计算线程中的射 线。

每一步包括标示另一个散射点，做终止决定，如果光子在碰撞后仍存活则继 续抽样新的方向，反之在源区生成新的光子这种策略能够保证线程都在忙碌运 行中经过调整，线程仍可能分歧确定一个光子是否吸收或者抽样一个新的散射 方向要根据当前的媒介密度；因而，这些任务需要解一个仅含少量变量的方程， 这个可以从实际位置确定即使因为 CDF 不能计算或者用符号表示反演而变得 复杂在这种情况下，数值方法可以得到应用自由程抽样甚至更差因为它需要 探查媒介的大部分区域，相当于大部分当前的 GPU 内存通道变得相对缓慢迭 代数据的根和访问不同密度的区域体元导致线程分歧最终降低GPU的利用率因为不同光子的跟踪是随机的，它们访问不同模拟体元的部分时，导致数据 存储在不同并行通道缓冲中这种方式的内存使用使得高速缓存的利用率递降 为了解决这个问题，也许有人会建议将类似路径的光子分配到一个 warp 中因 为随机路径是从均匀分布的随机数中产生的，这可能意味着随机数向量的分类 一个完整的路径可能联合了许多随机数（比如，超过 20 个），所有分类操作需要 在一个很高的维度工作，这在实际中是不可行的所有光子的存储是依据第一个 随机数并且分配到已经分好序列的 warp 中。

这一步提供很少的初始化数据，但 是早晚所有的线程必然能同时访问体元的不同部分为了保证线程在媒介中追踪射线的连续性，我们发展了一个无条件说明的方 案来抽样吸收和散射方向为了解决自由程问题同样为了限制数据通道的分歧， 我们提出一项技术，即使用低分辨率的体元数组在除原始的高分辨率的体元数组 的基础上这项技术显著的减少了到达高分辨率数组的组织数以及发现新的散射 点需要迭代的数目17.2.1 自由程抽样模拟一个光子的当前位置x，方向®，以及频率V，就是抽样自由程s (就是一 个光子在不碰撞下可以移动的距离)，也就是决定沿着射线Q(s)二x+Es的下一个 交互作用点光子与粒子发生碰撞的概率密度就是消光系数从这里，自由程 CDF可以用指数衰减模式来表达；因而自由程抽样等于下面关于标量r即单位间 隔均匀分布的量的解：1 - exp -fo (V, p(l))dl 二 rtI o 丿消光系数的积分称为光学深度，从而可以导出下面的方程：t (0, s)二(V, ~p(l ))dl 二一log(l- r)t0自由程抽样，也就是这个方程的解，当媒介均匀时形式很简单；换句话说当消光系数在媒介中处处相同时：—log(1- r)s = o (V)t当媒介不均匀时，消光系数不是常数，但是可以用一个数组表示；这种情况下，通常是方法是射线形变程即沿着射线每次行进一小步As，发现步数n在黎曼求和逼近光学深度且比-log(1- r)稍大：艺o (V, p(iAs))As < - log(1- r) W 工o (V, p(iAs))Astti=0 i=0步长As设置为体元的边缘长度，以便在穿程过程中略过体元。

不幸地是， 这种算法需要提取许多体元数组，特别是当体元数组是很大的而且平均消光很小 时更糟糕的是，组织提取的数目变化剧烈，取决于随机的射线长度再次回到 自由程抽样是一个简单公式，当在均匀媒介中时，光学深度是一个简单是代数表 达式然而，如果光学深度只能通过许多数据来计算，这种情况发生在消光系数 是由一个高分辨率数组确定，因而取样过程变得很慢幸运的是，如果我们至少 有一部分消光系数数据，可以用来减少获取数据的数目，这样我们就不需要付出 所有代价自由程抽样等于方程17.1关于路径长s的解，这里我们想用一个简单的方程 来替代光学深度，此方程可以用很少的变量计算为达到这个目的，我们添加一 些虚拟的物质或粒子来改变材料，从某种程度上可以使总密度服从一个简单的方 程有人可能会认为将材料混入原始物可能会改变体元内部的辐射场强度，导致 一个错误的结果，这当然不是我们期望的幸运的是，如果我们虚拟物质的另外 两个自由属性，即反照率和散射角概率密度，可以适当的定义，这样就完全不必 担心这种情况发生了虚拟物质不能改变媒介内部的场密度；就是说，决不能改变能力以及光子散 射时的方向这些要求的满足条件是，虚拟物质的反照率为 1，即 100%会反射 光子，散射角概率密度是一个为0的Dirac-delta 换句话说，就是完全不会改变光子方向。

更正式一些，我们通过混合不同虚拟物质掌握不同种类的体元来增加消光系 数，使其达到以及简单的上界函数a (v,X)对于原始的消光系数Q (v,X)，我 max t们确定虚拟物质的消光系数a (v,X)，这样在混合了实际物质和虚拟物质粒子的v消光系数就是a (v,X) =a (v,X) +a (v,X)max t v当碰撞发生在这种混合物质时，我们要确定到底发生在实际粒子上还是在虚 拟的粒子上抽样需要依据指定的概率密度来生成随机点，因而，足以用合适的 概率密度随机的解决这个问题因为消光参量定义了散射概率密度，比率a (v,X)/a (v,X)和a (v,X)/a (v,X)，通过该值我们有可能分别确定散射是发t max v max生在一个实际粒子还是虚拟粒子上通过加入虚拟粒子，我们可以通过下面的步骤来执行自由程抽样：(图 17.3)Higl> Low-dEnsity rtensiTy Constant densityregion region图 17.3虚拟粒子运用左图为原始实际粒子的自由程抽样(颜色为浅灰)右图描述了整个体元包 括虚拟粒子(黑色)和实际粒子的情况与虚拟粒子碰撞不改变光子方向。

1. 产生试验性的路径长s和试验性的碰撞位置Q(s)，利用上界消光系数函数c (v,X)这需要抽样最大消光的方程式的解：maxt (0, s) = (V, 0(l))dl = —log(1— r)max max02. 当一个试验性的碰撞点确定时，我们用概率c (v,X)/c (v,X)随机决定散射t max发生在实际粒子还是虚拟粒子上要是虚拟散射发生，那这个粒子的散射方 向不会改变，之后相同的抽样步骤将会从原点再次重复当一个实际粒子散 射发生时，这层循环结束注意到作为实际算法的抽样表明，通过原始消光系数 CDF 产生的，无论抽 样中建立了什么样的上界消光系数然而，抽样的效率受合适选择的上界系数影 响巨大一方面，上界系数需要一个针对方程 17.2 简单的解另一方面，如果 这个系数与实际系数，即虚拟粒子密度相比，差异大，虚拟散射发生的概率也会 很大，因此需要许。