WO2015150320A1

WO2015150320A1 - Segmentation of tubular organ structures

Info

Publication number: WO2015150320A1
Application number: PCT/EP2015/056882
Authority: WO
Inventors: Rahul Prasanna KUMAR
Original assignee: Universitetet I Oslo
Priority date: 2014-04-01
Filing date: 2015-03-30
Publication date: 2015-10-08
Also published as: GB201405837D0

Abstract

A method for segmentation of tubular organ structures in three dimensional (3D) image data comprises: extracting a two dimensional (2D) representation of a cross-section of the structure at the seed point; estimating a radius of the structure at the seed point; applying a circle structure enhancement filter to the two dimensional representation; thereby identifying a centre point of the structure at the seed point; extrapolating along a direction of the structure to find another point on the structure, which is defined as a new seed point; and repeating the method to thereby trace the direction and radius of the structure. The method may be used for segmentation of blood vessels using MRI, CT or ultrasound image data.

Description

SEGMENTATION OF

TUBULAR ORGAN STRUCTURES

The present invention relates to a method for segmentation of tubular organ structures, such as fluid carrying vessels within the body, for example segmentation of blood vessels. The method may be used to visualise or map the vessels within the body. The invention also extends to an apparatus and a computer programme product for performing the method.

The human body contains several different types of blood vessels that constitute a network of arteries and veins. Visualization of these blood vessels is important for improving the planning and navigation in several interventional procedures. It is specifically relevant in catheter based procedures like stent graft positioning and valve replacement and also, in planning and navigation of laparoscopic liver resection. Detailed images of these blood vessels are initially formed using computerized tomography (CT) and/or magnetic resonance (MR) depending on the part of the body and equipment available. These blood vessel contrast enhanced images provide a basis for visualizing the blood vessels in three dimensions, and this involves segmentation of blood vessels.

Along with a segmented model of the blood vessels, centreline and radius information at different cross-sections of the blood vessel are also crucial for navigating the catheter. The centreline of the blood vessels can be used as a guiding path for the catheter and radius information will help the interventionists in selecting the right valve or stent for the specific procedure. The centreline extraction may also be useful for fast registration and update of the blood vessel models generated during the intervention. Even though other methods are available for registration, centreline based methods have proven to be fast.

The techniques used for blood vessel segmentation can also be used for visualisation and mapping of other tubular organ structures in the body, airways, nerve axons and so on. The most common methods used today for blood vessel segmentation, and visualisation of other tubular organ structures in the body, are based on either multi scale - space or model detection approaches. Both of these techniques require a large amount of processing time and/or put a large load on computer processing resources.

Typically the so-called Frangi vesselness method is used for blood vessel segmentation ("Multiscale vessel enhancement fitting", Frangi et al, MICCAI, Lecture Notes in Computer Science, 1496, 130-137 (1998)). The Frangi vesselness method involves fitting a 3-D shape to the blood vessel based on a filtering process. The Hessian matrix is used in an analysis of image data to determine the degree to which a tubular structure in the image matches an ideal tubular structure. This allows for an estimation of, and hence a

visualisation of, the structure that has been imaged, which is typically a blood vessel. The Frangi method is often combined with thinning as described in T.-C. Lee, R. L. Kashyap, and C.-N. Chu, "Building skeleton models via 3-d medial surface axis thinning algorithms," CVGIP: Graphical Models and Image Processing, vol. 56, no. 6, pp. 462-478, 1994, and this combined method is typically referred to as Frangi and Lee.

US 2003/0056799 discloses one example of segmentation that uses the Frangi method. Cylinders are fitted in a cylindrical structure in order to understand the complete structure. The fitting is done by multiscale vesselness according to Frangi, with region growing. The method of US 2003/0056799 is considered to be relatively time consuming.

US 2007/0249912 is another example disclosure of the use of Frangi vesselness for vessel detection and vessel cross-section detection. Seed points are used to identify image regions that are known to be within the vessel of interest. The seed points are starting points for mapping of the vessel. The vessel boundaries are detected using a watershed algorithm in 2D cross-section and the visualisation of the vessel is propagated along the course of the imaged vessel using statistical front propagation of the previous found output.

Viewed from a first aspect, the present invention provides: a method for segmentation of tubular organ structures in three dimensional (3D) image data, the method comprising: extracting a two dimensional (2D) representation of a cross-section of the structure at the seed point; estimating a radius of the structure at the seed point; applying a circle structure enhancement filter to the two dimensional representation; thereby identifying a centre point of the structure at the seed point; extrapolating along a direction of the structure to find another point on the structure, which is defined as a new seed point; and repeating the method to thereby trace the direction and radius of the tubular organ structure.

In this method the inventors have made the realisation that accurate tracing of a tubular organ structures can be achieved using only 2D data at the point of interest. This greatly reduces the data processing requirement of the method. The 2D analysis can be repeated along the length of the structure without the need for any 3D data to be analysed. The use of a circle structure enhancement filter, also referenced herein as a 'circleness' filter, provides an efficient way to find a centre point of the 3D structure of the tubular organ based on a 2D representation of the cross-section.

The starting seed point may be provided by the user or it may be obtained via computer image recognition. Subsequent seed points, for propagation of the trace along the course of the structure in the body, are determined based on the previously calculated centre point of the structure and on the direction of the structure. The direction of the structure for the starting seed point may be provided by the user or it may be obtained via computer image recognition. For subsequent seed points it is preferred for the direction to be based on the previous results.

In a simple example the direction to find the next seed point could be based on extrapolation of a direction from the previous two or more seed points, for example based on a line fitted to the co-ordinates of prior seed points and extended. In one preferred embodiment the direction of the structure is determined in accordance with equation (1 ) below. Other techniques could also be used. What is important for the preferred method is that the previous centre points, determined by the circle structure enhancement filter, are used to estimate the position of the next seed point, and at this next seed point the circle structure enhancement filter is again used to find a new centre point. Thus, the next seed point does not necessarily need to be close to the centre of the structure, provided that it remains within the extent of the structure in the next two dimensional representation.

