threshold_segmentation.m

Simple threshold segmentation of medical CT data in DICOM data representation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%
%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all

[ filesource, pathname ] = uigetfile( '.ima', 'Wybierz obraz do Wyswietlenia' );  

X = dicomread(filesource);


figure(1)
subplot (2,2,1)
imagesc(X); title( 'Obraz oryginalny');
colormap(gray)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  filtracja medianowa
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X_bin = X;
X_med = medfilt2(X, [4,4]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


X= X_med;


[M N]=size(X);
[R C]=find((X>=1095));         % binaryzacja
I_bw=zeros(M,N);
for i=1:numel(R)
    I_bw(R(i),C(i))=1;                    % binaryzacja
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  skasowanie gantry
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=370:470
    for k=1:M
            X(i,k) = 0;
   
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  wysegmentowanie obszaru do transformacji dystansowej ??
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max(X);
for i=1:M
    for k=1:N
        if I_bw(i,k) == 1
            X(i,k) = X(i,k);
        else
            X(i,k) = 0;
    end
    end
end

subplot (2,2,2)
imagesc(X); title( 'Obraz po nalozeniu maski pierwszego progowania');
colormap(gray)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X = medfilt2(X, [8,8]);
%X = medfilt2(X, [8,8]);





[M N]=size(X);
[R C]=find((X>=1095)); % & (X<1600));        % binaryzacja
I_bw=zeros(M,N);
for i=1:numel(R)
    I_bw(R(i),C(i))=1;                    % binaryzacja
end;


subplot (2,2,3)
imagesc(I_bw); title( 'Obraz wysegmentowany');
colormap(gray)


se = strel('disk',2);
IM2 =  bwareaopen(I_bw,6);
IM2 =  imerode(IM2,se);

se = strel('disk',2);
IM2 = imclose(IM2,se);
% se_2 = strel ('disk', 1);
% IM2 = imerode(IM2,se_2);
subplot (2,2,4)
imagesc(IM2); title( 'Obraz po operacjach morfologicznych');
%colormap(gray)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

D = bwdist(~X);
%figure, imshow(D,[],'n'), title('Distance transform of ~bw')



L = watershed(D);
figure, imshow(L)
rgb = label2rgb(L,'jet'); % ,[.5 .5 .5]);
%figure, imshow(rgb), title('Watershed transform of D');
I = rgb2gray(rgb);

[M N]=size(I);
[R C]=find((I>10)); % & (X<1600));        % binaryzacja

I_bw_3=zeros(M,N);
for i=1:numel(R)
    I_bw_3(R(i),C(i))=1;                    % binaryzacja
end;



[M N]=size(D);
[R C]=find((D>2)); % & (X<1600));        % binaryzacja

I_bw_2=zeros(M,N);
for i=1:numel(R)
    I_bw_2(R(i),C(i))=1;                    % binaryzacja
end;

I_bw_2 =  bwareaopen(I_bw_2,8);
%figure, imshow(I_bw_2); title( 'Obraz po progowaniu po transformacji dystansowej');
%figure, imshow(I_bw & I_bw_2);

D_1 = bwdist(~IM2);
figure, imshow(D_1,[],'n'), title('Distance transform of ~bw')

D_1 = -D_1;
D_1(~I_bw) = -Inf;



L_1 = watershed(D_1);
rgb_1 = label2rgb(L_1,'jet',[.5 .5 .5]);
figure, imshow(rgb_1,'n'), title('Watershed transform of D');