1.INTRODUCTION

Magnetic Resonance Angiography (MRA) is increasingly used to provide volumetric information of vascular system. Vessel segmentation is one of demanding applications that has received a considerable attention Accurate assessment of MRA images requires that the vessel structures to be extracted from MRA data sets. Patient specific computerized 3D model of vascular has important applications in augmented realty based navigation for intervention and biopsy. While evolutionary schemes based on the level set theory have proved to be effective tools for vessel segmentation in high field MRA images, with current developments of lower field interventional MR scanners, they need to be robust in the presence of higher noise levels. One typical observation in level set based methods concerns with damaged intensity information of thin vessels that causes the evolution to stop before extracting the whole vessel. This is illustrated in Fig.1. Panel (a) is the maximum intensity projection from a portion of a 1.5 Tesla Time-Of- Flight MRA data set. The arrow indicates the pinching of a thin vessel caused by the image noise. However, in a higher magnetic field (3 Tesla) a smooth intensity pattern from the same vessel is obtained that clearly shows its extension. Such, noise "speckles" may stop the front evolution towards the thiner part, and consequently multiple distinct vessel fragments may be obtained. Masutani et al. [1] addressed this issue, using a mathematical morphology region-growing method. However the vessel structures remained un-interpolated between discontinuities.

In order to reduce the risk of front leakage to the background, a number of methods have been using shape constraints. For example topology constrained surface evolution has been proposed in [2], a method that uses 3D skeletons of the front to refine the spurious branches for iterative bifurcation and vessel segmentation. Also soft priors has been introduced in [3] for minimizing the leakage from noisy edges using a ball-filter which penalizes the deviations from tubular structures. Obviously, any geometric regularization of the segmentation process, with the ability to minimize the leakage is important for vessel segmentation.

2. RELATED WORK AND CONTRIBUTION

Currently, a number of techniques have been developed for vessel segmentation based on the advanced level set evolutionary methods. Lorgio et al [4] has proposed CURVES using image gradient strength information, and the surface minimum curvature as the smoothing term. A level set method is introduced by Vasilevsky [5] that integrates the directions of gradient vectors into the evolution equation so that the gradient flux through the evolving curve is maximized. Also capillary active contours is invented by Yan et al. [6], a method that is based on the capillary force acting on the free fluid surface through a capillary tube. In this research, it is assumed that vessels are curve-linear structures, this is the basic assumption as in [7] for enhancing the vascular structures. Based on this general assumption, a shape functional for effective geometrical regularization of tube-like structures is proposed. Our regularization method preserves cylindrical structures by imposing anisotropic front constraint in the level set variational framework. Our regularization method has three basic important properties:

1) Unlike the previous curvature based smoothing methods, the solution of our geometric energy minimizing scheme is not always a shrinking surface, but instead it can extend toward local "meaningful" surface features.

2) Since leakage develops isotropic structures, by applying anisotropic front constraint, it also reduces the risk of leakage. 3) It is a curvature dependent flow and therefore has smoothing effect. Therefore it can basically overcome some limitations encountered in previous methods.

We also utilize the idea of the ball-filter to extract information about the local segmented structures, however, our methodology is basically different from [3] since it can provide a robust image-independent estimation of the vessel direction for expansion. In fact, as the evolution may stop to extract the entire vessel, extracted structures show some elongations. One option is to utilize this "shape-induced" information to propagate the surface. This can be useful when the image content information is not reliable due to noise or intensity ambiguities. In that case our model enforces the surface to expand anisotropically in the main surface orientation so that it can pass over small noise speckels. The outlined framework can be combined with existing level set based vessel segmentation functionals.

3. METHOD

3.1 Basic assumptions

Assume for a given open region D, the evolving surface is represented as the zero level of the level set function Φ(x)where Φ(x) < 0 for inside of the object, and Φ(x) > 0 for outside. H(x) is representing Heaviside function such that H(x) = 1 if x > 0 otherwise H(x) = 0. Also δ(x)=dH(x)/dx is the Dirac delta function. Throughout this research the phrases such as: front, surface (3D) or contour (2D) are used interchangeably with implied dimensions.

3.2 Local shape structure

