Skip to content
BY 4.0 license Open Access Published by Oldenbourg Wissenschaftsverlag July 14, 2021

Segmentation of elongated structures processed on breast MRI for the detection of vessels

Segmentierung länglicher Strukturen zur Detektion von Blutgefäßen in Brust MRI
Marco Gierlinger

Marco Gierlinger will complete his bachelor’s degree in computer science in 2021. Subsequently, a master’s degree in electrical engineering will be started. Currently, he works as student assistant at Johannes Kepler University at Institute for Measurement Technology.

ORCID logo EMAIL logo
, Dinah M. Brandner

Dinah Brandner received her BSc degree in Mechatronics from the Johannes Kepler University Linz, Austria in 2019 and is currently enrolled there in the master’s program in Mechatronics specializing in Measurement Technology and Signal Processing. In July 2020 she finished her master’s thesis with the title ‘Sound Field Simulation in Breast Cancer Research’. Part of the work was done during her research stay as a Visiting Student Researcher at the Ferrara Lab as part of the Radiology Department of Stanford University, Palo Alto, USA while working on Quantitative Ultrasound (QUS) imaging using numerical simulations and experimental studies.

and Bernhard G. Zagar

Bernhard G. Zagar has been the head of the Institute for Measurement Technology of the Johannes Kepler University Linz since 2001. He has been working in this field and specializing in digital signal and image processing, sensors, laser-optical systems and magnetic tomography.

From the journal tm - Technisches Messen

Abstract

The multi-seed region growing (MSRG) algorithm from previous work is extended to extract elongated segments from breast Magnetic Resonance Imaging (MRI) stacks. A model is created to adjust the MSRG parameters such that the elongated segments may reveal vessels that can support clinicians in their diagnosis of diseases or provide them with useful information before surgery during e. g. a neoadjuvant therapy. The model is a pipeline of tasks and contains user-defined parameters that influence the segmentation result. A crucial task of the model relies on a skeletonization-like algorithm that collects useful information about the segments’ thickness, length, etc. Length, thickness, and gradient information of the pixel intensity along the segment helps to determine whether the extracted segments have a tubular structure, which is assumed to be the case for vessels. In this work, we show how the results are derived and that the MSRG algorithm is capable of extracting vessel-like segments even from noisy MR images.

Zusammenfassung

Wir bauen auf der bestehende Arbeit des Multi-Seed Region-Growing (MSRG) Algorithmus auf und erweitern ihn um damit längliche Segmente aus Brust-MR (Magnet Resonanz) Bildern zu extrahieren. Ein Modell wird von uns erstellt, das dazu dient die Eingabe-Parameter des MSRG so einzustellen, dass die länglichen Segmente potentielle Blutgefäße in der Brust darstellen um u. a. Chirurgen vor einem operativen Eingriff bei der Vorbereitung zu unterstützen. Das Modell stellt eine Aneinanderreihung von Aufgaben dar, die das Ergebnis beeinflussen sollen. Eine essentielle Aufgabe besteht darin Informationen wie Breite, Länge der Segmente und Gradient der Intensität der einzelnen Pixel entlang eines Segments zu erfassen. Mit dieser a priori Information werden Segmente extrahiert, die auf eine zylindrische Struktur hinweisen, da dies auf die zu detektierenden Blutgefäße zutrifft. In dieser Arbeit wird gezeigt, dass der MSRG Algorithmus in der Lage ist, auch bei erhöhtem Bildrauschen Segmente, die auf potentielle Gefäße hinweisen, extrahieren zu können.

1 Introduction

Seeded Region Growing (SRG) by Adams and Bischof [1] is an effective and well-known image segmentation algorithm. SRG grows one or more regions, initially called seeds, that can be single-pixel-sized or a set of adjacent pixels. The algorithm grows these distinct regions due to some homogeneous criterion until all pixels are assigned a region. The number of iterations equals the number of pixels, i. e. the algorithm halts when all pixels are partitioned.

In previous work, we presented a fast multi-seed region growing (MSRG) algorithm [2] that employs independent seeds, which can grow fully in parallel without a rendezvous before every pixel has been visited. Also, the mean intensity of a region is neglected (other than with SRG) and growth is promoted where two directly neighboring pixels meet the condition of some homogeneous criterion only.

In many region growing algorithms, k seeds typically grow k regions, and selecting a proper set of seed positions is a non-trivial task and crucial to the outcome. In our approach, we use one or more seeds but positions are homogeneously dispersed in pixel space and the number of seeds does not necessarily equal the number of extracted regions.

