How to Obtain Blood Vessels Segmentation in Retinal Images Using Matlab

Segmentation is the process of Retinal partitioning a digital image into multiple regions and extracting the significant ones. Retinal images are the digital images of the eyeball that show the retina, the optical nerves, and the blood vessels that send information to the brain. <!--more--> This practice helps the optometrist to find diseases and also give a health check on the eyes.

This segmentation process is essential in the field of medicine. When the vessels are segmented and viewed closely, the solution or cause of a given problem may be defined.

Therefore, it makes it a good application in this field. This article will look at how we can obtain the blood vessel segmentation in the retinal image.

Prerequisites #

To follow along with this tutorial, you will need to have:

• MATLAB installed.
• Proper understanding of MATLAB basics.
• Basic understanding of image processing using Matlab.

Matlab code for obtaining the segmentation #

To make this whole process easier, download your image and store it in Matlab's current folder. The retinal images can be downloaded directly from the internet.

Thereafter, we will open Matlab and create a new script. The first step is reading the images.

This is made possible by the imread function:

test_image = imread('retinal_image.jpg');      %Reading the image

Our image looks big at this point. We should therefore resize it using the imresize function. This function uses the image name and the preferred dimensions in vector form as the arguments.

The dimensions are in pixels. This means that decimal dimensions are not accepted:

resized_image = imresize(test_image, [584 565]);     %resizing the image

Since we will be segmenting very tiny blood vessels, we will need to convert the resized image into double data time using the im2double function.

This function only uses the image name as the argument:

converted_image = im2double(resized_image);       %converting the image to double data time

At this point, the image is an RGB image. We need to convert the image to CIE lab colour space. CIE Lab colour space was defined by the international commission of illumination.

Converting the image to CIE lab color space #

According to this colour space, the l represents the likeness of the colour. The likeliness ranges from 0 to 100, where 0 is black and 100 is white.

The a represents green to red, ranging from negative to positive. The higher the negative value, the brighter the green color, while the higher the positive value, the brighter the red colour.

The b represents blue to yellow. It also ranges from negative to positive.

The image will be converted as shown below:

lab_image = rgb2lab(converted_image);    %Converting the image from rgb color space to lab color space

Let us now use the cat function to concatenate 1, 0, and 0. Concatenation is the process of merging two or more variables together:

fill = cat(3, 1, 0, 0);              %cancating the image
filled_image = bsxfun(@times, fill, lab_image);

From the code above, the 3 represents the dimensions of the concatenated areas. Our image is in the CIE Lab colour space, which has 3 channels.

Then, we used the bsx function to perform an element-wise binary operation between the filled and lab images.

Reshaping the output image #

Next, we will reshape the filled image. The dimension arguments will be blank since we do not need any here.

Instead, we will use 3 since it is the existing dimension of the filled image, as shown below:

reshaped_lab_image = reshape(filled_image, [], 3);    %reshaping the image

Now we apply the principal component analysis (PCA) function. Also, we are using the reshaped_lab_image as the argument.

The function returns the coefficients and the score of the principal component. Variables C and S will store coefficients of the principal component.

You can then resize the scores based on the size of the lab image as shown:

[C, S] = pca(reshaped_lab_image);    %finding the coefficients of the pca
S = reshape(S, size(lab_image));

We need all the rows and the columns of the first channel:

S = S(:, :, 1);

Performing contrast enhancement and filter on the grayscale image. #

Now, let us convert the S into a grayscale image. First, we subtract S's minimum value from S and then divide it by the maximum and minimum of value of S.

The division is going to be an element-wise division, and that is why we use the dot before the division sign as shown in the code below:

gray_image = (S-min(S(:)))./(max(S(:)));    %converting image S into a grayscale

We need the contrast enhancement of this gray image. This will be done using the adaptive histogram equalization function adapthisteq, as shown:

enhanced_image = adapthisteq(gray_image, 'numTiles', [8 8], 'nBins', 128);     %enhancing the contrast

The numTiles, which stands for the number of tiles, means the enhancement is done block by block. The size of the block is going to be 8x8.

nBins indicates the number of beams will be 128.

The next step is to perform an average filter on the image. To do it, we define the filter using the special function. We then apply the imfilter function to compute the value of the output image in pixels as shown below:

avg_filter = special('average', [9 9]);    %filtering the image
filtered_image = imfilter(enhanced_image, avg_filter);

We will now view the filtered image and subtract it from the enhanced image using the imsubtract function. To view the image, we will use the imshow function as shown below:

figure, imshow(filtered_image)             %showing the resulting image
title('filtered image')                    %adding the title
subtracted_image = imsubtract(filtered_image, enhanced_image);

Calculating the threshold level #

Now we need to calculate the threshold level to segment the blood vessels. Let's create a function script named threshold_level.m for this.

This function will take the image as the argument as shown below:

function level = Threshold_level(image)

We then convert the image into uint8. This is a data type of 8-bits. This conversion is done using the im2uint8 function.

We then use the imhist function and the converted image as the argument to get the histogram count and beam number as shown:

image = im2uint8(image(:));     %converting the image to uint8
[Histogram_count, Bin_number] = imhist(image);    %calculating the histogram count and beam number
i = 1;     % Initializing the variable

Calculate the cumulative sum of histogram count using cumsum function as shown below:

cumulative_sum = cumsum(Histogram_count);    %calculating the cumulative sum

Let's now find the mean below and above T. T is the ratio of the sum of the multiplication of bin number and the histogram count to the cumulative sum indexed at the end.

T(i) = (sum(Bin_number.*Histogram_count))/cumulative_sum(end);  %calculating the ratio of the sum of the multiplication of bin number and the histogram count to the cumulative sum indexed at the end.

