A New Approach for Evaluating and Improving the Performance of Segmentation Algorithms on Hard-to-Detect Blood Vessels

We aim to characterize the local salience of blood vessel segments for each pixel belonging to the vasculature. We assume that each image has a respective ground truth manual annotation mask indicating the pixels that belong to the vessel. The first step of the methodology is to represent the ground truth annotation as a graph, where nodes indicate bifurcations and terminations of the vasculature and edges represent vessel segments. For this task, we use the Pyvane framework²²2https://github.com/chcomin/pyvane. Briefly, the framework first applies a skeletonization algorithm to calculate the medial axes of the vessels. Next, skeleton pixels having one neighboring pixel or three or more neighboring pixels become a node in the graph. A pruning procedure is applied to remove small spurious branches caused by the skeletonization algorithm. A merging strategy is also applied to merge neighboring graph nodes into a single node. The final graph contains the connectivity of the vasculature as well as the positions of the nodes. More importantly for our method, each edge of the graph contains the pixels of the medial axis of the respective blood vessel segment. We henceforth use MAS to refer to the medial axis segment of a vessel. Please refer to Freitas-Andrade et al. (2022) for a complete description of the method.

The aim is to calculate the salience of the vessel at each MAS pixel and expand the calculated values to nearby pixels. The salience is calculated as the relative difference in intensity between the vessel and the background around the pixel. More details are given below. For the calculation, it is first necessary to define the cross-section of the vessel at the position of the pixel. A possible approach is to define a normal vector with respect to the MAS. We found that this approach is unstable at sharp changes of direction of the vessel as well as close to bifurcations and terminations. Thus, a simple and robust method, that is guaranteed to find relevant vessel and background values for all MAS pixels, was developed.

First, the borders of the vessels are obtained using the ground truth mask. A parametric contour tracing algorithm is used Suzuki et al. (1985). Given a MAS pixel $p$, the closest vessel contour pixel $p_{c1}$ is identified using the Euclidean distance. A vector representing the direction from $p$ to $p_{c1}$ is calculated as $v_{c1}=p_{c1}-p$. Next, the second closest contour pixel $p_{ci}$ is identified and a respective vector $v_{ci}=p_{ci}-p$ is defined. The dot product between vectors $v_{c1}$ and $v_{ci}$ is calculated. If it is negative, it means that the vectors have opposite directions, and thus two opposite contour points were found. A positive dot product means that $p_{ci}$ is at the same side of the contour as $p_{c1}$, and thus should be discarded. In such a case, the third closest contour pixel is identified and the dot product calculated, and so on until a contour pixel having a negative dot product with $v_{c1}$ is found. The final opposite pixel is represented as $p_{c2}$.

Given $p$, $p_{c1}$ and $p_{c2}$, two straight lines representing a cross-section of the vessel can be defined. The first line, $l_{1}$, goes from $p_{c1}$ to $p$ while the second line, $l_{2}$, goes from $p$ to $p_{c2}$. Figure 2 shows an example of a MAS pixel and the respective shortest lines to nearby contour pixels. The pixels belonging to both lines are henceforth represented as $S_{v}$. The values of these pixels represent the local vessel signal that will be compared with the background.

Points $p_{c1}$ and $p_{c2}$ are used to identify nearby background pixels. All pixels with a Euclidean distance smaller than or equal to $r_{b}$ from $p_{c1}$ or from $p_{c2}$ are identified. The ground truth mask is used to discard pixels belonging to the blood vessel. Figure 2 shows an illustration of the background pixels identified at this step. These pixels are represented as $S_{b}$.

Pixels $S_{v}$ and $S_{b}$ tend to form a dumbbell-like shape representing local intensities associated with $p$. $r_{b}$ is a free parameter of the method defining the scale where background pixels will be searched, but using a value leading to a similar number of pixels in $S_{v}$ and $S_{b}$ tends to work well in practice.

Different approaches can be used to calculate the local vessel salience from pixels $S_{v}$ and $S_{b}$. We chose a simple calculation to define a value that can be easily interpreted. First, the average intensity values $I_{v}$ and $I_{b}$ of, respectively, $S_{v}$ and $S_{b}$ are calculated. Next, the relative difference in intensity is obtained as

We found that $\Delta I$ can have sharp changes in value, which hinders the definition of continuous regions for the recall measure presented in the next section. Thus, it is beneficial to smooth the values along the MAS. Since the pixels of the MAS are ordered parametrically along the segment, the final LVS index of a pixel $p_{i}$ is defined as

where $\Delta I(p_{j})$ represents the local difference for pixel $p_{j}$ of the MAS. Parameter $k$ sets the degree of smoothing.