The distance between seed points may be set by the user or adjusted automatically based on parameters of the structure. For example, the distance between seed points could be set based on the radius of the structure. This helps ensure that propagation of the structure proceeds at an appropriate 'resolution' give the size of the structure. In a preferred embodiment the distance between first and second seed points is the same as the radius at the first seed point.

It is preferred for the starting seed point to be placed at the main trunk in a network of structures, i.e. before any branching of the structures. In the case that the original image contains multiple separate networks of structures then the user may place multiple starting seed points, so that each network of structures can be traced.

The method may comprise determining eigenvalues and eigenvectors for the structure at the seed point by eigenanalysis to thereby find principle directions of curvature of the structure; and extracting the two dimensional (2D) representation of the cross-section of the structure at the seed point based on the results of the eigenanalysis. The eigenvectors provided by this eigenanalysis may be used for finding the radius of the structure in the two dimensional representation.

The eigenanalysis at the seed point is preferably based on the three dimensional Hessian matrix combined with a Gaussian convolution for a scale varying with the width of the structure. In one example this is carried out based on equation (3) below.

For finding the radius, an edge-detection filter, such as a Canny edge detection filter, may be used in order to identify the edge of the structure. It is preferred to measure the width of the structure in at least two cross-directions, preferably more than two, and to take the largest dimension as the radius. This allows for elliptical shapes. Preferably the circle structure enhancement filter is applied only to data identified as being within the structure. For example, this may be the data for the two dimensional representation that is within the edge of the structure. This further reduces the processing burden without any loss in accuracy.

The circle structure enhancement (circleness) filter may make use of an

eigenanalysis of the two dimensional representation. A two dimensional Hessian matrix may be used. If an eigenanalysis has already been performed at the seed point then it is particularly efficient to use the eigenvalues to analyse the structural information of the cross- section image. In one example a circle enhancement ratio is defined based on the two eigenvalues for the two dimensional representation, wherein the circle enhancement ratio is the modulus of the ratio of the sum of the eigenvalues divided by the difference between the two eigenvalues. This circle enhancement ratio will be larger for pixels at the centre of the cross-section of the structure than for pixels closer to the edge of the structure. The peak value of the circle enhancement ratio hence indicates a centre point for the structure.

In one example the circle enhancement ratio is as shown in equation (5) below. This may be used in a circleness filter as shown in equation (6) below, with a first component of the equation providing a measurement of how circular the image is and a second component of the equation providing noise reduction by means of a 'structureness' term as in equation (7) below.

It will be appreciated that there are alternative techniques for circle structure enhancement filtering of the two dimensional representation. Whilst the use of eigenvalues is preferred, any other suitable technique could also be used and would provide advantages compared to prior art 3D vesselness techniques. For example, one possibility would be to use Hough circle detection or similar techniques.

The method described above is highly effective to trace a trunk of a tubular structure or a tubular structure with just a single main path. When the structure splits, typically in a bifurcation, then it is necessary to detect this and to identify multiple new seed points if it is required to track all parts of the structure after a split. In the discussion below a bifurcation analysis is described. This is termed bifurcation analysis to reflect the fact that in a large proportion of cases bifurcations are the only form of split for a tubular organ structure, especially structures in the form of vessels such as blood vessels and airways. However, the bifurcation analysis is not necessarily only for bifurcations and may also identify trifurcations and even divisions of a trunk of the tubular structure into greater numbers of branches.

In preferred embodiments, when the radius increases significantly then the method determines that a splitting of the structure such as a bifurcation may be present and a bifurcation analysis is performed. A significant increase in radius may be defined as a radius increase above a set threshold for a given distance between seed points, for example a radius increase of greater than the distance between the current seed point and the previous seed point, a radius increase of greater than 50% of the distance between the current seed point and the previous seed point, or an increase in radius of above 150% of the previous radius. The method may alternatively or additionally determine that a bifurcation may be present by detecting an increase and decrease in radius in consecutive cross-sections.

Another example method for detecting a splitting of the structure, which may be used instead of or in addition to the use of radius changes, is to measure changes in circularity and/or radius variance of the structure. Circularity may be defined as the ratio between the square of the perimeter length of the structure cross-section and it's area. A perfect circle has the maximum area for a given perimeter and hence gives the smallest ratio. A less circular shape will give a larger value. The method may include calculating a circularity ratio and determining that a splitting of the structure is present when the circularity ratio is above a threshold value. Radius variance may be defined as the variance in distance from a centre point to the perimeter of the cross-section. The method may include determining radius variance and determining that a splitting of the structure is present when the radius variance exceeds a threshold value. Detection of a splitting of the structure may use circularity and radius variance in parallel, with a splitting being identified when one or both exceed a respective threshold value. As with the method above, when a splitting of the structure is detected then a bifurcation analysis is performed.

The bifurcation analysis may comprise determining multiple centre points for the structure at the point of bifurcation, such that subsequent propagation of the trace of the structure will have multiple seed points and can hence follow multiple paths. The method may include reverting to non-bifurcation analysis, with only a single centre point being found for each path, when it has been determined that the bifurcation is complete. This may be, for example, when the two dimensional cross-section shows multiple distinct structures.

