Contour evolution method for precise boundary delineation of medical images

ABSTRACT


INTRODUCTION
A basic problem in computer vision is to understand the high level of information contained in an image or a series of images.Many of the tasks in computer vision may look trivial when performed by a human.Tasks such as detecting and distinguishing faces and objects in images, or classifying hand gestures in videos are some of the examples of a problem that humans can do with relative ease but could be a challenge for a computer to perform.More complex tasks such as detection of anomalies in medical images pose an even greater challenge due to its high accuracy requirement and high-risk failure consequences.
Understanding medical images requires a specific set of skills and many years of experience working as a medical practitioner such as a radiologist.Thus, in building a computer-aided diagnosis (CAD) system, it is necessary to adopt the working procedure of the relevant radiologist and model it as an algorithm that computers can execute.This process is a complex task and may require an application of several image processing and computer vision techniques that as a whole complement each other.Image segmentation and boundary delineation are two important tasks in many modern image analysis processes including computer-assisted diagnosis through medical image analysis.Image segmentation has been successfully applied in many CAD systems working on many different types of medical images.Traditionally, image segmentation is achieved using clustering methods, such as genetic algorithm [1] and fuzzy-means [2], of a number of manually designed low-level features such as pixel values distribution and histograms of gradients.Stochastic methods are also a popular approach in image segmentation.A regression segmentation framework to detect vascular anomalies in cardiac magnetic resonance imaging (MRI) images by delineating boundaries of bi-ventricle from was proposed [3].In this framework, a regression model has been trained automatically on a deep belief network by using extracted DAISY feature [4] as input and using automatically generated boundary points as labels.The method was reported to yield high performance when tested on 2,900 images taken from 145 clinical subjects.
One of the major challenges in applying automatic image segmentation in medical images such as CT and MRI is the imperfection in the imaging process which can often result in inconsistent contrast and brightness levels, and low image sharpness and vanishing boundaries.Recent advances in deep learning produce vast improvements in the quality of image segmentation.As a result, research in this field has shifted its focus to using more contemporary techniques such as convolution neural networks [5,6] and random forest [7].In one of our earlier publications [5], we proposed a novel method for image segmentation and boundary delineation to assist diagnosis of lumbar spinal stenosis in MRI images using SegNet [8].
It is noted the accuracy of deep learning classifiers may also depend on the accuracy of the ground truth labels that are used to train the models.Producing label images manually is a laborious and lengthy process.It is very unlikely that medical experts' time is used to carry out this activity for hundreds or even thousands of images.Even if their time is used to label a few, they may not have the necessary skills in using digital applications to produce accurate segmentation.A practical solution to this is to have the medical experts to train several skilled image labellers on the important principles and what features to look for and assign the labellers to produce the ground truth label images for the classifier training.However, even with this approach, there is still a high level of probabilities of inaccuracies, especially along region boundaries [9].
To improve the segmentation accuracy along region boundaries, an active contour model can be applied to the label images.However, in our study of the important principles on image labeling provided by medical experts, we found that the existing active contour models have a number of drawbacks when are used to label medical images [10].Thus, we develop a new approach to contour evolution that is more intuitive but shares some common principles with the active contour model method.Through analysis of our experimental results, we show that our method has several properties that help overcome the problems in the existing active contour models.
This paper starts by discussing the principle of active contour models which is followed by the description of the problems in applying these to medical images.The design of our solution is provided in section 3 followed by the analysis of the experimental results in section 4. We use the lumbar spine images that we have collected [11] and their associated label images [12] in our experiment.This dataset contains axial MRI scans of the last three intervertebral discs of 515 patients.