The final step of the calculation is to expand the values calculated for the MAS to other vessel pixels. For each vessel pixel, the closest MAS pixel is identified, and the respective LVS index is attributed to the pixel. The LVS index of each vessel pixel is stored as an image so that it can be easily used on downstream calculations.

Given the LVS index of each vessel pixel, it is possible to identify regions with different vessel salience. Pixels with large LVS are easier to segment. Pixels with low LVS are more challenging. But notice that it is possible to define a third class of pixels having zero or negative LVS. These pixels represent regions of vessel discontinuities in the image. Such pixels are very challenging to segment. Interestingly, recent results show that humans tend to have a strong shape bias while CNNs tend to have a strong texture bias Geirhos et al. (2019). Thus, it is expected that humans perform better than CNNs at identifying vessel discontinuities, which can be located using the shape prior that two nearby and aligned, but discontinuous, segments are likely the same segment. This means that ground truth annotations might have discontinuities correctly annotated as being blood vessels, but CNNs will struggle to identify such regions.

Given the result of a segmentation algorithm to be compared with a ground truth annotation, we define as $G$ the set of ground truth vessel pixels and as $R$ the set of vessel pixels identified by the algorithm. From these two sets the number of true positive (TP) and false negative (FN) pixels can be calculated as

Then, the traditional recall metric can be defined as

Since the number of annotated blood vessel pixels is $P=TP+FN$, the recall metric quantifies the fraction of blood vessel pixels that were successfully identified by the algorithm.

The LSRecall is defined in a similar way to the usual recall metric. Given the LVS index of all vessel pixels, a threshold value $t$ is defined, and only pixels having an LVS index equal to or smaller than the threshold are considered. This defines a binary image containing challenging vessel pixels that can be used in a usual recall calculation. We represent as $G_{t}$ the set of such pixels. Notice that $G_{t}$ is a subset of $G$. Then, the respective metrics are calculated as

The respective recall metric can be defined as

The LSRecall quantifies how successful the algorithm was at identifying challenging blood vessel pixels. Notice that $G_{t}$ can be small. Thus, large changes in LSRecall might have respective insignificant changes on the traditional recall metric. The threshold $t$ can be set to a specific value to quantify if a segmentation method is identifying pre-defined challenging regions or a family of LSRecall values can be calculated using different values of $t$. When $t=1$ the LSRecall becomes the usual recall metric.

It is relevant to point out that, similarly to the usual recall metric, the LSRecall quantifies one aspect of an algorithm’s performance related to false negatives. Other metrics can be used to quantify the amount of false positives in the result.

It will be shown in the results section that segmentation methods reaching recall values of $\approx 0.9$ can have respective LSRecall scores as low as $0.2$. Thus, it is clear that, as expected, many methods tend to struggle in regions having low salience. Thus, we also propose a data augmentation methodology to improve the robustness of methods regarding blood vessels with low salience or discontinuities. In general terms, given the MAS pixels of a blood vessel, as defined in the previous section, the salience of a subset of the pixels is systematically reduced, with an optional discontinuity region added. That is, from a given initial point of a blood vessel segment, the intensity of the vessel is slowly reduced along the segment until it becomes equal to the local background intensity. The details of the method follow below.

Two parameters of the method define the length of the salience modification for a given blood vessel segment. Parameter $l$ sets the total salience modification length, that is, only a region of length $l$ of a segment is modified. Parameter $l_{d}<l$ sets the length of the discontinuity region. Thus, a region of length $l_{d}$ will have pixels with intensity identical to the local background intensity of the segment. We note that $l$ and $l_{d}$ as well as all distances of the procedure described below are calculated along the MAS, that is, they are calculated as path lengths, not as shortest distances between two points.

Given the MAS pixels of a specific blood vessel, a reference pixel $p_{c}$ is randomly selected. This pixel will be the center of the salience reduction region. Thus, the random selection only considers pixels having a path-length distance larger than $l/2$ from both endpoints of the segment. From $p_{c}$, the two MAS pixels that are a distance $l/2$ from $p_{c}$ are selected. These pixels, represented as $p_{1}$ and $p_{2}$, define the beginning and end of the salience reduction region. Two other pixels, $p_{d1}$ and $p_{d2}$, with path-length distance $l_{d}/2$ from $p_{c}$ are also identified. These pixels define the beginning and end of the discontinuity. Figure 3(a) illustrates the positions and lengths defined.

An intensity preservation factor $f_{1}$ is defined from $p_{1}$ to $p_{d1}$ as

