%?（醫學三維重建）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