CN103295259A - 一种自适应采样的最小梯度夹角预积分光照方法 - Google Patents

Abstract

一种自适应采样的最小梯度夹角预积分光照方法，包括以下步骤：1）导入医学体数据；2）首先建立多张不同步长的预积分表，表的数量可以根据实际应用的精度需求来决定，在采样过程中，根据标量场的变化趋势进行自适应采样，在标量值变化剧烈的区域，遇到极值，可以采用较小的采样步长；而在标量值变化相对平缓的区域，采样步长可以相应增大；3）通过插值重构的方式对预积分片段内的梯度进行重构，并计算预积分采样段中与光照方向夹角最小的梯度，并将该梯度用于预积分光照渲染中，从而加强了局部细节的光照效果；4）根据自适应采样策略得到的不同步长的预积分片段，通过查找不同步长的预积分表获得相应的颜色值和不透明度，然后进行体绘制积分。

Description

一种自适应采样的最小梯度夹角预积分光照方法

技术领域

本发明涉及一种自适应采样的最小梯度夹角预积分光照方法

背景技术

医学可视化是科学可视化中最为活跃的研究领域之一，在临床医学上有着重要的实际价值，标志着一个国家的医学水平的先进程度，是一个国家综合实力的一种体现。随着核磁共振成像技术（MRI），计算机断层投影技术（CT）等医学设备的快速发展，通过医学影像设备已经可以得到人体内部某一组织比较清晰的二维数字图像序列，这些图像序列包含了人体内部各个器官的三维信息。然而医生只能通过这些二维序列图像，凭借自己的临床经验对病人的病情进行诊断，如此做法可能会造成诊断治疗的偏差。为了提高诊断的准确性和科学性，旨在提供人体组织三维空间信息以及三维观察手段的医学可视化技术应运而生。医学可视化就是运用计算机图形学和图像处理技术，将一系列的二维数字图像构造出人体组织的三维图形，从而在屏幕上形象逼真地显示人体组织内部的复杂结构，医护人员通过旋转、缩放、平移以及图像分割的操作，可以详细的了解病灶组织的空间信息，帮助医生做出准确的诊断和制定正确的手术方案，提高医学治疗的准确性和科学性。

医学可视化技术中一种重要的方法是直接体绘制技术，相比面绘制技术，直接体绘制技术可以清晰的显示人体内部的组织信息，使得绘制结果具有极强的透视效果，具有较高的保真性。近年来，随着医学成像技术的进步和仪器设备精度的提高，直接体绘制技术面临的医学体数据量越来越大，尤其在增加了光照渲染之后，采样过程的计算开销过大，梯度和光照强度值的计算都会增加计算的复杂度，因此，在实际应用中很难获得良好的交互性。

另一方面，在实际采样过程中，由于采样率过低造成采样不足，导致体数据中的某些细节特征的丢失。预积分方法能够在一定程度上解决上述的问题，该方法能够有效的克服奈奎斯特采样定律的限制，并且对于非线性变化的传输函数也有良好的绘制效果。然而，为了获得高质量的绘制结果，使用预积分方法仍然需要较小的采样步长，无法达到良好的交互性能。

综合上述，如何获得高质量的体绘制结果，同时保证用户交互的实时性，成为目前直接体绘制中一个亟待解决的重要问题。

发明内容

为了在高质量的绘制结果和良好的交互性之间得到较好的统一，本发明提供一种交互性好，准确性高，可靠性强的自适应采样的最小梯度夹角预积分光照算法。

本发明解决其技术问题所采用的技术方案是：

一种自适应采样的最小梯度夹角预积分光照方法，所述的自适应采样的预积分光照方法包括以下步骤：

1）、导入医学体数据；

2）、首先建立多张不同步长的预积分表，表的数量可以根据实际应用的精度需求来决定，在采样过程中，根据标量场的变化趋势进行自适应采样，在标量值变化剧烈的区域，遇到极值（极小值和极大值），可以采用较小的采样步长；而在标量值变化相对平缓的区域，采样步长可以相应增大。

