医学成像原理实验报告

实验一DICOM图像的读取和显示
DICOM（Digital Imaging and Communications in Medicine）即医学数字成像和通信，是医学图像和相关信息的国际标准。

它定义了质量能满足临床需要的可用于数据交换的医学图像格式，利用不同的灰度值实现成像。

DICOM被广泛应用于放射医疗，心血管成像以及放射诊疗诊断设备（X射线，CT，核磁共振，超声等），并且在眼科和牙科等其它医学领域得到越来越深入广泛的应用。

当前大约有百亿级符合DICOM标准的医学图像用于临床使用。

I = dicomread('CT-MONO2-16-ankle.dcm');
info = dicominfo('CT-MONO2-16-ankle.dcm');
I = dicomread(info);
imshow(I,'DisplayRange',[]);
dicomwrite(I,'ankle.dcm');
info = dicominfo('CT-MONO2-16-ankle.dcm');
I = dicomread(info);
dicomwrite(I,'ankle.dcm',info);
（2）图像读取
程序如下：
I = dicomread('CT-MONO2-16-ankle.dcm');
imtool(I,'DisplayRange',[])
info = dicominfo('CT-MONO2-16-ankle.dcm');
info.SeriesInstanceUID
max(I(:))
min(I(:))
Imodified = I;
Imodified(Imodified == 4080) = 32;
imshow(Imodified,[])
20
40
60
80
100
120
20
40
60
80
100
120
20
40
60
80
100
120
20
40
60
80
100
120
20
40
6080100120
20
40
6080
100
120
20406080100120
20
40
60
80
100
120
实验二 MRI 图像显示和读取
MRI 可获得人体横面、冠状面、矢状面及任何方向断面的图像，实现三维定位图像。

完整程序如下：
load mri ; D = squeeze(D); figure('Colormap',map) image_num = 8;
image(D(:,:,image_num)) axis image ; x = xlim; y = ylim
Image num=8
Image num=2 Image num=12 Imagenum=18
20
40
60
80
100
120
2-D Contour Slice
figure('Colormap',cm)
contourslice(D,[],[],image_num) axis ij xlim(x) ylim(y)
daspect([1,1,1])
Displaying 3-D Contour Slices
igure('Colormap',cm)
contourslice(D,[],[],[1,12,19,27],8); view(3); axis tight
Applying an Isosurface to the MRI Data
figure('Colormap',map) Ds = smooth3(D);
hiso = patch(isosurface(Ds,5),... FaceColor',[1,.75,.65],... EdgeColor','none'); isonormals(Ds,hiso)
hcap = patch(isocaps(D,5),... FaceColor','interp',... EdgeColor','none'); view(35,30) ; axis tight ;
20
40
60
80
100
120
daspect([1,1,.4]); lightangle(45,30);
set(gcf,'Renderer','zbuffer'); lighting phong set(hcap,'AmbientStrength',.6)
set(hiso,'SpecularColorReflectance',0,'SpecularExponent',50)
程序：
figure('Colormap',map)
Ds = smooth3(D);
hiso = patch(isosurface(Ds,5),... 'FaceColor',[1,.75,.65],... 'EdgeColor','none');
isonormals(Ds,hiso)
实验三平行束投影仿真
实验原理:sheep-logan 模型用来模拟头部断层图像，通过得到sheep-logan 头型，得到不同方向上投影数据，利用模拟X 光平行束重建图像。

程序如下：
clc; clear all ; close all ; N=256; I=phantom(N); figure;
projection data
20
40
60
80
100
120
140
160
180
50100
150200250
300350
10
20
30
40
50
60
imshow(I) figure;
imshow(I,[0.8 1.0])
clc; clear all ; close all ; N=256; I=phantom(N); theta=0:179; P=radon(I,theta); figure; imshow(I);
title('original head model'); figure;
imagesc(P),colormap(gray),colorbar title('projection data');
实验四 利用Radon 函数直接反投影重建图像
程序如下：
clc;
clear all;
close all;
N=256;
I=phantom(N);
delta=pi/180;
theta=0:1:179;
theta_num=length(theta);
P=radon(I,theta);
[mm,nn]=size(P);
e=floor((mm-N-1)/2+1)+1;
P=P(e: N+e-1,:);
P1=reshape(P,N,theta_num);
rec=medfuncBackprojection(theta_num,N,P1,delta);
figure;
imshow(I,[]);
figure;
imshow(rec,[]);
function rec=medfuncBackprojection(theta_num,N,R1,delta)
rec=zeros(N);
for m=1:theta_num
pm=R1(:,m);
Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
for k1=1:N
for k2=1:N
Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
n=floor(Xrm);
t=Xrm-floor(Xrm);
n=max(1,n);
n=min(n,N-1);
p=(1-t)*pm(n)+t*pm(n+1);
rec(N+1-k1,k2) = rec(N+1-k1,k2)+p;
end
end
end
实验五滤波反投影算法重建实验
实验原理:
目前CT图像重建算法多采用滤波反投影算法.利用滤波反投影算法的基本原理,对R-L,S-L滤波函数分别进行了计算机仿真对比实验.即先对图像进行反投影，再利用滤波函数的卷积累加求和实现图像重建。

实验结果表明利用滤波反投影较好地重建图像,关键是滤波函数的选择.
主程序：
%filtered backprojection reconstruction
clc;
clear all;
close all;
N=256;
I=phantom(N);
delta=pi/180;
theta=0:1:179;
theta_num=length(theta);
d=1;
P=radon(I,theta);
[mm,nn]=size(P);
e=floor((mm-N-1)/2+1)+1;
P=P(e: N+e-1,:);
P1=reshape(P,N,theta_num);
fh_RL=medfuncRLfilterfunction(N,d);
rec= medfuncRLfilteredbackprojrction(theta_num,N,P1,delta,fh_RL);
figure;
imshow(I,[]);
figure;
imshow(rec,[]);
function fh_RL=medfuncRLfilterfunction(N,d)
fh_RL=zeros(1,N)
for k1=1:N
fh_RL(k1)=-1/(pi*pi*((k1-N/2-1)*d)^2)
if mod(k1-N/2-1,2)==0
fh_RL(k1)=0;
end
end
fh_RL(N/2+1)=1/(4*d^2);
function rec_RL = medfuncRLfilteredbackprojrction(theta_num,N,R1,delta,fh_RL) rec_RL=zeros(N);
for m=1:theta_num
pm=R1(:,m);
pm_RL=conv(fh_RL,pm,'same');
Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
for k1=1:N
for k2=1:N
Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
n=floor(Xrm);
t=Xrm-floor(Xrm);
n=max(1,n);
n=min(n,N-1);
m p_RL = (1-t)*pm_RL(n)+t*pm_RL(n+1);
rec_RL(N+1-k1,k2)=rec_RL(N+1-k1,k2)+p_RL;
end
end
end
end
R—L函数滤波反投影图像S—L滤波反投影成像。