CN105118083A

CN105118083A - Unbiased photon mapping drawing method

Info

Publication number: CN105118083A
Application number: CN201510489335.9A
Authority: CN
Inventors: 侯启明; 秦昊; 孙鑫; 周昆
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2015-08-11
Filing date: 2015-08-11
Publication date: 2015-12-02
Anticipated expiration: 2035-08-11
Also published as: CN105118083B

Abstract

本发明公开了一种无偏的光子映射绘制方法。光子映射被公认为是最为有效的绘制方法之一，然而之前的光子映射算法在收集光子的计算上是有偏的，导致最终绘制的结果与真实的结果之间不同。本发明通过将视线与光子所在的光路直接连接，准确地计算该光路传递的光能，更无偏地估计了视线采集光子的概率，从而以很小的额外代价得到了一个无偏的光子映射绘制方法。此外，进一步开发了一套多重重要性采样的权重，从而将无偏光子映射方法与双向光线跟踪技术相结合。本发明方法效率以及可靠性高。本发明方法的广泛应用，有望降低动漫、电影、广告以及建筑行业中图像渲染绘制的成本。The invention discloses an unbiased photon mapping drawing method. Photon mapping is recognized as one of the most effective rendering methods. However, previous photon mapping algorithms are biased in the calculation of collected photons, resulting in differences between the final rendering results and the real results. The invention directly connects the line of sight with the optical path where the photons are located, accurately calculates the light energy transmitted by the light path, and estimates the probability of collecting photons by the line of sight more unbiasedly, thereby obtaining an unbiased photon mapping at a small extra cost drawing method. Furthermore, a set of weights for multiple importance sampling is further developed to combine the unbiased photon mapping method with bidirectional ray tracing techniques. The method of the invention has high efficiency and reliability. The wide application of the method of the invention is expected to reduce the cost of image rendering and drawing in animation, film, advertisement and construction industries.

Description

An Unbiased Photon Mapping Drawing Method

技术领域technical field

本发明涉及图形绘制技术领域，尤其涉及一种光子映射的绘制方法。The invention relates to the technical field of graphic drawing, in particular to a photon mapping drawing method.

背景技术Background technique

本发明相关的研究背景简述如下：The relevant research background of the present invention is briefly described as follows:

一、光子映射1. Photon mapping

光子映射的核心思想是利用光子的密度来估计局部的光能，具体可参考JENSEN,H.W.2001.RealisticImageSynthesisUsingPhotonMapping.A.K.Peters,Ltd.,Natick,MA,USA。光子收集的局部性使得该方法对于具有挑战性的焦散以及多次反射的间接光照特别有效。此外，光子映射可以重用光子所对应的光路，从而减小了采样的代价。密度估计可以有效地降低噪声，但是同时也会给最终的结果引入偏差，因此之前所有的光子映射方法都是有偏的。虽然最新的一系列渐进光子映射算法可以在光子数量无限的时候一致收敛到正确的结果，具体方法可以参考HACHISUKA,T.,OGAKI,S.,ANDJENSEN,H.W.2008.Progressivephotonmapping.ACMTrans.Graph.27,5(Dec.),130:1–130:8；KNAUS,C.,ANDZWICKER,M.2011.Progressivephotonmapping:Aprobabilisticapproach.ACMTrans.Graph.30,3(May),25:1–25:13；HACHISUKA,T.,ANDJENSEN,H.W.2009.Stochasticprogressivephotonmapping.ACMTrans.Graph.28,5(Dec.),141:1–141:8；KAPLANYAN,A.S.,ANDDACHSBACHER,C.2013.Adaptiveprogressivephotonmapping.ACMTrans.Graph.32,2(Apr.),16:1–16:13。但是无一例外的这些方法在任意有限时间内的结果都是有偏的，而我们的方法自始至终都能得到无偏的绘制结果。也有工作通过多重重要性采样将有偏的光子映射与双向光线跟踪相结合，参考VORBA,J.2011.Bidirectionalphotonmapping.InProc.oftheCentralEuropeanSeminaronComputerGraphics(CESCG11)；TOKUYOSHI,Y.2009.Photondensityestimationusingmultipleimportancesampling.InACMSIGGRAPHASIA2009Posters,ACM,NewYork,NY,USA,SIGGRAPHASIA’09,37:1–37:1。而我们的工作改进了之前的多重重要性采样，首次使之将无偏的光子映射与其他的采样技术相结合。The core idea of photon mapping is to use the density of photons to estimate local light energy. For details, please refer to JENSEN, H.W.2001. Realistic Image Synthesis Using PhotonMapping. A.K.Peters, Ltd., Natick, MA, USA. The locality of photon collection makes this method particularly effective for challenging caustics and indirect lighting with multiple bounces. In addition, photon mapping can reuse the optical path corresponding to the photon, thereby reducing the cost of sampling. Density estimation can effectively reduce noise, but it also introduces bias to the final result, so all previous photon mapping methods are biased. Although the latest series of progressive photon mapping algorithms can consistently converge to the correct result when the number of photons is infinite, the specific method can refer to HACHISUKA, T., OGAKI, S., ANDJENSEN, H.W.2008.Progressivephotonmapping.ACMTrans.Graph.27, 5 (Dec.), 130:1–130:8; KNAUS, C., ANDZWICKER, M. 2011. Progressive photon mapping: Aprobabilistic approach. ACMTrans. Graph. 30, 3 (May), 25:1–25:13; HACHISUKA, T., ANDJENSEN, H.W.2009.Stochasticprogressivephotonmapping.ACMTrans.Graph.28,5(Dec.),141:1–141:8; KAPLANYAN,A.S.,ANDDACHSBACHER,C.2013.Adaptiveprogressivephotonmapping.ACMTrans.Graph.32,2( April.), 16:1–16:13. But without exception, the results of these methods are biased in any finite time, while our method can get unbiased drawing results from beginning to end. There is also work combining biased photon mapping with bidirectional ray tracing through multiple importance sampling, refer to VORBA, J.2011. Bidirectional photon mapping. InProc. of the Central European Seminaron Computer Graphics (CESCG11); TOKUYOSHI, Y.2009. NY, USA, SIGGRAPHASIA '09, 37:1–37:1. Our work improves upon previous multiple importance sampling to combine unbiased photon mapping with other sampling techniques for the first time.

二、双向光线跟踪以及多重重要性采样2. Two-way ray tracing and multiple importance sampling

Kajiya在1986年提出的渲染方程是对于光能传播的物理规律的数学描述，该方程是一个多重积分，可以通过蒙特卡洛的方法求解，该渲染方程可以参见KAJIYA,J.T.1986.Therenderingequation.SIGGRAPHComput.Graph.20,4(Aug.),143–150。双向光线跟踪则是通过分别采样视线以及光线，并将两条光路相连的方法来对完整的光路进行采样。该方法对于室内场景以及有限大小的光源和相机特别有效，具体方法可以参见LAFORTUNE,E.P.,ANDWILLEMS,Y.D.1993.Bi-directionalpathtracing.InProceedingsofCompuGraphics,vol.93,145–153；VEACH,E.,ANDGUIBAS,L.J.1994.Bidirectionalestimatorsforlighttransport.InProceedingsoftheFifthEurographicsWorkshoponRendering,Eurographics,147–162。The rendering equation proposed by Kajiya in 1986 is a mathematical description of the physical law of light energy propagation. This equation is a multiple integral that can be solved by the Monte Carlo method. The rendering equation can be found in KAJIYA, J.T.1986.Therenderingequation.SIGGRAPHComput. Graph. 20, 4 (Aug.), 143–150. Two-way ray tracing is to sample the complete light path by sampling the line of sight and light separately, and connecting the two light paths. This method is particularly effective for indoor scenes and light sources and cameras of limited size. For details, see LAFORTUNE, E.P., ANDWILLEMS, Y.D.1993. Bi-directional path tracing. InProceedings of CompuGraphics, vol.93, 145–153; Bidirectional estimators for light transport. In Proceeding of the Fifth Eurographics Workshop on Rendering, Eurographics, 147–162.

