Friday, June 24, 2016

matlab code to identify organs of a CT image

Contents

this matlab code provides some techniques and methods to view and identify organs of a CT image

The first method is to use histogram to segment the second method is to use color a CT image third method is filter.
% read a CT image
I=imread('./Abdominal wall normal anat (8).png');

plot histogram of the image

subplot(311); % plot to ensure a correct image is loaded.
imhist(I)

subplot(3,1,2:3);
imshow(I);
colorbar


perform the color map of a CT image

colormap(jet);
caxis
caxis([0,255]);
ans =

     0   255


colormap( [repmat([0,0,0],[50,1]); ... % 0-50
          repmat([0,1,0],[50,1]);...  % 51-100
          repmat([1,0,0],[75,1]);...  % 100-175
          repmat([0,0,1],[50,1]);...  % 176-225
          repmat([1,1,1],[56-25,1])]);
colorbar

HSV space

convert to single-precision floating point
I = single(I);
[ny,nx] = size(I);
Ihsv = zeros(ny,nx,3);
Ihsv(:,:,2) = 0.5;
Ihsv(:,:,3) = I/max(I(:));

% add color
th1 = 50;
th2 = 100;
th3 = 170;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I<th1) * 0;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I>th1 & I<th2) * 0.1;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I>th2 & I<th3) * 0.2;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I>th3 & I<225) * 0.5;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I>225) * 0.8;

imshow(hsv2rgb(Ihsv));



filtering demo on the CT data

clf;
imagesc(I);colormap(gray);
title('unfiltered');


average filter by hand

n = 5;
b = ones(n,n)/ (n*n);
I2 = conv2(I,b,'same');
imagesc(I2);


high-pass by hand

b = [-1, -1, -1; -1, 7.9866, -1; -1, -1, -1];
I2 = conv2(I,b,'same');
imagesc(I2);
colorbar
mean(I2(:))
caxis([-200,300])
ans =

   -0.0026



average filter

I=single(I);
b = fspecial('average',[9,9]);
I2 = conv2(I,b,'same');
imagesc(I2);


gaussian filter

note: pick filter size appropriately imagausfilt avoids this problem
I=single(I);
b = fspecial('gaussian',32,2);
I2 = conv2(I,b,'same');
imagesc(I2);


unsharp mask

b = fspecial('unsharp');
subplot(221);
imagesc(b);
colorbar;
subplot(222);
imagesc(I);
subplot(224);
I2 = conv2(I,b,'same');
imagesc(I2);


laplacian

clf
b = fspecial('laplacian');
I2 = conv2(I,b,'same');
imagesc(I2);


laplacian of gaussian

b = fspecial('log');
I2 = conv2(I,b,'same');
imagesc(I2);
colorbar


threshold

imagesc(abs(I2)>60);


prewitt edge filter