The multiple centre points may be determined by applying a circle structure enhancement (circleness) filter to identify multiple possible circle centres. When the cross- section of the structure in the two dimensional representation is at a bifurcation or other split then the increase in radius will be an increase in maximum radius and, for a bifurcation typically the result will be an elliptical shape at the start of the bifurcation, progressing to a two-lobed shape (for example, a figure eight) and then eventually becoming two separate structures. A split into more than two branches will form a more complex three lobed shape. A suitable circle structure enhancement (circleness) filter will produce multiple peak values indicative of multiple possible centre points. This then allows multiple seed points to be set for the next stage of propagation of the structure. For example, the two dimensional circle enhancement ratio described above may be applied, with the result being a mapping with multiple peaks of the ratio at multiple points spaced apart from the lobes of the structure cross-section.

Thus, by permitting multiple circleness with multiple centre points in response to a possible bifurcation then it is possible to continue the two dimensional calculation through a bifurcation.

As an alternative, it can be advantageous to switch to a vesselness analysis for the bifurcation analysis. That is to say, in response to a significant increase in radius detected during propagation of the structure, the method may switch to a bifurcation analysis using three dimensional data and fitting three dimensional shapes to the data. This increases the processing required, but it can provide advantages when the structure undergoes complex bifurcations (for example involving a reversal of direction for one branch after splitting) and/or where the structure splits into more than two branches. This can occur, for example, in the blood vessels of the human liver. A 3D method can more effectively trace complicated splitting patterns since the 3D samples of the image can overlap. When 2D representations are used then inevitably there will be 'gaps' between the 2D 'slices' of the structure and these gaps contain data that cannot be taken into account using the 2D method. By the use of a 3D method when a possible bifurcation is detected it is possible to keep the benefits of faster processing whilst tracing the majority of the structure, but also to access the benefits of 3D processes when mapping the more complex changes in the form of the tubular structure that occur during bifurcations.

The vesselness analysis could be any known technique, for example the commonly used Frangi method. It is preferred to use an enhanced 3D method developed by the current inventors, as described in detail in R. P. Kumar, F. Albregtsen, M. Reimers, T. Lang0, B. Edwin, and O. J. Elle, "3d multiscale vessel enhancement based centerline extraction of blood vessels," in SPIE Medical Imaging. International Society for Optics and Photonics, 2013, pp. 86 691X-86 691 X. It is of course important to understand that the current proposal differs significantly from the prior art in that the 3D analysis is applied only for bifurcations and ceases once the bifurcation has passed. In the prior art the 3D analysis is used for the whole structure, leading to significant disadvantages in processing time.

In one example the method comprises a vesselness equation determined from an eigenanalysis, preferably based on the Hessian matrix and a Gaussian convolution, for example a vesselness equation as set out in equation (8) below. It is preferred for the vesselness equation to be modified by addition of a term for noise correction, for example as set out in equation (10) below. By allowing for the vesselness analysis to return multiple values then multiple centre points can be determined for the structure at the seed point in a bifurcation. The method can then proceed to a new set of seed points, which are determined based on extrapolation from the preceding seed points. The vesselness analysis can be repeated if it is determined that a bifurcation is still in progress, for example if there are not yet multiple separate structures in cross-section. If the bifurcation is completed then the method preferably reverts to a two dimensional analysis using a circle structure enhancement (circleness) filter, as described above.

The image data that is utilised in the above discussed method can be any suitable three dimensional image data, for example MR images, Computer Tomography (CT) images, Rotational angiography, Confocal Microscopy, Focal Modulation Microscopy, Optical Coherence Tomography (OCT) or three dimensional ultrasound images. What is important is that the images can show the structures with sufficient contrast to the surrounding tissue, in order that the extent of the structures can be determined by edge detection or similar techniques.

The method is generally applicable to tubular organ structures that can be viewed in this type of image. It is envisaged that an important application will be segmentation of blood vessels, but the method could also be used for other fluid vessels, such as airways, urinary tracts and so on. It could also be used for other tubular organ structures, for example nerve axons. It is preferred for the method to be used as part of imaging of three dimensional networks of structures within the human or animal body.

Viewed from a second aspect the invention provides a computer programme product containing instructions that, when executed, will configure a data processing apparatus to perform the method of the first aspect and optionally any or all of the further method steps described above.

In a third aspect, and preferred embodiments thereof, the invention provides a data processing apparatus configured to perform the method of the first aspect and optionally any or all of the further method steps described above.

Certain preferred embodiments of the invention will now be described by way of example only and with reference to the accompanying drawings in which:

Figure 1 is a flowchart showing a segmentation process;

Figure 2 is a diagram illustrating a sequence of 2D cross-sections of a blood vessel as used for segmentation in the method described herein;

Figure 3 illustrates the sudden change in radius that occurs at a bifurcation;

Figure 4 shows the eigenvalues along principle directions of a tubular structure using the notation adopted herein; Figure 5 shows (a) a cross-section of a vessel at a seed point and (b) Canny output from the cross-section image, with width d1 and height d2;

Figure 6 shows (a) a 2D cross-section input image of a vessel and (b) output from a 'circleness' filter applied to the 2D cross-section as described herein;

Figure 7 shows (a) a 2D cross-section input image of a vessel at a bifurcation and (b)

2D output from a 'vesselness' volume filter applied at the bifurcation as described herein;

Figure 8 shows: (a)(f)(k)(p)(u) Maximum intensity projection (MIP) images of input medical image, with circles representing the seed points, (b)(g)(l)(q)(v) MIP images of centrelines extracted using Frangi and Lee, (c)(h)(m)(r)(w) MIP images of centrelines extracted using an enhanced prior art 3D/vesselness method, (d)(i)(n)(s)(x) MIP images of centrelines extracted using the method described herein, and (e)(j)(o)(t)(y) 3D model visualizations of blood vessels segmented using this method;

Figure 9 is a plot showing mean and standard deviation of centre estimation error of the propose method for various vessel radii;

Figure 10 shows the percentage of bifurcating vessels not detected at various bifurcation radii using the 3D bifurcation analysis described herein;

