一、數據來源
- 頭部組織的數據。此處直接引用了matlab自帶的mri數據。實際場景中，可以通過CT得到的數據進行轉換得到
- 插入異物的數據。此處我假設插入異物為一根細鐵絲。模擬為空間中的一條曲線。這個曲線的坐標我們可以通過一定的辦法獲取到。我已經將數據放入百度雲盤 曲線數據下載
二、技術實現要點
- 要根據CT數據繪制頭部的立體圖
- 異物曲線的任意兩點連線的切面的計算
- 切面切割CT立體圖的效果圖繪制
三、實現偽代碼
%% 數據准備
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 的大小，可以得到一系列切割圖。
四、實現代碼第二部分（此段代碼不注重切割方向）
如下代碼直接和上述代碼寫在一起即可，第二部分代碼有部分為參考他人的寫法
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)
效果圖：