之后多重重要性采样也被Veach等人提出，当时主要用于在双向光线跟踪之中结合各种光路不同的连接方式，参见VEACH,E.1998.RobustMonteCarloMethodsforLightTransportSimulation.PhDthesis,Stanford,CA,USA.AAI9837162。光子映射和双向光线跟踪都是通过连接视线和光线形成完成的光路采样。相比于双向光线跟踪而言，光子映射由于视线终点和光子非常接近，因此在采样比较复杂的焦散和多次反射的间接光照比较有优势；然而在直接光照等情况下，双向光线跟踪则更加有效。因此这两种方法在不同的类型的光路采样中是互补的，具体可以参见HASANˇ,M.,KRIVˇANEK′,J.,WALTER,B.,ANDBALA,K.2009.Virtualsphericallightsformany-lightrenderingofglossyscenes.ACMTrans.Graph.28,5(Dec.),143:1–143:6；VORBA,J.2011.Bidirectionalphotonmapping.InProc.oftheCentralEuropeanSeminaronComputerGraphics(CESCG11)。Later, multiple importance sampling was also proposed by Veach et al. At that time, it was mainly used to combine different connection methods of various optical paths in bidirectional ray tracing. See VEACH, E.1998.RobustMonteCarloMethodsforLightTransportSimulation.PhDthesis, Stanford, CA, USA.AAI9837162. Both photon mapping and bidirectional ray tracing are done by linking sight lines and rays to form light path sampling. Compared with two-way ray tracing, photon mapping has advantages in sampling complex caustics and indirect lighting with multiple reflections because the end of the line of sight is very close to photons; however, in situations such as direct lighting, two-way ray tracing is more effective. Therefore, these two methods are complementary in different types of light path sampling. For details, please refer to HASANˇ,M.,KRIVˇANEK′,J.,WALTER,B.,ANDBALA,K.2009.Virtualspherical lights form any-light rendering of glossy scenes.ACMTrans.Graph. 28, 5 (Dec.), 143:1–143:6; VORBA, J. 2011. Bidirectional photon mapping. In Proc. of the Central European Seminaron Computer Graphics (CESCG11).

三、统一采样3. Unified sampling

顶点合并(VertexMerging)将光子映射重新定义为一个光路以一定概率进行双向连接的过程，其中连接的概率由被连接的两点之间的距离所决定，GEORGIEV,I.,KRIVANEK′,J.,DAVIDOVIC,T.,ANDSLUSALLEK,P.2012.Lighttransportsimulationwithvertexconnectionandmerging.ACMTrans.Graph.31,6(Nov.),192:1–192:10。统一光路采样(UnifiedPathSampling)则是同一个问题的另一种角度，该工作将双向光线跟踪认为是一种光子和视线终点完全重合的光子映射Vertex Merging redefines photon mapping as a process of two-way connection of a light path with a certain probability, where the probability of connection is determined by the distance between the two connected points, GEORGIEV, I., KRIVANEK′, J., DAVIDOVIC, T., ANDSLUSALLEK, P. 2012. Light transport simulation with vertex connection and merging. ACM Trans. Graph. 31, 6 (Nov.), 192:1–192:10. Unified Path Sampling (UnifiedPathSampling) is another angle of the same problem. This work regards bidirectional ray tracing as a photon mapping in which photons and the end of the line of sight completely coincide.

发明内容Contents of the invention

本发明的目的在于针对现有技术的不足，提供一种无偏的光子映射绘制方法，在充分发挥光子映射绘制方法高效率的同时，消除了困扰该方法已久的误差问题。不但在理论上突破了原有光子映射方法的局限，并且具有很高的实用价值。The purpose of the present invention is to provide an unbiased photon mapping drawing method for the deficiencies of the prior art, which eliminates the long-standing error problem of the photon mapping drawing method while giving full play to the high efficiency of the photon mapping drawing method. It not only breaks through the limitations of the original photon mapping method in theory, but also has high practical value.

本发明的目的是通过以下技术方案来实现的：一种无偏的光子映射绘制方法，包括以下步骤：The object of the present invention is achieved through the following technical solutions: a kind of unbiased photon mapping drawing method, comprises the following steps:

(1)输入三维场景文件，并对三维场景文件进行解析；所述三维场景文件包括物体的几何信息、材质以及贴图、灯光信息和相机设置；(1) input three-dimensional scene file, and three-dimensional scene file is analyzed; Described three-dimensional scene file comprises the geometric information of object, material and decal, light information and camera setting;

(2)绘制初始化，对三维场景建立空间加速结构；(2) Rendering initialization, establishing a spatial acceleration structure for the 3D scene;

(3)开始绘制图像，根据用户指定的采样数量或者绘制时间开始执行采样循环；每个采样循环内对所有像素分别进行一个采样的绘制；每个采样循环内的计算过程包括以下子步骤：(3) Start to draw the image, and start to execute the sampling cycle according to the number of samples specified by the user or the drawing time; in each sampling cycle, a sample is drawn for all pixels; the calculation process in each sampling cycle includes the following sub-steps:

(3.1)进行光子的采样：光线从光源发射出来之后，在场景中进行一系列的反射折射；光线每一次的反射折射发生时，在该位置创建一个光子，并将光线当时携带的能量，以及当前光线反射次数的信息记录在该光子中；当所有的光线在场景中遍历完毕，得到一批光子，然后对这些光子建立空间加速结构；(3.1) Sampling photons: After the light is emitted from the light source, it undergoes a series of reflection and refraction in the scene; when each reflection and refraction of the light occurs, a photon is created at the position, and the energy carried by the light at that time, and The information of the number of reflections of the current light is recorded in the photon; when all the light rays have traversed the scene, a batch of photons is obtained, and then a spatial acceleration structure is established for these photons;

(3.2)进行视线的采样：视线的采样数量与光线一致，等于像素的个数；视线发射后同样在场景中进行遍历，每次在物体表面进行反射折射时，在反射折射点周围半径为d的球体范围内收集光子；之后对于每个光子建立一条光路，并对其光通量进行无偏估计，无偏估计的过程如下：(3.2) Sampling the line of sight: the number of samples of the line of sight is consistent with the light, which is equal to the number of pixels; after the line of sight is emitted, it also traverses the scene, and every time the reflection and refraction is performed on the surface of the object, the radius around the reflection and refraction point is d Collect photons within the range of the sphere; then establish a light path for each photon, and perform an unbiased estimation of its luminous flux. The unbiased estimation process is as follows:

(3.2.1)光子映射的误差分析(3.2.1) Error analysis of photon mapping

将光子映射中的光能传递函数修订为严格的光能传递函数公式如下：The radiosity transfer function in photon mapping Revised to strict radiosity transfer function The formula is as follows:

${f f}_{V V M m}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {f f}_{s the s} (({y the y}_{{s the s}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}} (({z z}_{{t t}^{' '}})) &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11})) G G (({z z}_{{t t}^{' '}} &LeftRightArrow; &LeftRightArrow; {z z}_{{t t}^{' '} - - 11})) {f f}_{s the s} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '}}))$

${f f}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {f f}_{s the s} (({y the y}_{{s the s}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11})) G G (({y the y}_{{s the s}^{' '}} &LeftRightArrow; &LeftRightArrow; {z z}_{{t t}^{' '} - - 11})) {f f}_{s the s} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}}))$

将光子映射中的连接概率修订为严格的连接概率公式如下：Linking probabilities in photon maps Revised to strict connection probabilities The formula is as follows:

${p p}_{V V M m}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {πd πd}^{22} {p p}_{x x} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}}))$

${p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {&Integral; &Integral;}_{S S (({y the y}_{{s the s}^{' '}},, d d))} {p p}_{x x} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; z z)) d d z z$