Figure 1 1 shows a workflow for another proposed segmentation process;

Figure 12 is a vessel cross-section image indicating vessel cross vectors v₂ and v₃ as used in the process of Figure 1 1 ;

Figure 13 illustrates shape descriptor values at various cross-section images along a blood vessel; and

Figures 14(a) and 14(b) show intensity plots as follows: (a) top row: vessel trunk cross-section image and its 3D intensity plot, (a) bottom row: single scale circleness image of trunk cross-section image and its 3D intensity plot; (b) top row: vessel bifurcation cross- section image and its 3D intensity plot, and (b) bottom row: multi-scale circleness image of bifurcation cross-section image and its 3D intensity plot.

The methods described herein are based on local structure analysis within the connected regions of the tubular organ structures, which are blood vessels in the current examples. It will of course be appreciated that the same techniques are applicable to other tubular organ structures in the body or elsewhere: segmentation of a blood vessel is just one preferred use for the methods described herein. The methods require that the structures can be shown with sufficient contrast in images such as MR, CT or ultrasound images.

One process is shown as a flow chart in Figure 1. This method is mainly divided into trunk analysis and bifurcation analysis. In the trunk analysis, the vessel is segmented using 2D cross-section analysis and in bifurcation analysis, the vessel is segmented using modified vesselness. The examples use angiogram images. On comparing the processing time for finding the blood vessel centreline, the proposed method is found to be on average more than twenty times faster than prior art multiscale (Frangi) vesselness with thinning and more than seven times faster than the inventors' own earlier method of blood vessel centreline extraction. Also, the centreline extraction was found to be highly accurate with a mean error of less than 1 mm on comparison to the corresponding geometric centres.

Overview

The preferred embodiments of the new method include blood vessel segmentation, centreline extraction and radius detection method, with the main aim of reducing the processing time by processing only the interested blood vessel regions of the image. The method lets the user select the vessels that need to be segmented and the method is performed only within the limits of these vessels.

The basic strategy of the method is to trace the vessel starting from a user provided seed point and direction. The tracing can continue to the end points of the continuous vessel or to the end of the image data, or it can of course be stopped by the user at some earlier point. Figure 1 shows an overall process flow diagram of an algorithm for the method.

Figure 2 illustrates a visualisation of the method as applied to a blood vessel with a bifurcation

Initially, the user initiates a seed point, a direction seed and an approximate blood vessel radius at the seed point. The seed is added into a seed list where the future seeds created by the algorithm will also be added. A direction seed is a point that the user initiates in the intended direction of tracking. The vector connecting the two points gives the tracking direction,

T = (x_s - x_d)l + (y_s - y_d)j + (z_s - z_d)k (1)

where T is the tracking direction, (x_s, y_s, z_s) is the seed point, (x_d, y_d, z_d) is direction seed and t is the unit vector along the tracking direction.

Using these inputs, the algorithm finds the radius and the local threshold at the 2D cross-section of blood vessel at the seed point (see below). Next, the radius and local threshold is used in the trunk analysis for finding 2D cross-section of the blood vessel. On analysis of the 2D cross-section, the centre, circleness output and next possible centre along the unit tracking direction of the blood vessel are calculated (see below). The next possible centre is then set as the new seed and added into the seed list. The method can detect bifurcations working on the principle that there is a significant change in radius at the blood vessel bifurcation compared to the neighbouring blood vessel trunks, as shown in Figure 3, which is an illustration on sudden change in radius at the blood vessel bifurcation compared to radius at the blood vessel trunks, where D_T1, D_T2 and D_T3 are the diameters at different vessel trunks, and D_B is the diameter at vessel bifurcation.

On detecting a significant change in radius, the process is passed into bifurcation analysis. In bifurcation analysis, a 3D modified multi-scale vesselness is applied to the 3D volume around the bifurcation and the output is used to find new centres at the bifurcation (see below). These new centres are then passed into the seed list as new seeds to continue the loop and the algorithm loop is terminated when the new seed found, is outside the blood vessel.

Finally, the whole 3D segmented volume of the blood vessel is found by combining all the previously found circleness and vesselness outputs. The circleness outputs include all the 2D circleness filtered images at each and every cross-section of the blood vessel trunk. The vesselness outputs include 3D multi-scale vesselness filtered images at all the blood vessel bifurcations. A neighbourhood around each 3D voxel position of the circleness output is also checked to see if the neighbouring positions are inside the vessel region and added into the final 3D blood vessel output. This eliminates the possible gaps on combining 2D cross-section outputs to form a 3D vessel output.

Radius Estimation

In the radius estimation process, the 2D cross-section image of the blood vessel is determined at a seed point and the radius at that cross-section. For determining the 2D cross-section image, first the blood vessel direction and vessel cross directions must be calculated. At the initial step, the tracking direction t is set as the vessel direction and vessel cross directions are set as cross vectors of t. However, for more precise structure analysis, calculation of vessel direction and vessel cross directions are needed.

Eigen analysis of the Hessian matrix at the seed point gives the principal directions of curvature of the vessel, as shown in Figure 4. The Hessian H of the image I describes the second-order structure of local intensity variations at each voxel. (See, for example, A. F.

Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, "Multiscale vessel enhancement filtering," in Medical Image Computing and Computer-Assisted Interventation— MICCAI'98. Springer, 1998, pp. 130-137; and Y. Sato, S. Nakajima, N. Shiraga, H. Atsumi, S. Yoshida, T. Koller, G. Gerig, and R. Kikinis, "Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images," Medical image analysis, vol. 2, no. 2, pp. 143-168, 1998, both describing the use of the Hessian matrix in 3D analysis.) The derivative computation of the Hessian matrix is combined with Gaussian G(o) convolution, for a scale σ which varies accordingly to the width of the blood vessel. For this purposes, previous radius is set as σ. The Hessian matrix H can be formulated as

d²9a d²9a a² gai

I * I * I *

d²x dxdy dxdz

H(x, y, z; σ) / * / * d²9o / * d²g_a

dydx d²y dydz (3) / * d²9a / * d²9a / * a²g_a

dzdx dzdy d²z \

where, g_a(x, y, z) = +y² +z²)/2a²