3）、通过插值重构的方式对预积分片段内的梯度进行重构，并计算预积分采样段中与光照方向夹角最小的梯度，并将该梯度用于预积分光照渲染中，从而加强了局部细节的光照效果。4）、根据自适应采样策略得到的不同步长的预积分片段，通过查找不同步长的预积分表获得相应的颜色值和不透明度，然后进行体绘制积分。

作为优选的一种方案：所述的步骤2）包括以下步骤：

Step1  根据实际应用的需求，确定多张不同采样步长的预积分表，这些表的步长按指数级变化，并且满足如下的约束条件：

Step_n=2^-n*Step_n-1（1）

其中Step₀表示第n级别的预积分表的采样步长，公式（1）初始条件为Step₀=L（L为采样过程中的步长），并且满足n=0,1,2,3...

Step2  构建不同采样步长的预积分表，积分公式如下：

$C_i={\∫}_0^1\mathrm{\τ}((1-\mathrm{\ω})s_f+{\mathrm{\ωs}}_b)c((1-\mathrm{\ω})s_f+{\mathrm{\ωs}}_b)$ （2）

$\×\exp (-{\∫}_0^{\mathrm{\ω}}\mathrm{\τ}((1-{\mathrm{\ω}}^{\′})s_f+{\mathrm{\ω}}^{\′}{s}_b)\mathrm{Ld}{\mathrm{\ω}}^{\′})\mathrm{Ld\ω}$

其中s_f和s_b分别表示采样段Slab的入点Pa和出点Pb的标量值，L表示采样段Slab的长度；

Step3 计算采样光线s(x)上d_i处的方向导数D_v(g(x))，我们根据公式（3）来计算方向导数。

$D_v\left(g\left(x\right)\right)=\▿g\left(x\right)\·v---\left(3\right)$

其中梯度值

可以通过常规的梯度计算方法得到，即中心作差法，向前作差法等。采样点的梯度值可以预先存放在3D纹理中，或是在绘制时实时计算得到。v是光线的采样方向。

Step4根据采样段[d_i,d_i+1]的入点和出点的方向导数来判断该采样段中是否含有标量极值，从而决定是否需要对该段采样进一步细分，当入点d_i和出点d_i+1的方向导数乘积小于0时，即满足公式（4），采样段[d_i,d_i+1]需要进一步细分。

D_v(g(d_i))D_v(g(d_i+1))<0（4）

其中，D_v(g(d_i))和D_v(g(d_i+1))分别是采样光线s(x)上d_i和d_i+1处的方向导数。

Step5若公式（4）成立，则对该采样段进行自适应的划分。假定方向导数在较小的范围内是线性变化的，即在一个采样段内可以认为是线性变化的，则在方向导数为0的位置d_x取得标量值的极值。d_x的位置可以通过公式（5）计算得到：

$d_x=\frac{\left|D_v\left(g\left(d_i\right)\right)\right|}{\left|D_v\left(g\left(d_i\right)\right)\right|+|D_v\left(g\left(d_{i+1}\right)\right)}---\left(5\right)$

Step6判定d_x的大小，若d_x小于1/2，则极值点位于采样段[d_i,d_i+1]的前半段，否则位于采样段[d_i,d_i+1]的后半段；若位于前（后）半段，则计算前（后）半段两个端点的方向导数，然后计算前（后）半段内方向导数为0的位置d_x；

Step7重复step6的过程，直到采样段内最小划分的采样步长小于或等于预先设定的最小步长。

作为优选的另一种方案：所述的步骤3）还包括以下步骤：

Step1本文采用Catmull-Rom曲线来重构位于采样段Slab内的梯度值，插值公式如下：

N(t)=N₀*B₀(t)+N₁*B₁(t)+N₂*B₂(t)+N₃*B₃(t)（6）