其中，表示一条光子映射的光路，由长度为s′的光线与长度为t′-1的视线连接而成；y和z分别表示光线和视线上的顶点，下标代表了该点在各自子光路中的位置，y_s′代表光线上的最后一个点，也就是光子，而y_s′-1代表光子的上一个点，z_t′和z_t′-1分别代表视线的终点和上一个点，z_t′-2是视线终点上一个点的再前一个点；f_s是物体表面的材质反射系数，表示了物体上具有某种材质的一点，给定光线入射角度和视线观察角度之后的光能反射量，f_s(z_t′-2→z_t′-1→z_t′)表示了在z_t′-1点位置，光能由z_t′-2点而来经过该点后反射到z_t′点的反射系数，f_s(y_s′-1→y_s′(z_t′)→z_t′-1)表达的是在y_s′点位置，光能从y_s′-1点入射之后沿z_t′到z_t′-1的方向出射的反射系数；G代表了两个点之间的几何系数，p_x(z_t′-2→z_t′-1→z)代表了从z_t′-2点打到z_t′-1的光线，在z_t′-1点反射后打到z点的概率；积分范围S(y_s′,d)指的是y_s′点周围半径为d的采集范围；in, Indicates a photon-mapped optical path, which is formed by connecting a ray of length s′ to a line of sight of length t′-1; y and z represent the vertex on the ray and the line of sight respectively, and the subscript represents the point in each sub-optical path , y _s' represents the last point on the ray, that is, the photon, and y _s'-1 represents the previous point of the photon, z _t' and z _t'-1 represent the end point and the previous point of the line of sight respectively, z _t′-2 is the point before the point on the end point of the line of sight; f _s is the material reflection coefficient on the surface of the object, which indicates the light after a given point of light incident angle and line of sight observation angle at a point on the object with a certain material The amount of reflection, f _s (z _t′-2 → z _t′-1 → z _t′ ) indicates that at point z _t′-1 , light energy is reflected from point z _t′-2 after passing through this point The reflection coefficient to the point z _t′ , f _s (y _s′-1 →y _s′ ₍ z _t′ )→z _t′-1 ) expresses the light energy from y _s′- Reflection coefficient after incident at point ₁ along the direction from z _t′ to z _t′-1 ; G represents the geometric coefficient between two points, p _x (z _t′-2 →z _t′-1 →z) Represents the probability of the light hitting z _t′-1 from point z _t′-2 and hitting point z after being reflected at point z _t′-1 ; the integral range S(y _s′ ,d) refers to y _{s The} collection range with a radius of d around the point;

(3.2.2)无偏的连接概率积分倒数的估计：对于每条光线都从光子的上一个点发射出一条试探光线，该光线的分布与实际的光线采样完全一致，该光线采样能够被视线终点的采集范围所接受的概率为在同样的配置下反复生成该试探光线形成一个伯努利试样过程；对于连接概率积分倒数的无偏估计通过伯努利试样中的第一个成功的试探光线的编号N来估计，也就是第一次击中视线终点所在的采集范围的试探光线的编号；(3.2.2) Estimation of the reciprocal of the unbiased connection probability integral: For each ray, a trial ray is emitted from the previous point of the photon, the distribution of the ray is exactly the same as the actual ray sampling, and the ray sampling can be detected by the line of sight The accepted probability of the acquisition range of the end point is Repeatedly generating the trial ray under the same configuration forms a Bernoulli sample process; the unbiased estimation of the reciprocal of the connection probability integral is estimated by the number N of the first successful trial ray in the Bernoulli sample, and also It is the number of the test ray that first hits the acquisition range where the line of sight end point is located;

(4)将视线采集到的每一个光子的光能都累加到图像的相应像素上，在指定的采样循环运行完毕后，输出图像，绘制完毕。(4) Accumulate the light energy of each photon collected by the line of sight to the corresponding pixel of the image, and output the image after the specified sampling cycle runs, and the drawing is completed.

进一步地，所述步骤(3)还包括对步骤(3.2.2)中的伯努利试样中的试探光线进行角度限制的步骤，具体如下：Further, the step (3) also includes the step of limiting the angle of the test light in the Bernoulli sample in the step (3.2.2), specifically as follows:

将采集范围投影到试探光线的起点所在的单位球上，只生成在投影后的采集范围内的试探光线；将投影后的采集范围转化为一个定义在采样的随机数空间中的与坐标轴对齐的包围盒；之后通过限制随机数生成的范围得到该积分范围内的一个采样，角度限制后的伯努利试样的个数从N降为N^b，并且连接概率积分的倒数相应发生改变：Project the acquisition range onto the unit sphere where the starting point of the trial ray is located, and only generate the trial rays within the projected acquisition range; convert the projected acquisition range into a random number space defined in the sampling and aligned with the coordinate axis The bounding box of ; then a sample within the integration range is obtained by limiting the range of random number generation, the number of Bernoulli samples after angle limitation is reduced from N to N ^b , and the reciprocal of the connection probability integral changes accordingly:

$\frac{11}{{p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))} = = \frac{E E. [[{N N}^{b b} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))]]}{{p p}^{b b} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))}$

${p p}^{b b} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {&Integral; &Integral;}_{{Ω Ω}_{b b}} {p p}_{x x} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; z z)) d d z z$

其中积分范围Ω_b从整个半球面缩减为采集范围投影的区域。Wherein the integration range Ω _b is reduced from the entire hemisphere to the area projected by the collection range.

进一步地，通过步骤(3)得到一条光路的无偏估计后，通过多重重要性采样将该估计与其他的光路估计相结合，采用指数权重，权重中不同采样技术的采样概率如下：Further, after obtaining an unbiased estimate of a light path through step (3), this estimate is combined with other light path estimates through multiple importance sampling, using exponential weights, and the sampling probabilities of different sampling techniques in the weights are as follows:

双向光线跟踪的方法生成一条光路的概率为：The method of bidirectional ray tracing generates the probability of a ray path for:

${p p}_{s the s,, t t}^{B B D D. P P T T} ((\overset{&OverBar; &OverBar;}{x x})) = = {p p}^{L L} (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t})) {p p}^{E E.} (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t}))$

其中和分别是光线的生成概率以及视线的生成概率；in and are the generation probability of light and the generation probability of line of sight;

无偏光子映射生成一条光路的概率为：Probability of unpolarized photon mapping to generate a light path for:

${p p}_{{s the s}^{' '},, {t t}^{' '}}^{U u P P G G} ((\overset{&OverBar; &OverBar;}{x x})) = = {p p}^{L L} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {p p}^{E E.} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))$

将连接概率近似为：Approximate the join probability as:

${p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) \approx \approx m m i i n no (({πd πd}^{22} {p p}_{x x} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}})),, 11)) . .$

本发明的有益效果是：本发明在理论上分析了光子映射的误差来源，通过创新性的伯努利试探光线来彻底消除了该误差，并通过试探光线的角度限制将额外的代价降至最低。由于将光子映射变为了无偏的绘制技术，因此可以利用多重重要性采样与其他互补的绘制技术相结合。本发明全面提高了光子映射这一重要绘制技术的效率，不但减少了用户调节参数的负担，并且使用户可以在更短的时间内得到满意的图像。该技术的广泛应用可以提高动漫、电影、广告以及建筑业的渲染绘制效率，从而降低生产成本。The beneficial effects of the present invention are: the present invention theoretically analyzes the error source of photon mapping, completely eliminates the error through the innovative Bernoulli test ray, and minimizes the additional cost by the angle limit of the test ray . Since photon mapping becomes an unbiased rendering technique, multiple importance sampling can be exploited in combination with other complementary rendering techniques. The invention comprehensively improves the efficiency of photon mapping, which is an important drawing technology, not only reduces the user's burden of adjusting parameters, but also enables users to obtain satisfactory images in a shorter time. The wide application of this technology can improve the rendering and drawing efficiency of animation, film, advertisement and architectural industries, thereby reducing production costs.

附图说明Description of drawings

图1是本发明中关于伯努利试探光线以及角度限制的示意图说明，(a)图表示了不含角度限制的伯努利试探光线采样，(b)图表示了加上角度限制之后的伯努利试探光线采样。Fig. 1 is the schematic diagram explanation about Bernoulli's trial ray and angle limitation among the present invention, (a) figure has represented the Bernoulli's trial ray sampling without angle restriction, (b) figure has represented Bernoulli's trial ray sampling after adding angle restriction Play around with ray sampling.