By eigen analysis of the hessian matrix H , the eigenvalues Λ_ΐ5 λ₂ and λ₃ are calculated. On sorting lAJ < \λ₂ \≤ \λ₃ \ , λ₁ is the eigenvalue along the direction of the vessel, while λ₂ and λ₃ are the eigenvalues along the vessel cross section. Their corresponding eigenvectors are v , v₂, and v₃. Accordingly, v is set as the new tracking direction or vessel direction, and v₂ and v₃ are set as the vessel cross directions. Vessel cross directions, v₂ and v₃ , are used to find an interpolated 2D blood vessel cross-section image at the seed point, similar to the one shown in Figure 5a. An interpolated cross-section image is a more precise representation of the cross-section, rather than cross-section images from discrete locations of original image.

Calculating radius of the blood vessel at the cross-section requires determination the borders of the blood vessel. As the blood vessel intensity is different from that of its bordering pixels, the edges can be detected by a Canny edge detection filter, as shown in Figure 5b. Along the vessel cross directions, two lengths of the cross-section are found by propagating outwards from the seed point towards the boundaries. The radius is estimated from the largest length, so as to consider more elliptical shapes. An average of all the edge intensities provides the local threshold at that blood vessel cross-section.

Trunk Analysis

Here, the blood vessel trunk cross-sections are analysed to find the blood vessel cross-section filtered image and the centre of the cross-section for finding the next possible centre along the blood vessel. The blood vessel cross-section filtered image is an image where the cross-section is enhanced, resulting in a Gaussian profile where the peak is the centre of the vessel cross-section.

A new circle structure enhancement filter or 'circleness' filter is formulated for getting a blood vessel cross-section enhanced image. The circleness method is based on 2D eigenanalysis of 2D Hessian matrix obtained from the cross-section image. Similar to Eq. 3,

The eigenvalues 1₂DI ^and ^-202 calculated from eigenanalysis of the two dimensional Hessian matrix H_2D are used for analysing the structural information of the cross-section image. The circle enhancement ratio using the sorted eigenvalues, |A₂DI I ≤ 2D2 I is formulated as

The circle enhancement ratio ^L R gives higher values at the centre of the cross- section compared to pixels near the border. The final equation for circleness filter using the circle enhancement te

where k affects the rate of increase of Gaussian profile of the circleness C, and the structureness term S,

reduces the effect of noise in the filter.

Figure 5a shows a 2D cross-section input image with blood vessel connected region; and Figure 5b shows the 'circleness' filter output of 2D cross-section image. For these calculations, the circleness C is only applied at connected region inside the blood vessel cross-section. This eliminates the enhancement of possible neighbouring blood vessel cross- sections or other circle-like objects that fall within the 2D cross-section image. Figure 6b shows an example of circleness filter applied to the blood vessel connected region of the 2D cross-section as shown in Figure 6a. The circleness output gives a Gaussian-like profile, where the peak of the output is the centre of the blood vessel cross-section.

From the 3D position of the centre in the original input image, the next possible centre or the seed point for the next cross-section along the vessel direction of blood vessel is calculated using the local vessel direction. The new seed is then saved into seed list to continue processing.

Bifurcation Analysis

On detecting a significant change in radius between neighbouring blood vessel cross- sections, the algorithm moves to a bifurcation analysis. In this example a 3D structural analysis is performed to understand the different vessel tracks being bifurcated. The use of a 3D analysis will of course require more data processing time than a 2D bifurcation analysis, as described below but this may be balanced by increases in accuracy. It has been found that at a bifurcation a 3D analysis can produce results that are sufficiently more accurate to justify the extra processing burden, in some situations. For example, when there is a trifurcation or several bifurcations close together then the 2D approach may not have sufficient 'resolution' to correctly follow the vessel paths without user input. This is because there is inevitably a volume between each 2D 'slice' of the vessel that is not taken into account. With a 3D analysis each volume that is analysed can overlap, which means that all available data can be used.

In the 3D bifurcation analysis, initially a sub-volume is created around the bifurcation and an eigenanalysis of the Hessian matrix H from Eq. 3 at the sub-volume gives the directions of curvature. After sorting eigenvalues lAJ < \λ₂ \≤ \λ₃ \ , the modified multiscale vesselness equation from previous work [1 1] is used for enhancing all the vessel-like structures in the image. The modified vesselness at a single-scale is expressed as,

A structureness term S 3D '

^S3D = ^ -ι² + -2² + -3²

is then added to the vessel enhancement equation, Eq. 8, to reduce noise in the enhanced image. The new equation using the structureness term is,

^ (' - S )E'. - ¾ - ¾) (i - ^) (10)

A second scale dependent weighing term is multiplied to the vesselness equation, as Eq. 8 gives exponentially reduced vessel peak values as the scale is increased. So the new vesselness term for single-scale is