ACTIVE CONTOUR MODEL
Active contour model has been used to segment and delineate boundaries in different types of medical images such as echocardiographs [13], MRI [14], and Computed Tomography (CT) [15].It is a popular class of image segmentation and boundary delineation method due to its ability to fit a curve to an object boundary by iteratively expanding or contracting its boundary estimate.The iterative process starts by defining an initial contour  0 and at each iteration , the contour   evolves from the previous iteration  −1 such that   =  −1 +   × ∆, where   is the contour evolution speed.The contour evolution speed is not a singular entity but rather a multi-value one as it is defined adaptively for each point  along the contour position at each iteration.This is illustrated in Figure 1.
The contour evolution speed, i.e., the red arrow illustrated in Figure 1, (which is often also referred to as Force denoted as ), is the movement of each point along the contour.The force has direction perpendicular to the tangent of the curve.The value of  depends on local properties such as curvature, global properties such as shape, and other independent properties [16].One of the earliest and most popular active contour models is Snake [17].It can arguably be considered as the basis of almost all subsequent active contour models to date.In this model, the force is calculated to minimize an energy function  that is a function of the curve's internal factors e.g., the first and second-order derivatives of the curve as well as external factors (e.g., image dependent features such as image gradient or edges).The value of  is the integral of three components expressed in the (1):  The implementation of the Snake active contour algorithm uses the zero level set (ZLS) framework [18].This framework propagates a shape boundary perpendicularly to the contour by means of an Eulerian formulation.Instead of representing the boundary and solve the evolution parametrically, ZLS represents 2D contours implicitly as a set of points of where the value of a function crosses zero (either from positive to negative, or vice versa).This function is an everywhere-differentiable image function :  →  with domain    2 .The function  is a special image function that is derived from a label image.Given a binary image denoting a closed set of labelled pixels with two classes namely  and ¬, where  is the class of the object of interest in the image and ¬ is its complement class, the function  has a negative value at pixels that belong to  and a positive value at pixels which do not.Since  is differentiable everywhere, its value should be zero, or close to zero, at boundary locations between  and ¬ regions.The function is often illustrated in 3D with the  and  axes making up the image plane and the  for the value of the function.The intersection of the function with the  = 0 plane will create an outline of the boundary of the two regions in a segmented image.This set of points is called zero level set and is denoted as .Contour evolution using the ZLS method is achieved by evolving the  function, hence the state of the function at iteration  is denoted specifically as   .Therefore, at any iteration , the zero-level set   = {(, )|  = 0} contains the intermediate solution to the boundary delineation problem.The process to construct  from a binary label image is illustrated in Figure 2.
The initial value of , i.e., the state of the function at the start of the iteration, denoted as  0 can be calculated from the initial segmented image.The initial segmented image is a binary image with pixel values of 1 for every pixel that is inside the shape and 0 otherwise as shown in Figure 2 (a).One popular method to construct  from a segmented image is by calculating the signed distance function (sdf) [19].For each pixel,  is calculated as the distance between that pixel to the closest point on .This may sound like an egg-and-chicken problem with regards to  and , but in practice, the  can be calculated using the Euclidian distance transform (Edt) of the binary segmented image.More precisely, the  is calculated as the difference of the Euclidian distance transform of the binary segmented image and the Euclidian distance transform of the inverse of the binary segmented image, i.e.,