图2是宫殿场景用顶点合并算法绘制10分钟的结果。Figure 2 is the result of drawing the palace scene for 10 minutes with the vertex merging algorithm.

图3是宫殿场景用本发明绘制10分钟的结果。Fig. 3 is the result of rendering the palace scene for 10 minutes with the present invention.

图4是宫殿场景用双向光线跟踪绘制24小时的结果，作为正确结果参考。Figure 4 is the result of drawing the palace scene for 24 hours with two-way ray tracing, as a reference for the correct result.

图5是楼梯场景用双向光线跟踪绘制1小时的结果。Figure 5 is the result of rendering the staircase scene with bi-directional ray tracing for 1 hour.

图6是楼梯场景用本发明中的光子映射部分绘制1小时的结果。Fig. 6 is the result of one hour rendering of the staircase scene using the photon mapping part of the present invention.

图7是楼梯场景用本发明中光子映射与双向光线跟踪相结合的方法绘制1小时的结果。Fig. 7 is the result of rendering the stair scene for 1 hour with the method of combining photon mapping and bidirectional ray tracing in the present invention.

图2-图7的效果图中，为了更好的比较绘制效果的差异，均在原图下方加入了局部放大图。被放大的局部区域用方框标示在了原图中。In the effect diagrams of Fig. 2-Fig. 7, in order to better compare the differences in rendering effects, partial enlarged images are added below the original images. The magnified local area is marked with a box in the original image.

具体实施方式Detailed ways

本发明的核心技术是对光子映射中的每个光子所对应的光路所携带光能的无偏估计，为此，需要首先将视线采集到的光子连接成为一条完整的光路，然后对该光路的贡献进行无偏的估计，最后通过多重重要性采样将该种采样技术与其他的采样技术相结合。The core technology of the present invention is the unbiased estimation of the light energy carried by the optical path corresponding to each photon in the photon mapping. The contribution is estimated unbiasedly, and finally this sampling technique is combined with other sampling techniques through multiple importance sampling.

1.输入三维场景，所述三维场景包括物体的几何信息、材质以及贴图、灯光信息和相机设置。在本方法中采用了Mitsuba格式的场景文件(参考https://www.mitsuba-renderer.org/)，通过读取该场景文件并进行解析，来获取场景相关的所有信息。1. Input a 3D scene, the 3D scene includes the geometric information, material and texture of the object, lighting information and camera settings. In this method, a scene file in Mitsuba format is used (refer to https://www.mitsuba-renderer.org/ ), and all information related to the scene can be obtained by reading the scene file and parsing it.

2.绘制初始化，对场景建立空间加速结构，我们选用了SBVH的加速结构(参照SpatialSplitsinBoundingVolumeHierarchies，MartinStichetal.,HPG'09ProceedingsoftheConferenceonHighPerformanceGraphics2009,Pages7-13)。该加速结构用于为几何信息在空间中建立索引，加速之后的光线跟踪计算。2. Draw initialization, and establish a spatial acceleration structure for the scene. We choose the acceleration structure of SBVH (refer to Spatial Splits in Bounding Volume Hierarchies, Martin Stiche et al., HPG'09 Proceeding of the Conference on High Performance Graphics 2009, Pages7-13). This acceleration structure is used to index geometric information in space to accelerate subsequent ray tracing calculations.

3.开始绘制图像，根据用户指定的采样数量或者绘制时间开始执行采样循环。每个采样循环内对所有像素分别进行一个采样的绘制。每个采样循环内部的计算过程包括以下子步骤：3. Start to draw the image, and start to execute the sampling cycle according to the number of samples specified by the user or the drawing time. All pixels are drawn for one sample in each sampling cycle. The calculation process inside each sampling cycle consists of the following sub-steps:

3.1进行光子的采样。在每个循环中采样的光线数量等于像素的个数，每个循环中采样的光线数量与视线数量相等。光线从光源发射出来之后，会在场景中进行一系列的反射折射。光线每一次的反射折射发生时，在该位置创建一个光子，并将光线当时携带的能量，以及当前光线反射次数的信息记录在该光子中。当所有的光线在场景中遍历完毕，得到一批光子，然后会对这些光子建立空间加速结构，以加速之后采集光子的过程。3.1 Sampling of photons. The number of rays sampled in each cycle is equal to the number of pixels, and the number of rays sampled in each cycle is equal to the number of views. After the light is emitted from the light source, it will undergo a series of reflection and refraction in the scene. Every time the reflection and refraction of light occurs, a photon is created at the position, and the energy carried by the light at that time and the information of the number of reflections of the current light are recorded in the photon. When all the light rays have traversed the scene, a batch of photons will be obtained, and then a spatial acceleration structure will be established for these photons to accelerate the process of collecting photons later.

3.2进行视线的采样，视线的采样数量与光线一致，等于像素的个数。视线发射后同样在场景中进行遍历，每次在物体表面进行反射折射，就会在反射折射点附近收集光子。之后对于每个光子建立一条光路，并对其光通量进行无偏的估计，。无偏光路的具体估计方法如下：3.2 Sampling the line of sight, the sampling number of the line of sight is consistent with the light, which is equal to the number of pixels. After the line of sight is launched, it also traverses the scene. Every time it performs reflection and refraction on the surface of the object, it will collect photons near the reflection and refraction point. Afterwards, a light path is established for each photon, and its luminous flux is estimated unbiasedly. The specific estimation method of the unpolarized optical path is as follows:

3.2.1光子映射的误差分析3.2.1 Error analysis of photon mapping

首先我们回顾一下光路积分以及顶点合并，这是本发明无偏估计的理论基础。First, let's review the light path integration and vertex merging, which are the theoretical basis of the unbiased estimation of the present invention.

蒙特卡洛光路积分是对整个光路空间的积分，光路的一端是照相机，另一端是光源，中间所有的光路形成了光能所有可能的传播方式，因此对于这所有的光路的积分也就是对于当前的场景配置下完全准确的绘制结果。由于该积分不可能通过解析的方式得到准确解，因此一般通过蒙特卡洛方法进行无偏估计：Monte Carlo optical path integration is the integration of the entire optical path space. One end of the optical path is the camera, the other end is the light source, and all the optical paths in the middle form all possible propagation modes of light energy. Therefore, the integral of all the optical paths is also for the current Completely accurate rendering results under the scene configuration. Since it is impossible to obtain an accurate solution of the integral analytically, it is generally estimated unbiasedly by the Monte Carlo method:

$I I = = &Integral; &Integral; f f ((\overset{&OverBar; &OverBar;}{x x})) d d \overset{&OverBar; &OverBar;}{x x} \approx \approx Σ Σ \frac{f f ((\overset{&OverBar; &OverBar;}{x x}))}{p p ((\overset{&OverBar; &OverBar;}{x x}))}$

其中x是一条特定的光路，f(x)是该光路携带的光能，p(x)是该光路生成的概率，I是准确的图像绘制结果。传统的双向光线跟踪是通过分别采样一条视线z和一条光线y，然后将两条光线连接形成一条完整的光路，其贡献函数如下：Where x is a specific light path, f(x) is the light energy carried by the light path, p(x) is the probability generated by the light path, and I is the accurate image rendering result. The traditional two-way ray tracing is to sample a line of sight z and a ray y respectively, and then connect the two rays to form a complete light path. The contribution function as follows:

${C C}^{* *} (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t})) \frac{f f (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t}))}{p p (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t}))} = = {α α}^{L L} (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t})) {f f}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t})) {α α}^{E E.} (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t}))$

其中和分别是和光线以及视线相关的项，由于这些项和我们的方法无关，这里我们主要关注中间的项，也就是连接部分的贡献：in and They are items related to light and line of sight. Since these items have nothing to do with our method, here we mainly focus on the middle term, which is the contribution of the connected part:

${f f}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t})) = = {f f}_{s the s} (({y the y}_{s the s - - 11} &RightArrow; &Right Arrow; {y the y}_{s the s} &RightArrow; &Right Arrow; {z z}_{t t})) G G (({y the y}_{s the s} &LeftRightArrow; &LeftRightArrow; {z z}_{t t})) {f f}_{s the s} (({z z}_{t t - - 11} &RightArrow; &Right Arrow; {z z}_{t t} &RightArrow; &Right Arrow; {y the y}_{s the s}))$

其中，表示一条光子映射的光路，由长度为s′的光线与长度为t′-1的视线连接而成；y和z分别表示光线和视线上的顶点，下标代表了该点在各自子光路中的位置，如y_s′代表光线上的最后一个点，也就是光子，而y_s′-1代表光子的上一个点，以此类推，z_t′和z_t′-1分别代表视线的终点和上一个点；f_s是物体表面的材质反射系数，表示了物体上具有某种材质的一点，给定光线入射角度和视线观察角度之后的光能反射量，如f_s(z_t′-2→z_t′-1→z_t′)表示了在z_t′-1点位置，光能由z_t′-2点而来经过该点后反射到z_t′点的反射系数，一个特例是f_s(y_s′-1→y_s′(z_t′)→z_t′-1)，该函数表达的是在y_s′点位置，光能从y_s′-1点入射之后沿z_t′到z_t′-1的方向出射的反射系数；G代表了两个点之间的几何系数。由于双向光线跟踪是以1概率进行连接，因此对于双向光线跟踪的连接部分而言，这里的贡献就只有光能函数所包含的项，而没有与概率相关的项。in, Indicates a photon-mapped optical path, which is formed by connecting a ray of length s′ to a line of sight of length t′-1; y and z represent the vertex on the ray and the line of sight respectively, and the subscript represents the point in each sub-optical path For example, y _s′ represents the last point on the ray, that is, the photon, and y _s′-1 represents the previous point of the photon, and so on, z _t′ and z _t′-1 represent the end points of the line of sight respectively and the previous point; f _s is the material reflection coefficient on the surface of the object, which indicates the amount of light energy reflected by a point with a certain material on the object, given the incident angle of light and the viewing angle of the line of sight, such as f _s (z _{t′- 2} →z _t′-1 →z _t′ ) represents the reflection coefficient of the light energy from point z _t′-2 at point z t′ _-1 and then reflected to point z _t′ after passing through this point, a special case is f _s (y _s′-1 →y _s′ (z _t′ )→z t′ _- ₁ ), this function expresses that at the position of point y _s′ , light energy along the The reflection coefficient of the outgoing direction from z _t′ to z _t′-1 ; G represents the geometric coefficient between two points. Since the two-way ray tracing is connected with a probability of 1, for the connection part of the two-way ray tracing, the contribution here is only the item contained in the light energy function, and there is no item related to the probability.

在光子映射中，视线终点处的入射光能是通过密度估计来进行计算的，密度估计则是通过在一定距离范围内搜索所有的光子进行的。搜索光子的距离一般是一个设定值d。根据顶点合并对于光子映射的定义，将视线与每一个采集到的光子配对起来，在视线和光线进行连接的终点处，将光线的最后一个点移除，从而形成一条光路，并且其贡献函数如下：In photon mapping, the incident light energy at the end of the line of sight is calculated by density estimation, which is performed by searching all photons within a certain distance range. The distance to search for photons is generally a set value d. According to the definition of photon mapping by vertex merging, the line of sight is paired with each collected photon, and at the end point where the line of sight and light are connected, the last point of the light is removed to form a light path, and its contribution function as follows:

${C C}_{V V M m}^{* *} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {α α}^{L L} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {α α}_{V V M m}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {α α}^{E E.} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))$

${α α}_{V V M m}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = \frac{{f f}_{V V M m}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))}{{p p}_{V V M m}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))}$

其中上一个公式中的连接项贡献变成了其中包含了光能传递函数以及连接概率 where the connected term in the previous formula contributes became which contains the light energy transfer function and the connection probability

顶点合并将光子采集描述成一个虚拟的俄罗斯轮盘(RussianRoulette)的过程。光子采集的概率对应到从光子的上一个点向视线终点的采集范围内进行独立的采样的接受概率。所谓接受，也就是光子距离视线的终点的距离小于等于d(用户设定的采集半径)。在此基础上，可以将光子映射描述成另一个双向的光路采样方法，并找出其估计的偏差所在。Vertex merging describes photon collection as a virtual Russian Roulette process. Probability of photon collection Corresponds to the acceptance probability of independent sampling in the collection range from the last point of the photon to the end of the line of sight. The so-called acceptance means that the distance between the photon and the end point of the line of sight is less than or equal to d (the collection radius set by the user). On this basis, photon mapping can be described as another two-way optical path sampling method, and the deviation of its estimation can be found.

为了得到无偏的光子映射，将与进行比较，可以看到首先要保证光能传递函数与f^c相等，然后需要是光子被采集到的真正的概率。然而在传统的光子映射或者顶点合并中，只是通过密度估计进行计算的，因此这两个要求都不能满足。In order to get unbiased photon mapping, the and For comparison, it can be seen that the light energy transfer function must first be guaranteed is equal to f ^c , then requires is the true probability that a photon is collected. However, in traditional photon mapping or vertex merging, is only computed through density estimation, so neither requirement is met.

具体来说，在顶点合并中，光能传递函数可以表示为：Specifically, in vertex merging, the radiosity transfer function can be expressed as:

然而严格的光能传递函数应该是：However, the strict light energy transfer function should be:

相比之下包括BSDF函数以及几何项在内的部分都存在着近似。In contrast, there are approximations in parts including BSDF functions and geometric terms.

而对于光子采集的概率，真正的采集概率应该是一个积分：And for the probability of photon collection, the true collection probability Should be an integral:

其中p_x(z_t′-2→z_t′-1→z)代表了从z_t′-2点打到z_t′-1的光线，在z_t′-1点反射后打到z点的概率；积分范围S(y_s′,d)指的是y_s′点周围半径为d的采集范围。Among them, p _x (z _t′-2 → z _t′-1 → z) represents the light from point z _t′-2 to z _t′-1 , which hits point z after being reflected at point z _t′-1 The probability of ; the integration range S(y _s′ ,d) refers to the acquisition range around the point y _s′ with a radius of d.

然而顶点合并中近似了这个概率，也就是：However, this probability is approximated in vertex merging, that is:

可以看到一方面顶点合并用了一个采样的蒙特卡洛方法来估计这个积分，同时将积分区域近似成一个半径为d的圆盘。It can be seen that on the one hand, vertex merging uses a sampling Monte Carlo method to estimate the integral, and at the same time approximates the integral area as a disk with a radius of d.

如果我们同时对于光能传递函数以及概率的计算采用严格的公式，则我们就能得到一个无偏的光子映射算法，该算法的贡献函数为：If we use strict formulas for both the light energy transfer function and the calculation of the probability, we can get an unbiased photon mapping algorithm whose contribution function for:

${C C}_{R R R R}^{* *} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {α α}^{L L} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {α α}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {α α}^{E E.} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))$

${α α}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = \frac{{f f}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))}{{p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))}$

其中的阿尔法函数描述了一个假象的俄罗斯轮盘赌过程，也就是在光子采集的时候，从光线的上一个点发出的光子落在视线的采集半径内的贡献的期望。z就是一个独立于视线的采样，其采样过程和普通的光线采样相同，是从光子的上一个点发射出来的。如果这个采样落在视线终点的采集半径之内，则接受这个采样，否则则拒绝这个采样，这样一来该采样落在采集范围内的概率的期望就是注意到该过程实际上就是视线和光线两段光路可以通过光子相连的过程，因此这个概率也就是一条通过光子采集生成的光路所产生的概率。The alpha function describes an imaginary Russian roulette process, that is, the expectation of the contribution of the photon emitted from the previous point of the ray falling within the collection radius of the line of sight when the photon is collected. z is a sampling independent of the line of sight. Its sampling process is the same as that of ordinary light sampling, and it is emitted from the previous point of the photon. If the sample falls within the collection radius of the end point of the line of sight, the sample is accepted, otherwise the sample is rejected, so that the expectation of the probability that the sample falls within the collection range is Note that this process is actually a process in which the two optical paths of line of sight and light can be connected through photons, so this probability is also the probability generated by an optical path generated by photon collection.