B₀(t)=0.5(-t³+2t²-t)

B₁(t)=0.5(3t³-5t²+2)

B₂(t)=0.5(-3t³+4t²+t)

B₃(t)=0.5(t³-t²)

0≤t≤1

其中N₀,N₁,N₂,N₃分别是一条光线上连续的四个采样点各自的梯度值，并且已经归一化处理。

公式（6）经过整理之后，得到公式（7）

N(t)=0.5*(-N₀+3.0*N₁-3.0*N₂+N₃)*t³

+0.5*(2.0*N₀-5.0*N₁+4.0*N₂-N₃)*t²（7）

+0.5*(-N₀+N₂)*t

+0.5*(2.0*N₁)

Step2 最小梯度夹角预积分算法的关键在于求L·N(t)和H·N(t)的最大值，即L、H与N的夹角最小。根据公式（7）得到公式（9）

L·N(t)=0.5*(-L·N₀+3.0*L·N₁-3.0*L·N₂+L·N₃)*t³

+0.5*(2.0*L·N₀-5.0*L·N₁+4.0*L·N₂-L·N₃)*t²（8）

+0.5*(-L·N₀+L·N₂)*t

+0.5*(2.0*L·N₁)

0≤t≤1

令a=0.5*(-L·N₀+3.0*L·N₁-3.0*L·N₂+L·N₃)

b=0.5*(2.0*L·N₀-5.0*L·N₁+4.0*L·N₂-L·N₃)

c=0.5*(-L·N₀+L·N₂)

d=0.5*(2.0*L·N₁)，并且a，b，c，d都是常数

上述问题转为求F(t)=at³+bt²+ct+d在区间0≤t≤1的最大值F(t)_max。

H·N(t)求最大值的过程同理可得。

Step3 在得到F(t)_max的函数之后，接下来计算体绘制的光照渲染，在引入Blinn-Phong局部光照模型之后，体绘制的积分公式如下：

$I={\&Integral;}_0^D[c\left(s\left(\mathrm{\&lambda;}\right)\right)\mathrm{\&tau;}\left(s\left(\mathrm{\&lambda;}\right)\right)(K_a+K_d(L\&CenterDot;\&dtri;g\left(\mathrm{\&lambda;}\right)))$ （9）

$+c_s\left(s\left(\mathrm{\&lambda;}\right)\right)K_s{(H\&CenterDot;\&dtri;g\left(\mathrm{\&lambda;}\right))}^{n_{\mathrm{\&delta;}}}]e^{-{\&Integral;}_0^{\mathrm{\&lambda;}}\mathrm{\&tau;}\left(s\left(\mathrm{\&lambda;}\right)\right)\mathrm{d\&lambda;}}\mathrm{d\&lambda;}$

其中K_a，K_d，K_s分别代表的是环境光ambient，漫反射光diffuse，镜面高光specular的系数，n_δ表示光亮规模系数，c_s表示镜面高光的颜色，L表示光源到采样点方向的单位向量，H表示L与视线方向夹角的一半，

表示采样点的梯度值。

作为优选的另一种方案：所述的步骤4）还包括以下步骤：

Step1根据步骤2）中的自适应采样策略，原始采样段Slab的长度为l，假设被细分成两个采样小段Slab₁和Slab₂，长度分别为l₁和l₂，通过查找预积分表可以得到Slab₁和Slab₂颜色值分别为

和

，以及得到不透明度分别为

和

。

Step2通过将

和进行融合，可以得到采样段Slab的颜色值

，如下：

$C_{[d_i,d_{i+1}]}=C_{i,i+1,l_1}+(1+{\mathrm{\&alpha;}}_{i,i+1,l_1})C_{i,i+1,l_2}---\left(10\right)$

Step3通过将

和

进行融合，可以得到采样段Slab的颜色值如下：