𝑠𝑑𝑓 = 𝐸𝑑𝑡(𝐵) − 𝐸𝑑𝑡(¬𝐵)
(2 where  is the binary segmented image and ¬ is its inverse.A distance transform is an operator that is applied to every pixel in its input binary image to calculate the distance between the pixel to the nearest non-zero pixel [20].Using ZLS, the evolution of the contour is achieved through an evolution of this signed distance function.The first ,  0 is required as the initial value to start with.The relationship between the new  and the contour evolution speed  is given as the level set equation [16]: since its introduction, there are a number of variations of the ZLS framework.Some methods apply different energy minimization procedures whereas others use different image features when calculating . One of the most popular is the edge-based Geodesic Active Contours (GAC) [21].In this approach, the force  acting on the signed distance function is calculated as: where As the name suggests, the GAC method calculates the feature   using the gradient magnitude information of the input image.In particular, it is calculated per pixel as: where the feature used is this feature is more popularly known as the normalized gradient magnitude of the input image .
The operation (  * ) is a convolution of the image  with a Gaussian kernel   with width .

Drawbacks in using active contour models to segment medical images
Accurate segmentation of images, including medical images, with significant intensity inhomogeneity within the same region, is often difficult to achieve [22].In our study of several medical image types, topology, and geometry, and while guided by the relevant medical experts, we found that the existing active contour models have a number of drawbacks when they are used to segment medical images.Although in this paper we only focus on the segmentation of lumbar spine MRI images as a case study, our rationale and findings should apply generally.
Borrowing from the approach we took in our previous study when delineating this type of images [5] to detect lumbar spinal stenosis [23], there are five regions of interests in each image.They are the intervertebral discs (IVD), posterior element (PE), thecal sac (TS), area between anterior and posterior (AAP) vertebrae elements, and Other (OT).The drawbacks that we found are summarized as follow: − The true boundary of an object region may not be found by just using a single type of image features because it could depend on several different ones.For example, the boundary between the IVD with AAP regions is characterized by strong pixel intensity differences however the boundary between IVD and TS is not.Therefore, if the image gradient feature alone is used to improve the boundary of IVD it may result in its deterioration instead.While many different models of active contour algorithms have been proposed in the literature which use different types of image features, there is none, to the best of our knowledge, which allows multiple different ones to be adaptively applied along the curvature.
− The inclusion of all three components of speed functions in one calculation makes it harder to find the right combination to yield the best results.The effect of adjusting one control parameter may counter the effect of adjusting the others.In GAC, for example, generally using a large γ value pushes the contour harder towards an edge resulting in a contour evolution that moves strongly towards hard edges which will make it fits these edges rigidly.However, increasing the value of β could cancel that effect out because it penalizes high curvature.As a result, finding the right combination of parameter values is very tricky and can be application dependent.This phenomenon is illustrated in Figure 3, that shows the delineation of an unmodified label image (Figure 3 (a)) and the effects of applying the GAC algorithm using different combinations of γ and β values to the label image (Figure 3 (b)).The boundary of each region pair is color-coded to allow easier inspection.The figure shows that having low values of γ and β produces the worst results (top row and first four columns).The results are better when a large value of β is used (top row-fifth column).While fixing this β value, the effect of changing γ values does not produce significantly different results (last column).However, when the largest γ value is used the effect of changing β appears to reverse.This time lowering β value shows an improvement in the quality of the delineation.
On a further note, we notice that most active contour models are applied mainly to an image in its original size.This limits the precision of boundary points to half a pixel.Increasing the size by up sampling the image and its label  number of times, where  ∈ ℕ, can increase the precision to the nearest 2 −(+1) .However, this approach significantly increases the number of pixels to consider and our experience shows that active contour models can scale poorly with the number of pixels.Based on these factors, we specify three requirements that must be met for a suitable contour evolution method.They are 1) it must allow an application of evolution on a subset of contour as opposed to the entire contour.This will make it possible to tailor the evolution process using different image feature type for different parts of the contour and 2) it should decouple the different components of evolution function so that the boundary fitting can be optimized independently in each component space.Lastly, as the third requirement 3) the method must scale well with the increase in pixel numbers due to up sampling to increase the delineation precision.

THE PROPOSED CONTOUR EVOLUTION METHOD
Image segmentation is a process that clusters pixels into multiple regions.The most straightforward way to represent the regions in an image is by using a matrix of labels with the same size as the original image.As an example, an instance where a 4x4 image that is segmented into three regions is illustrated in Figure 4.The figure depicts the region boundaries as red lines between pixels that have different labels.

