http://blog.csdn.net/taigw/article/details/28627697

一，光线穿过体数据的建模

体绘制描述了在某一密度条件下，光线穿越物体时，每个体素对光线强度的影响，包括吸收，发射，散射等。对该过程有不同的建立模型的方法：

1. 吸收模型（ Absorption only ）：体素对光线只是吸收，本身既不发射光线，也不反射、透射光线；

2. 发射模型（ Emission only ）：体数据中的体素只是发射光线，不吸收光线；

3. 吸收和发射模型（ Absorption plus emission ）：这种光学模型使用最为广泛，体数据中的体素本身发射光线，并且可以吸收光线，但不对光线进行反射和散射。

4. 散射和阴影模型（ Scattering and Shading/shadowing ）：体素可以散射（反射和折射）外部光源的光线，并且由于体素之间的遮挡关系，可以产生阴影；

5. 多散射模型（ Multiple Scattering ）：光线在被眼睛观察之前，可以被多个体素散射。

通常我们使用吸收和发射模型。为了增强真实感，也可以加上阴影（包括自阴影）计算。

二，光线合成公式

吸收和发射模型的基本公式为：

公式1

其中I是光照强度，k表示体素对光的吸收系数，整个指数部分表示物体的透明度，q为体素发出的光的强度。

透明度定义为：

公式2

公式1可改写为：

公式3

将其离散化，得到：

公式4

用Ti表示第i段的透明度，Ci表示第i段所发的光：

公式5

当对光线方向上的采样点依次计算时，计算每个采样点对该条光线最终的颜色的贡献：可看成是每个采样点处发出的光Ci穿过体数据被吸收而到达观测点时剩余的部分。计算方向可以使用自前向后或者自后向前的方法，以自前向后为例，即从离观测点最近的采样点出发，依次迭代采样点直到离观测点最远的采样点，在每个采样点处用Csrc表示该点处发出的光，Cdst表示计算到该采样点时该条光线上已计算过的所有采样点合成的颜色值。依据上述公式，每当迭代到一个新的采样点时，Cdst的更新公式为：

公式6

这里用alpha=1-T表示不透明度。考虑到每个采样点处发出的光本身的不透明度，公式6进一步完善为：

公式7

三，不透明度校正

光线离散化的过程通常假设相邻采样点之间的距离是均等的。当采样率发生改变时就会产生一个问题：离散化的不透明度和颜色值需要校正，因为它们的值取决于采样间距（见公式5）。

根据公式2，当光线经过长度为delta_x的一段物体时，若假设该处的吸收系数k为常数，其透明度为T：

公式8

同理，长度为delta_x~的物体其透明度为T~=exp(-k*delta_x~), T 和T~的关系为：

公式9

用不透明度表示，得到：

公式10

四，判断采样点是否在体数据范围内

假设采样点是以纹理坐标表示的，纹理坐标的范围是texMin=vec3(0.0,0.0,0.0)到texMax=vec3(1.0,1.0,1.0)之间，判断一个点是否在体数据范围内科使用符号函数。符号函数的特性是：如果输入小于0，则返回结果-1；如果输入等于0，则返回结果0；如果输入大于0，则返回结果1. 假设采样点的纹理坐标是tPos，若采样点在体数据范围之内，则sign（tPos-texMin）与sign（texmax-tPos）的结果都是vec3(1,1,1),对这两个向量做点积则得到结果3.如果采样点不在体数据范围内，则该点积值小于3.

五，GLSL例子代码

1，使用相机位置和起点确定光线方向

[cpp] view plaincopy

void dvr_raycasting()
{
//get the 3D world coordinate for start point，assuming vUV
//is the texture coordinate of current frament
vec3 vPos=FromtPosTovPos(vUV);
//Getting the ray marching direction. camPos is a unifrom
//variable in this shader, which is the camera position
vec3 geomDir = normalize(vPos - camPos);
//multiply the raymarching direction with the step size to get the
//sub-step size we need to take at each raymarching step
float samplingStepSize_=1.0/288.0/sample_rate;
vec3 dirStep=geomDir * samplingStepSize_;
bool stop = false;
//for all samples along the ray
for (int i = 0; i < MAX_SAMPLES; i++) {
//get the position of one new sampling point
vPos=vPos+dirStep;
// convert world coordinate to texture coordinate
vec3 tPos=FromvPosTotPos(vPos);
stop = dot(sign(tPos-texMin),sign(texMax-tPos)) < 3.0;
//if the stopping condition is true we brek out of the ray marching loop
if (stop)
break;
// data fetching from the red channel of volume texture
float sample = texture(volume, tPos).r;
// get the gradient at the sample point,we use sample value and gradient
//to construct a 2D transfer function
float gradient=GetGradientMagnitude(tPos);
float src_alpha=0;
vec3 src_composition=vec3(0);
// color resulting from transfer function
vec4 lut_value;
lut_value=TransferFunction(sample,gradient);
src_alpha=lut_value.a;
src_composition=lut_value.rgb;
//opasity correlation
src_alpha=1.0-pow(1.0-src_alpha,samplingStepSize_ * SAMPLING_BASE_INTERVAL_RCP);
//Opacity calculation using compositing:
vFragColor.rgb+=src_alpha*src_composition*(1-vFragColor.a);
vFragColor.a+=src_alpha*(1-vFragColor.a);
//early ray termination
//if the currently composited colour alpha is already fully saturated
//we terminated the loop
if( vFragColor.a>0.99)
break;
}
}