注意到该连接过程与顶点合并算法不同，因为严格遵守了整个光路无偏估计的条件。为了实现这样的算法，我们需要将以上数学表达中的每一项都严格进行计算从而不影响整体的无偏性质。其中最主要的挑战来自于连接概率的计算。该概率中不止包含了光线前一个点上的材质反射函数，还包括了两点之间的几何信息，包括相对的距离，两个点上各自出射的角度，以及可见性，从而使得解析地计算这个连接概率是不可能的。Note that this connection process is different from the vertex merging algorithm because the condition of unbiased estimation of the entire optical path is strictly observed. In order to implement such an algorithm, we need to strictly calculate each of the above mathematical expressions so as not to affect the unbiased property of the whole. One of the most important challenges comes from the connection probability calculation. This probability not only includes the material reflection function at the previous point of the ray, but also includes the geometric information between the two points, including the relative distance, the angles at which the two points are emitted, and the visibility, so that it can be calculated analytically This connection probability is impossible.

因此我们引入一个独立的蒙特卡洛过程来估计这个连接概率。然而，直接对这个概率的积分函数应用蒙特卡洛估计会在最终的结果中引入偏差，因为这个概率积分在最终的估计中是出现在分母上的。Jensen不等式说明了这一点：We therefore introduce an independent Monte Carlo process to estimate this connection probability. However, applying Monte Carlo estimation directly to the integral function of this probability introduces bias in the final result, since this probability integral appears in the denominator in the final estimate. Jensen's inequality illustrates this:

$\frac{11}{E E. [[{p p}^{c c}]]} &GreaterEqual; &Greater Equal; E E. [[\frac{11}{{p p}^{c c}}]]$

该不等式只有在估计的噪声为零时可以实现相等，然而在我们的问题中噪声为零是不可能实现的。因此我们需要一个方法来无偏地估计概率积分的倒数。This inequality can only be achieved when the estimated noise is zero, which is not possible in our problem. So we need a way to unbiasedly estimate the reciprocal of the probability integral.

3.2.2无偏的连接概率积分倒数的估计3.2.2 Estimation of the reciprocal of the unbiased connection probability integral

如附图1所示，我们按照之前显式的光子采集过程，对于每条光线都从光子的上一个点发射出一条试探光线，该光线的分布与实际的光线采样完全一致。根据之前的描述，该光线采样能够被视线终点的采集范围所接受的概率为在同样的配置下反复生成该试探光线就形成了一个伯努利试样过程。假设第一次被接受的试探光线是第N条，则N符合如下的几何分布：As shown in Figure 1, according to the previous explicit photon collection process, we emit a trial ray from the previous point of the photon for each ray, and the distribution of this ray is completely consistent with the actual ray sampling. According to the previous description, the probability that the ray sampling can be accepted by the collection range of the line of sight end point is Generating the test light repeatedly under the same configuration forms a Bernoulli sample process. Assuming that the first tentative ray accepted is the Nth one, then N conforms to the following geometric distribution:

$Pr PR ((N N (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = i i)) = = {p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {((11 - - {p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))))}^{i i - - 11}$

然后N的期望恰好是我们所需要的概率积分的倒数：Then the expectation of N is exactly the inverse of the probability integral we need:

$\begin{matrix} E E. [[N N (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))]] = = {Σ Σ}_{i i = = 11}^{+ + ∞ ∞} {ip ip}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {((11 - - {p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))))}^{i i - - 11} \\ \frac{11}{{p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))} \end{matrix}$

因此，对于概率积分的倒数的无偏估计可以通过伯努利试样中的第一个成功的试探光线的编号来估计，也就是第一次击中视线终点所在的采集范围的试探光线的编号。Therefore, an unbiased estimate of the reciprocal of the probability integral can be estimated from the number of the first successful trial ray in the Bernoulli sample, that is, the number of the first trial ray that hits the acquisition range where the end of the line of sight is located .

3.2.3试探光线的角度限制3.2.3 Test the Angle Limitation of Rays

之前我们提到的概率积分倒数的无偏估计的算法存在一个显著的问题，就是伯努利试样可能要经过任意长的时间才能得到第一个成功的尝试，然而这个代价是无法承受的。为了使我们的算法实际可用，我们必须尽可能避免明显不可能的连接。There is a significant problem in the unbiased estimation algorithm of the reciprocal of the probability integral that we mentioned before, that is, the Bernoulli sample may take an arbitrarily long time to get the first successful attempt, but this cost is unbearable. In order for our algorithm to be practical, we must avoid obviously impossible joins as much as possible.

为此，我们利用了采集范围对于试探光线起点的相对大小。具体来说，我们将采集范围(也就是中心在视线终点的一个球体)投影到试探光线的起点所在的单位球上，然后只生成在投影后的采集范围内的试探光线，如附图1(b)所示。我们将投影后的采集范围转化为一个定义在采样的随机数空间中的AABB(坐标轴对齐的包围盒)。之后我们就可以直接通过限制随机数生成的范围来得到该积分范围内的一个采样。To do this, we take advantage of the relative size of the acquisition range to the origin of the trial ray. Specifically, we project the collection range (that is, a sphere whose center is at the end of the line of sight) onto the unit sphere where the starting point of the trial ray is located, and then only generate trial rays within the projected collection range, as shown in Figure 1 ( b) as shown. We convert the projected acquisition range into an AABB (Axis Aligned Bounding Box) defined in the sampled random number space. Then we can directly get a sample within the integral range by limiting the range of random number generation.

由此，我们将伯努利试样的个数有效降低了，降低的效果与试探光线起点与视线终点的距离有关，此时伯努利试样的个数满足以下的关系：Therefore, we have effectively reduced the number of Bernoulli samples, and the reduction effect is related to the distance between the starting point of the test light and the end point of the line of sight. At this time, the number of Bernoulli samples Satisfy the following relationship:

$\frac{E E. [[{N N}^{b b} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))]]}{E E. [[N N (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))]]} \approx \approx \frac{{πd πd}^{22}}{| | | | {y the y}_{{s the s}^{' '}} - - {z z}_{{t t}^{' '} - - 11} | | {| |}^{22}}$

基于我们的实验，该角度限制可以显著降低伯努利采样的数量，将其复杂度从O(n^2)降到O(n)，n是视线或者光线的数量。Based on our experiments, this angle limit can significantly reduce the number of Bernoulli samples, reducing its complexity from O(n^2) to O(n), where n is the number of sight lines or rays.

加上角度限制之后，概率积分的倒数也需要相应发生改变：After adding the angle limit, the reciprocal of the probability integral also needs to be changed accordingly:

其中积分范围Ω_b从整个半球面缩减为采集范围投影的区域。值得注意的是，这里的由于不存在可见性的问题，因此我们必须解析地计算这个概率积分。Wherein the integration range Ω _b is reduced from the entire hemisphere to the area projected by the collection range. It is worth noting that here Since there is no visibility problem, we must compute this probability integral analytically.

在这一部分中，我们将光子映射的绘制问题放到了路径积分的框架下进行分析，并与现有的无偏以及有偏方法进行比较分析，从而得到了无偏光子映射的方法。更进一步的，我们首创性地提出了试探光线的方法来计算其中一个概率积分的倒数，成为了我们这个方法得以实现的决定性因素。然而简单地实现该试探性光线的方法会带来性能上的巨大损失，我们又进一步地提出了对试探光线进行角度限制的方法，将该部分的额外代价降至最低。综合以上的步骤，我们最终能够无偏地估计视线所收集到的光子所包含的正确的能量。In this part, we analyze the drawing problem of photon mapping under the framework of path integral, and compare it with the existing unbiased and biased methods, so as to obtain the method of unbiased photon mapping. Furthermore, we pioneered the method of testing light to calculate the reciprocal of one of the probability integrals, which became the decisive factor for the realization of our method. However, the method of simply implementing the tentative ray will bring about a huge loss in performance. We further proposed a method to limit the angle of the tentative ray to minimize the additional cost of this part. Combining the above steps, we can finally unbiasedly estimate the correct energy contained in the photons collected by the line of sight.