Boundary information representations
Boundary locations are naturally sub-pixels because it exists between adjacent pixels.They can be stored as a boundary grid,  which is a matrix whose size is twice that of the label image minus one.A version of the Boundary Grid that is used in [24] using the label image shown in Figure 4 (b) is shown in Figure 5.In the above example, cells with boundary information are shaded in red.Not all cells in a boundary grid can be boundary cells and in general only cells that have either 1) an even row number and an odd column number or 2) an even column number and an odd row number can be boundary cells.Furthermore, boundary cells in category (1) mark only horizontal boundaries whereas those in category (2) mark only vertical boundaries.Any of the cells that meet the criteria but are not boundary cells are illustrated as blank squares.
We made a small modification to this idea by including image feature values into these cells.The type of image feature chosen is a design decision and it is important to the accuracy and suitability of the algorithm to a given problem.However, in this paper, we use a similar type of image feature as GAC to ensure its direct relevance when comparing to the latter.We refer to this modification as the modified boundary grid, denoted as ′.The image function of ′ is given as: where   is the  ℎ element of  and ℒ = {1, … , } where  is the maximum number of regions in the image.The algorithm to construct our ′ is as follows: Let :  →  be the input image function and :  → ℒ be its associated label image function, both with domain    2 .At locations where ( 1 mod 2) ≠ 0 ∧ ( 2 mod 2) ≠ 0,  ′ contains the label information as prescribed in .At locations where ( 1 mod 2) = 0 ∧ ( 2 mod 2) ≠ 0,  ′ contains vertical gradients and at locations where ( 1 mod 2) ≠ 0 ∧ ( 2 mod 2) = 0,  ′ contains horizontal gradients.
In practice, the gradient values are calculated using the convolution of the input image  with a Gaussian kernel   of width  to minimize high-frequency elements in the derivative operation.Hence the horizontal and vertical gradient image functions  ℎ :  →  and   :  →  are defined as  ℎ = ∇ ℎ (  * ) and   = ∇  (  * ), respectively.The size of the horizontal and vertical gradient images are  × ( − 1) and ( − 1) × , respectively.To illustrate the process to construct our ′ let's assume that the horizontal and vertical gradient images of the input images we used previously are shown in Figure 6 (a) and Figure 6 (b), respectively.The resulting ′ is shown in Figure 6 (c).This representation is, however, neither efficient nor compact if we want to use it in an iteration to search for the coordinates of a specific region's boundary.To compensate this shortcoming, we also use a sparse boundary representation, denoted as   , as a look-up table to index pairs of neighbouring regions.This is a much more compact and efficient representation of the boundary information because any search or information retrieval operations that are carried out in this representation are carried out faster due to them being dependent only on the number of boundary pixels rather than the number of pixels.We can construct   from ′ by parsing the latter once, to gather all edge coordinates for each neighboring region pairs.Table 1 shows the sparse boundary representation of the Boundary Grid example we used earlier.Note that matrix subscripts are expressed as row and column order.

The contour evolution algorithm
In this section, we will discuss the algorithm to evolve the region boundaries using the two boundary representations described previously.Let  be a set containing the matrix subscript pairs of all the boundary points that we want to evolve using our chosen image feature.We construct  directly by querying   .For each element in , we find its location in ′, and based on the values of the coordinate we can ascertain the type of edge it is.A horizontal edge is evolved horizontally by evaluating the feature vectors in the horizontal direction.Likewise, vertical edges are evolved vertically.We need to decide the value of a parameter , where  ∈ ℕ and  ≠ 0, which is the search width (in pixel) of the evolution.The value of  affects the speed and accuracy of the evolution and typically we want to use a small number.Large  value may result in a label that is very different from the original estimate.The evolution of the boundary contour will be based directly on  and is done by altering the label contents of ′ as described in Algorithm 1 shown in Figure 7.Our new contour is then obtained directly by parsing the evolved  ′ while generating the new   .The algorithm is executed iteratively until the number of pixels in the Boundary Grid that change from one iteration to the next, averaged out over  number of iterations, converges below a specified threshold .