Detection of blood vessels is a challenging task in Magnetic Resonance Imaging (MRI) due to noise and low resolution, and particularly in breast MRI, connective and fat tissues sometimes resemble elongated structures similar to blood vessels [3].

Although region growing is known for its weakness in noisy environments, the proposed novel algorithm is capable of extracting homogeneous regions even in images exhibiting noise artefacts. We make the assumption that breast vessels have turbular structure and process spatial information of adjacent MR images of given MR image stacks to distinguish elongated fat tissue and blurry-edged structures like the breast contour from allegedly segmented vessels.

While the algorithm is still in active development, it has shown promising results on which we rely on in this current contribution. The algorithm generates a finite two-dimensional solution space of extracted homogeneous regions from a given input image. We extend the previous work with a skeletonization-like algorithm that is capable of filtering elongated segments from this solution space that may help to detect vessel-like structures in breast MRI. Information about location, size and morphology of vessels may provide clinicians with useful information before surgery during e. g. a neoadjuvant therapy.

Figure 1 
Region based segmentation applied to a breast MRI dataset [5].
Figure 1

Region based segmentation applied to a breast MRI dataset [5].

In Section 2, the algorithm of the previous work is explained in detail. The subsequent section will describe the skeletonization-like algorithm and show the results of allegedly extracted vessels. We evaluate the results intrinsically: The proposed method is applied to two noise exhibiting input MRI stacks from the same patient (and same breast), that are acquired during two sessions to show the likelihood of the occurrence of artefacts. As research shows [3], [4], blood vessels in breast MRI are determined subjectively by radiologists, which is considered as the gold standard. In [3], authors claim to be the first who perform 3D blood vessel segmentation in breast MRI. They employ a Hessian method to enhance tubular structures.

2 Previous work

The following algorithm is described for n-bit grayscale images in N 0 w × h , however, adapting it for higher dimensions should be straightforward. Multiple seeds are employed and aligned as grid with a user-defined width u and height v with u , v N, u w and v h. For every k N 0 , with k < u v, a seed position vector ( x , y ) N 0 2 is defined as:

(2.1) x = w u ( k mod u ) + w 2 u , y = h v k u + h 2 v ,

as schematically depicted in red ( u = v = 3) on the blue pane of Fig. 1. The algorithm walks through pixel space in an 8-adjacency flood-fill manner independently for each of the k seeds as follows: Let p N 0 3 = ( p 0 , p 1 , p 2 ) be a pixel, where ( p 0 , p 1 ) is its position in pixel space and let I ( p ) = p 2 be its intensity value. A bucket queue is used to hold data objects ( p , q ), where q N 0 with q < 2 n is the bucket’s index. A bucket is prioritized in the queue by its index and contains only elements with the same index q. Initially, the queue is filled with the seed only, which is a single pixel only. An iteration i is initiated by polling a bucket B i . a B i b N ( a ), a cost function δ f : ( a , b ) { r N 0 r < 2 n } is applied, where N ( a ) denotes all non-visited neighbors of a . Cost function δ f is defined as:

(2.2) δ f ( a , b ) = I ( a ) I ( b ) .

At the end of an iteration, each result is added to the queue as ( b , δ f ( a , b ) ). A pixel counts as ‘visited’ when it is polled from (and not added to) the bucket queue.

Figure 2 
Left: Cumulative frequency of number of assigned pixels 


m


i

{m_{i}} on ‘pliers’ 
(
u
,
v
)
=
(
640
,
480
)(u,v)=(640,480) and largest step 


m


max

{m_{\max }} is found at highlighted iteration i. Right: Image with initial seed positions indicated.
Figure 2

Left: Cumulative frequency of number of assigned pixels m i on ‘pliers’ ( u , v ) = ( 640 , 480 ) and largest step m max is found at highlighted iteration i. Right: Image with initial seed positions indicated.

2.1 Heuristic

Additionally, the number of newly assigned pixels m i is tracked for each iteration i. A map M s N 0 w × h , drawn as an orange pane in Fig. 1, is used for the s t h of k seeds and every r a 0 , a 1 M s , where ( a 0 , a 1 ) is the position of a , is set to m max ( i ) = max ( m 0 , m 1 , ... , m i ) at iteration i. Finally, when the queue is empty, i. e. every pixel was visited, M s is converted into a binary map by

(2.3) r M s = 0 r < m max ( i ) 1 otherwise

Then, all k binary maps are added up elementwise to s = 1 k M s , i. e. additive image fusion is applied, and normalized to the range from zero to 2 n 1 to suppress regions that occur less often than others. An example is highlighted in Fig. 1 (right) with random colors assigned for regions containing the same numerical value. The more binary maps are created the more details are revealed in the segmentation result in form of sub segments.