where $p_{i}$ represents a MAS pixel between $p_{1}$ and $p_{d1}$ and $d(p_{x},p_{y})$ represents the path-length distance between two points. Intuitively, the preservation factor is one for pixel $p_{1}$ and linearly decreases along the segment until reaching a value of zero at $p_{d1}$. A preservation factor $f_{2}$ is similarly defined for pixels between $p_{2}$ and $p_{d2}$. Between pixels $p_{d1}$ and $p_{d2}$, a preservation factor of zero is set. A final preservation factor $f$ combining all calculated values is then defined. Figure 3(b) shows the preservation factor for the MAS illustrated in Figure 3(a).

The preservation factor is calculated only for the MAS pixels. Thus, they need to be expanded to the remaining vessel pixels. For each vessel pixel, the closest MAS pixel is identified (using the Euclidean distance) and the respective factor is associated with the vessel pixel. The final result is a preservation factor for every vessel pixel between points $p_{1}$ and $p_{2}$. A preservation factor of one is set to all other pixels of the vessel when necessary.

Next, the intensities of the vessel segment are transformed. The minimum value needs to be similar to the local background intensity of the sample. This is done by taking the average value of all background pixels that are inside a rectangular region defined by the upper-left point $(p_{c}(x)-l,p_{c}(y)-l)$ and the lower-right point $(p_{c}(x)+l,p_{c}(y)+l)$ on image coordinates. $(p_{c}(x),p_{c}(y))$ are the coordinates of point $p_{c}$. The calculated value is represented as $I_{m}$. Next, the intensity $I(p)$ of each vessel pixel is transformed as

Thus, pixels close to $p_{1}$ and $p_{2}$ have little change in intensity, while pixels close to $p_{d1}$ and $p_{d2}$ become very similar to the background. Pixels between $p_{d1}$ and $p_{d2}$ become equal to $I_{m}$ since the preservation factor is zero in this region.

A specific procedure was developed to make the discontinuous region as similar as possible to the background. Vessel pixels having a preservation factor of zero define a continuous image region $S_{a}$ that must have the appearance of the image background. A procedure was developed to replace the intensities at $S_{a}$ with a background region having the same shape as $S_{a}$. This was done by inverting the ground truth annotation and applying a binary erosion using $S_{a}$ as a structuring element with central point $p_{a}$. The remaining pixels are all centers of candidate background regions that can replace the intensities in $S_{a}$. However, the new background region might be significantly distinct from the local background of the discontinuous region $S_{a}$. Thus, a criterion was used to find good background candidates for replacement.

Given $S_{a}$, the outer and inner contours of the set are obtained, as shown in Figure 4. The outer contour $C_{o}$ is defined as all background pixels having at least one 8-neighbor belonging to $S_{a}$. The inner contour $C_{i}$ is defined as all pixels in $S_{a}$ having at least one 8-neighbor in $C_{o}$. Thus, $C_{o}$ and $C_{i}$ define the pixels at the interface between the vessel and the background at the region where a vessel discontinuity is desired. The intensities $I_{o}$ of the outer contour are obtained and stored. The positions in $C_{i}$ are translated to a candidate background region found with the procedure described above, and the respective intensities $I_{i}$ are calculated. The average absolute difference between intensities of the outer and translated inner contour is obtained. The procedure is repeated to all candidate background regions.

A randomly selected background region having an absolute intensity difference smaller than a given threshold $t_{b}$ is identified, and the intensities are copied to the respective pixels of $S_{a}$. Thus, vessel pixels between $p_{d1}$ and $p_{d2}$ acquire new intensities copied from a background region having similar intensity to the local background intensity of pixels between $p_{d1}$ and $p_{d2}$. Threshold $t_{b}$ sets the largest acceptable difference between the local background and the background region selected to replace $S_{a}$.

For augmenting the images using the procedure described, the values $l$ and $l_{d}$ are randomly selected for each segment from a range of possible values that are parameters of the augmentation. The number of vessel segments to be augmented in an image is also randomly selected from a pre-specified range.

A dataset containing 2641 fluorescence microscopy images of the mouse brain vasculature is used to analyze the potential of the proposed methodologies. The samples were acquired under different experimental conditions and from different animals. Details about the image acquisition procedure and the characteristics of the samples can be found in Freitas-Andrade et al. (2022); da Silva et al. (2024). The dataset is interesting because the blood vessels contained in the samples have many different characteristics. Interesting variations include samples with different amounts of noise, contrast, and vessel caliber as well as samples containing imaging artifacts and significant intensity variation along blood vessel segments. Two examples of the dataset are shown in Figure 1.

To investigate if the LVS index can successfully identify challenging regions to segment, the LVS indices of all blood vessels in the dataset were calculated. For all experiments, a value of $r_{b}=4$ was used for identifying background pixels for the LVS calculation, and $k=15$ was used for smoothing the values (Equation 2). Figure 5 shows example results for six images. The index tends to correlate with the intensity of the blood vessels, but in regions where the intensity of the background is higher, or the background is noisier, the LVS index tends to result in smaller values, as expected. The Pearson correlation coefficient between the intensities of the pixels and respective LVS values was calculated for each sample and averaged among all samples of the dataset. A value of 0.43 was obtained, which indicates that the two quantities contain distinct information about the vessels. The regions indicated in red on the last row of plots in Figure 5 are challenging regions to segment identified using the index.