${\mathrm{\&alpha;}}_{[d_i,d_{i+1}]}={\mathrm{\&alpha;}}_{i,i+1,l_1}+(1-{\mathrm{\&alpha;}}_{i,i+1,l_1}){\mathrm{\&alpha;}}_{i,i+1,l_2}$ （11）

Step4将

和

代入到原始的体绘制积分公式中，得到如下公式：

$C_{i+1}=C_i+(1-{\mathrm{\&alpha;}}_i)C_{[d_i,d_{i+1}]}$ （12）

$A_{i+1}=A_i+(1-A_i){\mathrm{\&alpha;}}_{[d_i,d_{i+1}]}$ （13）

其中C_i+1和A_i+1分别表示积累颜色值和积累不透明度值。

Step5指定积累不透明度阈值A_acc，当A_i+1≥A_acc时，结束当前的光线的采样。

本发明的技术构思为：导入需要绘制的医学体数据；首先建立多张不同步长的预积分表，表的数量可以根据实际应用的精度需求来决定，在采样过程中，根据标量场的变化趋势进行自适应采样，在标量值变化剧烈的区域，遇到极值（极小值和极大值），可以采用较小的采样步长；而在标量值变化相对平缓的区域，采样步长可以相应增大。通过插值重构的方式对预积分片段内的梯度进行重构，并计算预积分采样段中与光照方向夹角最小的梯度，并将该梯度用于预积分光照渲染中，从而加强了局部细节的光照效果。根据自适应采样策略得到的不同步长的预积分片段，通过查找不同步长的预积分表获得相应的颜色值和不透明度，然后进行体绘制积分。

本发明的有益效果是：实现一种自适应采样的最小梯度夹角预积分光照方法，可以在临床应用中将高质量的绘制结果和良好的交互性之间得到较好的统一，同时这个过程不需要或是尽可能少的要求用户去干预。

附图说明

图1是自适应采样的最小梯度夹角预积分光照方法的系统总体结构图。

图2是自适应采样的预积分光照方法的程序简要流程。

图3是自适应采样过程的简要流程图。

图4是计算最小梯度夹角的流程图。

图5是自适应采样的预积分体绘制积分的简要流程图。

具体实施方式

下面结合附图对本发明作进一步描述：

参照图1～图5，一种自适应采样的最小梯度夹角预积分光照方法，所述自适应采样的预积分光照算法包括以下步骤：

1）、导入医学体数据；

2）、首先建立多张不同步长的预积分表，表的数量可以根据实际应用的精度需求来决定，在采样过程中，根据标量场的变化趋势进行自适应采样，在标量值变化剧烈的区域，遇到极值（极小值和极大值），可以采用较小的采样步长；而在标量值变化相对平缓的区域，采样步长可以相应增大。

3）、通过插值重构的方式对预积分片段内的梯度进行重构，并计算预积分采样段中与光照方向夹角最小的梯度，并将该梯度用于预积分光照渲染中，从而加强了局部细节的光照效果。4）、根据自适应采样策略得到的不同步长的预积分片段，通过查找不同步长的预积分表获得相应的颜色值和不透明度，然后进行体绘制积分。

图1中，自适应采样的最小梯度夹角预积分光照方法系统总体结构图是本方法的总体方案，数据导入以及处理部分使用C++语言实现，自适应采样及光照渲染部分利用GLSL语言编写，整个程序在VS2008平台上实现。

图2中，自适应采样的最小梯度夹角预积分光照方法的步骤是：首先导入医学体数据，然后根据标量场的变化趋势进行自适应采样，将预积分采样段进一步细分，通过插值重构的方式对预积分片段内的梯度进行重构，计算得到预积分片段中与光照方向夹角最小的梯度，最后查找不同部长的预积分表获得相应的颜色值和不透明度，接着进行体绘制积分。

在上述的自适应采样的预积分光照方法的步骤中，如何确定预积分采样段内极值点的位置，对后续进一步细分采样步长起着至关重要的作用，也即直接影响着绘制结果的质量，是该方法的核心步骤。

