This tutorial will teach you how to utilize MatLab's image processing features to take an MRI scan of a brain with a tumor and isolate the image to show just the tumor as well as give some anatomical details about it. Before starting it is recommended to have MatLab updated as well as some prior basic knowledge in programming or image processing.

Step 1: Finding Proper MRI Images

Sure anyone can take a lame picture off of google images... but you're better than that! While these kinds of scans can be quite difficult to find at times, it is important to make sure that you are using proper images without any watermarks or additional noise added because of the poor image handling from third party websites. There are many online databases that post some images for research purposes that you can use. Some databases include Aylward and OpenfMRI. Here's a link to get you started: http://www.aylward.org/notes/open-access-medical-.... A clear and precise MRI image is the first step to becoming an image processing pro.

Step 2: Converting to Binary Image

In order to properly process your image, you want to first convert it into a binary image. A binary image is simply an image that consists of 1's and 0's that correspond to either black or white pixels. For your image, we are going to set a threshold on our image so that it takes pixel values over a set threshold (you may have to play with number a little yourself. Maybe start off with 150) and converts them to 1's. This way we can immediately get rid of any parts of the image that are too dark to be the tumor.

In matlab, we can do this by entering in the code:

I = imread('image_file_name');

figure(2);
thresha = I>150 & I < 260;

imagesc(thresha);colormap(gray);

This will display the binary image that removes pixels that are either too dim or too bright, in this case it will set pixel values between 150 and 260 to 1 and other values to 0.

Step 3: Threshold Finder (optional)

Although it's not entirely necessary, we felt that making the threshold an automatic process would add a little more accuracy to our code and also allow for dynamic analysis of a larger variety of images. This automatic threshold function relied on analyzing the pixel intensities in each image to obtain a spread of the pixel values. In MatLab, this can be done by simply using a series of if/elseif statements nested inside two for loops to count the number of times each pixel entry occurs. The below code illustrates a snippet of this process.

[iRow, iCol] = size(ImageFile);

for iLoop = 1:iRow
for iLoop0 = 1:iCol

if (I(iLoop,iLoop0) >= 250)

iPixelList(25) = iPixelList(25) + 1;

elseif (I(iLoop,iLoop0) >= 240)

iPixelList(24) = iPixelList(24) + 1;

From subsequent groupings of these intensities values into bins with a width of 10 pixels (as shown above), pixel counts that were too small in value were ignored (1), the largest consecutive group of these bins were isolated using a for loop (2), the index of the minimum and maximum values of the total counts for these pixel entries were then removed (3), and the average was calculated for the remaining counts of the pixel entries (4). A brief snippet of this algorithm is displayed below.

(1)

switch iLoop
case 1

iFilter = 6000; % set pixel count value to be ignored

case 2

iFilter = 5000;

case 3

iFilter = 4000;

end

% loop to remove unwanted pixels

for iLoop = 2:length(iPixelList)

if (iCopyPixelList(iLoop) < iFilter) % check whether is unwanted pixel

iCopyPixelList(iLoop) = 0; % remove the pixel value from list

end

(2)

iMax = 0;
% loop to find the bigest pixel group

for iLoop = 1:length(iPixelList)

if (iPixelGroup(iLoop) > iMax)

iMax = iPixelGroup(iLoop);

iBiggestPixelGroupIndex = iLoop;

end

end;

if (iBiggestPixelGroupIndex == 0) % check whether found biggest pixel group

iThreshold = 0; % assign no threshold

return % nothing found, then exit function

end

(3)

iPixelList(iMinIndex) = 0;

iPixelList(iMaxIndex) = 0;

remainingPixelGroup = iPixelGroup(iBiggestPixelGroupIndex) - 2;

(4)

iAvg = iSum / remainingPixelGroup; % calculate the average of the biggest pixel group

The function then employs three methods to determine the best threshold value. The first method utilizes the pixel intensity value closest to the index of the calculated average as the threshold, the second method uses the index of the minimum value of the total counts for the grouped pixel entries as the threshold, and the third method uses the index of the maximum value of the total counts for the grouped pixel entries as the threshold. The function returns the threshold value that produces the greatest number of pixels to avoid cropping too much of the original image. This is done by the following:

% find the average/min/max threshold which returns the most pixels
if (iAvgThreshold > 0 && iAvgPixelCtr > iMinPixelCtr && iAvgPixelCtr > iMaxPixelCtr)

iThreshold = iAvgThreshold; % set the final threshold

elseif (iMinThreshold > 0 && iMinPixelCtr > iAvgPixelCtr && iMinPixelCtr > iMaxPixelCtr)

iThreshold = iMinThreshold * 10; % set the final threshold

elseif (iMaxThreshold > 0)

iThreshold = iMaxThreshold * 10; % set the final threshold

end

Now you've successfully developed a function to dynamically determine an appropriate threshold value without having to use trial and error for conversion of your MRI to a binary image.

Step 4: Filter Out Noise From the Image