T(i) = round(T(i));    %Rounding the T(i)

We then need to find the cumulative sum of the histogram count from 1 to T(i). We use this to find the mean above T(MAT) and MBT.

cumulative_sum_2 = cumsum(Histogram_count(1:T(i)));         %finding the cumulative sum at the second index.
MBT = sum(Bin_number(T(i)).*Histogram_count(1:T(i)))/cumulative_sum_2(end); %finding the MBT 
cumulative_sum_3 = cumsum(Histogram_count(T(i):end));      %finding the cumulative sum at the second index.
MAT = sum(Bin_number(T(i):end).*Histogram_count(T(i):end))/cumulative_sum_3(end);    %finding the MBT 

Finding the average of MAT, MBT and the obsolete value of T above one #

We will now find the threshold. This is achieved by the code below:

i = i+1;
T(i) = round((MAT+MBT)/2);      %Finding the average of MBT and MAT

Making the obsolete value great than one #

Let's now introduce a while loop to the conditions to make the absolute value of T to be greater than one:

while abs(T(i)-T(i-1))>=1
cumulative_sum_2 = cumsum(Histogram_count(1:T(i)));    %finding the histogram count at the second index
MBT = sum(Bin_number(T(i)).*Histogram_count(1:T(i)))/cumulative_sum_2(end);  %finding the histogram count at MBT
cumulative_sum_3 = cumsum(Histogram_count(T(i):end));    %finding the histogram count at the third index
MAT = sum(Bin_number(T(i):end).*Histogram_count(T(i):end))/cumulative_sum_3(end);
i = i+1;                      %looping i
T(i) = round((MAT+MBT)/2);    %rounding off the average mat and mbt.
Threshold = T(i);             % making T(i) the threshold
end

Finally we normalize the threshold as follows:

level = (Threshold-1)/(Bin_number(end)-1);
end

Going back to our initial script, we are going to call this function threshold_level. This function uses the subtracted_image we found before and returns its threshold level:

level = Threshold_level(subtracted_image);

Now we can convert this subtracted_image to binary image using the im2bw function.

binary_image = im2bw(subtracted_image, level-0.008);    %converting to binary

We then use the imshow function to display this binary image in a figure as shown below:

figure, imshow(binary_image)
title('binary image')

The resulting outputs and the further modifications #

This is what we have when running the code at this point:

This is not the result required. We need to remove the noise from the binary_image we have displayed above. To do that, we use the bwareaopen function and the binary image as the argument as shown below:

clean_image = bwareaopen(binary_image, 100);
figure, imshow(clean_image)
title('clean_image')

We will then see this new image:

Since this is a binary image, we need to convert it into a colour image. We will take the complement of the binary image first, using the incomplement function, which returns the complemented image.

A complement image is an image in which the pixels are subtracted from the maximum pixel that the class can support. It is mainly used to improve the contrast.

To convert, add the following block of code:

complemented_image = imcomplement(clean_image);   %complemented image
figure, imshow(complemented_image)
title('complemented image')

This is how the complemented image will look like:

Colorizing the image #

We can now colorize this image. In order to do that, we will create a function named colorized_image.m.

The function will return the color image. To avoid any difficulty, we will use the default color defined by DEFAULT COLOR, as shown:

function color_image = colorize_image(resized_image, complemented_image, colorspace_defination)
DEFAULT_COLOR = [1 1 1];     % default color is white.

Let's introduce the condition for the color channels. We will be using the nargin function; such that, if nargin is less than 3, we will use the default color.

nargin is a function that returns the number of the input argument of the function.

The code will be as follows:

if nargin<3
colorspace_defination = DEFAULT_COLOR;   %calling the default color function
end

We now make the complemented image to be a logical image, like this:

complemented_image = (complemented_image~=0);      %colored image

Then convert the resized image and the color space into uint8:

resized_uint8 = im2uint8(resized_image);       %converting the resized image and the color space into uint8.
color_uint8 = im2uint8(colorspace_defination);

If the dimensions of the resized_uint8 is 2, then it is a grayscale. We use the ndims to return the dimension of the input array.

If this is the case, we have to initialize all our input channel with same values, as follows:

if ndims(resized_uint8) == 2
red_channel = resized_uint8;            %initializing the colors
green_channel = resized_uint8;
blue_channel = resized_uint8;

However, in other cases, we have to define the channel, as shown below:

%defining the colors
else
red_channel = resized_uint8(:, :, 1);
green_channel = resized_uint8(:, :, 2);
blue_channel = resized_uint8(:, :, 3);
end

From the code above:

• 1 represents red, 2 represents green and 3 represents blue.

Now apply the colors to the complemented image:

%applying color to the images
red_channel(complemented_image) = color_uint8(1);
green_channel(complemented_image) = color_uint8(2);
blue_channel(complemented_image) = color_uint8(3);

Finally, we concatenate the red, green and the blue channel into a 3-D array with the following code which will give us the coloured image:

color_image = cat(3, red_channel, green_channel, blue_channel); %color image
end

We now go back to the script and call this colorize_image function, after which we will be able to view the final image:

final_results = colorize_image(resized_image, complemented_image, [0 0 0]);
figure, imshow(final_results)
title('final_result')

Finally, this is how our final image will look like:

Conclusion #

Through this article we have seen that using Matlab for segmentation is simple. This is because it uses the in-built functions that you do not need to define.

These in-built functions make the work more accessible, and the code does not seem to be bulky.

Also, Matlab allows the use of functions, and this makes the code appear organized. The output given by Matlab is reliable. It makes it even more efficient for segmentation.

I hope you find this tutorial beneficial.

Happy coding!

---

Peer Review Contributions by: Monica Masae