The following pseudo code in Algorithm 1 gives an additional overview of the algorithm for a single seed.

Other than described above, with slight adaption, it is not necessary to visit every pixel. The iteration can halt when the number of non-visited pixels becomes smaller than m max , since m max cannot change its value anymore. m max can then be assigned in the map M s to all positions of non-visited pixels without further examination. For the sake of simplicity, not every considered optimization is noted.

Algorithm 1 
Generating binary Map 


M


s

{M_{s}} for a single seed.
Algorithm 1

Generating binary Map M s for a single seed.

Growth is depicted in Fig. 2 and Fig. 3 for two arbitrary seed positions. Figure 2 shows a cumulative histogram of the pixel assimilation process and highlights the iteration for each seed where the largest peak is detected and Fig. 3 shows the corresponding growing process.

Figure 3 
Growth of blue seed at (60,60) and orange seed at (320,240) on ‘pliers’.
Figure 3

Growth of blue seed at (60,60) and orange seed at (320,240) on ‘pliers’.

2.2 Preprocessing

The algorithm’s result is influenced by an input image and the selected grid size so far. However, it is desired to dynamically ‘scan’ for multiple acceptable solutions. Image quantization is used as a preprocessing step to decrease the image’s bit depth, which has shown to be very useful to find practical solutions. Scaling parameter g R with 0 < g 1 is used to scale the image range to result in a lower bit depth. A user can tune the grid size u v, i. e. the density of dispersed seeds, and scaling parameter g as shown in Fig. 4 until a suitable solution is found. We want to refer to [6], where image quantization is investigated when applied as a preprocessing step for dimensionality reduction in image classification pipelines.

Figure 4 
Tuning through the algorithm’s solution space by adjusting the seed grid size 
u
vuv and scaling parameter g. Random colors are assigned to distinctive regions. The fourth column from left contains the same regions as the third one and regions are colored by its mean intensity of the pixels of the original image.
Figure 4

Tuning through the algorithm’s solution space by adjusting the seed grid size u v and scaling parameter g. Random colors are assigned to distinctive regions. The fourth column from left contains the same regions as the third one and regions are colored by its mean intensity of the pixels of the original image.

3 Filter

In this section we initially detail an algorithm that filters elongated segments from the result set of the multi-seed region growing (MSRG) algorithm. Furthermore, in the second subsection, a model is depicted to show how these two algorithms interact to obtain vessel-like structures from an input breast MRI stack. The stack is defined as follows: Let S i be the i th image on the breast MRI stack, that is an 8-bit grayscale image, we define C i as an RGB image, that is composed of S i as red channel, S i + 1 as green channel and S i + 2 as blue channel. For example, the top of the second stack from left ( C 24 ) in Fig. 5 shows such an 8-bit 256×256 RGB image as one of 34 possible slices.

Figure 5 
From left to right: Five of a set of 36 Breast MRI slices of a single session → RGB image composited by three grayscale images from the left stack → Filtered MSRG algorithm 
fi
(
C
,


g


off


+
n
·


g


inc


)
→\operatorname{fi}(C,{g_{\text{off}}}+n\cdot {g_{\text{inc}}})\to  Result stack showing a single result slice.
Figure 5

From left to right: Five of a set of 36 Breast MRI slices of a single session → RGB image composited by three grayscale images from the left stack → Filtered MSRG algorithm fi ( C , g off + n · g inc ) Result stack showing a single result slice.

3.1 Algorithm

Elongated vessel-like segments are filtered out from the set of segments derived by the MSRG. We apply a skeletonization-like algorithm sk ( msrg ( C , g , u , v ) ) = M sk , where C is an arbitrary RGB image and g, u, v the other MSRG input parameters as described in the previous section. Result M sk is a binary map containing the filtered segments. The algorithm sk ( msrg ( C , g , u , v ) ) is described for a single segment of the MSRG result set as follows.

Initially, we select a random pixel p from given segment A and 8-adjacently flood-fill until every pixel has been ‘filled’. Let N ( p ) be the set of adjacent filled pixels of p in pixel space of A. Each iteration, we prioritize from the set of potential new pixels N ( A ) = { r A ¬ isfilled ( r ) N ( r ) } the subset N max ( A ) N ( A ), defined as the largest group of adjacent pixels as highlighted in Fig. 6 left.

Figure 6 
A Flood-fill algorithm that promotes nearly homogeneous flooding by prioritizing the largest adjacent group 


N


max


(
A
)
⊆
N
(
A
){N_{\max }}(A)\subseteq N(A).
Figure 6