b = fspecial('prewitt');
I2 = conv2(I,b.','same');
imagesc(abs(I2));
colorbar
caxis([0,300])


imagesc(abs(I2)>150)


sobel edge filter

b = fspecial('sobel');
I2 = conv2(I,b,'same');
imagesc(I2);
colorbar


imagesc(abs(I2)>200)


standard deviation edge filter

I2 = nlfilter(I, [3,3], 'std2');
imagesc(I2);
colorbar;


imagesc(I2>10);


skew filter

f = @(X) skewness(X(:));
I2 = nlfilter(I, [3,3], f);
imagesc(I2);
colorbar;


kurtosis filter

f = @(X) kurtosis(X(:));
I2 = nlfilter(I, [16,16], f);
imagesc(I2);
colorbar;


Matlab code for comparing two image files

Contents

Matlab Code:

Author: Minh Anh Nguyen (minhanhnguyen@q.com)
%This script will perform the followign tasks:
%1. read in two ".png" file, display both images on the computer screen.
%2. histogram of the images
% 3. compare both image files.
% 4. display a message on the computer screen if both image are not the
% same.


clear all; % Erase all existing variables.
close all; % Close all figures (except those of imtool.)
clc;% Clear the command window.
imtool close all; % Close all imtool figures.
workspace; % Make sure the workspace panel is showing.
fontSize = 20;


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

team: edit your file here

set the path

folder = 'C:\Users\Minh anh\Desktop\Image_demo';

%reading images as array to variable 'a' & 'b'.
image1 = imread('.\Loaded-Mid Channel.PNG');
% image2= imread('.\Loaded-Mid Channel.PNG');

image2= imread('.\Gross Misloaded-Mid Channel.PNG');
% image2= imread('.\Mild Misloaded-Mid Channel.PNG');

info

imfinfo('.\Loaded-Mid Channel.PNG')
imfinfo('.\Gross Misloaded-Mid Channel.PNG')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do not edit anything below here!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a = image1;
b = image2;

size

 a1= size(a);
 b1= size (b);
    figure;
% [rows columns numberOfColorBands] = size(a);
subplot(2, 2, 1); imshow(a, []);
title('Correct loading image', 'Fontsize', fontSize);
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
redPlane = a(:, :, 1);
greenPlane = a(:, :, 2);
bluePlane = a(:, :, 3);


% [rows columns numberOfColorBands] = size(b);
subplot(2, 2, 2); imshow(b, []);
title('User loading image', 'Fontsize', fontSize);
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
redPlane1 = b(:, :, 1);
greenPlane1 = b(:, :, 2);
bluePlane1 = b(:, :, 3);

% Let's get its histograms.
[pixelCountR, grayLevelsR] = imhist(redPlane);
subplot(2, 2, 3); plot(pixelCountR, 'r');
title('Histogram of correct loading image', 'Fontsize', fontSize);
% xlim([0 grayLevelsR(end)]); % Scale x axis manually.
xlim([0,256]); ylim([0,82000]);
xlabel('Intensity values','Fontsize', fontSize);
ylabel('Number of pixels','Fontsize', fontSize);

indexmin1 = find(min(pixelCountR) == pixelCountR);
xmin = grayLevelsR(indexmin1);
ymin = pixelCountR(indexmin1);
% strmax2a = ['Minimum = ',num2str(xmin2),' ','    ', num2str(ymin2),' '];
% text(xmin2,ymin2,strmax2a,'HorizontalAlignment','Right');



indexmax1 = find(max(pixelCountR) == pixelCountR);
xmax1 = grayLevelsR(indexmax1);
ymax1 = pixelCountR(indexmax1);
strmax1 = ['Number of intensity =',num2str(xmax1),'  ','Number of pixels = ', num2str(ymax1),''];
text(xmax1,ymax1,strmax1,'HorizontalAlignment','Left');





[pixelCountR2, grayLevelsR1] = imhist(redPlane1);
subplot(2, 2, 4);
plot(pixelCountR2, 'r');
title('Histogram of user loading image', 'Fontsize', fontSize);
% xlim([0 grayLevelsR1(end)]); % Scale x axis manually.
 xlim([0,256]); ylim([0,82000]);
 xlabel('Intensity values','Fontsize', fontSize);
ylabel('Number of pixels','Fontsize', fontSize);
indexmin2 = find(min(pixelCountR2) == pixelCountR2);
xmin2 = grayLevelsR1(indexmin2);
ymin2 = pixelCountR2(indexmin2);
% strmax2a = ['Minimum = ',num2str(xmin2),' ','    ', num2str(ymin2),' '];
% text(xmin2,ymin2,strmax2a,'HorizontalAlignment','Right');



indexmax2 = find(max(pixelCountR2) == pixelCountR2);
xmax2 = grayLevelsR1(indexmax2);
ymax2 = pixelCountR2(indexmax2);
strmax2 = ['Number of intensity = ',num2str(xmax2),'    ','Number of pixels =', num2str(ymax2),''];
text(xmax2,ymax2,strmax2,'HorizontalAlignment','Left');


%check image for different
 different = a1- b1;
 if different==0
     disp('The images are same')%output display
 else
    disp('the images are not same')
    msgbox('Centrifuge loading error. Please check channel is seated properly')
    end;
the images are not same