Subpixel boundary location
To further improve the precision of our method, we upsample our input and label images at the end of each complete evolution.We can repeat this for  number of times, where  ∈ ℕ, to increase the boundary precision to the nearest 2 −(+1) of a pixel as described in algorithm 2 shown in Figure 7.

Adaptive curve smoothing
We design our method to allow the decoupling of the different components of the contour evolution speed.To this effect, we apply an adaptive curve smoothing function at the end of the above contour evolution process.The rationale for this is to allow us to adjust the tightness of the curve to the feature and TELKOMNIKA Telecommun Comput El Control  Contour evolution method for precise boundary delineation of medical images (Friska Natalia) 1629 only concern about smoothing the results afterward.The curve smoothing is adaptive, in such a way that in areas where strong image features are present it applies a weaker smoothing function to preserve the edge than in areas where the image features are weak.We achieve this by means of variable width moving average method.Let  = ⟨  | ∈ ℕ and  < card()⟩ where   ∈  2 be a sorted sequence of set .We define the sort operation such that ‖  −  −1 ‖ for 1 ≤  < card() is minimized.We then apply a moving average on  at every point  along its curvature with variable half-width   ∈ ℕ such that: for ∀ that meets the   <  < card() −   requirement.The value of   is set between  min and  max and tied to the image feature at the location of the boundary points.In our experiment, we use the second derivative of the image function ′′ as the image feature and set   to  min when the ′′ is at its lowest value, to  max when the ′′ is at its highest value and linearly interpolated and rounded to the nearest integer in between.

EXPERIMENTAL RESULTS AND ANALYSIS
We tested our method on a dataset of lumbar spine MRI images that we have collected [11] and their associated label images [12].The dataset contains axial MRI scans of the last three intervertebral discs of 515 patients.We present in this section, one-in our view the best-example case because it contains all of the problems to solve as described in the last part of section 2. This example case is described and illustrated in Figure 8.The figure shows the color-coded boundaries of a manually created label image.The figure shows a number of problems with the boundaries.Areas labelled (a)-(d) are examples of gaps between regions that were meant to be neighbours.Area labelled (e) illustrates inaccurate boundaries that are a few pixels away from the true edge.And lastly, we note that these boundaries in general have high curvature which can be seen from its jaggedness such as those labelled (f).
The parameters that we use are as follows: image features  = ∇(  * ), the search width  = 2, the number of upsampling levels  = 2, the iteration averaging window  = 10, the standard deviation threshold  = 2, the min and max smoothing window half-widths  min = 5 and  max = 20, respectively.We applied our method to only evolve boundaries that we have prior knowledge to be gradient dependent while also improving the smoothness of the others.For comparison, we also apply the GAC algorithm [21] to the label image, sequentially for all regions except OT.We experimented with several combinations of ,  and  values to produce an acceptable compromise between accuracy and smoothness, and we decided to use  = 1,  = 160 and  = 500 to be the best.We also apply a morphological closing on the input label image prior to applying each technique to remove any small holes and gaps in the image.On average, the execution speed of our method compared to GAC at up sampling levels  = 0, 1 and 2 are 10.2, 8.6, and 5.2 times faster, respectively.We observed that the main reason why our method is significantly faster is that it carries out a shorter number of iteration before it converges.As can be seen in Figure 9, our method solves the bulk of the pixel changes in the first iteration resulting in a much quicker convergence than GAC.The contour evolution results of the two methods are shown in Figure 10.There are two notable improvements of our method compared to GAC.Our method primarily produces boundaries that are more tightly aligned with the true edge-marked as (e) in Figure 8 while at the same time produces a generally smoother contour than the latter.Note also the fact that GAC worsens the boundary between IVD and TS since this type of boundary is not characterized by strong intensity differences.To quantitatively measure the method's accuracy, the boundary-improved manually-labeled images were used to train a SegNet image segmentation model [5] which is then used to automatically label lumbar spine MRI images.The method was then subsequently applied to the segmented images to further improve their boundary.The method's performance in improving region boundaries can be quantitatively measured by BF-score semantic contour-based metric [25] of the boundary-improved automatically segmented label images with respect to its boundary-improved manually-labeled images.We used the with distance threshold values of 1, 2, and 3 pixels for this purpose.The results are shown in Table 2.For comparison, the same metric is calculated when GAC was used instead and when no modification was applied.The result clearly shows that our method is better at improving the boundary accuracy of the label images than GAC and when no modification was applied.It is worth noting that the implementation of our method can be parallelized hence its effectiveness can be further optimized by utilizing GPU computation.Using the criteria set in [26], our method is suited for GPU computation due to its high data parallelism and low memory usage and branch divergence.The only drawback in this respect is the iterative nature of the process which requires synchronization at the end of each iteration.