2，使用光线起点和终点确定光线方向（要用到frame buffer，参见上一篇博文）

[cpp] view plaincopy

vec4 directRendering(vec3 first,vec3 last) {
vec4 result = vec4(0.0);
float samplingStepSize_=1.0/288.0/sample_rate;
vec3 direction = last - first;
vec3 normalized_dir=normalize(direction);
int steps = int(floor(length(direction)/samplingStepSize_));
vec3 diff1 = direction / float(steps);
for (int i=0; i<steps; i++) {
first += diff1;
float sample = texture(volume, first).r;
float gradient=GetGradientMagnitude(first);
float src_alpha=0;
vec3 src_composition=vec3(0);
lut_value=TransferFunction(sample,gradient);
src_alpha=lut_value.a;
src_composition=lut_value.rgb;
//opasity correlation
src_alpha=1.0-pow(1.0-src_alpha,samplingStepSize_ * SAMPLING_BASE_INTERVAL_RCP);
result.rgb+=src_alpha*src_composition*(1-result.a);
result.a+=src_alpha*(1-result.a);
//early ray termination
//if the currently composited colour alpha is already fully saturated
//we terminated the loop
if( result.a>0.99)
break;
}
return result;
}
void dvr_raycasting_front_back_position()
{
vec2 fragCoord = vec2(gl_FragCoord.xy/viewportSize);
vec3 frontPos_t = texture2D(frontTex, fragCoord.xy).rgb;
vec3 backPos_t = texture2D(backTex, fragCoord.xy).rgb;
vFragColor=directRendering(frontPos_t,backPos_t);
}

参考阅读：Ray casting的原理： http://blog.csdn.net/timzc/article/details/6020260

Real time volume graphics http://download.csdn.net/download/dragonag/3566275

体绘制的原理和Raycasting的实现相关推荐

Mybatis插件原理和PageHelper结合实战分页插件（七）
今天和大家分享下mybatis的一个分页插件PageHelper,在讲解PageHelper之前我们需要先了解下mybatis的插件原理.PageHelper 的官方网站:https://github ...
HBase学习指南之HBase原理和Shell使用
HBase学习指南之HBase原理和Shell使用参考资料: 1.https://www.cnblogs.com/nexiyi/p/hbase_shell.html,hbase shell 转载于: ...
IAP的原理和stm8的IAP
一.引出(IAP的原理和stm8上实现IAP的问题) 具有IAP功能的单片机,程序可以分为两部分:IAP和APP.APP是用来实现真正功能的程序,而IAP是用来远程重新编程APP的程序.单片机上电时会 ...
单链表反转的原理和python代码实现
链表是一种基础的数据结构,也是算法学习的重中之重.其中单链表反转是一个经常会被考察到的知识点. 单链表反转是将一个给定顺序的单链表通过算法转为逆序排列,尽管听起来很简单,但要通过算法实现也并不是非常容 ...
计算机网络原理和OSI模型与TCP模型
计算机网络原理和OSI模型与TCP模型作者:尹正杰版权声明:原创作品,谢绝转载!否则将追究法律责任. 一.计算机网络的概述 1.计算机网络的定义计算机网络是一组自治计算机的互连的集合 2.计算机 ...
HTTPS原理和CA证书申请
转载自:HTTPS原理和CA证书申请(满满的干货) 众所周知,WEB服务存在http和https两种通信方式,http默认采用80作为通讯端口,对于传输采用不加密的方式,https默认采用443,对于 ...
class加载原理和Dex加载的原理-----android插件化技术
2019独角兽企业重金招聘Python工程师标准>>> class加载原理和Dex加载的原理转载于:https://my.oschina.net/quguangle/blog/15 ...
转载椭圆曲线原理和openssl命令操作
原文地址:https://www.johannes-bauer.com/compsci/ecc/ 椭圆曲线原理和openssl命令操作
mvc原理和mvc模式的优缺点
mvc原理和mvc模式的优缺点一.mvc原理 mvc是一种程序开发设计模式,它实现了显示模块与功能模块的分离.提高了程序的可维护性.可移植性.可扩展性与可重用性,降低了程序的开发难度.它主要分 ...

体绘制的原理和Raycasting的实现