mtracking.m

运用双向轮廓跟踪方法实现医学图象外部论文提取的算法

a=dicomread('F:\CT数据\数据1\IM10'); %读入CT图像
b=uint16(a);  %数据类型转换
figure;imshow(b,[]);title('初始图像');
c(1:512,1:512)=0;  %消除头架矩阵
for i=70:470
for j=60:450
c(i,j)=1;
end
end
for i=70:100  %除去左上角边缘
for j=60:80
c(i,j)=0;
end
end
for i=70:100 %除去右上角边缘
for j=430:450
c(i,j)=0;
end
end
for i=400:470 %除去左下角边缘
for j=60:120
c(i,j)=0;
end
end
for i=400:470 %除去右下角边缘
for j=390:450
c(i,j)=0;
end
end
for i=350:400 %除去左侧头架
for j=60:90
c(i,j)=0;
end
end
for i=440:470%去除左下侧头架
for j=120:160
c(i,j)=0;
end
end
for i=440:470    %去除右下侧头架
for j=350:390
c(i,j)=0;
end
end
for i=350:400   %去除右侧头架
for j=420:450
c(i,j)=0;
end
end
c1=uint16(c);  %矩阵类型转变
g=immultiply(b,c1);
for i=1:512      %处理图像头部的阴影部分
for j=1:512
if (g(i,j)<760)
g(i,j)=0;
end
end
end
figure;imshow(g,[]);title('初始化后的图像');
g1=double(g);
e=edge(g1,'sobel',-177);
[m n]=size(e); %确定要检测图像的大小；矩阵e是图像所对应的矩阵；
f=zeros(m,n);%设置边缘点矩阵f,
i=1;j=m/2;
while (e(i,j)<1)
i=i+1;
end
m0=i;n0=j;
f(m0,n0)=1;         %初始边缘点
%确定右跟踪的起始点（mr,nr）;
if (e(m0+1,n0+1)==1) %检测邻域4区域中的点；
mr=m0+1;nr=n0+1;
f(mr,nr)=1;
end
if (e(m0,n0+1)==1)   %检测邻域5区域中的点
mr=m0;nr=n0+1;
f(mr,nr)=1;
end
if (e(m0-1,n0+1)==1) %检测邻域6区域中的点
mr=m0-1;nr=n0+1;
f(mr,nr)=1;
end
%确定左跟踪的起始点（ml,nl）;
if (e(m0+1,n0-1)==1)  %检测邻域2区域中的点；
   ml=m0+1;nl=n0-1;
   f(ml,nl)=1;
end
if (e(m0,n0-1)==1)  %检测邻域1区域中的点；
   ml=m0;nl=n0-1;
   f(ml,nl)=1;
end
if (e(m0-1,n0-1)==1)  %检测邻域0区域中的点；
   ml=m0-1;nl=n0-1;
   f(ml,nl)=1;
end
%设置一个标识矩阵d，标识该像素点是否被检查过，如果检查过就不再检查；
d=zeros(m,n);
d(m0,n0)=1;  %初始边缘点已检测过；
d(mr,nr)=1;  %右起始点已被检测过；
d(ml,nl)=1;  %左起始点已被检测过； 
nextl=1;
nextr=1;
while ((nextr<9)&(nextl<9))
t1=zeros(1,9);
t2=zeros(1,9);
t1(1,9)=1;
t2(1,9)=1;
if (d(mr-1,nr-1)==0)  %右跟踪的点检测邻域0区域的点
   d(mr-1,nr-1)=1;
   if (e(mr-1,nr-1)==1)
       f(mr-1,nr-1)=1;
       t1(1,1)=1;
   end 
end   
if (d(ml-1,nl-1)==0)  %左跟踪的点检测邻域0区域的点
   d(ml-1,nl-1)=1;
   if (e(ml-1,nl-1)==1)
       f(ml-1,nl-1)=1;
       t2(1,1)=1;
   end 
end   
if (d(mr,nr-1)==0)  %右跟踪的点检测邻域1区域的点
  d(mr,nr-1)=1; 
  if (e(mr,nr-1)==1)
       f(mr,nr-1)=1;
       t1(1,2)=1
   end 
end 
if (d(ml,nl-1)==0)  %左跟踪的点检测邻域1区域的点
  d(ml,nl-1)=1; 
  if (e(ml,nl-1)==1)
       f(ml,nl-1)=1;
       t2(1,2)=1
   end 