Different threshold values can be used to identify specific hard-to-segment blood vessels to focus on when evaluating segmentation algorithms. Figure 6 shows examples of thresholded indices using different threshold values. The threshold sets the expected difficulty in identifying vessel segments. Very low thresholds can be used to define regions where the blood vessels are almost indistinguishable from the background. In some cases, continuity criteria need to be used to identify regions connecting pairs of aligned, but discontinuous, vessel segments.

For the sample shown in Figure 6, segments having an LVS larger than 0.5 are easy to identify. Most segmentation methods should be able to identify these segments. Thus, as discussed above, it is more relevant to focus the performance analysis on vessels having low LVS.

To verify if the LSRecall metric can aid in measuring the performance of segmentation algorithms at difficult-to-segment regions, a CNN was trained on the blood vessel dataset and the LSRecall metric was calculated on the results. The neural network used is based on the U-net architecture Ronneberger et al. (2015) using residual blocks. Figure 7 illustrates the architecture. The network was trained on a subset of 45 randomly selected images from the dataset and validated on another 5 randomly selected images. The trained network was then applied to all images of the dataset, including those used for training and validation. This split using a small set of images for training aims at reproducing real-world scenarios where only a few samples are annotated. The training set was included in the final results because we are not interested in absolute performance values, but only in comparing different performance metrics.

The network was trained for 500 epochs using a learning rate of 0.01 with a polynomial learning rate scheduler and a batch size of 10 images. The AdamW optimizer Loshchilov and Hutter (2019) was used with a momentum of 0.9. The model with the lowest validation loss found during training was used. To avoid class imbalances, the background and vessel classes were weighted according to the inverse frequency of the pixels for loss function calculation. The only augmentation used consisted of random crops with a size of $256\times 256$ pixels. The images were normalized using z-score for each image independently (i.e., each image was normalized by its mean and standard deviation values). The network reached an average Dice score of 0.82 on the whole dataset.

Figure 8(a) shows a comparison between the LSRecall and the recall metrics calculated for all images. A threshold of 0.2 was used for the LSRecall. It is clear that the segmentation is much worse for vessels having low salience, with some images displaying LSRecall of zero. Interestingly, many samples having recall as high as 0.9 display a low LSRecall of around 0.2.

The LSRecall threshold parameter sets the largest salience to be considered when measuring the segmentation performance. Figure 8(b) shows the average LSRecall of all samples in the dataset for different threshold values. There is a clear trend of the CNN missing more and more vessel pixels as the threshold gets lower. When the threshold is one, the LSRecall is equal to the recall metric.

Having a metric that quantifies segmentation performance at low salience regions allows studying the influence of hyperparameter changes for improving the results at such regions. One possible change is the inclusion of specific data augmentations to make the method more robust to blood vessel intensity changes. This approach is studied next.

The augmentation procedure described in Section III.3 was applied to verify if the CNN can be made more robust to low-salience blood vessels. The training protocol was the same as presented in the previous section, with the addition of the salience augmentation. Specifically, for each image the number of blood vessel segments to be augmented, $n$, is randomly selected with uniform probability in the range $[50,100]$. Next, parameters $l$ and $l_{d}$ are randomly selected with uniform probability in, respectively, the ranges $[20,100]$ and $[0,30]$. A segment with length larger than $l$ is then randomly selected and the procedure described in Section III.3 is applied. The process is repeated until $n$ segments have been augmented. A value of $t_{b}=5$ was used for finding suitable background regions to generate vessel discontinuities.

Figure 9(a) shows a comparison between the LSRecall of the CNNs trained with and without augmentation. A threshold of 0.5 was used for the LSRecall calculation. The result shows that the augmentation procedure indeed tends to make the CNN more robust to blood vessels having low salience. The absolute difference between LSRecall values tends to be small, but this difference tends to be more relevant for samples having low LSRecall. This can be verified using an alternative approach to compare the values. Figure 9(b) shows the relative improvement in LSRecall using the augmentation plotted as a function of LSRecall with no augmentation. The relative improvement is calculated as the difference between LSRecall values with and without augmentation divided by the LSRecall without augmentation. The figure shows that the augmentation procedure significantly improved the accuracy for samples having low LSRecall.

The average precision (positive predictive value) of the CNNs trained with and without augmentation were, respectively, 0.72 and 0.73. Thus, the precision did not significantly change between the two training procedures.