A Flood-fill algorithm that promotes nearly homogeneous flooding by prioritizing the largest adjacent group N max ( A ) N ( A ).

This rule leads to a fast method such that the segment is flood-filled nearly homogeneously along an elongated segment as shown in Fig. 6 right. Each iteration and until every pixel has been examined, we derive an updated set of N max ( A ). We use this set for statistical analysis: E. g., N max ( A ) gives a rough estimation of the thickness of an elongated shape each iteration due to the homogeneous growth. This statistical information is accumulated to filter out segments that are too short, too thick, or have a too large variance in thickness or intensity along the potential vessel or along its cross section. With this filter, the number of parameters increases but the parameters could be adjusted by medical experts to fit the realistic requirements of the vessels in question.

3.2 Model

We generate with the MSRG a one-dimensional solution space: The seed grid size ( u , v ) is set to (4,4), i. e. we select a column as shown in Fig. 4. These values are selected based on user experience with the MSRG. We define within this subsection the filtered MSRG algorithm fi ( C , g ) = sk ( msrg ( C , g , 4 , 4 ) ), where sk ( msrg ( C , g , 4 , 4 ) ) denotes the filter algorithm from the previous subsection. The algorithm fi ( C , g ) is applied to RGB (3 channels) image C N 0 w × h × 3 with image width w and image height h and results in a binary map M fi N 0 w × h . The image quantization factor g R with 0 < g 1 is now the only parameter that influences the result.

To keep the computation time low, we do not iterate over all possible values of the image quantization factor g. We limit the range from g off to g off + g inc · g res , where g off denotes the lower bound, g res the resolution, and g off + g inc · g res the upper bound of the dimension. Again, the values are selected based on user experience of the MSRG and should be selected such that the desired features are found within these bounds.

The model is schematically depicted in Fig. 5. The results are binary images with elongated structures that potentially reveal vessels. Formally, we generate a result stack (see Fig. 5 right) for all j N with 1 < j ( MRI stack size 3 ), defined as:

(3.1) Y j = n = 0 g res co ( C j 1 , C j , C j + 1 , g off + n · g inc )

where co ( R 1 , R 2 , R 3 , g ) returns a single binary map that is composed of the union of the two adjacent slice pairs. Formally,

(3.2) co ( R 1 , R 2 , R 3 , g ) = fi ( R 1 , g ) fi ( R 2 , g ) fi ( R 2 , g ) fi ( R 3 , g )

where ∧ and ∨ denote element-wise binary operators. Finally, as a post processing step, we set

(3.3) r Y j = 0 r < t 1 otherwise

where t is a user-defined threshold.

Figure 7 
Left: A potential vessel with color gradient varying lengthwise. Right: Breast contour with color gradient varying along the cross section.
Figure 7

Left: A potential vessel with color gradient varying lengthwise. Right: Breast contour with color gradient varying along the cross section.

Figure 8 
Each row contains four adjacent slices of the result stack as shown in Fig. 5 right. Each row’s input stack (Fig. 5 on the left) was acquired from the same patient (Patient:1010 of I-SPY dataset [5]) but several weeks elapsed between these sessions.
Figure 8

Each row contains four adjacent slices of the result stack as shown in Fig. 5 right. Each row’s input stack (Fig. 5 on the left) was acquired from the same patient (Patient:1010 of I-SPY dataset [5]) but several weeks elapsed between these sessions.

With this technique, it is possible to differ the contours of the skin from elongated segments that are within the region of interest. We assume that the color gradient more likely varies along the cross section as shown in Fig. 7. This can be detected by the algorithm of the previous subsection by statistical analysis of N max ( A ).

4 Results and evaluation

We apply the proposed method of the previous section to two different breast MRI stacks that were acquired from the same patient (and same breast) during two sessions. Several weeks elapsed between these two sessions. Algorithm parameters are set for both stacks equally as follows: g off = 0.1, g res = 100, g inc = 10 4 , seed-grid-size = 4 × 4. Each slice has a256×256×8 bit resolution and a stack has 36 slices, obtained from the I-SPY dataset [5]. The result is shown in Fig. 8 for four adjacent slices of each result stack. Slices of these two stacks are shifted such that they roughly match each other from column to column in Fig. 8. We derive similar results for both stacks, although the region of interest is shaped differently as seen in the two RGB images of Fig. 8. However, these binary results do not cover all elongated structures as visible in the RGB images. A better strategy is required to efficiently look for elongated segments in the solution space of the MSRG algorithm to derive more results. Also, a slice represents 3 mm thickness and obviously, a much higher resolution would improve the results for thinner vessel-like structures.