CONCLUSION
We have proposed a new method of contour evolution that is suited to improving manually segmented medical images.Our method solves two main problems in Active Contour Models namely 1) their inability to use different image features at different segment of the contour and 2) the interdependency of the parameters of the contour evolution speed which makes finding good and suitable combination parameter values a difficult task.We tested our method on lumbar spine MRI images and our experimental results show that it can, not only, improve the accuracy of the boundary delineation of the manually segmented images, by 10 percentage points as measured using the BF-Score semantic contour-based metric, but also produce visually accurate and yet smoother contour than the traditional GAC method.Ours finding should apply generally even though we only focus on the segmentation of lumbar spine MRI images as a case study.


ISSN: 1693-6930 TELKOMNIKA Telecommun Comput El Control, Vol. 18, No. 3, June 2020: 1621 -1632 1622 and feature() is the image feature value calculated at the contour location.The model uses three parameters to control the effect of each component on the contour evolution.Parameters  and  control the contribution of the internal factors whereas  controls the contribution of the external factors.Almost all of the subsequently proposed active contour models in the literature are based on this principle.Differences amongst them are mainly on the different ways of calculating and determining each of the three components and their parameters.

Figure 1 .
Figure 1.The evolution of the contour   (brown line) at time  from the previous contour  −1 (blue line) with the red arrows marking the contour evolution speed at each position on the contour

Figure 2 .
Figure 2.An illustration of the construction of the signed distance function from a binary label image; (a) input label image is used to develop it, (b) corresponding , (c) shows the same  as viewed in 3D with its zero level set marked in red

Figure 3 .
Figure 3. (a) Delineation of the original label image of a lumbar spine MRI image, (b) The result of applying the GAC algorithm using different combinations of  and  values; From left to right, the  values used are 0.1, 0.5, 1.0, 1.2, and 5.0.From top to bottom the  values used are 0.1, 0.5, 1.0, and 1.2.The value of  is fixed at 1 throughout

Figure 4 .
Figure 4. (a) An example of a4x4 image, (b) The image segmented into three regions

Figure 5 .
Figure 5. Boundary Grid of the label image shown in Figure 4 (b)

Figure 6 .
Figure 6.Gradient images of the input image shown in Figure 5 (a), (a) horizontal, (b) vertical, (c) the result of modified boundary grid

Figure 7 .
Figure 7. (a) Algorithm 1 to evolve the boundaries stored in a Boundary Grid and a Sparse Boundary Representation and (b) Algorithm 2 that iteratively up-samples the input and label images and makes calls to Algorithm 1 to achieve sub-pixel boundary evolution.

Figure 8 .
Figure 8. Color-coded boundaries of a manually created label image illustrating (a-d) gaps between supposedly neighbouring regions, (e) inaccurate boundary and (f) high curvature

Figure 9 .
Figure 9.A typical number of pixel changes between successive iterations, (a) proposed, (b) geodesic active contours

TELKOMNIKAFigure 10 .
Figure 10.(a) The contour evolution results using our method, (b) The contour evolution results using GAC

Table 2 .
BF-Score of the automatically segmented label images