matlab三维地质体绘制,（医学三维重建）MATLAB体绘制算法：光线投射（RC）

(医学三维重建)MATLAB体绘制算法：光线投射(RC)

算法原理

代码实现

实验结果

在这里插入图片描述

其他

by HPC_ZY

算法原理

其原理与体绘制算法MIP类似，只是在灰度计算上有区别。

核心就是 “按照透明度对光线穿过的像素值进行累加, 同时更新总不透明度”

(偷个懒，就不上公式了，可参考各种毕业论文……)

代码实现

注：为更好的展示算法原理，方便大家理解，以下代码未做任何优化，基本按照C语言的格式书写。

数据准备

关于model怎么生成，请查看三维体数据的生成

% 不得不说的事：

% 这个model是从医院拿的DICOM数据，然后事先读取进来的。

% 也不一定非得医学图像，你可以用任意你自己准备好的三维数据。

% 再简单一点讲，model就是一个M*M*D的数组，里面很很多数值，就能用来重建。

load model % 我的数据

model = im2double(model); % 模型

[mRows,mCols,mDims] = size(model); % 模型尺寸

mCenter = ([mRows,mCols,mDims]+1)/2; % 模型中心(因为下标从1开始)

参数设置

% 相机成像范围(默认等于尺寸中最大值)

cRows = max([mRows,mCols,mDims]);

cCols = max([mRows,mCols,mDims]);

cDims = max([mRows,mCols,mDims]);

% 相机成像分辨率(输出尺寸)

rRows = 200;

rCols = 200;

rDims = 100;

% 精度(采样间隔)

precison = ([cRows,cCols,cDims]-1)./([rRows,rCols,rDims]-1);

% 相机角度

alpha = 90;

beta = 90;

网格生成

我们有(rRows，rCols)条射线，每条射线采样rDims个点。所以建立一个(rRows，rCols，rDims)的网格，存储每条射线上每个采样点的坐标位置。

% 转为弧度

alpha = alpha/180*pi;

beta = beta/180*pi;

% 初始化网格

mesh.x = zeros(rRows,rCols,rDims);

mesh.y = zeros(rRows,rCols,rDims);

mesh.z = zeros(rRows,rCols,rDims);

% 计算网格

for d = 1:rDims

for r = 1:rRows

for c = 1:rCols

% 赋初值

x = 1+(r-1)*precison(1);

y = 1+(c-1)*precison(2);

z = 1+(d-1)*precison(3);

% {中心旋转}

% 移至原点

x = x-mCenter(1);

y = y-mCenter(2);

z = z-mCenter(3);

% x轴逆时针旋转

tmp = [x,y,z]; % 避免旋转中覆盖原值

y = tmp(2)*cos(alpha)-tmp(3)*sin(alpha);

z = tmp(2)*sin(alpha)+tmp(3)*cos(alpha);

% y轴逆时针旋转

tmp = [x,y,z]; % 避免旋转中覆盖原值

x = tmp(1)*cos(beta)+tmp(3)*sin(beta);

z = -tmp(1)*sin(beta)+tmp(3)*cos(beta);

% 移回中心

mesh.x(r,c,d) = x+mCenter(1);

mesh.y(r,c,d) = y+mCenter(2);

mesh.z(r,c,d) = z+mCenter(3);

end

end

end

注:由于xy轴配合旋转，就能实现z轴旋转的效果。所以代码中省去z轴旋转，节省计算时间。

4. 数据采样

V = zeros(rRows,rCols,rDims);

% 注：为演示原理没有使用interp3()函数

for d = 1:rDims

for r = 1:rRows

for c = 1:rCols

% 此处使用最近邻插值

x = round(mesh.x(r,c,d));

y = round(mesh.y(r,c,d));

z = round(mesh.z(r,c,d));

if x>=1&&x<=mRows && y>=1&&y<=mCols && z>=1&&z<=mDims

V(r,c,d) = model(x,y,z);

end

end

end

end

RC绘制

im = zeros(rRows,rCols);

for r = 1:rRows

for c = 1:rCols

a = 1;

d = 1;

while a>0.01 && d

v = V(r,c,d);

atmp = V(r,c,d);

im(r,c) = im(r,c)+a*(atmp*v); % 灰度按权累加

a = a*(1 - atmp); % 更新累积透明度

d = d+1;

end

end

end

imshow(im)

实验结果

模型尺寸为512x512x400

左上(rRows = 256，rCols = 256，rDims = 200，alpha = 0，beta = 0)

右上(rRows = 256，rCols = 256，rDims = 100，alpha = 0，beta = 0)

左下(rRows = 128，rCols = 128，rDims = 100，alpha = 0，beta = 0)

右下(rRows = 256，rCols = 256，rDims = 200，alpha = 90，beta = 90)

其他

以上代码在网格、采样、绘制上都有很大的优化空间，可以大大减少内存消耗和时间成本。

关于我文中的model，是利用手里的dicom数据生成的，具体方法可参考文献一，文献二。如果你还没有dicom数据，可以用MATLAB自带的核磁图像进行测试，调用方法如下：

load mri;

有任何问题欢迎讨论，最后还是把测试代码上传

(包含上述完整代码)

https://download.csdn.net/download/xsz591541060/11388472

如果你买了年VIP，请果断下载。

提示：上述代码数据过大，要运行很久！记得自己优化。

基于MATLAB的优化可参考另一篇文章结尾

https://blog.csdn.net/xsz591541060/article/details/92792627