3.3通过多重重要性采样与其他互补的采样技术相结合缩小光路估计噪声3.3 Reduce the optical path estimation noise by combining multiple importance sampling with other complementary sampling techniques

得到一条光路的无偏估计后，通过多重重要性采样将该估计与其他的光路估计相结合。Once an unbiased estimate of a lightpath is obtained, this estimate is combined with other lightpath estimates through multiple importance sampling.

加入多重重要性采样之后的光路积分为：The optical path integral after adding multiple importance sampling is:

$I I \approx \approx {Σ Σ}_{i i = = 11}^{m m} \frac{11}{{n no}_{i i}} {Σ Σ}_{j j = = 11}^{{n no}_{i i}} {w w}_{i i} (({\overset{&OverBar; &OverBar;}{x x}}_{i i,, j j})) \frac{f f (({\overset{&OverBar; &OverBar;}{x x}}_{i i,, j j}))}{{p p}_{i i} (({\overset{&OverBar; &OverBar;}{x x}}_{i i,, j j}))}$

其中m表示采样技术的个数，n_i代表对于第i个采样技术所采样的个数。Among them, m represents the number of sampling techniques, and n _i represents the number of sampling for the i-th sampling technique.

此外，我们采用被广泛接受的指数权重 Furthermore, we employ the widely accepted index weights

${w w}_{i i} ((\overset{&OverBar; &OverBar;}{x x})) = = \frac{{(({n no}_{i i} {p p}_{i i} ((\overset{&OverBar; &OverBar;}{x x}))))}^{β β}}{{Σ Σ}_{j j = = 11}^{m m} (({n no}_{j j} {p p}_{j j} ((\overset{&OverBar; &OverBar;}{x x})))) β β}$

我们采用指数β为2。由于该权重满足之前提到的无偏估计的条件，因此多重重要性采样的主要挑战在于计算同一条光路对于不同采样技术的相对概率密度。We use an exponent β of 2. Since this weight satisfies the condition of unbiased estimation mentioned earlier, the main challenge of multiple importance sampling is to calculate the relative probability density of the same light path for different sampling techniques.

给定一条长度为k的光路，同时有k+1个双向光线跟踪的采样技术以及k-1个无偏光子映射的采样技术可以生成完全相同的光路，其区别主要在于连接位置的不同。双向光线跟踪的采样技术较多是由于它可以处理长度为零的子光路。双向光线跟踪的方法生成一条光路的概率为：Given an optical path of length k, there are k+1 sampling techniques of bidirectional ray tracing and k-1 sampling techniques of unpolarized photon mapping can generate exactly the same optical path, the difference mainly lies in the different connection positions. The sampling technique of bidirectional ray tracing is more because it can handle sub-raypaths of zero length. The method of bidirectional ray tracing generates the probability of a ray path for:

其中和分别是光线的生成概率以及视线的生成概率，由于双向光线跟踪在连接处是以1概率连接，因此最终光路的生成概率就是两段概率的乘积。in and They are the generation probability of light and the generation probability of line of sight. Since two-way ray tracing is connected with a probability of 1 at the connection, the final generation probability of the light path is the product of the two probabilities.

另一方面，无偏光子映射生成一条光路的概率为：On the other hand, unpolarized photon mapping generates the probability of a light path for:

相比于双向光线跟踪多了一项p^c，也就是连接概率。但是这个连接概率无法解析计算，而我们之前提到的对于这个概率的无偏估计则会破坏多重重要性采样的一致性。因此我们在这里对于连接概率采取了近似。注意到虽然我们对于多重重要性采样的权重函数引入了近似，但只要该权重函数满足之前提到的无偏估计的条件，则不会破坏最终估计的无偏性质。我们将该连接概率近似为：Compared with two-way ray tracing, there is one more item p ^c , which is the connection probability. But this connection probability cannot be calculated analytically, and the unbiased estimation of this probability we mentioned before will destroy the consistency of multiple importance sampling. We therefore make approximations here for the join probabilities. Note that although we introduce an approximation to the weight function of multiple importance sampling, as long as the weight function satisfies the aforementioned unbiased estimation conditions, it will not destroy the unbiased nature of the final estimate. We approximate this connection probability as:

${p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) \approx \approx m m i i n no (({πd πd}^{22} {p p}_{x x} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}})),, 11))$

这个最小值截断的原因非常简单，因为p^c是概率密度函数的积分，也就是一个概率，而概率是不可能大于1的，因此我们在这里加一个截断可以提高近似的质量。如果不加这个截断，在实际计算中会出现远大于1的概率，从而在权重函数中获得过大的权重，从而引入不必要的噪声。The reason for this minimum truncation is very simple, because p ^c is the integral of the probability density function, that is, a probability, and the probability cannot be greater than 1, so adding a truncation here can improve the quality of the approximation. If this truncation is not added, there will be a probability much greater than 1 in the actual calculation, so that an excessive weight will be obtained in the weight function, thereby introducing unnecessary noise.

4.将视线采集到的每一个光子的光能都累加到图像的相应像素上，在指定的采样循环运行完毕后，输出图像，绘制完毕。4. Accumulate the light energy of each photon collected by the line of sight to the corresponding pixel of the image. After the specified sampling cycle is completed, the image is output and the drawing is completed.

实施实例Implementation example

发明人在一台配备Inteli73.40GHz的四核中央处理器，以及16GB内存的计算机上实现了上文所描述的算法，并实现了与我们算法紧密相关的其他绘制方法，生成了附录中的渲染结果。实验结果表明相比于传统的有偏的光子映射，我们的方法可以有效的消除图像错误，同时由于不需要缩小采集半径，也可以更加充分地发挥光子映射的优势，使图像的噪声更小。将我们的方法与双向光线跟踪相结合后，我们的方法与顶点合并(VCM)方法相比可以在避免图像错误的情况下更加充分地利用光子映射采样技术的优势。在将我们的方法与马尔科夫链蒙特卡洛方法(MCMC)结合后，我们的方法可以显著帮助原MCMC方法提高在多次反射的光路以及金属材质反射上的效率。The inventor implemented the algorithm described above on a computer equipped with an Inteli7 3.40GHz quad-core CPU and 16GB of memory, and implemented other rendering methods closely related to our algorithm, and generated the rendering in the appendix result. The experimental results show that compared with the traditional biased photon mapping, our method can effectively eliminate image errors, and at the same time, because it does not need to reduce the collection radius, it can also take full advantage of the advantages of photon mapping, making the image noise smaller. Combining our method with bidirectional ray tracing, our method can more fully exploit the advantages of photon-mapping sampling techniques while avoiding image errors compared to vertex merging (VCM) methods. After combining our method with the Markov chain Monte Carlo method (MCMC), our method can significantly help the original MCMC method to improve the efficiency of the multiple reflection optical path and the reflection of metal materials.

图2、图3、图4展示了我们的方法与顶点合并方法绘制结果的比较。图2的顶点合并算法由于光子映射有偏的问题，严重模糊了模型表面的细节。图3中我们的方法显著修正了这一问题，同时噪声水平也和图2相当。图4中是无偏方法通过极长时间的绘制完全收敛后的结果，可见我们的结果与该结果是完全一致的。Figure 2, Figure 3, and Figure 4 show the comparison of our method and the rendering results of the vertex pooling method. The vertex merging algorithm in Figure 2 seriously blurs the details of the model surface due to the problem of biased photon mapping. Our method in Fig. 3 corrects this problem significantly, and the noise level is also comparable to Fig. 2. Figure 4 is the result of the unbiased method after a very long time of drawing to complete convergence, it can be seen that our results are completely consistent with this result.