图3中描述了如何对预积分采样段进行自适应采样的划分，首先计算采样光线s(x)上d_i处的方向导数D_v(g(x))

$D_v\left(g\left(x\right)\right)=\&dtri;g\left(x\right)\&CenterDot;v$

其中梯度值

可以通过常规的梯度计算方法得到，即中心作差法，向前作差法等。v是光线的采样方向。

根据采样段[d_i,d_i+1]的入点和出点的方向导数来判断该采样段中是否含有标量极值，从而决定是否需要对该段采样进一步细分，当入点d_i和出点d_i+1的方向导数乘积小于0时，采样段[d_i,d_i+1]需要进一步细分。

D_v(g(d_i))D_v(g(d_i+1))<0

其中，D_v(g(d_i))和D_v(g(d_i+1))分别是采样光线s(x)上d_i和d_i+1处的方向导数。

假定方向导数在较小的范围内是线性变化的，即在一个采样段内可以认为是线性变化的，则在方向导数为0的位置d_x取得标量值的极值。d_x的位置可以通过如下公式计算得到：

$d_x=\frac{\left|D_v\left(g\left(d_i\right)\right)\right|}{\left|D_v\left(g\left(d_i\right)\right)\right|+\left|D_v\left(g\left(d_{i+1}\right)\right)\right|}$

判定d_x的大小，若d_x小于1/2，则极值点位于采样段[d_i,d_i+1]的前半段，否则位于采样段[d_i,d_i+1]的后半段；若位于前（后）半段，则计算前（后）半段两个端点的方向导数，然后计算前（后）半段内方向导数为0的位置d_x，直到采样段内最小划分的采样步长小于或等于预先设定的最小步长。

图4描述了计算最小梯度夹角的过程，首先采用Catmull-Rom曲线来重构位于采样段Slab内的梯度值：

N(t)=N₀*B₀(t)+N₁*B₁(t)+N₂*B₂(t)+N₃*B₃(t)

B₀(t)=0.5(-t³+2t²-t)

B₁(t)=0.5(3t³-5t²+2)

B₂(t)=0.5(-3t³+4t²+t)

B₃(t)=0.5(t³-t²)

0≤t≤1

其中N₀,N₁,N₂,N₃分别是一条光线上连续的四个采样点各自的梯度值，并且已经归一化处理。

整理之后，得到如下公式：

N(t)=0.5*(-N₀+3.0*N₁-3.0*N₂+N₃)*t³

+0.5*(2.0*N₀-5.0*N₁+4.0*N₂-N₃)*t²

+0.5*(-N₀+N₂)*t

+0.5*(2.0*N₁)

接下来计算L·N(t)和H·N(t)的最大值，即L、H与N的夹角最小。

L·N(t)=0.5*(-L·N₀+3.0*L·N₁-3.0*L·N₂+L·N₃)*t³

+0.5*(2.0*L·N₀-5.0*L·N₁+4.0*L·N₂-L·N₃)*t²

+0.5*(-L·N₀+L·N₂)*t

+0.5*(2.0*L·N₁)

0≤t≤1

令a=0.5*(-L·N₀+3.0*L·N₁-3.0*L·N₂+L·N₃)

b=0.5*(2.0*L·N₀-5.0*L·N₁+4.0*L·N₂-L·N₃)

c=0.5*(-L·N₀+L·N₂)

d=0.5*(2.0*L·N₁)，并且a，b，c，d都是常数

上述问题转为求F(t)=at³+bt²+ct+d在区间0≤t≤1的最大值F(t)_max。

H·N(t)求最大值的过程同理可得。

在得到F(t)_max的函数之后，接下来计算体绘制的光照渲染，在引入Blinn-Phong局部光照模型之后，体绘制的积分公式如下：