V = \°' ^{if V}° ^{≤ 0}' (1 1 ) [ν_σε^ασ, othervise

where a is a constant set to 0.5, based on experiments and σ is the chosen scale.

A multiscale analysis is performed using a range of scales between o_min and o_max for detecting vessels of varying width. The multiscale modified vesselness equation is

V_mui_ti = max_{amin≤(J≤(Jmax} V{a) (12) where a_max is half of the maximum length of bifurcation cross-section and a_min is set less than one third of o_max , to consider all possible vessel ranges.

By applying modified multiscale vesselness, two distinct peaks representing different bifurcating parts of the vessel are found. Figure 7 shows (a) the 2D cross-section images at bifurcation input and (b) the corresponding modified vesselness output at the bifurcation cross-section.

At each of the peak positions, an eigen analysis is performed using the Hessian matrix Eq. 3, to find the directions of the bifurcating vessels. The vessel directions are then used to find the next possible centre or seed locations at each of the bifurcating vessels. The new seeds are saved into the seed list and the algorithm loop continues until seeds go out of the connected blood vessels or outside the image region.

Results

The proposed method, using a combination of 2D for vessel trunks and 3D when bifurcations are detected, was evaluated on five contrast enhanced MR angiogram images and three synthetic blood vessel-like images. In the medical images, three were obtained from publicly available DICOM medical image database of OSIRIX software and the other two were obtained from the Institute of Medical Science and Technology, University of Dundee. In Figure 8, images (a), (f), (k), (p) and (u) are the maximum intensity projects along the anterior-posterior direction of the five contrast enhanced MR angiogram images. These evaluated medical images are very relevant for catheter based interventions. The synthetic images were made from an extension to the traditional Lindenmayer system (L- system) that generates synthetic 3D blood vessels by adding stochastic rules (for example, as discussed in M. A. Galarreta-Valverde, M. M. Macedo, C. Mekkaoui, and M. P. Jackowski, "Three-dimensional synthetic blood vessel generation using stochastic l-systems," in SPIE Medical Imaging. International Society for Optics and Photonics, 2013, pp. 86 6911-86 6911; and X. Liu, H. Liu, A. Hao, and Q. Zhao, "Simulation of blood vessels for surgery simulators," in Machine Vision and Human-Machine Interface (MVHI), 2010 International Conference on. IEEE, 2010, pp. 377-380.

All processing was performed on an Intel(R) Core(TM) i7CPU @2.93GHz with 8GB RAM for consistency. The method is compared with multiscale vesselness (as described in A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, "Multiscale vessel enhancement filtering," in Medical Image Computing and Computer-Assisted Interventation— MICCAI'98. Springer, 1998, pp. 130-137) combined with 3D thinning (as described in T.-C. Lee, R. L. Kashyap, and C.-N. Chu, "Building skeleton models via 3-d medial surface axis thinning algorithms," CVGIP: Graphical Models and Image Processing, vol. 56, no. 6, pp. 462-478, 1994) and also with an earlier enhanced 3D method for blood vessel centreline extraction that preceded the currently described method (see R. P. Kumar, F. Albregtsen, M. Reimers, T. Lang0, B. Edwin, and O. J. Elle, "3d multiscale vessel enhancement based centerline extraction of blood vessels," in SPIE Medical Imaging. International Society for Optics and Photonics, 2013, pp. 86 691X-86 691X). Figure 8 shows a comparison of results obtained with Frangi and Lee, the earlier enhanced 3D method and the proposed method on the medical images. Figure 8 also shows the blood vessel segmentation output from the proposed method for the different MR angiogram images. All methods were coded using C++ and ITK libraries (as described in L. Ibanez, W. Schroeder, L. Ng, and J. Cates, "The itk software guide," 2003).

Table 1 shows the processing time for centreline extraction of blood vessels by different methods. All the images discussed have different image sizes and number of bifurcations, but have the same voxel size of 1 mm x 1 mm x 1 mm. In Table 1 , images 1-5 are medical images of blood vessels and images 6-8 are synthetic blood vessel-like images. On average, the proposed method is shown to be more than twenty times faster than the Frangi and Lee method and more than seven times faster than the earlier enhanced 3D method.

Table 1 - comparison of processing time

Data Image Size Frangi and Lee Prior art enhanced Circleness and 3D vesselness 3D vesselness vesselness method

Image 1 384x384x82 386s 1 1 1 s 21 s

Image 2 352x384x95 240s 36s 9s

Image 3 320x384x96 371 s 48s 13s

Image 4 356x330x124 403s 149s 13s

Image 5 272x499x88 197s 88s 24s

Image 6 300x300x300 581 s 273s 15s

Image 7 150x150x150 72s 34s 7s

Image 8 225x225x225 255s 108s 16s

The evaluation of the centreline extracted using the proposed method was done against the geometric centre. The geometric centre at each blood vessel cross-section was estimated by applying Hough circle detection and calculating the centre of the detected circle, which is the blood vessel cross-section (see, for an example, E. Davies, "A modified Hough scheme for general circle location," Pattern Recognition Letters, vol. 7, no. 1 , pp. 37- 43, 1988). Figure 9 shows the mean and standard deviation of centre estimation error of the proposed method. Figure 10 shows the percentage of bifurcating vessels not detected at various bifurcating vessel radius. The figure clearly shows that the accuracy of the algorithm in detecting bifurcation increase with increase in radius.

As an alternative to the process of Figure 1 , Figure 1 1 shows another possible method. This is similar but in place of 3D vesselness for the bifurcation analysis the method uses 2D circleness fitting multiple circles. There are some other differences as well, including an alternative way to detect the bifurcation. It will be appreciated that various elements of this second method, such as the bifurcation detection steps, could be used in place of the equivalent parts of the first method described above.

The method of Figure 1 1 focuses on segmenting the 3D connected blood vessels by tracking its cross-sections from a user-initialized seed to the ends of the blood vessel. The user initially initializes a seed point, a direction seed point and an approximate blood vessel cross-section radius at the seed point, for the proposed method to begin. Figure 1 1 shows the workflow of the proposed method for segmenting the blood vessels by analysis and tracking of the blood vessel cross-section images.

The method, is dived into four parts, namely, cross-section image, preprocessing for cross-section image analysis, bifurcation detection and circleness filter. The first part describes how the cross-section image is calculated at a seed point and its corresponding cross-sections. The second part describes all the preprocessing steps required for further cross-section analysis. The third part describes how a cross-section is classified as a bifurcation. The fourth and final part describes the novel circleness filter, which is used for circle enhancement and subsequent centerline extraction.

Cross-section Image

At the beginning, the tracking direction of the connected blood vessel of interest is estimated using the user-provided seed point and direction seed point. The direction seed point used to specify the direction in which the blood vessel is to be segmented. The tracking direction is the vector connecting the seed point and the direction seed point,

T = x_s - x_d)l + y_s - y_d)j + (z_s - z_d)k (13) i⁼ W\ ⁽¹⁴⁾ where T is the tracking direction, (x_s, y_s, z_s) is the seed point, (x_d, y_d, z_d) is direction seed and t is the unit vector along the tracking direction.

The tracking direction is used only at the seed point to understand the direction of segmentation. At the seed, the tracking direction is set as the vessel direction of flow and its cross vectors are set as the approximate vectors representing the cross-section of the blood vessel. However, an eigenanalysis of the Hessian matrix is calculated for precise information of the vessel cross-section vectors.

Eigen analysis can geometrically interpret the second order derivatives of an image at each point. The second order differential quantity for a volume / (x, y, z) with a Gaussian convolution g_a(x,y,z), is given by the indefinite Hessian matrix,

d²9a , d²9a d²9a

/ * I * I *

d²x dxdy dxdz

d²g_a d²g_a d²g_a (15)

H(x,y,z σ) = / * / * I *

dydx d²y dydz

d²9a d²9a d² 9σ

/ * / * I *

dzdx dzdy d²z

where,

g_a{x,y,z) =

and σ is the scale parameter set according to the radius of the blood vessel. At the initial seed point, the user defined initial radius is set as the σ, and for the subsequent cross- sections analysis, the previous cross-section radius is set as the σ.

Let eigenvalues of Hessian matrix H, be λ₁,λ₂ and λ₃, and their respective normalized eigenvectors be represented by v , v₂ and v₃. On sorting the eigenvalues as lA < \λ₂\ ≤ \λ₃\, the eigenvector v₁ represents the vessel direction and the eigenvectors v₂ and v₃ represent the vessel cross-section vectors. These vessel cross-section vectors are used to create the vessel cross-section image as shown in Figure 12. The vessel cross-section image is found by interpolating at discrete set of points along the 2D plane represented by the two cross vectors, v₂ and v₃.

Preprocessing for Cross-Section Image Analysis

Before moving into the cross-section analysis, some preprocessing steps have to be performed for calculating vessel cross-section border, radius and local threshold. These preprocessing steps help later on understand the structure at the cross-section and thus in segmenting the blood vessel itself.

The first and foremost step is to find the border of the vessel cross-section. This is found by applying a Canny edge detection filter on the cross-section image, as the blood vessel intensity is different from that of its surroundings. 2D Gaussian smoothing filter is applied at the beginning of the Canny filter to smooth out noise, where the variance is proportional to the square of previous radius or the initial user-set radius. Figure 5 shows the border of a blood vessel cross-section after applying the canny edge detection method. Each border pixel obtained through canny represents the discrete border contour of the vessel cross-section.

Depending on the cross-section being perfectly circular or a bit elliptical, there might be difference in the diameters calculated by moving along 0° or 90° vectors from the seed or center candidate. The center candidate is a possible center pixel at a cross-section

determined by moving along the vessel direction from the previous center position. The radius of cross section is found by taking the half of the highest of the diameters. Figure 5 also shows diameters found by moving along the different angles. Finally, the local threshold is found by averaging all the blood vessel cross-section intensities along the border found through canny.

Bifurcation Detection

After performing the preprocessing steps, each of the vessel cross-section images is checked to determine if it is a bifurcation. The contour of vessel bifurcation cross-section is very different from the contour of the vessel trunk cross-section. Figure 13 shows the difference in the shape of the vessel cross-sections between trunk and bifurcation. By calculating the shape descriptors (circularity and radius variance) bifurcating vessel cross- sections can easily be identified.

Circularity is a measure of compactness of a structure or how circular a given contour is. It can be defined as

Circularity = P² /A (17) where P is the Perimeter of the vessel cross-section shape and A is the area of vessel cross-section.

Radius variance is the variance in distance from the center candidate to the border points or vessel cross-section contour. Let B (n) represent the boundary contour points and cp be the candidate center point. Then, the set of distances from the center to the boundary contour points is D = {dt₁, dt₂, dt₃, ... , dt_n}, where dt_{ = \\B(i)— cp \\ . Radius variance is defined as,

where a is a constant and D is the average distance.

Whenever there is a sudden increase above a certain threshold or if there is a sudden change in the shape descriptors: circularity and radius variance, the corresponding cross- sections are analyzed as a bifurcation cross-section. Fig. 13 illustrates how the shape descriptor values change as the cross-section changes from vessel trunk to the vessel bifurcation. The sudden change in the shape descriptors resembles the change in vessel cross-section border, signifying bifurcation.

Circleness Filter

On determining vessel cross-section to be part of a vessel truck or a vessel bifurcation, the process goes forward in enhancing the structure accordingly and finding the center of cross-section. Here a novel circleness filter is implemented enhancing circular-near images in the image and provide a good Gaussian profile for the output intensity. The filter adopts the scale space approach with the possibility of performing the operation in single or multiscale approaches.

In the proposed method, single scale approach is used to enhance the cross-section at vessel trunk and multiscale space approach is used to enhance the cross-section at vessel bifurcation. This makes the bifurcation detection step a very crucial step to determine whether to apply multiscale or single scale circleness to detect the centers or center accordingly.

The circleness filter or circle enhancement filter is based on 2D Eigen analysis on the 2D Hessian matrix computed at each pixel of the vessel cross-section image. Similar to Eq. 15, 2D Hessian matrix is,

where I_2D is the vessel cross-section image, g_2Da is the 2D Gaussian filter and σ is the scale parameter set according to the vessel cross-section radius, which is found at the preprocessing step.

Let λ_2Ό1 and _2D2 be the eigenvalues from the 2D Hessian matrix. These eigenvalues are used to understand the structural information of at pixel of the cross-section image. The circleness filter is formulated with the knowledge that both the eigenvalues will be high at the center of the cross-section as the cross-section is near circular in nature. After sorting the eigenvalues (|/L_2D1 | < U_2D2 I)_> a circleness ratio is coined (2°)

The novel circleness equation at a single scale using the circleness ratio C_R is

C = (1 - e-^(c«^fe)2) M - e 2cA (₂1) where k affects the rate of increase of the Gaussian profile of the circleness filter and the structureness term

reduces the effect of noise in the filtered output.

For more efficient processing, the circleness filter is applied only inside the vessel cross-section region. The region of interest is determined by applying the local threshold found at the preprocessing step. Figure 14(a) shows the single scale circleness filter output on vessel trunk cross-section image. The single peak obtained from the single scale circleness filter applied to vessel trunk cross-section is the center of that cross-section. Moving along the vessel direction v from the current center, the next possible center or the next center candidate for the next cross-section is found at the vessel trunk.

At the bifurcation cross-section, multiscale circleness filter is applied. This allows determination of multiple peaks, where each corresponds to a bifurcating vessel. The multiscale circleness is formulated as,

Cmuiti = max C{&) (23) στηϊη-^σ-^στηαχ

where is o_min is the minimum radius and o_max is the maximum radius. o_max is set as the radius of the current cross-section calculated at the preprocessing step and a_min is set as one-third value of a_max.

Figure 14(b) shows the multiscale circleness filter output with multiple centers at the vessel bifurcation cross-section image, where each center corresponds to a different bifurcating vessel. Each of the centers found at the bifurcation cross-section is set as a new seed for the whole process to start again.

At vessel cross-section, the immediate surrounding 3D neighborhood is also checked and added as part of the vessel, if they fall within the local vessel threshold. This helps in reducing the gaps that might be caused by processing 2D slices along a 3D blood vessel. The vessel cross-section tracking finally stops, when the newly found center candidate falls outside the connected blood vessel region.

Claims

CLAIMS:

1. A method for segmentation of tubular organ structures in three dimensional (3D) image data, the method comprising:

extracting a two dimensional (2D) representation of a cross-section of the structure at the seed point;

estimating a radius of the structure at the seed point;

applying a circle structure enhancement filter to the two dimensional representation; thereby identifying a centre point of the structure at the seed point;

extrapolating along a direction of the structure to find another point on the structure, which is defined as a new seed point; and

repeating the method to thereby trace the direction and radius of the tubular organ structure.

2. A method as claimed in claim 1 , comprising: determining eigenvalues and eigenvectors for the structure at the seed point by eigenanalysis to thereby find principle directions of curvature of the structure; and extracting the two dimensional (2D)

representation of the cross-section of the structure at the seed point based on the results of the eigenanalysis.

3. A method as claimed in claim 2, wherein the eigenvectors provided by the eigenanalysis are used for finding the radius of the structure in the two dimensional representation.

4. A method as claimed in claim 2 or 3, wherein the eigenanalysis at the seed point is based on the three dimensional Hessian matrix combined with a Gaussian convolution for a scale varying with the width of the structure.

5. A method as claimed in any preceding claim, wherein the circle structure enhancement filter is applied only to image data identified as being within the structure cross- section.

6. A method as claimed in any preceding claim, wherein the circle structure enhancement filter makes use of an eigenanalysis of the two dimensional representation.

7. A method as claimed in claim 6, wherein a circle enhancement ratio is defined based on the two eigenvalues for the two dimensional representation, wherein the circle enhancement ratio is the modulus of the ratio of the sum of the eigenvalues divided by the difference between the two eigenvalues such that a peak value of the circle enhancement ratio hence indicates a centre point for the structure.

8. A method as claimed in any preceding claim, comprising detecting a splitting of the structure and performing a bifurcation analysis when it is determined that the structure is splitting.

9. A method as claimed in claim 8, wherein a splitting of the structure is detected by measurement of one or more of: radius of the structure cross-section, circularity of the cross-section, or radius variance.

10. A method as claimed in claim 8 or 9, wherein the bifurcation analysis comprises determining multiple centre points for the structure at the point of bifurcation, such that subsequent propagation of the trace of the structure will have multiple seed points and can hence follow multiple paths.

1 1. A method as claimed in claim 10, comprising reverting to non-bifurcation analysis for each of the multiple paths, with only a single centre point being used for subsequent seed points in each path, when it has been determined that the splitting of the structure is complete.

12. A method as claimed in claim 1 1 , wherein the multiple centre points determined by applying a circle structure enhancement filter to identify multiple possible circle centres.

13. A method as claimed in claim 1 1 , wherein the multiple centre points are determined by a vesselness analysis, and the vesselness analysis is repeated until it is determined that the splitting of the structure is completed at which point tracking of the structure reverts to a two dimensional analysis using a circle structure enhancement filter.

14. A method as claimed in any preceding claim when used as part of imaging of three dimensional networks of tubular organ structures within the human or animal body, such as for segmentation of blood vessels.

15. A computer programme product containing instructions that, when executed, will configure a data processing apparatus to perform the method of any preceding claim.

16. A data processing apparatus configured to perform the method of any of claims 1 to 14.