医学影像分割

1.迭代阈值法

%迭代阈值法
clc;clear all;close all
I=imread('spine.tif');
% I=dicomread('CT-MONO2-16-ankle.dcm');
% I=rgb2gray(I);
subplot(1,2,1);imshow(I);title('(a) 原图');
I=double(I);T=(min(I(:))+max(I(:)))/2;
d=false;i=0;
while ~ d
    u1=find(I<=T);
    u2=find(I>T);
    Tn=(mean(I(u1))+mean(I(u2)))/2;
    d=abs(Tn-T)<1;
    T=Tn;
    i=i+1;
end
I(u1)=0;I(u2)=1;
subplot(1,2,2);imshow(I);title('(b) 处理结果');

2.分水岭法

%分水岭法
clc;clear all;close all
I=imread('mri.tif');
%I=dicomread('C:\MR-MONO2-12-angio-an1\MR-MONO2-12-angio-an1');
%I=dicomread('C:\MR-MONO2-16-knee\MR-MONO2-16-knee');
% I=rgb2gray(I);
subplot(2,2,1);imshow(I);title('(a) 原图');
I=double(I);hv=fspecial('prewitt');hh=hv.'
gv=abs(imfilter(I,hv,'replicate'));
gh=abs(imfilter(I,hh,'replicate'));
g=sqrt(gv.^2+gh.^2);
L=watershed(g);wr=L==0;
subplot(2,2,2);imshow(wr);title('(b) 分水岭');
I(wr)=255;
subplot(2,2,3);imshow(uint8(I));title('(c) 分割结果');
rm=imregionalmin(g);
subplot(2,2,4);imshow(rm);title('(d) 局部极小值');

3.边缘探测法

%边缘探测法
clc;clear all;close all
% I=imread('mri.tif');
I=dicomread('C:\MR-MONO2-12-angio-an1\MR-MONO2-12-angio-an1');
% I=dicomread('C:\MR-MONO2-16-knee\MR-MONO2-16-knee');
% I=rgb2gray(I);
figure;subplot(2,3,1);imshow(I);title('(a) 原图');
I1=im2bw(I);
subplot(2,3,2);imshow(I1);title('(b) 二值图像');
I2=edge(I,'roberts');
subplot(2,3,3);imshow(I2);title('(c) roberts算子分割');
I3=edge(I,'sobel');
subplot(2,3,4);imshow(I3);title('(d) sobel算子分割');
I4=edge(I,'Prewitt');
subplot(2,3,5);imshow(I4);title('(e) Prewitt算子分割');
subplot(2,3,6);imhist(I);title('(f) 直方图');

4.区域生长法

%区域生长法
clc;clear all;close all
%I=imread('mri.tif');
I=dicomread('C:\MR-MONO2-12-angio-an1\MR-MONO2-12-angio-an1');
%I=dicomread('C:\MR-MONO2-16-knee\MR-MONO2-16-knee');
% I=rgb2gray(I);
figure;
% seedx=[256, 128, 480];
% seedy=[128, 256, 384];
seedx=[128, 64, 255];
seedy=[64, 128, 192];
hold on
plot(seedx,seedy,'gs','linewidth',1);title('(a) 原图图像及种子质量');
figure;subplot(1,2,1);imshow(I);title('(a) 原图');
% seedx=[256, 128, 480];
% seedy=[128, 256, 384];
% hold on
% plot(seedx,seedy,'gs','linewidth',1);
% title('(a) 原图图像及种子质量');
I=double(I);
markerim=I==I(seedy(1),seedx(1));
for i=2:length(seedx)
    markerim=markerim | (I==I(seedy(i),seedx(i)));
end
thresh=[15, 10, 15];
maskim=zeros(size(I));
for i=1:length(seedx)
    g=abs(I-I(seedy(i), seedx(i)))<=thresh(i);
    maskim=maskim | g;
end
[g, nr]=bwlabel(imreconstruct(markerim, maskim),8);
g=mat2gray(g);
subplot(1,2,2);imshow(g)title('(b) 三个种子点区域生长结果');

5.直方图阈值分割法

% 直方图阈值分割法
clc;clear all;close all
I=imread('cell.tif');
% I=rgb2gray(I);
figure;subplot(1,3,1);imshow(I);title('(a) 原图');
subplot(1,3,2);imhist(I);title('(b) 直方图');
[m, n]=size(I);
for i=1:m
    for j=1:n
        if I(i,j)>140    %阈值根据实际情况设置
            I(i,j)=255;
        end
    end
end
subplot(1,3,3);imshow(I);title('(c) 分割结果');

6.Otsu分割法

% Otsu分割法
clc;clear all;close all
I=imread('cell.tif');
% I=rgb2gray(I);
figure;subplot(1,2,1);imshow(I);title('(a) 原图');
count=imhist(I);
[r,t]=size(I);
N=r*t;
L=256;
count=count/N;
for i=2:L
    if count(i)~=0
        st=i-1;
        break
    end
end
for i=L:-1:1
    if count(i)~=0
        nd=i-1;
        break
    end
end
f=count(st+1:nd+1);
p=st;
q=nd-st;
u=0;
for i=1:q
    u=u+f(i)*(p+i-1);
    ua(i)=u;
end
for i=1:q
    w(i)=sum(f(1:i));
end
d=(u*w-ua).^2/(w.*(1-w));
[y,tp]=max(d);
th=tp+p;
if th<=140
    th=tp+p;
else
    th=140;
end
y1=zeros(r,t);
for i=1:r
    for j=1:t
        x1(i,j)=double(I(i,j));
    end
end
for i=1:r
    for j=1:t
        if x1(i,j)>=th
            y1(i,j)=x1(i,j);
        else
            y1(i,j)=0;
        end
    end
end
subplot(1,2,2);imshow(y1);title('(b) 分割结果');