Now we want to get rid of all the extra bits such as the outline of the skull or the some of the other specks that still reside in the image. There are many different kinds of techniques to go about filtering out these extra sections that you can use and it is recommended to use a few to get the best results (more information here: https://www.mathworks.com/help/ images/linear-filtering.html.)

Our first step was to perform a series of morphological functions on the code (ie. filtImage = bwmorph(filtImage, 'close')) to perform some very basic filtering and get rid of some of the easy filters on the image. After this, we used an algorithm to locate all the blobs in the filtered image and sort them by area. Next, we checked to see which blobs in the head had a significantly large enough area to be a tumor. We then removed all the blobs that were too small to be a tumor and were left with the final tumor in the end. (See second image above for example code.) Lastly, we created a masked image to further isolate the intended region of interest by increasing the contrast between the surrounding structures and the tumor region of the MRI scan.

Step 5: Highlighting the Tumor

Our final operations performed on the image allowed us to impose a colored shape around the tumor in the original image in order to highlight it to the viewer. In order to perform this step we used the code:

BWoutline = bwperim(BW); % find perimeter of tumor in binary image
SE3 = strel('square', 3); % use structuring element to increase thickness of outline by a width of 3 pixels

outlineImg = II;

thickOutlines = imdilate(BWoutline, SE3); % increase thickness of outline

outlineImg(thickOutlines) = 50; % set color of outline

subplot(1,4,4);

imshow(outlineImg); % display the outlined tumor on the original MRI

title('Outlined Tumor');

To do this we first used another morphological function to find the outer perimeter of the isolated tumor. Once we attained the perimeter, we set the corresponding pixels in the original image to 50 to change their color and make them brighter. Finally, after this, we displayed the image onto a subplot with the other desired images attained previously.

Step 6: Processing Image and Estimating Anatomical Information

Now that we have the final tumor extracted, we can perform a few functions on the image in order to give a few more specific characteristics about it. For our program, we wanted to calculate the size, area, and location of the tumor in relation to the overall image. In order to do this we were able to utilize some of the built in functions from MatLab to make this process very simple. The regionprops function from matlab enables us to calculate most of the desired information with only a few lines of code (more information here: https://www.mathworks.com/help/images/ref/regionp...) In order to calculate the location and size in relation the image and save it to a text file to save the information, we used this code:

%Saving the results as a text file
FileName = 'MRI Brain Tumor Results';

fid = fopen(FileName, 'wt'); %opens the file

stats2 = regionprops(BW, 'Centroid','Area');

center = stats2.Centroid;

[rows,cols] = find(BW);

Hrow = max(rows); Lrow = min(rows); rowBounds = [Hrow, Lrow];

Hcol = max(cols); Lcol = min(cols); colbounds = [Hcol,Lcol];

After using regionprops to give us the location of the center of the tumor and the amount of pixels it fills, we wanted to find the region that the tumor spanned. To do this we found all the rows and columns which had 1's, then searched for the highest and lowest row and column to give us the total space spanned by the tumor. In order to print all this information to the command window we used:

fprintf(fid,'The center of the tumor is located at coordinates (%0.f , %0.f)\n', center(1),center(2));
fprintf(fid,'The tumor tumor spans from %0.d to %0.d on the x-axis and %0.d and %0.d on the y-axis.\n',colbounds(2),colbounds(1),rowBounds(2),rowBounds(1));

fprintf(fid,'The tumor fills an area of %0.d pixels.\n', stats2.Area)

If you really want to go above and beyond, you can add an extra few lines of code to tell you what percent of the brain your tumor fills. The following code below can be used as a guideline to help you move along:

thresh2 = II >= 45; %new threshold to get rid of the black background
thresh = 1/3*(thresh2(:,:,1)+ thresh2(:,:,2)+ thresh2(:,:,3));

numberExtract = 1;

[labeledImage, numberOfBlobs] = bwlabel(thresh);

stats3 = regionprops(labeledImage, 'Area');

allAreas = [stats3.Area];

[sortedAreas, sortIndices] = sort(allAreas, 'descend');

biggestBlob = ismember(labeledImage, sortIndices(1:numberExtract)); %filter out everything except the head as a whole

binaryImage = biggestBlob > 0;

BW = bwmorph(binaryImage, 'close');

BW = bwmorph(BW, 'fill');

BW = imfill(BW,'holes');

Perc_Tumor = (stats2.Area/stats_head.Area)*100;

In order to compare the tumor to the overall head, we first recreated another binary image, but this time with a very low threshold so that it would convert the entire head structure to a binary image. After this, we then found the biggest blob in the image and removed any other noise in the background (because we know that sometimes you just need to use a watermarked image...) in order to isolate the head. We are then able to find the entire area of the head which we then use the area of the tumor to simply calculate the percentage.

more information of fprintf here: (https://www.mathworks.com/help/matlab/ref/fprintf.html)

Step 7: Conclusion: Final Product and Further Goals That Can Be Achieved

Congratulations! You're now finished with your first (maybe) MatLab project regarding image processing! The data achieved by this code is a starting point in regards to brain tumor diagnostics. While we were able to get most of the information of the tumor regarding its physical presence in the brain, a major step further comes to the actual classification of the tumor as benign or malignant. We can save this for a further project! In the meantime, you can now show off to your friends what an amazing and cool gadget you just built on your computer using nothing but your cunning wits and keen mind, and then watch as they slowly stop texting you and find new and more interesting friends instead. But hey, that just means more free time to go and apply these skills elsewhere! So go try something else with what you've just learned. There are many different projects you can do, and you can process pretty much any image you find, so have fun!

Introduction: Isolating MRI Brain Tumor Using Matlab