$I={\&Integral;}_0^D[c\left(s\left(\mathrm{\&lambda;}\right)\right)\mathrm{\&tau;}\left(s\left(\mathrm{\&lambda;}\right)\right)(K_a+K_d(L\&CenterDot;\&dtri;g\left(\mathrm{\&lambda;}\right)))$

$+c_s\left(s\left(\mathrm{\&lambda;}\right)\right)K_s{(H\&CenterDot;\&dtri;g\left(\mathrm{\&lambda;}\right))}^{n_{\mathrm{\&delta;}}}]e^{-{\&Integral;}_0^{\mathrm{\&lambda;}}\mathrm{\&tau;}\left(s\left(\mathrm{\&lambda;}\right)\right)\mathrm{d\&lambda;}}\mathrm{d\&lambda;}$

其中K_a，K_d，K_s分别代表的是环境光ambient，漫反射光diffuse，镜面高光specular的系数，n_δ表示光亮规模系数，c_s表示镜面高光的颜色，L表示光源到采样点方向的单位向量，H表示L与视线方向夹角的一半，

表示采样点的梯度值。

图5描述了自适应采样的预积分体绘制积分的过程。首先根据自适应采样策略，原始采样段Slab的长度为l，假设被细分成两个采样小段Slab₁和Slab₂，长度分别为l₁和l₂，通过查找预积分表可以得到Slab₁和Slab₂颜色值分别为

和，以及得到不透明度分别为

和

。

通过将

和

进行融合，可以得到采样段Slab的颜色值

如下： $C_{[d_i,d_{i+1}]}=C_{i,i+1,l_1}+(1+{\mathrm{\&alpha;}}_{i,i+1,l_1})C_{i,i+1,l_2}$

通过将

和进行融合，可以得到采样段Slab的颜色值

如下：

${\mathrm{\&alpha;}}_{[d_i,d_{i+1}]}={\mathrm{\&alpha;}}_{i,i+1,l_1}+(1-{\mathrm{\&alpha;}}_{i,i+1,l_1}){\mathrm{\&alpha;}}_{i,i+1,l_2}$

将和

代入到原始的体绘制积分公式中，得到如下公式：

$C_{i+1}=C_i+(1-{\mathrm{\&alpha;}}_i)C_{[d_i,d_{i+1}]}$

$A_{i+1}=A_i+(1-A_i){\mathrm{\&alpha;}}_{[d_i,d_{i+1}]}$

其中C_i+1和A_i+1分别表示积累颜色值和积累不透明度值。

最后根据设定的积累不透明度阈值A_acc，当A_i+1≥A_acc时，结束当前的光线的采样。

Claims (4)

1.一种自适应采样的最小梯度夹角预积分光照方法，所述的自适应采样的预积分光照方法包括以下步骤：

1）、导入医学体数据；

2）、首先建立多张不同步长的预积分表，表的数量可以根据实际应用的精度需求来决定，在采样过程中，根据标量场的变化趋势进行自适应采样，在标量值变化剧烈的区域，遇到极值（极小值和极大值），可以采用较小的采样步长；而在标量值变化相对平缓的区域，采样步长可以相应增大；

3）、通过插值重构的方式对预积分片段内的梯度进行重构，并计算预积分采样段中与光照方向夹角最小的梯度，并将该梯度用于预积分光照渲染中，从而加强了局部细节的光照效果；

4）、根据自适应采样策略得到的不同步长的预积分片段，通过查找不同步长的预积分表获得相应的颜色值和不透明度，然后进行体绘制积分。

2.如权利要求1所述的方法，其特征在于：所述的步骤2）包括以下步骤：

Step1根据实际应用的需求，确定多张不同采样步长的预积分表，这些表的步长按指数级变化，并且满足如下的约束条件：

Step_n=2^-n*Step_n-1（1）

其中Step₀表示第n级别的预积分表的采样步长，公式（1）初始条件为Step₀=L（L为采样过程中的步长），并且满足n=0,1,2,3...

Step2构建不同采样步长的预积分表，积分公式如下：

