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

文章：（医学三维重建）MATLAB体绘制算法：光线投射（RC）的demo

% （医学三维重建）MATLAB体绘制算法：光线投射（RC）
%  demo
%
% by HPC_ZY 20190719
clear; close all; clc

load model % 我的数据

tic
%% 数据准备
model = im2double（model）; % 模型
[mRowsmColsmDims] = size（model）; % 模型尺寸
mCenter = （[mRowsmColsmDims]+1）/2; % 模型中心（因为下标从1开始）


%% 设置相机参数
% 相机成像范围（默认等于尺寸中最大值）
cRows = max（[mRowsmColsmDims]）;
cCols = max（[mRowsmColsmDims]）;
cDims = max（[mRowsmColsmDims]）;
% 相机成像分辨率（输出尺寸）
rRows = 256;
rCols = 256;
rDims = 200;
% 精度（采样间隔）
precison = （[cRowscColscDims]-1）./（[rRowsrColsrDims]-1）;
% 相机角度
alpha = 90;
beta = 90;


%% 网格生成
% 转为弧度
alpha = alpha/180*pi;
beta = beta/180*pi;
% 初始化网格
mesh.x = zeros（rRowsrColsrDims）;
mesh.y = zeros（rRowsrColsrDims）;
mesh.z = zeros（rRowsrColsrDims）;
% 计算网格
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 = [xyz]; % 避免旋转中覆盖原值
            y = tmp（2）*cos（alpha）-tmp（3）*sin（alpha）;
            z = tmp（2）*sin（alpha）+tmp（3）*cos（alpha）;
            % y轴逆时针旋转
            tmp = [xyz]; % 避免旋转中覆盖原值
            x = tmp（1）*cos（beta）+tmp（3）*sin（beta）;
            z = -tmp（1）*sin（beta）+tmp（3）*cos（beta）;
            % 移回中心
            mesh.x（rcd） = x+mCenter（1）;
            mesh.y（rcd） = y+mCenter（2）;
            mesh.z（rcd） = z+mCenter（3）;
        end
    end
end


%% 数据采样
V = zeros（rRowsrColsrDims）;
% 注：为演示原理没有使用interp3（）函数
for d = 1:rDims
    for r = 1:rRows
        for c = 1:rCols
            % 此处使用最近邻插值
            x = round（mesh.x（rcd））;
            y = round（mesh.y（rcd））;
            z = round（mesh.z（rcd））;
            if x>=1&&x<=mRows && y>=1&&y<=mCols && z>=1&&z<=mDims
                V（rcd） = model（xyz）;
            end
        end
    end
end


%% RC绘制
im = zeros（rRowsrCols）;
for r = 1:rRows
    for c = 1:rCols
        a = 1;
        d = 1;
        while a>0.01 && d            v = V（rcd）;
            atmp = V（rcd）;
            im（rc） = im（rc）+a*（atmp*v）; % 灰度按权累加
            a = a*（1 - atmp）; % 更新累积透明度
            d = d+1;
        end
    end
end

toc

imshow（im）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2585  2019-07-19 16:07  demo.m

----------- ---------  ---------- -----  ----

                 2585                    1