图5、图6、图7展示了我们的无偏光子映射方法与双向光线跟踪算法相结合的效果。图5中的双向光线跟踪算法在建筑内部噪声水平很高，然而在光源直接照射的区域效果不错。图6中我们的方法在室内的噪声水平显著降低，但是在光源直接照射的区域绘制效果却有所下降。因此通过将这两种互补的采样技术相结合，可以得到更加强健的绘制技术，结果如图7所示。结合后的方法继承了两种方法各自的有点，并进行了有机的结合，得到超过任一单一技术的绘制效果。Figures 5, 6, and 7 demonstrate the effect of our unbiased photon mapping method combined with a bidirectional ray tracing algorithm. The bi-directional ray tracing algorithm in Figure 5 has a high noise level inside buildings, but works well in areas directly illuminated by light sources. In Figure 6 our method achieves a significant reduction in the noise level indoors, but suffers from poor rendering performance in areas directly illuminated by the light source. Therefore, by combining these two complementary sampling techniques, a more robust rendering technique can be obtained, and the result is shown in Figure 7. The combined method inherits the respective points of the two methods, and organically combines them to obtain a drawing effect that exceeds that of any single technique.

Claims

1. an unbiased photon mapping drawing method, is characterized in that, comprises the following steps:

(1) input three-dimensional scene file, and three-dimensional scene file is analyzed; Described three-dimensional scene file comprises the geometric information of object, material and decal, light information and camera setting;

(2) Rendering initialization, establishing a spatial acceleration structure for the 3D scene;

(3) Start to draw the image, and start to execute the sampling cycle according to the number of samples specified by the user or the drawing time; in each sampling cycle, a sample is drawn for all pixels; the calculation process in each sampling cycle includes the following sub-steps:

(3.1) Sampling photons: After the light is emitted from the light source, it undergoes a series of reflection and refraction in the scene; when each reflection and refraction of the light occurs, a photon is created at the position, and the energy carried by the light at that time, and The information of the number of reflections of the current light is recorded in the photon; when all the light rays have traversed the scene, a batch of photons is obtained, and then a spatial acceleration structure is established for these photons;

(3.2) Sampling the line of sight: the number of samples of the line of sight is consistent with the light, which is equal to the number of pixels; after the line of sight is emitted, it also traverses the scene, and every time the reflection and refraction is performed on the surface of the object, the radius around the reflection and refraction point is d Collect photons within the range of the sphere; then establish a light path for each photon, and perform an unbiased estimation of its luminous flux. The unbiased estimation process is as follows:

(3.2.1) Error analysis of photon mapping

The radiosity transfer function in photon mapping Revised to strict radiosity transfer function The formula is as follows:

{f f}_{V V M m}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {f f}_{s the s} (({y the y}_{{s the s}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}} (({z z}_{{t t}^{' '}})) &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11})) G G (({z z}_{{t t}^{' '}} &LeftRightArrow; &LeftRightArrow; {z z}_{{t t}^{' '} - - 11})) {f f}_{s the s} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '}}))

{f f}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {f f}_{s the s} (({y the y}_{{s the s}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11})) G G (({y the y}_{{s the s}^{' '}} &LeftRightArrow; &LeftRightArrow; {z z}_{{t t}^{' '} - - 11})) {f f}_{s the s} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}}))

Linking probabilities in photon maps Revised to strict connection probabilities The formula is as follows:

{p p}_{V V M m}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {πd πd}^{22} {p p}_{x x} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}}))

{p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {&Integral; &Integral;}_{S S (({y the y}_{{s the s}^{' '}},, d d))} {p p}_{x x} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; z z)) d d z z

in, Indicates a photon-mapped optical path, which is formed by connecting a ray of length s′ to a line of sight of length t′-1; y and z represent the vertex on the ray and the line of sight respectively, and the subscript represents the point in each sub-optical path , y _s' represents the last point on the ray, that is, the photon, and y _s'-1 represents the previous point of the photon, z _t' and z _t'-1 represent the end point and the previous point of the line of sight respectively, z _t′-2 is the point before the point on the end point of the line of sight; f _s is the material reflection coefficient on the surface of the object, which indicates the light after a given point of light incident angle and line of sight observation angle at a point on the object with a certain material The amount of reflection energy, f _s (z _t′-2 →z _t′-1 →z _t′ ) indicates that at the position of point z _t′-1 , the light energy is reflected from point z _t′-2 after passing through this point The reflection coefficient to point z _t′ , f _s (y _s′-1 →y _s′ (z _t′ )→z _t′-1 ) expresses the light energy at point y _s′ from y _s′- Reflection coefficient after incident at point ₁ along the direction from z _t′ to z _t′-1 ; G represents the geometric coefficient between two points, p _x (z _t′-2 →z _t′-1 →z) Represents the probability of the light hitting z _t′-1 from point z _t′-2 and hitting point z after being reflected at point z _t′-1 ; the integral range S(y _s′ ,d) refers to y _{s The} collection range with a radius of d around the point;

(3.2.2) Estimation of the reciprocal of the unbiased connection probability integral: For each ray, a trial ray is emitted from the previous point of the photon, the distribution of the ray is exactly the same as the actual ray sampling, and the ray sampling can be detected by the line of sight The accepted probability of the acquisition range of the end point is Repeatedly generating the trial ray under the same configuration forms a Bernoulli sample process; the unbiased estimation of the reciprocal of the connection probability integral is estimated by the number N of the first successful trial ray in the Bernoulli sample, and also It is the number of the test ray that first hits the acquisition range where the line of sight end point is located;

(4) Accumulate the light energy of each photon collected by the line of sight to the corresponding pixel of the image, and output the image after the specified sampling cycle runs, and the drawing is completed.

2. a kind of unbiased photon mapping drawing method according to claim 1, it is characterized in that, described step (3) also comprises to the test ray in the Bernoulli sample in step (3.2.2) carries out angle The restricted steps are as follows:

Project the acquisition range onto the unit sphere where the starting point of the trial ray is located, and only generate the trial rays within the projected acquisition range; convert the projected acquisition range into a random number space defined in the sampling and aligned with the coordinate axis The bounding box of ; then a sample within the integration range is obtained by limiting the range of random number generation, the number of Bernoulli samples after angle limitation is reduced from N to N ^b , and the reciprocal of the connection probability integral changes accordingly:

\frac{11}{{p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))} = = \frac{E E. [[{N N}^{b b} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))]]}{{p p}^{b b} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))}

{p p}^{b b} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) = = {&Integral; &Integral;}_{{Ω Ω}_{b b}} {p p}_{x x} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; z z)) d d z z

Wherein the integration range Ω _b is reduced from the entire hemisphere to the area projected by the collection range.

3. A kind of unbiased photon mapping drawing method according to claim 1, it is characterized in that, after obtaining the unbiased estimate of an optical path by step (3), compare this estimate with other optical path estimates by multiple importance sampling Combined, using exponential weights, the sampling probabilities of different sampling techniques in the weights are as follows:

The method of bidirectional ray tracing generates the probability of a ray path for:

{p p}_{s the s,, t t}^{B B D D. P P T T} ((\overset{&OverBar; &OverBar;}{x x})) = = {p p}^{L L} (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t})) {p p}^{E E.} (({\overset{&OverBar; &OverBar;}{x x}}_{s the s,, t t}))

in and are the generation probability of light and the generation probability of line of sight;

Probability of unpolarized photon mapping to generate a light path for:

{p p}_{{s the s}^{' '},, {t t}^{' '}}^{U u P P G G} ((\overset{&OverBar; &OverBar;}{x x})) = = {p p}^{L L} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) {p p}^{E E.} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11}))

Approximate the join probability as:

{p p}^{c c} (({\overset{&OverBar; &OverBar;}{x x}}_{{s the s}^{' '},, {t t}^{' '} - - 11})) \approx \approx m m i i n no (({πd πd}^{22} {p p}_{x x} (({z z}_{{t t}^{' '} - - 22} &RightArrow; &Right Arrow; {z z}_{{t t}^{' '} - - 11} &RightArrow; &Right Arrow; {y the y}_{{s the s}^{' '}})),, 11)) . .