Due to the breast varying in shapes in different postures, it is very challenging to detect common features, however, this method might be used for automatic image registration.

5 Conclusions

With our proposed filter and combined with the MSRG algorithm from the previous work, it seems possible to assist clinicians in detecting vessel-like segments. The employed algorithm filtered elongated structures from the solution space of the MSRG algorithm. The results seem promising, however, they could not be evaluated for the correct detection of vessels without medical expertise. Instead, the proposed method was applied to two input MRI stacks from the same patient (and same breast), that were acquired during two sessions. The evaluation showed how barely the results deviate from each other. Although vessels will be positioned slightly differently between two sessions due to shaping the breast dependent on the actual posture in the MRI scanner, this method seemingly extracted similar elongated segments for both input MRI stacks. This suggests the assumption that the extracted segments are not artefacts. Since noise is a common problem in MRI, our approach may help to circumvent this issue.

Based on user experience, the solution space of the MSRG algorithm was searched through. Due to the complexity and computation time, it is impractical to search through the whole solution space. It is required to adjust the MSRG parameters accordingly, however, the MSRG algorithm is not guaranteed to be complete.

Further investigations are required for more efficient search strategies. Future work may focus much more on the automatic adjustment of the MSRG algorithm parameters via state of the art machine learning approaches.

Funding statement: This work has been supported by the LCM K2 Center within the framework of the Austrian COMET-K2 Programme.

About the authors

Marco Gierlinger

Marco Gierlinger will complete his bachelor’s degree in computer science in 2021. Subsequently, a master’s degree in electrical engineering will be started. Currently, he works as student assistant at Johannes Kepler University at Institute for Measurement Technology.

Dinah M. Brandner

Dinah Brandner received her BSc degree in Mechatronics from the Johannes Kepler University Linz, Austria in 2019 and is currently enrolled there in the master’s program in Mechatronics specializing in Measurement Technology and Signal Processing. In July 2020 she finished her master’s thesis with the title ‘Sound Field Simulation in Breast Cancer Research’. Part of the work was done during her research stay as a Visiting Student Researcher at the Ferrara Lab as part of the Radiology Department of Stanford University, Palo Alto, USA while working on Quantitative Ultrasound (QUS) imaging using numerical simulations and experimental studies.

Bernhard G. Zagar

Bernhard G. Zagar has been the head of the Institute for Measurement Technology of the Johannes Kepler University Linz since 2001. He has been working in this field and specializing in digital signal and image processing, sensors, laser-optical systems and magnetic tomography.

References

1. Adams R., Bischof L. Seeded region growing. IEEE Transactions on Pattern Analysis and Machine Intelligence, 16(6), pp. 641–647, 1994.10.1109/34.295913Search in Google Scholar

2. Gierlinger M., Brandner D., Zagar B. Multi-seed region growing algorithm for medical image segmentation. M. Heizmann and T. Laengle: Proceedings of FORUM BILDVERARBEITUNG 2020, KIT Scientific Publishing, pp. 267–278, 2020.Search in Google Scholar

3. Kahala G., Sklair M., Spitzer H. Multi-scale blood vessel detection and segmentation in breast mris. Journal of Medical and Biological Engineering, 39, pp. 424–430, 2019.10.1007/s40846-018-0418-6Search in Google Scholar

4. Lin M., Chen J.H., Nie K., Chang D., Nalcioglu O., Su M.Y. Algorithm-based method for detection of blood vessels in breast mri for development of computer-aided diagnosis. Journal of magnetic resonance imaging, 30, pp. 817–824, 2009.10.1002/jmri.21915Search in Google Scholar PubMed PubMed Central

5. Newitt D., Hylton N., on behalf of the I-SPY 1 Network and ACRIN 6657 Trial Team. Multi-center breast dce-mri data and segmentations from patients in the i-spy 1/acrin 6657 trials. The Cancer Imaging Archive. http://doi.org/10.7937/K9/TCIA.2016.HdHpgJLK, 2016.Search in Google Scholar

6. Ponti M., Nazare T.S., Thume G.S. Image quantization as a dimensionality reduction procedure in color and texture feature extraction. Neurocomputing, 173, pp. 385–396, 2016.10.1016/j.neucom.2015.04.114Search in Google Scholar

Received: 2021-02-13
Accepted: 2021-07-30
Published Online: 2021-07-14
Published in Print: 2021-08-27

© 2021 Gierlinger et al., published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.

Downloaded on 6.12.2022 from frontend.live.degruyter.dgbricks.com/document/doi/10.1515/teme-2021-0017/html
Scroll Up Arrow