Matlab 沿三維任意方向切割CT圖的仿真計算

1、數據來源

1. 頭部組織的數據。此處直接引用了matlab自帶的mri數據。實際場景中，能夠經過CT獲得的數據進行轉換獲得
2. 插入異物的數據。此處我假設插入異物爲一根細鐵絲。模擬爲空間中的一條曲線。這個曲線的座標咱們能夠經過必定的辦法獲取到。我已經將數據放入百度雲盤 曲線數據下載

2、技術實現要點

1. 要根據CT數據繪製頭部的立體圖
2. 異物曲線的任意兩點連線的切面的計算
3. 切面切割CT立體圖的效果圖繪製

3、實現僞代碼

%% 數據準備
clc,clear
load mri;
D = squeeze(D);
D=double(D);
% img = importdata('image00.mat');
% c_start=1;
% D=img(:,:,:);
[m,n,q]=size(D);
%% 準備曲線數據
Curve = importdata('reference0.txt');  
dataS=1000;                         %這兩個參數 是我處理數據用的，數據太多，截取了一段曲線
dataLen=2000;
k=1; % 控制係數，讓血管和腦補靠近
x_cur = (Curve(dataS:dataLen,1)/k);
y_cur = (Curve(dataS:dataLen,2)/k);
z_cur =-35+Curve(dataS:dataLen,3)/k;
figure
plot3(x_cur,y_cur,z_cur,'r*','linewidth',2)
title('單獨曲線圖')
grid on
xlabel('x');
ylabel('y');
zlabel('z');
grid on;
%爲了畫圖方便
%畫圖
%另外還常常用到點法式方程：A(x-x0)+B(y-y0)+C(z-z0)=0
% x=x0+kt  y=y0+mt  z=z0+nt（x0，y0，z0）表示直線通過的一個點，t爲任意實數，
% N_num=[100,200,300,400,500,600,700,800,900,1000];
% i=5;
% 而向量（k，m，n）表示直線的方向。
N=200;         % 第多少個點；最好設置爲100 200 300 ...1000.此處N 爲控制切片位置的參數
x_start=x_cur(N-1);
y_start=y_cur(N-1);
z_start=z_cur(N-1);
x_end=x_cur(N+1);
y_end=y_cur(N+1);
z_end=z_cur(N+1);
x_mid=x_cur(N);
y_mid=y_cur(N);
z_mid=z_cur(N);
A=x_end-x_start;
B=y_end-y_start;
C=z_end-z_start;
x0=x_mid;
y0=y_mid;
z0=z_mid;
%%生成與法線垂直的切片平面
[xs, ys] = meshgrid(0:1:128);
zs=-1*(((A*(xs-x_mid)+B*(ys-y_mid))./C)-z_mid);
[xxx,yyy,zzz]=meshgrid(1:m,1:n,1:q);  %構造這個平面
figure
plot3(x_cur,y_cur,z_cur,'r','linewidth',2)
xlabel('x');
ylabel('y');
zlabel('z');
grid on;
hold on;
slice(xxx,yyy,zzz,D,xs,ys,zs);
view(3)
shading interp
hold on

調節N 的大小，能夠獲得一系列切割圖。

4、實現代碼第二部分（此段代碼不注重切割方向）

以下代碼直接和上述代碼寫在一塊兒便可，第二部分代碼有部分爲參考他人的寫法

 　

image_num=8;
image(D(:,:,image_num))%將矩陣D顯示成圖，D中每個元素表明圖像中一個長方形塊的顏色
axis image%等同於axis equal設置寬高比使3個方向數據單位相同
x=xlim;%xlim返回當前x軸的界限
y=ylim;
 
cm=brighten(jet(length(map)),-.5);%使句柄爲jet（length（map））的圖形子對象變得更亮/暗 負爲變暗
%jet, by itself, is the same length as the current figure's colormap. If no figure exists, MATLAB uses the length of the default colormap.
figure('Colormap',cm)
plot3(x_cur,y_cur,z_cur,'r*','linewidth',2)
hold on
contourslice(D,[],[],image_num)%在體積切平面中繪製等高線
axis ij%將座標系的原點放在左上角
xlim(x)
ylim(y)
daspect([1,1,1])%設置軸數據的寬高比，此處設置x:y:z=1:1:1
 
figure('Colormap',cm)
plot3(x_cur,y_cur,z_cur,'r*','linewidth',2)
hold on
contourslice(D,[1,2],[],[1,12,19,27],8);
view(3);%設置視角爲默認的三維視圖
axis tight%設置軸的限度爲數據的範圍

效果圖以下：

立體包絡面展現

figure('Colormap',map)
plot3(x_cur,y_cur,z_cur,'r*','linewidth',2)
hold on
Ds=smooth3(D);%W = smooth3(V) smooths the input data V and returns the smoothed data in W.平滑輸入數據D，輸出Ds 變成double型數據
hiso = patch(isosurface(Ds,5),...%返回patch建立的塊對象句柄；從塊體數據中提取等值表面數據，返回等表面的面和頂點，可直接傳遞給patch；fv = isosurface(V,isovalue) assumes the arrays X, Y, and Z are defined as [X,Y,Z] = meshgrid(1:n,1:m,1:p) where [m,n,p] = size(V)
     'FaceColor',[1,.75,.65],...
     'EdgeColor','none');
isonormals(Ds,hiso)%計算等值表面頂點的法向 isonormals(V,p) and isonormals(X,Y,Z,V,p) set the VertexNormals property of the patch identified by the handle p to the computed normals rather than returning the values.
 
hcap=patch(isocaps(D,5),...%計算帽端等表面幾何 D爲塊體數據 輸出帽端的面、頂點和顏色數據給patch
    'FaceColor','interp',...
    'EdgeColor','none');
 
view(35,30)%兩值設定了視角
axis tight
daspect([1,1,.4])
 
lightangle(45,30);%在球座標系中建立並放置一個發光體，兩值設定了視角lightangle(az,el) creates a light at the position specified by azimuth and elevation. az is the azimuthal (horizontal) rotation and el is the vertical elevation (both in degrees). The interpretation of azimuth and elevation is the same as that of the view command.
set(gcf,'Renderer','zbuffer');
lighting phong%指定馮氏明暗處理算法
set(hcap,'AmbientStrength',.6)
set(hiso,'SpecularColorReflectance',0,'SpecularExponent',50)

效果圖：