我正在实现一种体积渲染算法“GPU 光线投射单通道”。为此，我使用强度值的浮点数组作为 3d 纹理（此 3d 纹理描述了球坐标中的常规 3d 网格）。
这里有数组值的示例：
75.839354473071637,
64.083049468866022,
65.253933716444365,
79.992431196592577,
84.411485976957096,
0.0000000000000000,
82.020319431382831,
76.808403454586994,
79.974774618246158,
0.0000000000000000,
91.127273013466336,
84.009956557448433,
90.221356094672814,
87.567422484025627,
71.940263118478072,
0.0000000000000000,
0.0000000000000000,
74.487058398181944,
..................,
..................
（完整数据如下：[链接]（https://drive.google.com/file/d/1lbXzRucUseF-ITzFgxqeLTd0WglJJOoz/view?usp=sharing))
球形网格的尺寸为(r,theta,phi)=(384,15,768)，这是加载纹理的输入格式：
glTexImage3D(GL_TEXTURE_3D, 0, GL_R16F, 384, 15, 768, 0, GL_RED, GL_FLOAT, dataArray)
这是我的可视化图像：
问题是可视化应该是一个磁盘，或者至少是一个类似的形式。
我认为问题是我没有正确指定纹理的坐标（在球面坐标中）。
这是顶点着色器代码：
#version 330 core
layout(location = 0) in vec3 vVertex; //object space vertex position
//uniform
uniform mat4 MVP; //combined modelview projection matrix
smooth out vec3 vUV; //3D texture coordinates for texture lookup in the fragment shader
void main()
{
//get the clipspace position
gl_Position = MVP*vec4(vVertex.xyz,1);
//get the 3D texture coordinates by adding (0.5,0.5,0.5) to the object space
//vertex position. Since the unit cube is at origin (min: (-0.5,-0.5,-0.5) and max: (0.5,0.5,0.5))
//adding (0.5,0.5,0.5) to the unit cube object space position gives us values from (0,0,0) to
//(1,1,1)
vUV = vVertex + vec3(0.5);
}
这是片段着色器代码：
#version 330 core
layout(location = 0) out vec4 vFragColor; //fragment shader output
smooth in vec3 vUV; //3D texture coordinates form vertex shader
//interpolated by rasterizer
//uniforms
uniform sampler3D volume; //volume dataset
uniform vec3 camPos; //camera position
uniform vec3 step_size; //ray step size
//constants
const int MAX_SAMPLES = 300; //total samples for each ray march step
const vec3 texMin = vec3(0); //minimum texture access coordinate
const vec3 texMax = vec3(1); //maximum texture access coordinate
vec4 colour_transfer(float intensity)
{
vec3 high = vec3(100.0, 20.0, 10.0);
// vec3 low = vec3(0.0, 0.0, 0.0);
float alpha = (exp(intensity) - 1.0) / (exp(1.0) - 1.0);
return vec4(intensity * high, alpha);
}
void main()
{
//get the 3D texture coordinates for lookup into the volume dataset
vec3 dataPos = vUV;
//Getting the ray marching direction:
//get the object space position by subracting 0.5 from the
//3D texture coordinates. Then subtraact it from camera position
//and normalize to get the ray marching direction
vec3 geomDir = normalize((vUV-vec3(0.5)) - camPos);
//multiply the raymarching direction with the step size to get the
//sub-step size we need to take at each raymarching step
vec3 dirStep = geomDir * step_size;
//flag to indicate if the raymarch loop should terminate
bool stop = false;
//for all samples along the ray
for (int i = 0; i < MAX_SAMPLES; i++) {
// advance ray by dirstep
dataPos = dataPos + dirStep;
stop = dot(sign(dataPos-texMin),sign(texMax-dataPos)) < 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, dataPos).r;
vec4 c = colour_transfer(sample);
vFragColor.rgb = c.a * c.rgb + (1 - c.a) * vFragColor.a * vFragColor.rgb;
vFragColor.a = c.a + (1 - c.a) * 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;
}
}
我如何指定坐标以便在球坐标中可视化 3D 纹理中的信息？
UPDATE:
顶点着色器：
#version 330 core
layout(location = 0) in vec3 vVertex; //object space vertex position
//uniform
uniform mat4 MVP; //combined modelview projection matrix
smooth out vec3 vUV; //3D texture coordinates for texture lookup in the fragment shader
void main()
{
//get the clipspace position
gl_Position = MVP*vec4(vVertex.xyz,1);
//get the 3D texture coordinates by adding (0.5,0.5,0.5) to the object space
//vertex position. Since the unit cube is at origin (min: (-0.5,- 0.5,-0.5) and max: (0.5,0.5,0.5))
//adding (0.5,0.5,0.5) to the unit cube object space position gives us values from (0,0,0) to
//(1,1,1)
vUV = vVertex + vec3(0.5);
}
和片段着色器：
#version 330 core
#define Pi 3.1415926535897932384626433832795
layout(location = 0) out vec4 vFragColor; //fragment shader output
smooth in vec3 vUV; //3D texture coordinates form vertex shader
//interpolated by rasterizer
//uniforms
uniform sampler3D volume; //volume dataset
uniform vec3 camPos; //camera position
uniform vec3 step_size; //ray step size
//constants
const int MAX_SAMPLES = 200; //total samples for each ray march step
const vec3 texMin = vec3(0); //minimum texture access coordinate
const vec3 texMax = vec3(1); //maximum texture access coordinate
// transfer function that asigned a color and alpha from sample intensity
vec4 colour_transfer(float intensity)
{
vec3 high = vec3(100.0, 20.0, 10.0);
// vec3 low = vec3(0.0, 0.0, 0.0);
float alpha = (exp(intensity) - 1.0) / (exp(1.0) - 1.0);
return vec4(intensity * high, alpha);
}
// this function transform vector in spherical coordinates from cartesian
vec3 cart2Sphe(vec3 cart){
vec3 sphe;
sphe.x = sqrt(cart.x*cart.x+cart.y*cart.y+cart.z*cart.z);
sphe.z = atan(cart.y/cart.x);
sphe.y = atan(sqrt(cart.x*cart.x+cart.y*cart.y)/cart.z);
return sphe;
}
void main()
{
//get the 3D texture coordinates for lookup into the volume dataset
vec3 dataPos = vUV;
//Getting the ray marching direction:
//get the object space position by subracting 0.5 from the
//3D texture coordinates. Then subtraact it from camera position
//and normalize to get the ray marching direction
vec3 vec=(vUV-vec3(0.5));
vec3 spheVec=cart2Sphe(vec); // transform position to spherical
vec3 sphePos=cart2Sphe(camPos); //transform camPos to spherical
vec3 geomDir= normalize(spheVec-sphePos); // ray direction
//multiply the raymarching direction with the step size to get the
//sub-step size we need to take at each raymarching step
vec3 dirStep = geomDir * step_size ;
//flag to indicate if the raymarch loop should terminate
//for all samples along the ray
for (int i = 0; i < MAX_SAMPLES; i++) {
// advance ray by dirstep
dataPos = dataPos + dirStep;
float sample;
convert texture coordinates
vec3 spPos;
spPos.x=dataPos.x/384;
spPos.y=(dataPos.y+(Pi/2))/Pi;
spPos.z=dataPos.z/(2*Pi);
// get value from texture
sample = texture(volume,dataPos).r;
vec4 c = colour_transfer(sample)
// alpha blending function
vFragColor.rgb = c.a * c.rgb + (1 - c.a) * vFragColor.a * vFragColor.rgb;
vFragColor.a = c.a + (1 - c.a) * vFragColor.a;
if( vFragColor.a>1.0)
break;
}
// vFragColor.rgba = texture(volume,dataPos);
}
这些是生成边界立方体的点：
glm::vec3 vertices[8] = {glm::vec3(-0.5f, -0.5f, -0.5f),
glm::vec3(0.5f, -0.5f, -0.5f),
glm::vec3(0.5f, 0.5f, -0.5f),
glm::vec3(-0.5f, 0.5f, -0.5f),
glm::vec3(-0.5f, -0.5f, 0.5f),
glm::vec3(0.5f, -0.5f, 0.5f),
glm::vec3(0.5f, 0.5f, 0.5f),
glm::vec3(-0.5f, 0.5f, 0.5f)};
//unit cube indices
GLushort cubeIndices[36] = {0, 5, 4,
5, 0, 1,
3, 7, 6,
3, 6, 2,
7, 4, 6,
6, 4, 5,
2, 1, 3,
3, 1, 0,
3, 0, 7,
7, 0, 4,
6, 5, 2,
2, 5, 1};
这是它生成的可视化效果：
我不知道你渲染什么以及如何渲染。有许多技术和配置可以实现它们。我通常使用覆盖屏幕/视图的单通道单四边形渲染，而几何图形/场景作为纹理传递。当你的对象是 3D 纹理时，我认为你也应该这样做。这是它的完成方式（假设透视，均匀的球形体素网格作为 3D 纹理）：
CPU端代码
简单地渲染单个QUAD覆盖场景/视图。为了使这个更加简单和精确，我建议您使用球体局部坐标系作为传递给着色器的相机矩阵（它将大大简化光线/球体相交计算）。
Vertex
在这里，您应该投射/计算每个顶点的光线位置和方向，并将其传递给片段，以便为屏幕/视图上的每个像素进行插值。
所以相机是通过它的位置（焦点）和视角方向（通常是透视中的Z轴）来描述的。OpenGL）。光线从焦点投射(0,0,0)在相机局部坐标中znear plane (x,y,-znear)也在相机本地坐标中。在哪里x,y是在屏幕/视图不是正方形时应用纵横比校正的像素屏幕位置。
因此，您只需将这两个点转换为球体局部坐标（仍然是笛卡尔坐标）。
射线方向只是两点相减......
Fragment
首先规范化从顶点传递的光线方向（由于插值，它不会是单位向量）。之后，只需从外到内测试球体素网格每个半径的射线/球体相交，以便测试球体rmax to rmax/n where rmax是 3D 纹理可以拥有的最大半径，n是与半径对应的轴的 ids 分辨率r.
每次点击时将笛卡尔交点位置转换为球坐标。将它们转换为纹理坐标s,t,p并获取体素强度并将其应用到颜色（具体取决于您渲染的内容和方式）。
所以如果你的纹理坐标是(r,theta,phi)假设phi是经度，角度标准化为<-Pi,Pi> and <0,2*Pi> and rmax是 3D 纹理的最大半径：
s = r/rmax
t = (theta+(Pi/2))/Pi
p = phi/(2*PI)
如果您的球体不透明，则在体素强度不为空的第一次命中时停止。否则更新光线开始位置并再次执行整个项目符号，直到光线离开场景 BBOX 或没有发生相交。
您还可以通过在物体边界命中时分割光线来添加斯涅耳定律（添加反射折射）...
这里有一些相关的QA使用此技术或拥有有效信息可以帮助您实现此目标：
[Edit1] 示例（输入 3D 纹理最终发布后
因此，当我将上面（和评论中）的所有内容放在一起时，我想出了这个。
CPU端代码：
//---------------------------------------------------------------------------
//--- GLSL Raytrace system ver: 1.000 ---------------------------------------
//---------------------------------------------------------------------------
#ifndef _raytrace_spherical_volume_h
#define _raytrace_spherical_volume_h
//---------------------------------------------------------------------------
class SphericalVolume3D
{
public:
bool _init; // has been initiated ?
GLuint txrvol; // SphericalVolume3D texture at GPU side
int xs,ys,zs;
float eye[16]; // direct camera matrix
float aspect,focal_length;
SphericalVolume3D() { _init=false; txrvol=-1; xs=0; ys=0; zs=0; aspect=1.0; focal_length=1.0; }
SphericalVolume3D(SphericalVolume3D& a) { *this=a; }
~SphericalVolume3D() { gl_exit(); }
SphericalVolume3D* operator = (const SphericalVolume3D *a) { *this=*a; return this; }
//SphericalVolume3D* operator = (const SphericalVolume3D &a) { ...copy... return this; }
// init/exit
void gl_init();
void gl_exit();
// render
void glsl_draw(GLint prog_id);
};
//---------------------------------------------------------------------------
void SphericalVolume3D::gl_init()
{
if (_init) return; _init=true;
// load 3D texture from file into CPU side memory
int hnd,siz; BYTE *dat;
hnd=FileOpen("Texture3D_F32.dat",fmOpenRead);
siz=FileSeek(hnd,0,2);
FileSeek(hnd,0,0);
dat=new BYTE[siz];
FileRead(hnd,dat,siz);
FileClose(hnd);
if (0)
{
int i,n=siz/sizeof(GLfloat);
GLfloat *p=(GLfloat*)dat;
for (i=0;i<n;i++) p[i]=100.5;
}
// copy it to GPU as 3D texture
// glClampColorARB(GL_CLAMP_VERTEX_COLOR_ARB, GL_FALSE);
// glClampColorARB(GL_CLAMP_READ_COLOR_ARB, GL_FALSE);
// glClampColorARB(GL_CLAMP_FRAGMENT_COLOR_ARB, GL_FALSE);
glGenTextures(1,&txrvol);
glEnable(GL_TEXTURE_3D);
glBindTexture(GL_TEXTURE_3D,txrvol);
glPixelStorei(GL_UNPACK_ALIGNMENT, 4);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S,GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T,GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R,GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER,GL_NEAREST);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER,GL_NEAREST);
glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE,GL_MODULATE);
xs=384;
ys= 15;
zs=768;
glTexImage3D(GL_TEXTURE_3D, 0, GL_R16F, xs,ys,zs, 0, GL_RED, GL_FLOAT, dat);
glBindTexture(GL_TEXTURE_3D,0);
glDisable(GL_TEXTURE_3D);
delete[] dat;
}
//---------------------------------------------------------------------------
void SphericalVolume3D::gl_exit()
{
if (!_init) return; _init=false;
glDeleteTextures(1,&txrvol);
}
//---------------------------------------------------------------------------
void SphericalVolume3D::glsl_draw(GLint prog_id)
{
GLint ix;
const int txru_vol=0;
glUseProgram(prog_id);
// uniforms
ix=glGetUniformLocation(prog_id,"zoom" ); glUniform1f(ix,1.0);
ix=glGetUniformLocation(prog_id,"aspect" ); glUniform1f(ix,aspect);
ix=glGetUniformLocation(prog_id,"focal_length"); glUniform1f(ix,focal_length);
ix=glGetUniformLocation(prog_id,"vol_xs" ); glUniform1i(ix,xs);
ix=glGetUniformLocation(prog_id,"vol_ys" ); glUniform1i(ix,ys);
ix=glGetUniformLocation(prog_id,"vol_zs" ); glUniform1i(ix,zs);
ix=glGetUniformLocation(prog_id,"vol_txr" ); glUniform1i(ix,txru_vol);
ix=glGetUniformLocation(prog_id,"tm_eye" ); glUniformMatrix4fv(ix,1,false,eye);
glActiveTexture(GL_TEXTURE0+txru_vol);
glEnable(GL_TEXTURE_3D);
glBindTexture(GL_TEXTURE_3D,txrvol);
// this should be a VAO/VBO
glColor4f(1.0,1.0,1.0,1.0);
glBegin(GL_QUADS);
glVertex2f(-1.0,-1.0);
glVertex2f(-1.0,+1.0);
glVertex2f(+1.0,+1.0);
glVertex2f(+1.0,-1.0);
glEnd();
glActiveTexture(GL_TEXTURE0+txru_vol);
glBindTexture(GL_TEXTURE_3D,0);
glDisable(GL_TEXTURE_3D);
glUseProgram(0);
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
当 GL 已经初始化时，在应用程序启动时调用 init，在应用程序退出之前退出，同时 GL 仍然工作并在需要时绘制...代码是C++/VCL基于 so 端口到您的环境（文件访问、字符串等）。我还使用二进制形式的 3D 纹理，因为加载 85MByte ASCII 文件对我来说有点太多了。
Vertex:
//------------------------------------------------------------------
#version 420 core
//------------------------------------------------------------------
uniform float aspect;
uniform float focal_length;
uniform float zoom;
uniform mat4x4 tm_eye;
layout(location=0) in vec2 pos;
out smooth vec3 ray_pos; // ray start position
out smooth vec3 ray_dir; // ray start direction
//------------------------------------------------------------------
void main(void)
{
vec4 p;
// perspective projection
p=tm_eye*vec4(pos.x/(zoom*aspect),pos.y/zoom,0.0,1.0);
ray_pos=p.xyz;
p-=tm_eye*vec4(0.0,0.0,-focal_length,1.0);
ray_dir=normalize(p.xyz);
gl_Position=vec4(pos,0.0,1.0);
}
//------------------------------------------------------------------
它或多或少是体积射线追踪器链接的副本。
分段：
//------------------------------------------------------------------
#version 420 core
//------------------------------------------------------------------
// Ray tracer ver: 1.000
//------------------------------------------------------------------
in smooth vec3 ray_pos; // ray start position
in smooth vec3 ray_dir; // ray start direction
uniform int vol_xs, // texture resolution
vol_ys,
vol_zs;
uniform sampler3D vol_txr; // scene mesh data texture
out layout(location=0) vec4 frag_col;
//---------------------------------------------------------------------------
// compute length of ray(p0,dp) to intersection with ellipsoid((0,0,0),r) -> view_depth_l0,1
// where r.x is elipsoid rx^-2, r.y = ry^-2 and r.z=rz^-2
float view_depth_l0=-1.0,view_depth_l1=-1.0;
bool _view_depth(vec3 _p0,vec3 _dp,vec3 _r)
{
double a,b,c,d,l0,l1;
dvec3 p0,dp,r;
p0=dvec3(_p0);
dp=dvec3(_dp);
r =dvec3(_r );
view_depth_l0=-1.0;
view_depth_l1=-1.0;
a=(dp.x*dp.x*r.x)
+(dp.y*dp.y*r.y)
+(dp.z*dp.z*r.z); a*=2.0;
b=(p0.x*dp.x*r.x)
+(p0.y*dp.y*r.y)
+(p0.z*dp.z*r.z); b*=2.0;
c=(p0.x*p0.x*r.x)
+(p0.y*p0.y*r.y)
+(p0.z*p0.z*r.z)-1.0;
d=((b*b)-(2.0*a*c));
if (d<0.0) return false;
d=sqrt(d);
l0=(-b+d)/a;
l1=(-b-d)/a;
if (abs(l0)>abs(l1)) { a=l0; l0=l1; l1=a; }
if (l0<0.0) { a=l0; l0=l1; l1=a; }
if (l0<0.0) return false;
view_depth_l0=float(l0);
view_depth_l1=float(l1);
return true;
}
//---------------------------------------------------------------------------
const float pi =3.1415926535897932384626433832795;
const float pi2=6.2831853071795864769252867665590;
float atanxy(float x,float y) // atan2 return < 0 , 2.0*M_PI >
{
int sx,sy;
float a;
const float _zero=1.0e-30;
sx=0; if (x<-_zero) sx=-1; if (x>+_zero) sx=+1;
sy=0; if (y<-_zero) sy=-1; if (y>+_zero) sy=+1;
if ((sy==0)&&(sx==0)) return 0;
if ((sx==0)&&(sy> 0)) return 0.5*pi;
if ((sx==0)&&(sy< 0)) return 1.5*pi;
if ((sy==0)&&(sx> 0)) return 0;
if ((sy==0)&&(sx< 0)) return pi;
a=y/x; if (a<0) a=-a;
a=atan(a);
if ((x>0)&&(y>0)) a=a;
if ((x<0)&&(y>0)) a=pi-a;
if ((x<0)&&(y<0)) a=pi+a;
if ((x>0)&&(y<0)) a=pi2-a;
return a;
}
//---------------------------------------------------------------------------
void main(void)
{
float a,b,r,_rr,c;
const float dr=1.0/float(vol_ys); // r step
const float saturation=1000.0; // color saturation voxel value
vec3 rr,p=ray_pos,dp=normalize(ray_dir);
for (c=0.0,r=1.0;r>1e-10;r-=dr) // check all radiuses inwards
{
_rr=1.0/(r*r); rr=vec3(_rr,_rr,_rr);
if (_view_depth(p,dp,rr)) // if ray hits sphere
{
p+=view_depth_l0*dp; // shift ray start position to the hit
a=atanxy(p.x,p.y); // comvert to spherical a,b,r
b=asin(p.z/r);
if (a<0.0) a+=pi2; // correct ranges...
b+=0.5*pi;
a/=pi2;
b/=pi;
// here do your stuff
c=texture(vol_txr,vec3(b,r,a)).r;// fetch voxel
if (c>saturation){ c=saturation; break; }
break;
}
}
c/=saturation;
frag_col=vec4(c,c,c,1.0);
}
//---------------------------------------------------------------------------
它对体积射线追踪器链接进行了轻微修改。
请注意，我假设纹理内的轴是：
latitude,r,longitude
分辨率暗示（经度应该是纬度的双倍分辨率），所以如果它与您的数据不匹配，只需重新排序片段中的轴...我不知道体素单元的值意味着什么，所以我将它们相加为强度/最终颜色的密度，一旦饱和度和达到，光线跟踪就会停止，但你应该计算你想要的东西。
这里预览：
我用的是这个相机矩阵eye for it:
// globals
SphericalVolume3D vol;
// init (GL must be already working)
vol.gl_init();
// render
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glDisable(GL_CULL_FACE);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glTranslatef(0.0,0.0,-2.5);
glGetFloatv(GL_MODELVIEW_MATRIX,vol.eye);
vol.glsl_draw(prog_id);
glFlush();
SwapBuffers(hdc);
// exit (GL must be still working)
vol.gl_init();
射线/球体命中工作正常，球坐标中的命中位置也正常工作，所以唯一剩下的就是轴顺序和颜色算术......