如何使用Matlab获得视网膜图像中的血管分割结果

分割是指视网膜将数字图像分割成多个区域并提取重要区域的过程。视网膜图像是眼球的数字图像，显示了视网膜、光神经和向大脑发送信息的血管。

这种做法有助于验光师发现疾病，也能给眼睛做健康检查。

这种分割过程在医学领域至关重要。当血管被分割并被仔细观察时，可以确定某个问题的解决方案或原因。

因此，这使得它在这个领域有很好的应用。本文将探讨我们如何在视网膜图像中获得血管分割的问题。

前提条件

要跟上本教程，你需要具备以下条件。

• 安装了[MATLAB]。
• 对[MATLAB]基础知识的正确理解。
• 对使用Matlab进行[图像处理]的基本了解。

获得分割的Matlab代码

为了使整个过程更容易，请下载你的图像并将其存储在Matlab的当前文件夹中。视网膜图像可以直接从网上下载。

此后，我们将打开Matlab并创建一个新的脚本。第一步是读取图像。

这是由imread 函数实现的。

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

我们的图像在这时看起来很大。因此，我们应该使用imresize 函数来调整它的大小。这个函数使用图像名称和矢量形式的首选尺寸作为参数。

尺寸的单位是像素。这意味着不接受十进制的尺寸。

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

由于我们将对非常微小的血管进行分割，我们将需要使用im2double 函数将调整后的图像转换为双倍数据时间。

这个函数只使用图像名称作为参数。

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

在这一点上，图像是一个RGB 。我们需要将图像转换为CIE lab 颜色空间。CIE Lab色彩空间是由国际照明委员会定义的。

将图像转换为CIE实验室颜色空间

根据这个颜色空间，l 代表颜色的likeness 。likeliness 范围从0到100，其中0是黑色，100是白色。

a 代表绿色到红色，范围从负到正。负值越大，绿色越亮，而正值越大，红色越亮。

b 代表蓝色到黄色。它的范围也是由负到正。

图像将被转换为如下所示。

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

现在让我们使用cat 函数来连接1、0和0。连接是将两个或多个变量合并在一起的过程。

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

从上面的代码来看，3 代表了串联区域的尺寸。我们的图像是在CIE Lab色彩空间中，它有3个通道。

然后，我们用bsx 函数在filled 和lab 图像之间进行一个从元素到二进制的操作。

重塑输出图像

接下来，我们将对填充的图像进行重塑。尺寸参数将是空白的，因为我们在这里不需要任何参数。

相反，我们将使用3 ，因为它是filled 图像的现有尺寸，如下图所示。

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

现在我们应用主成分分析（PCA）函数。同时，我们使用reshaped_lab_image 作为参数。

该函数返回主成分的系数和分数。变量C 和S 将存储主成分的系数。

然后你可以根据lab 图像的大小来调整分数的大小，如图所示。

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

我们需要第一通道的所有行和列。

S = S(:, :, 1);

在灰度图像上进行对比度增强和过滤。

现在，让我们把S ，转换成灰度图像。首先，我们从S 中减去S's 的最小值，然后用S 的最大和最小值除以它。

这个除法将是一个从元素开始的除法，这就是为什么我们在除法符号前使用点，如下图所示。

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

我们需要对这个灰色图像进行对比度增强。这将使用自适应直方图均衡功能adapthisteq ，如图所示。

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

numTiles ，它代表的是瓦片的数量，意味着增强是逐块进行的。块的大小将是8x8 。

nBins 表示梁的数量将是 。128

下一步是对图像进行平均过滤。为了做到这一点，我们使用special 函数定义过滤器。然后我们应用imfilter 函数来计算输出图像的像素值，如下图所示。

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

现在我们将查看filtered image ，并使用imsubtract 函数从增强后的图像中减去它。为了查看图像，我们将使用imshow 函数，如下图所示。

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

计算阈值水平

现在我们需要计算阈值水平来分割血管。让我们为此创建一个名为threshold_level.m 的函数脚本。

这个函数将把图像作为参数，如下图所示。

function level = Threshold_level(image)

然后我们将图像转换为uint8 。这是一个8-bits 的数据类型。这个转换是通过im2uint8 函数完成的。

然后我们使用imhist 函数和转换后的图像作为参数，得到直方图计数和光束数，如图所示。

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

使用cumsum 函数计算直方图计数的累积和，如下图所示。

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

现在让我们找出低于和高于T 。T ，是bin数和直方图数的乘积与最后索引的累积和的比率。

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)

然后我们需要找到从1 到T(i) 的直方图计数的累积和。我们用它来找到T(MAT)和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 

找到MAT、MBT和T以上的过时值的平均值

我们现在要找到阈值。这可以通过下面的代码实现。

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

使过时的值大于1

现在让我们给条件引入一个while循环，使T 的绝对值大于1。

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

最后，我们将阈值规范化，如下所示。

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

回到我们的初始脚本，我们将调用这个函数threshold_level 。这个函数使用我们之前找到的subtracted_image ，并返回其阈值水平。

level = Threshold_level(subtracted_image);

现在我们可以使用im2bw 函数将这个subtracted_image 转换为二进制图像。

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

然后我们使用imshow 函数将这个二进制图像显示在一个图中，如下图所示。

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

得到的输出和进一步的修改

这就是我们在此时运行代码时得到的结果。

这不是所需的结果。我们需要从我们上面显示的binary_image ，去除噪音。要做到这一点，我们使用bwareaopen 函数和二进制图像作为参数，如下图所示。

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

然后我们将看到这个新的图像。

由于这是一个二进制图像，我们需要将其转换为彩色图像。我们将首先使用incomplement 函数来获取二进制图像的补码，并返回补码图像。

补色图像是一种图像，其中的像素被减去了该类图像所能支持的最大像素。它主要用于提高对比度。

要转换，请添加以下代码块。

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

这就是补充后的图像的模样。

给图像上色

我们现在可以给这个图像上色。为了做到这一点，我们将创建一个名为colorized_image.m 的函数。

该函数将返回彩色图像。为了避免任何困难，我们将使用DEFAULT COLOR 所定义的默认颜色，如图所示。

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

我们来介绍一下颜色通道的条件。我们将使用nargin 函数；这样，如果nargin 小于3，我们将使用默认颜色。

nargin 是一个函数，返回函数的输入参数的数字。

代码将如下。

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

我们现在把补全的图像做成逻辑图像，像这样。

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

然后将调整后的图像和颜色空间转换为uint8 。

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

如果resized_uint8 的尺寸是2，那么它就是一个灰度。我们用ndims 来返回输入数组的尺寸。

如果是这种情况，我们必须用相同的值初始化我们所有的输入通道，如下所示。

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

然而，在其他情况下，我们必须定义通道，如下所示。

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

从上面的代码来看。

• 1 代表红色， 代表绿色， 代表蓝色。2 3

现在将这些颜色应用于补充后的图像。

%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);

最后，我们用下面的代码将红、绿和蓝通道串联成一个三维数组，这样就得到了彩色的图像。

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

现在我们回到脚本中，调用这个colorize_image 函数，之后我们就可以看到最终的图像了。

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

最后，这就是我们最终图像的样子。

结论

通过这篇文章我们看到，使用Matlab进行分割是很简单的。这是因为它使用了内置的函数，你不需要定义。

这些内置的函数使工作更容易进行，而且代码似乎并不笨重。

另外，Matlab允许使用函数，这使代码显得有条理。Matlab所给出的输出是可靠的。这使得它在分割方面的效率更高。