$C_i={\&Integral;}_0^1\mathrm{\&tau;}((1-\mathrm{\&omega;})s_f+{\mathrm{\&omega;s}}_b)c((1-\mathrm{\&omega;})s_f+{\mathrm{\&omega;s}}_b)$ （2）

$\&times;\exp (-{\&Integral;}_0^{\mathrm{\&omega;}}\mathrm{\&tau;}((1-{\mathrm{\&omega;}}^{\&prime;})s_f+{\mathrm{\&omega;}}^{\&prime;}{s}_b)\mathrm{Ld}{\mathrm{\&omega;}}^{\&prime;})\mathrm{Ld\&omega;}$

其中s_f和s_b分别表示采样段Slab的入点Pa和出点Pb的标量值，L表示采样段Slab的长度；

Step3计算采样光线s(x)上d_i处的方向导数D_v(g(x))，我们根据公式（3）来计算方向导数。

$D_v\left(g\left(x\right)\right)=\&dtri;g\left(x\right)\&CenterDot;v---\left(3\right)$

其中梯度值

可以通过常规的梯度计算方法得到，即中心作差法，向前作差法等。采样点的梯度值可以预先存放在3D纹理中，或是在绘制时实时计算得到。v是光线的采样方向。

Step4根据采样段[d_i,d_i+1]的入点和出点的方向导数来判断该采样段中是否含有标量极值，从而决定是否需要对该段采样进一步细分，当入点d_i和出点d_i+1的方向导数乘积小于0时，即满足公式（4），采样段[d_i,d_i+1]需要进一步细分。

D_v(g(d_i))D_v(g(d_i+1))<0（4）

其中，D_v(g(d_i))和D_v(g(d_i+1))分别是采样光线s(x)上d_i和d_i+1处的方向导数。

Step5若公式（4）成立，则对该采样段进行自适应的划分。假定方向导数在较小的范围内是线性变化的，即在一个采样段内可以认为是线性变化的，则在方向导数为0的位置d_x取得标量值的极值。d_x的位置可以通过公式（5）计算得到：

$d_x=\frac{\left|D_v\left(g\left(d_i\right)\right)\right|}{\left|D_v\left(g\left(d_i\right)\right)\right|+|D_v\left(g\left(d_{i+1}\right)\right)}---\left(5\right)$

Step6判定d_x的大小，若d_x小于1/2，则极值点位于采样段[d_i,d_i+1]的前半段，否则位于采样段[d_i,d_i+1]的后半段；若位于前（后）半段，则计算前（后）半段两个端点的方向导数，然后计算前（后）半段内方向导数为0的位置d_x；

Step7重复step6的过程，直到采样段内最小划分的采样步长小于或等于预先设定的最小步长。

3.如权利要求2所述的方法，其特征在于：：所述的步骤3）还包括以下步骤：

Step1本文采用Catmull-Rom曲线来重构位于采样段Slab内的梯度值，插值公式如下：

N(t)=N₀*B₀(t)+N₁*B₁(t)+N₂*B₂(t)+N₃*B₃(t)（6）

B₀(t)=0.5(-t³+2t²-t)

B₁(t)=0.5(3t³-5t²+2)

B₂(t)=0.5(-3t³+4t²+t)

B₃(t)=0.5(t³-t²)

0≤t≤1

其中N₀,N₁,N₂,N₃分别是一条光线上连续的四个采样点各自的梯度值，并且已经归一化处理。

公式（6）经过整理之后，得到公式（7）

N(t)=0.5*(-N₀+3.0*N₁-3.0*N₂+N₃)*t³

+0.5*(2.0*N₀-5.0*N₁+4.0*N₂-N₃)*t²（7）

+0.5*(-N₀+N₂)*t

+0.5*(2.0*N₁)

Step2最小梯度夹角预积分算法的关键在于求L·N(t)和H·N(t)的最大值，即L、H与N的夹角最小。根据公式（7）得到公式（9）

