Universality in the merging dynamics of parametric active contours: a study in MRI-based lung segmentation

Measurement of lung ventilation is one of the most reliable techniques of diagnosing pulmonary diseases. The time consuming and bias prone traditional methods using hyperpolarized H${}^{3}$He and ${}^{1}$H magnetic resonance imageries have recently been improved by an automated technique based on multiple active contour evolution. Mapping results from an equivalent thermodynamic model, here we analyse the fundamental dynamics orchestrating the active contour (AC) method. We show that the numerical method is inherently connected to the universal scaling behavior of a classical nucleation-like dynamics. The favorable comparison of the exponent values with the theoretical model render further credentials to our claim.

Ventilation analysis is an authentic way of diagnosing lung airway diseases.The ratio of the volume of ventilated (functional) portions of lungs to the total lung volume is known as lung ventilation, which is used in validating pulmonary drugs [1].The process involves two complementary magnetic resonance (MR) imaging modalities, the hyperpolarized helium-3 (H 3 He) imagery and the proton ( 1 H) imagery.Lung functionality, including the volume of ventilated lungs, can be obtained from the former modality while lung anatomic details, including total lung volume, are accessed through the latter [1].Since manual investigation of the MR imagery to compute lung ventilation is extremely time consuming, an active contour (AC) or snake based automated method has been proposed [1,2] to compute the total lung volume from proton MR imagery on a 2-D sliceby-slice basis [1].A snake is defined as a massless 2-D thin string (closed or open) that can move on the image domain driven by two types of forces, internal elastic forces and external image forces [3].Under the influence of these two forces an initial contour (snake) clings to image edges and delineates an object.A snake always gives continuous edges unlike any traditional edge detector (such as the Canny method [4]), thereby eliminating any post-processing steps to connect the detected broken edges.These two properties are particularly useful when the object outline is broken and noisy as in most of the 1 H MR imagery.The snake method of finding the object boundary relies on the initial snake placement inside the image.If a small initial snake is placed inside a lung cavity on the MR image, while growing, the snake may be stopped by the associated numerical artifacts and may not capture the actual lung outline [1].Starting with a larger snake may result in missing the lung cavity completely.A possible solution is to start with multiple non-overlapping small snakes inside the lung cavity and evolve (grow) them until they merge with each other and capture the cavity outline [1].During such a process, the growing snakes merge with each other into a single contour.This automatic merging of non-overlapping snakes is characterized by certain attributes: a) during the evolution process no two snakes overlap with each other, b) every snake stops evolving at the object edge as a single snake does during its course of evolution, and c) growing convex shaped snakes (e.g., circular or rectangular snakes) inside a convex object recovers the object boundary.Although merging of multiple snakes is experimentally verified in a multitude of cases, a concrete theoretical understanding of this merging snake approach remains a challenge.
The aim of this letter is to bridge this gap by showing that the underlying principle of multiple snake-merging is governed by a universal power law behavior which originates from an inherent nucleating structure.From a biomedical engineering perspective, the study of the efficacy of the lung cavity delineation method is crucial for robust clinical application.The lung cavity segmentation by merging snake method involves three steps: a) initially small non-overlapping contours are placed inside the lung cavity, b) generalized gradient vector flow (GGVF) fields [2] are computed with a Dirichlet boundary condition on the initial circular snakes [1], and c) all the snakes are evolved simultaneously and independently of each other with the GGVF force field as the external force for the snakes.This automated lung cavity segmentation is attractive for a number of reasons.While other merging snake algorithms, such as the one proposed by McInerney and Terzopoulos [5] is computationally non-trivial compared to the original snake evolution algorithm of Kass et al. [3], the merging snake algorithm by Ray et al. maintains the same computational simplicity of Kass et al. algorithm.Also, based on the position of an initial snake, the rigidity parameters of the snakes can be varied, so that on one hand delicate high curvature features, such as costophrenic angles, can be accurately captured, while on the other hand, snakes can be made sufficiently stiff in order to avoid capturing artifacts.Fig. 1 shows multiple snake initialization, evolution, merging and delineation of lung boundary by the Ray et all method.
In order to map the merging snake scenario to statistical thermodynamical systems, we first provide the mathematical background for a parametric active contour or snake.A snake is a curve C(s) = (p(s), q(s)) defined by the parameter s ∈ [0, 1].The snake is evolved in such a way that it minimizes the energy functional [2] where the first two terms give the internal energy of the snake (α, β ≥ 0) and E ext represents the external energy added to the system.C ′ and C ′′ are the first and second derivatives of the snake with respect to s.An example of external force for the snake is the gradient force: , where I(x, y) is the image pixel intensity.In general, in the presence of an external force V(C(s, t)), the time dynamics of such a snake is given by [ Note that the snake in eqn.( 2) is now a function of time t as well.The stationary solution of eqn.( 2) corresponds to a snake that minimizes the energy functional.Ray et al proposed an external force V (x, y) obtained by solving the following partial differential equation (PDE) applying Dirichlet boundary condition [1]: where g(α) = exp(−Kα) and f (x, y) = −E ext (x, y), K being a tunable parameter controlling the smoothness of the external snake force field.Here we study the resultant dynamics due to the evolution of a number of such snakes, defined by the above system of forces.
When multiple snakes are evolved inside the desired closed boundary, the growth algorithm confirms that they all finally merge into a single snake after a finite time.This merging of snakes is basically a topological effect and naively the phenomenology reminds one of "nucleation" as seen in classical first-order thermodynamic systems [6].In our case, we allow a finite number of GGVF snakes, each with a finite starting radius, to evolve in a 2-D plane and then numerically evaluate a few measurables -the nucleation time (NT), the bounded area (BA) after nucleation has occurred, the critical radius (CR) at the time of nucleation and also the nucleation rate (NR), all as functions of snake evolution time.The respective quantities are defined as followsnucleation time is the time required for all the snakes to merge together as a single unit, bounded area is the sum of the areas of the initial snakes before complete merging and the area under the single snake after merging, critical radius is the radius of curvature of the merged structure and nucleation rate is the ratio of the number of snakes to the bounding area, before the merging has actually taken place.One should note that by complete merging we refer to the critical phase when all the initial snakes merge together for the first time.
We perform two numerical experiments with merging.In the first experiment, we start with a circular binary image of radius R containing N number of circular snakes, each of radius r < R, randomly distributed inside.As described before, the initial snakes are driven by GGVF forces and they maintain a non-overlapping dynamics.In another numerical experiment, we vary the radii of the smaller circles and later also vary their total number.Both numerics are repeated with varying sizes of the initial domain R. The enclosed figures 3 and 4 show the variation of the bounding area and the critical radius against time respectively, in non-dimensionalized units, in loglog plots.We find that all the graphs show excellent scaling in a broad domain, with the nonlinear zones referring to the saturation limits in each case.In the figures shown, the bounding image radius is 80 units while the corresponding starting snake radius is 5 units.We have simulated using different number of initial snakes, indicated in the legend of the figures.The plots show power law variations of each of the measured quantities with time, i. e. y ∼ x α and the two relevant exponents, that of critical radius and bounding area, have the values α BA ≈ 0.9 − 1.0, α CR = 0.25.We have used multiple combinations of the parameter values but the exponents remain universal, that is they invariably follow a power law statistics.
To analyse the results further, we begin with the hypothesis that the phenomenology of the whole process, that of various snakes merging together under certain boundary conditions in the presence of suitably defined external forces, is fundamentally similar to that of a domain growth process as encountered in the realms of classical nucleation.In what follows, we try to evaluate the time rate of evolution of a system of snakes before and after all the snakes have completely merged with each other.We start from a model of solutes trying to nucleate in a solution and then compare the results with the merging mechanism of snakes we find in the AC method.
Our starting description is that of the diffusive growth of two dimensional spherical droplets (circles, in 2-D) in a Lifshitz-Slyozov [6,7,8] type of continuum theory.We start with the thermodynamic definition of the work ∆W required to form a nucleus from an aggregate of solute particles where δ(E − T 0 S) is the increment in the free energy of the system and δ(P 0 V ) is the associated external work done.Optimizing this total work factor, one can show [6] that the critical radius of a nucleating system, this being the rate of merger of snakes in AC method, is given by Here µ ′ 0 and s ′ are the chemical potential and molecular surface area of the nucleus, β is the surface tension and µ ′ is the chemical potential of the solute in the solution.As more and more solutes start nucleating i. e. snakes merge, the solution approaches the critical limit of supersaturation and in the steady state, the critical radius grows as where c(R) represents the spherically symmetric concentration distribution of snakes around a nucleus of radius a, D being the diffusion coefficient.Following the treatment of Lifshitz-Slyozov [6,7], one can show that if R CR (0) is the critical radius at the beginning of the merging process, then for a diffusive nature of relaxation, the radius follows a dynamics We now define the dimensionless quantities x(t) = RCR(t) RCR(0) , u(t) = R(t) RCR(t) and τ = 2 log x(t).The last one increases monotonically from 0 to ∞ as the time t increases likewise.Combining ( 6) and (7) in 2D, the dynamics is given by where γ = dt RCR(0)dx > 0. A linear stability analysis of the above nonlinear equation gives a fixed point at γ 0 = 27 4 .We are interested in the dynamics around an ǫ-neighborhood of this point γ (ǫ → 0 as γ → ∞), then near the critical point (u 0 = (γ/2) 1/3 = 3/2) the merging snakes follow a dynamics defined by The time variation of the merging nuclei (snakes) is Thus the bounding area (∼ x 2 (t)) of the merging snakes grow at the rate of t which is roughly speaking our numerical estimate (0.9-1.0) also (Fig. 3).However, the above analysis is only valid before complete merging of the snakes has occurred.In the merged phase when the resultant asymmetrical structure continues growing finally to coalesce with the bounding image, the system dynamics is modified.To analyse the situation, we start from eqn. (7).The equation may alternatively be represented as T in the above equation is the temperature of the solution in the equivalent thermodynamic problem which in our case, is a measure of the average surface energy E, E ∝ T , of the system.Evidently, in 2D, E goes as 2πRβ.Plugging these values in the above equation, we see that after complete merging of the snakes has taken place, the effective radius of curvature of the resultant structure grows as R CR (t) ∼ t 1/4 (Fig. 4).Once again our numerical result is in exact harmony with the theory.Our above analysis, both numerical and analytical clearly suggests that below the apparently simplistic level of the GGVF application, the system dynamics has a more fundamental symmetry.This symmetry comes from the fact that the GGVF method lies in the same universality class as that of a classical nucleation model.This first-order critical response of the system is what provides the subtlety of the underlying physics.There is, however, a notable shortcoming with the GGVF technique, that it does not allow us the liberty of starting with an initial condition at an arbitrary location.If the snake starts at a position in which a major portion of the initial snake is outside the desired boundary or vice versa, then the snake driven by GGVF will not converge to the actual boundary.Although it can be shown by the Reed-Simon's theorem [11] that a convex set (a circle, say) growing within a larger convex set (a larger circle or rectangle, say) will always merge with the outer boundary under the action of isotropic driving forces, to the best of our knowledge, no such mathematical lemma exists for a convex set growing in a concave set or vice versa.
In summary, our achievement in the present paper has been to analyze the physical foundation of the GGVF merging technique and to provide answer to the rather puzzling question as to why it works so accurately.In the process, we have shown that the answer lies in the general scaling behavior of the underlying nucleation dynamics, defined by proper scaling laws.This clearly indicates that an active contour system is in the same universality class as the nucleation model we considered.Our results, in unison with probable biomedical applications, are expected to inspire further studies in the understanding of lung-based diseases.

FIG. 1 :
FIG. 1: Evolution of multiple snakes in a 2-D slice of the lungs by the Ray, et al [1] method.