end 
if (d(mr+1,nr-1)==0)  %右跟踪的点检测邻域2区域的点
   d(mr+1,nr-1)=1;
   if (e(mr+1,nr-1)==1)
       f(mr+1,nr-1)=1;
       t1(1,3)=1;
   end 
end   
if (d(ml+1,nl-1)==0)  %左跟踪的点检测邻域2区域的点
   d(ml+1,nl-1)=1;
   if (e(ml+1,nl-1)==1)
       f(ml+1,nl-1)=1;
       t2(1,3)=1;
   end 
end 
if (d(mr+1,nr)==0)  %右跟踪的点检测邻域3区域的点
   d(mr+1,nr)=1;
   if (e(mr+1,nr)==1)
       f(mr+1,nr)=1;
       t1(1,4)=1;
   end 
end  
if (d(ml+1,nl)==0)  %左跟踪的点检测邻域3区域的点
   d(ml+1,nl)=1;
   if (e(ml+1,nl)==1)
       f(ml+1,nl)=1;
       t2(1,4)=1;
   end 
end  
if (d(mr+1,nr+1)==0)  %右跟踪的点检测邻域4区域的点
   d(mr+1,nr+1)=1;
   if (e(mr+1,nr+1)==1)
       f(mr+1,nr+1)=1;
       t1(1,5)=1;
   end 
end 
if (d(ml+1,nl+1)==0)  %左跟踪的点检测邻域4区域的点
   d(ml+1,nl+1)=1;
   if (e(ml+1,nl+1)==1)
       f(ml+1,nl+1)=1;
       t2(1,5)=1;
   end 
end 

if (d(mr,nr+1)==0)  %右跟踪的点检测邻域5区域的点
   d(mr,nr+1)=1;
   if (e(mr,nr+1)==1)
       f(mr,nr+1)=1;
       t1(1,6)=1;
   end 
end  
if (d(ml,nl+1)==0)  %左跟踪的点检测邻域5区域的点
   d(ml,nl+1)=1;
   if (e(ml,nl+1)==1)
       f(ml,nl+1)=1;
       t2(1,6)=1;
   end 
end 
if (d(mr-1,nr+1)==0)  %右跟踪的点检测邻域6区域的点
   d(mr-1,nr+1)=1;
   if (e(mr-1,nr+1)==1)
       f(mr-1,nr+1)=1;
       t1(1,7)=1;
   end 
end  
if (d(ml-1,nl+1)==0)  %左跟踪的点检测邻域6区域的点
   d(ml-1,nl+1)=1;
   if (e(ml-1,nl+1)==1)
       f(ml-1,nl+1)=1;
       t2(1,7)=1;
   end 
end  
if (d(mr-1,nr)==0)  %右跟踪的点检测邻域7区域的点
   d(mr-1,nr)=1;
   if (e(mr-1,nr)==1)
       f(mr-1,nr)=1;
       t1(1,8)=1;
   end 
end  
if (d(ml-1,nl)==0)  %左跟踪的点检测邻域7区域的点
   d(ml-1,nl)=1;
   if (e(ml-1,nl)==1)
       f(ml-1,nl)=1;
       t2(1,8)=1;
   end 
end  
%确定下一个右检测点
i=1;
while (t1(1,i)<1)
    i=i+1;
end;
nextr=i;
if (nextr==1)
    pr=mr-1;qr=nr-1;
end
if (nextr==2)
    pr=mr;qr=nr-1;
end
if (nextr==3)
    pr=mr+1;qr=nr-1;
end
if(nextr==4)
    pr=mr+1;qr=nr;
end
if(nextr==5)
    pr=mr+1;qr=nr+1;
end
if(nextr==6)
    pr=mr;qr=nr+1;
end
if(nextr==7)
    pr=mr-1;qr=nr+1;
end
if(nextr==8)
    pr=mr-1;qr=nr;
end
mr=pr;nr=qr; 
%确定下一个左检测点；
i=1;
while (t2(1,i)<1)
    i=i+1;
end;
nextl=i;
if (nextl==1)
    pl=ml-1;ql=nl-1;
end
if (nextl==2)
    pl=ml;ql=nl-1;
end
if (nextl==3)
    pl=ml+1;ql=nl-1;
end
if(nextl==4)
    pl=ml+1;ql=nl;
end
if(nextl==5)
    pl=ml+1;ql=nl+1;
end
if(nextl==6)
    pl=ml;ql=nl+1;
end
if(nextl==7)
    pl=ml-1;ql=nl+1;
end
if(nextl==8)
    pl=ml-1;ql=nl;
end
ml=pl;nl=ql; 
end
figure;imshow(f,[]);title('皮肤组织轮廓线');