L·N(t)=0.5*(-L·N₀+3.0*L·N₁-3.0*L·N₂+L·N₃)*t³

+0.5*(2.0*L·N₀-5.0*L·N₁+4.0*L·N₂-L·N₃)*t²（8）

+0.5*(-L·N₀+L·N₂)*t

+0.5*(2.0*L·N₁)

0≤t≤1

令a=0.5*(-L·N₀+3.0*L·N₁-3.0*L·N₂+L·N₃)

b=0.5*(2.0*L·N₀-5.0*L·N₁+4.0*L·N₂-L·N₃)

c=0.5*(-L·N₀+L·N₂)

d=0.5*(2.0*L·N₁)，并且a，b，c，d都是常数

上述问题转为求F(t)=at³+bt²+ct+d在区间0≤t≤1的最大值F(t)_max。

H·N(t)求最大值的过程同理可得。

Step3  在得到F(t)_max的函数之后，接下来计算体绘制的光照渲染，在引入Blinn-Phong局部光照模型之后，体绘制的积分公式如下：

$I={\&Integral;}_0^D[c\left(s\left(\mathrm{\&lambda;}\right)\right)\mathrm{\&tau;}\left(s\left(\mathrm{\&lambda;}\right)\right)(K_a+K_d(L\&CenterDot;\&dtri;g\left(\mathrm{\&lambda;}\right)))$ （9）

$+c_s\left(s\left(\mathrm{\&lambda;}\right)\right)K_s{(H\&CenterDot;\&dtri;g\left(\mathrm{\&lambda;}\right))}^{n_{\mathrm{\&delta;}}}]e^{-{\&Integral;}_0^{\mathrm{\&lambda;}}\mathrm{\&tau;}\left(s\left(\mathrm{\&lambda;}\right)\right)\mathrm{d\&lambda;}}\mathrm{d\&lambda;}$

其中K_a，K_d，K_s分别代表的是环境光ambient，漫反射光diffuse，镜面高光specular的系数，n_δ表示光亮规模系数，c_s表示镜面高光的颜色，L表示光源到采样点方向的单位向量，H表示L与视线方向夹角的一半，

表示采样点的梯度值。

4.如权利要求3所述的方法，其特征在于：：所述的步骤4）还包括以下步骤：

Step1 根据步骤2）中的自适应采样策略，原始采样段Slab的长度为l，假设被细分成两个采样小段Slab₁和Slab₂，长度分别为l₁和l₂，通过查找预积分表可以得到Slab₁和Slab₂颜色值分别为

和

以及得到不透明度分别为

和

Step2 通过将和

进行融合，可以得到采样段Slab的颜色值

如下：

$C_{[d_i,d_{i+1}]}=C_{i,i+1,l_1}+(1-{\mathrm{\&alpha;}}_{i,i+1,l_1})C_{i,i+1,l_2}---\left(10\right)$

Step3 通过将

和

进行融合，可以得到采样段Slab的颜色值

如下：

${\mathrm{\&alpha;}}_{[d_i,d_{i+1}]}={\mathrm{\&alpha;}}_{i,i+1,l_1}+(1-{\mathrm{\&alpha;}}_{i,i+1,l_1}){\mathrm{\&alpha;}}_{i,i+1,l_2}---\left(11\right)$

Step4 将

和

代入到原始的体绘制积分公式中，得到如下公式：

$C_{i+1}=C_i+(1-{\mathrm{\&alpha;}}_i)C_{[d_i,d_{i+1}]}---\left(12\right)$

$A_{i+1}=A_i+(1-A_i){\mathrm{\&alpha;}}_{[d_i,d_{i+1}]}---\left(13\right)$ （13）

其中C_i+1和A_i+1分别表示积累颜色值和积累不透明度值。

Step5指定积累不透明度阈值A_acc，当A_i+1≥A_acc时，结束当前的光线的采样。

Images