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)

289c7d665690c4e72771a11f8b039905.png

其他

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

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

load mri;

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

(包含上述完整代码)

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

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

提示:上述代码数据过大,要运行很久!记得自己优化。

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

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

  • 0
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值