How to perform skull stripping using matlab?

I would like to do a project on skull stripping. What is the effective method to perform skull stripping/removal? Could anyone help me in getting code for that?

This is my image sir.

[Matz Johansson Bergström]

Maybe you should change the term "Skull stripping" to a more generic term like "3d slice visualization". You might get more help that way.

[manasa]

Thank you Mr.Johansson

Answer 1

[Image Analyst]

Try this code:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 24;

% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

%===============================================================================
% Read in a standard MATLAB gray scale demo image.
folder = 'D:\Temporary stuff';
baseFileName = 'Jones-54-1-jan10-f3.jpg';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- didn't find it there.  Check the search path for it.
  fullFileNameOnSearchPath = baseFileName; % No path this time.
  if ~exist(fullFileNameOnSearchPath, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.
% numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(grayImage);
if numberOfColorBands > 1
  % It's not really gray scale like we expected - it's color.
  % Convert it to gray scale by taking only the green channel.
  grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the original gray scale image.
subplot(2, 3, 1);
imshow(grayImage, []);
axis on;
title('Original Grayscale Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

% Let's compute and display the histogram.
[pixelCount, grayLevels] = imhist(grayImage);
subplot(2, 3, 2);
bar(grayLevels, pixelCount);
grid on;
title('Histogram of original image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

% Crop image to get rid of light box surrounding the image
grayImage = grayImage(3:end-3, 4:end-4);
% Threshold to create a binary image
binaryImage = grayImage > 20;
% Get rid of small specks of noise
binaryImage = bwareaopen(binaryImage, 10);
% Display the original gray scale image.
subplot(2, 3, 3);
imshow(binaryImage, []);
axis on;
title('Binary Image', 'FontSize', fontSize);

% Seal off the bottom of the head - make the last row white.
binaryImage(end,:) = true;
% Fill the image
binaryImage = imfill(binaryImage, 'holes');
subplot(2, 3, 4);
imshow(binaryImage, []);
axis on;
title('Cleaned Binary Image', 'FontSize', fontSize);

% Erode away 15 layers of pixels.
se = strel('disk', 15, 0);
binaryImage = imerode(binaryImage, se);
subplot(2, 3, 5);
imshow(binaryImage, []);
axis on;
title('Eroded Binary Image', 'FontSize', fontSize);

% Mask the gray image
finalImage = grayImage; % Initialize.
finalImage(~binaryImage) = 0;
subplot(2, 3, 6);
imshow(finalImage, []);
axis on;
title('Skull stripped Image', 'FontSize', fontSize);

msgbox('Done with demo');

[manasa]

Thank sir.I will try this.

[manasa]

Sir, We have done upto background removal..After that we have to draw histogram for only head part to extract brain.How can we draw histogram for only foreground in an image?

[Image Analyst]

headPixels = grayImage(binaryImage); % Head pixels only.
counts = imhist(headPixels);

[Erfan Ta'atizadeh]

Really Thanks! Really!!!

[Maria]

This is incredibly useful & very elegant. Thanks!

Maria

[MMSAAH]

if I change the original Image, should I change the threshold of the histogramme to create the binary image ? and how Can I do it ?

[Image Analyst]

If it's not working for you, you'd change this line:

% Threshold to create a binary image
binaryImage = grayImage > 20;

[MMSAAH]

Hello Image Analyst

Please, how to generate the skull striped image for all slices (j=1 :14). I changed the code as following but it doesn't work. Could you please Help me ?

% Mask the gray image
for j = 1: 14
 finalImage(:,:,j) = grayImage(:,:,j); % Initialize.
 finalImage(~binaryImage(:,:,j)) = 0;
end

[Kowsalya Nandhakumar]

This code is not working for other brain MRI images which are in sagittal view. Any suggestion ??

[Image Analyst]

Not sure why it would only work with certain views. It should work with any view as long as the slice has a skull completely surrounding the brain. Just make the adaptations so that it works with your view.

[Magdalena Gierczynska]

How you choose trashold point? In your code this is 20, why?

[Image Analyst]

I need a threshold that discriminates between body and background (air). Looking at the histogram it seems like 20 would be that threshold so I tried it and it worked.

Answer 2

[Alfonso Nieto-Castanon]

Assuming that you do not have simply that slice but the entire MRI scan (e.g. a filename.nii or filename.img file) I would suggest using SPM (this is Matlab software). Applying segmentation to your structural volume will generate a number of masks (files named c#filename.nii) indicating the different tissues of interest.

To perform skull-stripping you typically want to keep the grey-matter, white-matter, and CSF masks (masks c1filename.nii to c3filename.nii, respectively) and remove everything else. You can do this after segmentation simply using SPM 'imcalc' option and multiplying your original structural volume by the sum of these masks (c1- to c3).

You may try searching through or posting to SPM email list for more info

Answer 3

[Image Analyst]

To strip/remove the skull, look for my code and search for things like skull and brain in this Answers forum. I know I gave a demo for this some time back. A quick check didn't find my demo but I didn't look at a lot of questions on this, and there are a lot of them. I did find one answer I gave and it seems valid for your case.

To strip the skull you threshold it so that you get the skull separated - not connected to the brain. Then you label it with bwlabel() and then use ismember() to isolate the region with label 1, which should be the skull assuming it's the outermost object (no letters or annotation present in the image). Or you can use ismember to extract out everything with labels 2 and up to get everything except the skull.

Attach your image and attempt at coding that up for more help.

[Chandita Nath barbhuyan]

How to remove the white skull region outside?

[Image Analyst]

Adapt my attached skull stripping demo.

Answer 4

[Matz Johansson Bergström]

I found something that might be of help: link to slicing. Or better: MRI slicing.

Answer 5

[Magdalena Gierczynska]

Hi, I have a question. I understand everything to "%Erode away 15 layers of pixels.". You used erosion to reduce border pixels, but what exactly happen in "% Mask the gray image"? On the next image you show "skull stripped image", how? I don't understand this part of code. Can you explain me? I will be greatful.

[Magdalena Gierczynska]

What mean "~" in your code?

[Image Analyst]

In code, ~ means inverse. So if the expression immediately after the ~ is false, using ~ will produce a true value. And vice versa

% Example 1
t = true;
f = ~t;  % f will have a value of false
% Example 2
f = false;
t = ~f; % t will have a value of true
% Example 3
mask = grayImage < 100; % mask is true if pixels have values less than 100
invertedMask = ~mask;   % invertedMask is true if pixels have a value >= 100

[Ayush singhal]

@Image Analyst

Hello, I have a simikar doubt reagarding this.

I have two images. one with stripes and another object with stripes. I would like to extract the object from these two images. These two are uint8 images. Also I want to find the depth of the object by showing images in the surf plot.

Could you give some hint?

I have used imsubtract but it won't give the desire results.

Pprs3 = imread('P1000571-min.JPG'); % Colour Image

Pprs1 = rgb2gray(Pprs3); % Grayscale Image

x = 0:size(Pprs1,2)-1;

y = 0:size(Pprs1,1)-1;

[X,Y] = meshgrid(x,y); % Coordinate Matrices (Not Necessary)

figure(1)

meshc(X, Y, Pprs1) % Mesh Plot

grid on

xlabel('X')

ylabel('Y')

zlabel('Intensity')

colormap(jet)

hold on

a = imread('P1000578-min.JPG'); % Colour Image

b = rgb2gray(a); % Grayscale Image

x = 0:size(b,2)-1;

y = 0:size(b,1)-1;

[X,Y] = meshgrid(x,y); % Coordinate Matrices (Not Necessary)

figure(2)

meshc(X, Y, b) % Mesh Plot

grid on

xlabel('X')

ylabel('Y')

zlabel('Intensity')

colormap(jet)

hold on

c = imsubtract(Pprs1,b);

x = 0:size(c,2)-1;

y = 0:size(c,1)-1;

[X,Y] = meshgrid(x,y); % Coordinate Matrices (Not Necessary)

figure(3)

meshc(X, Y, c) % Mesh Plot

grid on

xlabel('X')

ylabel('Y')

zlabel('Intensity')

colormap(jet)

In the attachment, images are given.

[Image Analyst]

@Ayush singhal, can you please start a new thread and I'll answer it there?

[Ayush singhal]

@Image Analyst

Yes, I have posted it seperately.

Thanks.