The local structure of the evolving surface at point x, can be expressed by the correlation matrix of gradient vectors of surface signed distance transform (SDT) in its neighborhood. Assuming that re-initialization of level set function to SDT of the evolving surface enforces|▽Φ(x)|=1, we define the following matrix:

M(x)=∫B(x,x')H*(x')▽Φ(x')▽tΦ(x')dx' (1)

Where H*(x)=1-H(x), ▽tΦ(x') is the transpose of the gradient vector ▽2Φ(x') and B(x,x') is the neighborhood function with the general property of: B(x,x')= B(x',x) >0 Here we define: B(x,x')= 1 if x' lies inside the cubic neighborhood of size n3 pixels around x. Since the iso-level surfaces tend to be spherical as they get farther from the surface, inclusion of H*(x')is intended to reduce the directional ambiguity. At a given point estimation of local structure of the evolving surface, is defined by a modified correlation matrix of gradient vectors of its signed distance transform within a user selected scale, i.e. M(x).

Application of anisotropy constraint is achieved by evaluating the availability of a major local orientation and propagating the surface at that direction. This is accomplished by analyzing the eigenvalues and vectors of M(x).

3.3 Local shape measure

It is shown in [8] for a given point such as placed on the border, the smallest eigenvalue of the correlation matrix M(x) is zero if the local structure has a cylindrical form and the corresponding eigenvector specifies the cylinder axis. Fig.1 is a two dimensional illustration, we observe that fig. 1(A), the local structure is rather spherical compared to (B) whish has more "tapered" shape. It can be shown that in smallest eigenvalue of M(x)in (B) is smaller than (A). Since we are following a variational based optimization frame work, we chose a differentiable form of f(TrM(x)(-1)) where the function f:R+->R+ is non-increasing function and TrM(x)(-1) is the trace of inverse matrix. It is easy to see that f will take its minimum if the local structure has a zero eigenvalue.

3.4 Shape energy functional

Based on this definition, we define a geometric geodesic active contour model that minimizes a shape energy:

Es(Φ)=∫δ(Φ(x))f(TrM(x)(-1))|▽Φ(x)|dx. (2)

Note that because depends f(TrM(x)(-1)) on Φ(x), the corresponding minimization problem is different from conventional geodesic active contour model. In the next step this is minimized using gradient-decent method and with constraint|▽Φ(x)|=1, the corresponding Euler-Lagrange can be obtained as:

∂Φ/∂t=δ(Φ)[fk2+▽tf.▽Φ/|▽Φ|+▽tΦ.L▽Φ] (3)

where k2 is the surface minimal curvature the dependency on x is implied and the matrix L(x) is defined as below:

L(x)=∫δ(Φ)f'(TrM(-1)(x'))M(-2)(x'))|▽Φ|dx' (4)

where f' denotes the derivative of f. we note that for a given x, the matrix L(x) is the weighted average of positive definite matrices M(-2)(x') of its neighborhood determined by B(x,x'). It is interesting to look at the the properties of (8) in further etail. The right hand side of evolution consists of three different terms:

Smoothing: f.k2 is a mean curvature dependent smoothing term. Note that by multiplication of f, the curvature remains effective if local structure is isotropic, therefore annihilation of narrow structures with lower values of f is limited.

Advection: Second term is an advection term that attracts the object's border toward lower values of f. It reduces the orientation ambiguity, and also minimizes the leakage from spurious noisy edges.

Propagation: To analyse this term, we note that since f' < 0 we have▽tΦ.L▽Φ<0, i.e., this term propagates the surface outward. The important point is that depending on the orientation of ▽Φ(x) and eigenvectors of L(x) the expansion is anisotropic. For segmentation of tubular structures, provided that the size of neighborhood is large enough, L(x) maintains its main component in the axial orientation. As a result propagation may only appear at the endings of those structures. This can be useful when the image content information is not reliable due to noise or intensity ambiguities.

3.5 Implementation

Central-differencing scheme was used for computing of M(x) and L(x). Having the values of cost function f, advection term is calculated by simple up-winding scheme. The first order forward-Euler method is used for digitization over time and Heaviside and delta functions are evaluated using their smeared-out versions Hε(x) andδε(x) with ε= 1.5 pixel, as explained in [9]. In this research we have set f(x)=1/(x+σ), where σ is a small positive number to prevent singularity of division to zero. This setting is optional but, practically we obtained better results using this definition. The algorithm is implemented using a fast narrow band level set method [10]. We note that in (4) computation of matrix L(x) is expensive in terms of CPU cycles. As an optimization, at each iteration the evolving front is compared to its status in previous iteration and L(x) is computed at places where the front has displacements. By this means a large portion of the surface that remains intact does not require updating of structural matrix and the program executes much faster.

3.6 Sample evolutions

Fig.3 is the evolution of a cylinder. For this experiment, the neighborhood radius size is set to four pixels, i.e, the same radius of the cylinder. Initial shape is indicated on the top, and at every iteration the right hand side of (3) is color mapped onto the evolving front. Note that the expansion occurs exclusively in the axial orientation and as the evolution proceeds, the both ends turn into needle like structures with anisotropy. Basically, by extra iterations development of sharper heads is possible, however it is practically limited to the resolution of the implementation grid.

Fig.4 is a comparison between our proposed regularization method and surface minimum curvature method . For both methods the top left image is the initial surface. The left middle and bottom figures correspond to shape evolutions after 20 and 60 iterations under surface minimal curvature. Note the shrinkage starting from the ends. The right side of Fig.4 is the result of our proposed method. Note that the surface remains smooth and the anisotropic expansion appearing exclusively on axial orientations closes the ring. This shows that, using our method while preserving smoothness, improvement of vessel segmentation is possible even if the image information is noisy or not adequate.

3.7 Intracranial vascular segmentation

A combination of above defined shape dependent energy functional with the geodesic active Contours [11] Eg(Φ)=∫δ(Φ(x))g(|▽I(x)|)|▽Φ(x)|dx where g(|▽I(x)|) is a uniform decreasing function of the image gradients, is considered as the total quantity to be minimized and the Euler-Lagrange equation is derived to minimize the energy E=Eδ+αEg.The final level set equation (with implied dependency on x is given:

∂Φ/α=δ(Φ)[(1+αf/g)k2+p(▽2Φ.▽I)▽2g/g.▽Φ/|▽Φ|+α▽2ΦL▽Φ/g] (5)

here p has the same definition as in [4]. The MRA data sets are initially smoothed by a Gaussian filter and interpolated along side the axial axis. Evolution starts from a user defined threshold value for the available MRA data, and is forwarded in time using forward Euler method, and stops after a predefined number of iterations. Every few iterations the level set function is reinitialized to a signed distance transform of the zero level set.

4. RESULTS

We have applied our method to MRA data sets and consistently obtained robust extraction of vascular structures. A typical example is presented here. For qualitative assessment, results are compared with the pioneer CURVE algorithm proposed by Lorgio.et.al in [4]. Figure 5-A maximum projection intensity of a Phase Contrast MRA data set with the size of 256×256×60 voxels and spacing of 0.625×0.625×0.9 mm that is obtained from a 0.4 Tesla MR scanner. Figure 5-B is the result of our segmentation method specified equation (5). Segmentation achieved using CURVE algorithm is indicated in (C). Segmentation using our method is indicated in figure 2-D from a different angle of view. Comparison of these two figures, particularly the areas indicated by an arrow, reveals that anisotropic propagation term can be very useful to improve the continuity of the extracted vessels.

5. DISCUSSION AND CONCLUSION

Our method combines geometrical and image content information. The geometry is captured using the signed distance transform of the evolving surface and is useful to improve the continuity of vessels. Though accurate comparison with other segmentation methods requires access to the same data sets, initial results are very encouraging. Validation using images obtained from higher field MR scanners stays as our future research plan. Compared to the minimal surface curvature smoothing method that is widely used for level set based vessel segmentations, the proposed method is computationally expensive but it can provide anisotropic elongation which is particularly important for thin structure segmentation. Using our regularization method reasonable smoothing can also be obtained, nevertheless if an explicit smoothing is required it can be combined with minimal surface curvature smoothing scheme. One possible interesting direction for research is the development of a mechanism to suppress our regularization model at vessel branchings. Whether this can be achieved using geometrical features or image content features remains as the future research activities.

Fig. 1. (a) A portion of a TOFMRA 1.5 Tesla data set, the arrow indicates a damaged pattern of a thin vessel, (b) The same portion in 3 Tesla indicates an smooth vascular pattern.

Fig. 2. The contour is indicated in bold and other level sets in dot line, arrows indicate the eigenvector corresponding to smallest eigenvalue of structure matrix M(x). Not that the inner level sets around x in (A) are rather circular compared to (B). This property can be measured by Tr M(x)-1.

Fig. 3. Evolution of a cylinder; the speed of propagation on the surface is color coded by its correspondant color bar in right, Top: initial cylinder, Middle: after 100 iterations, Bottom: after 200 iterations, elongation occurs only in axial orientation,

Fig. 4. Evolution of a ring with a few embedded gaps; Left from top to bottom: initial surface, after 20 and 60 iterations under the surface minimum curvature respectively, the surface shrinks to zero by further iterations. Right from top to bottom: after 20, 30 and 50 iterations using our proposed shape regularization method, anisotropic expansion appears exclusively on the "tips" so that the ring closes itself.

Fig. 5: A): Maximum intensity projection image of a sample PCMRA data set, (B): Segmented vessels using our method, (C), (D): Segmented structures using CURVE and our method from different angle of view.

論文題目「VARIATIONAL LEVEL SET BASED METHODS FOR VASCULAR SEGMENTATION (和訳:変分レベルセット法による血管セグメンテーション)」の学位論文は,ノイズの多いMagnetic Resonance Angiography (MRA)などの画像情報に対して,セグメンテーションアルゴリズムに血管の幾何学的モデリングを組み込むことを提案し,それにより重要な血管形状情報をノイズとして処理せずに正しくセグメンテーションが行えることを示したものである.

三次元的な血管セグメンテーションは,コンピュータ支援画像診断やコンピュータ外科の分野には不可欠であるが,従来多くの血管セグメンテーション法が開発されてきたものの,それらは理想とは程遠いのが現状である.血管構造の複雑さ,微細構造のスペックトルノイズに対する脆弱性,画像のゆがみの存在などが,より高度な血管セグメンテーション手法の開発を必要とする主な理由である.これに対して,本研究のレベルセット変分法による医用画像セグメンテーション法では,コンピュータ支援による血管近傍の患部でのナビゲーションとして,2つのレベルセットによる血管セグメンテーション法を新たに提案した.主な対象は,Phase-Contrast and Time-Of-Flight MRアンギオグラフィ(MRA)のデータを用いた脳血管3次元モデルである.

第1番目のアルゴリズムは,幾何学的正則化法であり,細長い円筒構造を作るのに優れている.血管の管形状は誰でもが知っているように異方性構造を有している.すなわち,さまざまな大きさの管状パターンを有する構造的異方性は血管の基本的特性であるため,変分レベルセットによるセグメンテーション方法では,血管の異方性特性を効率的に強調するはできない.また,最小表面曲率に依存する円滑化展開は,十分に細径構造を引き伸ばすことが困難な場合が多い.さらに,強度パターンで一部破損している微細血管をセグメンテーションすると,分断された複数の領域を残したまま計算が終了してしまうことが多い.これらの問題点を解決するために,本研究では,血管の異方性先端を強調することで,抽出した血管を伸張させることが可能であることを示した.また,形状オフセットレベルセットによる先端部分の展開に,新たに形状エネルギ関数を定義した.定義した関数の制限付き最適化により,抽出構造の延長を改善できる拡張フローが,血管セグメンテーションに有用であることが確認できた.

第2番目のアルゴリズムとして,エッジ結合のための変分手法を新たに提案した.ステップ状で無いなだらかなエッジからなる領域の統計手段と方向に関する情報を,新しいエッジ概念に基づいてエネルギ最小化手法に取り込んだ.エッジ強度の効果をシミュレーションする領域駆動項は,この最小化手法から直接得られる.以上に述べた手法を複数の2次元網膜の血管画像およびMRA三次元データを用いて評価した.その結果を従来のセグメンテーション手法との比較を行い,開発した正則化手法が血管抽出に有効であることが示された.

以上により,本研究で開発した方法は,これまで困難であった血管情報とノイズとを区別して,血管画像情報に対してエッジ先端部の正しい延長が行え,血管のセグメンテーションをより正確に行うことが可能であることが示された.

よって本論文は博士(情報理工学)の学位請求論文として合格と認められる.