Abstract

I. INTRODUCTION

The recent advances in medical imaging (e.g. high resolution patient data), the availability of methods for image analysis and computer graphics, the technological resources provided by micro-devices and the design of minimal access surgical procedures have opened the road to new concepts. Virtual endoscopy represents one of these areas and points out the applicative potential of three dimensional (3D) imaging. It leads to less invasive diagnosis and therapeutic achievements, and addresses important issues in computer aided diagnosis, interventional planning, and medical education.

Volumes provided by CT or MRI allow the user to move virtually within the anatomical structures. Different approaches addressing various clinical applications, such as virtual colonoscopy [1]–[7], virtual bronchoscopy [8]–[12], virtual angioscopy [13]–[19] and others, have been proposed in the literature. Most of them assume that a pre-computed model (organ inner surface, sensor path) of the anatomical structures to be explored virtually is available. Beyond the interactive observation of hollow organs, they are now dealing with the problems of precomputing a central path and automatic path planning [20]–[23], controlling the position and orientation of the virtual camera to fully examine the interior organ surface [24] [25], quantitative analysis of the observed lesions [26], fusion of multi-modality information [27], simulation of the physical properties of the endoscope [28], matching virtual and real endoscope views [29] [30], and evaluation and clinical validation [2] [31]. Although these problems are generic, their difficulty can vary, according to the explored organ, whose size, geometry, topology, arrangement in the numerical dataset and pathology are subject to significant disparities.

We investigate a new approach which aims to explore automatically vascular structures using virtual angioscopy. Conventional fly-through approaches require a centerline which is used as the virtual angioscopy path. It may be calculated off-line or on-line. The off-line computation of the centerline in tubular objects requires the preprocessing of the volumetric dataset by skeletonization and segmentation algorithms using region growing [32], multi-scale filtering [33], or deformable models [34] [35]. The on-line extraction of the centerline was recently investigated by locally analyzing volumetric data within a plane shifted perpendicularly to the tracked centerline [36] [37] or with 3D geometrical moments [38]. These approaches, and more generally virtual path determination approaches, have still to face some difficulties. Due to the detection problems (low contrast, noisy data), the free-form shape of the structures (tortuousity), the presence of different types of lesions (thrombus, calcified plaques), the possible branching of the cavities (leading to multiple candidate paths), automatic path planning is a challenging issue in virtual angioscopy. Moreover, it is generally performed prior to the visualization process which constitutes one of the primary functions of virtual angioscopy and may provide meaningful information.

This paper focuses on scene analysis and automatic path planning. The proposed method is based on the “active sensor” notion and can be related to vision guided navigation for mobile robots. Active navigation combines the online computation of the scene view and its analysis, to automatically define the 3D virtual path and the structural description of the unknown environment explored. No pre-modeling of the anatomical structures is used: surface detection, computation of optical-like virtual images, scene analysis, and path estimation are performed throughout the navigation. Thanks to a scene analysis process, short range to long range visions are also available. The former provides the information about the immediate surroundings (walls or obstacle) and determines the displacement of the sensor and the potential for collision. The latter gives useful indications about the target area to be reached.

The general framework of this study, whose basic functions have already been reported in the literature [39], [40], is summarized in Section II. Section III presents a complete description of the scene analysis process, going from the computation of feature points to the generation of a graph showing the topology of the observed structures. In order to plan the virtual sensor path, a navigation strategy based on the use of the information derived from the scene analysis is proposed in Section IV. The effectiveness of our approach is highlighted in Section V by its application to several experimental datasets acquired by helical CT-scan modality under clinical conditions, and then discussed in Section VI.

II. VIRTUAL EXPLORATORY NAVIGATION: BASICS FOR ACTIVE NAVIGATION

An overview of the functions involved in active navigation is given in Fig. 1. The main functionalities a virtual endoscope device should provide are: image computation and display, scene analysis, path estimation, collision avoidance, and operator controllability. Different loops can then be identified to automate the navigation: they deal with the motion of the sensor to be applied, according to the detection or not of a collision, and, in case of multiple cavities, the storage of candidate paths.

Open in a separate window

Fig. 1

Overall flowchart of the active navigation.

The geometric model of the virtual camera introduced in the discrete volume is based on a classical perspective projection in which the focal length and the size of the imaging area determine the field of view. The rotations between the volume coordinate system and the local coordinate system attached to the virtual sensor are expressed in terms of Euler / Cardan angles. Yaw and pitch rotations are used to describe the virtual camera orientation. Roll rotation around the viewing axis is fixed at zero degree. Given a 3D image, the first task consists in interactively setting the virtual sensor parameters, which include its initial pose (position and orientation) within the volume and its field of view (FOV). These geometrical parameters must be completed by specifying the 3D discrete array resolution and the detection algorithm parameters. From this point on, two operational modes can be applied: (i) the semi-interactive mode, where some commands are entered interactively (stopping the navigation, deciding on the preferred path to track, …), and (ii) the fully automatic mode where a free (i.e. auto-constrained) navigation is carried out between an initial point and a target point, both defined by the user.

Virtual image computation is based on a ray casting technique integrating analysis functions. Several types of detection operators might be used [41, 42]. In this study we used a simple solution. Due to the large resolution difference between the local sensor image and the initial volume data (in the case of vascular structures notably), the trilinear interpolation of the voxel values and the surface normals is performed in order to avoid unrealistic discontinuities in the virtual image. Thus, a relatively good image quality can be obtained, even when considering small sized structures (1 mm of diameter), by seeking the transition surface between the contrast medium (or the cavity) and the vessel wall, with a coarse to fine progression rate along the ray. The surface point is found when the difference between the interpolated values at the origin and along the ray becomes higher than a user defined surface detection threshold.

Several solutions to path planning from the analysis of image sequences have been proposed in the computer vision literature [43]. They deal with problems of vision based-navigation for indoor and outdoor robots. In [44] a robot path planning approach [45] has been adapted to automatic path planning in virtual endoscopy. Concerning path generation for autonomous robot navigation in natural environments, one can generally distinguish reactive techniques based on local perception of the environment [46] [47], from planning based techniques which make use of a more global representation of the environment [48] [49]. Since the virtual image computation provides a set of cavity surface points whose coordinates are known within the camera reference system and also within the global frame of the volume image, the proposed method consists in reducing the problem to only one information: the depth map. Thus, this approach can be qualified as Map-Building-Based navigation. A number of studies dealing with depth map segmentation have already been carried out. Among them, one can mention: plane fitting and clustering [50], [51], iterative labeling of objects by means of Gaussian curvatures [52], or surface primitives (such as spheres and cylinders) fitting [53]. It must be emphasized that the environments we are dealing with are poorly featured when compared to indoor or outdoor environments. They cannot be described by means of simple primitives such as straight lines, corners or planes.

Assuming the piecewise linear approximation of the path, collision (an object must not pass through another) has to be examined. Based on the collision detection and the correction of its position, the camera-like sensor can be displaced and re-oriented. The last feedback implies the activation of a more complex decision process to choose the cavity to be explored when multiple candidates are detected. A backtracking procedure can be activated if the user-defined target is not reached (wrong previous path), as well as if the termination of a branch has been reached for instance. In the following sections we focus on the scene analysis and path planning process.

III. SCENE ANALYSIS

In our approach, the scene analysis process is based on the progressive construction of a global structure describing the hollow organ observed. At the local level (for each position of the sensor), a tree like representation is derived from the depth map processing. A node is given by a point featuring the frontier between different branches (border points), and a leaf corresponds to the maximum depth in a branch.

A. Depth map analysis

Thanks to the ray casting algorithm, the image computation process directly provides the distance between each detected surface point and the optical center of the virtual camera. This depth information can be coded by an image the same size as the virtual image (Fig. 2(a)–(c)). In this depth map, noted $\mathcal{M}$, the possible branches are associated to light gray regions. Thus, $\mathcal{M}$ can be processed as a gray level image to detect local maxima, to select the relevant ones characterizing the branches and to detect points separating the branches in the bifurcation areas (border points).

Open in a separate window

Fig 2

A typical depth map encountered in vascular structures: (a) the virtual endoscopy image, (b) the 2D gray level image of the depth and, (c) the 2.5D representation of the depth map. (d) An illustration of the detection of the depth maxima in a vascular bifurcation observed through a perspective geometry based virtual endoscope (M_i: 3D surface point corresponding to a depth maximum) and (e) the binarized depth map (mi: 2D projection of M_i on the image plane).

1) Detection of the depth maxima

A non linear filtering of the depth map is performed by imposing equality/inequality constraints which a depth maximum has to satisfy in relation to its neighbors. The maxima detection filters ${\mathcal{F}}_k$ are defined as follows:

$$\begin{array}{l}{\mathcal{F}}_k:{\mathbb{N}}^2\to \{0,1\}/{\mathcal{F}}_k(i,j)=1\\ \iff \{\begin{array}{l}\forall (i^{\prime },j^{\prime })\in {\text{V}}_1=\left[i-1,i+1\right]\times [j-1,j+1],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathcal{M}(i^{\prime },j^{\prime })\le \mathcal{M}(i,j)\\ \forall (i^{\prime },j^{\prime })\in {\text{V}}_2=([i-k,i+k]\times [j-k,j+k])\backslash {\text{V}}_1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathcal{M}(i^{\prime },j^{\prime })<\mathcal{M}(i,j)\end{array}\end{array}$$

(1)

A good compromise between over- and under-detection is reached by using a ${\mathcal{F}}_3$ filter. If no maximum is detected (or if the number of detected maxima is superior to the predefined depth selectivity parameter N_max), the filtering is performed with ${\mathcal{F}}_2$ (respectively with more selective filters ${\mathcal{F}}_4$ to ${\mathcal{F}}_7$). The set S of local depth maxima m_i (associated with 3D surface points M_i) can then be reduced by considering a distance criterion. If two depth maxima are sufficiently close together (< d_max), the smallest maximum is removed from S. For all the experiments presented in this paper, d_max was fixed at 10 pixels.

2) Labeling of the depth maxima

The objective of this step based on the segmentation of the binarized depth map is twofold: (i) to label the depth maxima according to the branch they belong to, and (ii) to characterize each branch by its deepest maximum. In order to point out the regions characterizing the branches, a binary map B_i is created for each depth maximum m_i. The binarization threshold Th is determined from two parameters: the value of the current depth maximum m_i considered in the list S, and the backward parameter, b, which defines the minimal height of the hills that can be separated in the depth map: $Th(m_i)=\mathcal{M}(m_i)-b(m_i)$. The backward parameter b is a sigmoid function of the depth which guaranties a fine separation of the branches close to the virtual sensor and limits the number of regions beyond the depth x₀ of its inflection point.

$$b(x)=\frac{b_{+\infty }-b_{-\infty }}{1+e^{-c(x-x_0)}}+b_{-\infty }\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{with}\hspace{0.17em}\{\begin{array}{l}{x}_0=40\hfill \\ c=0.1\hfill \\ b_{+\infty }=13\hfill \\ b_{-\infty }=1\hfill \end{array}$$

(2)

The labeling process consists in attributing a binary code to each depth maximum according to the region it belongs to. The coding scheme enables the distinction between regions of different natures as well as the simple handling of inclusion relationships. The labeling process works as follows:

• - Initialization of n to zero.

• - For each depth maximum m_i of S listed in decreasing order:

  • Binarization of the depth map with the threshold Th(m_i).

  • Segmentation of the binary map by means of a classical region growing algorithm using m_i as seed point.

  • Labeling of the depth maxima:

    • If the detected region does not contain already labeled maxima, labeling of the depth maxima of the region with the binary word L_n equal to 2ⁿ and incrementation of n.

    • If the detected region contains already labeled maxima, assignment of a binary word to the non labeled depth maxima, given by the logic OR operation between the binary words affected to the already labeled maxima.

For example, in Fig. 2(d)–(e) showing the binary map obtained when considering the depth maximum m₂, the region R₂₃ is segmented. The depth maximum m₁ (deeper than m₂) has been treated previously and labeled by L₀. The region R₂₃ corresponding to a single branch at the considered position of the virtual sensor – the backward parameter taking a high value –, the depth maxima m₂ and m₃ are coded by the binary word L₁. When the depth maximum m₄ is considered, a region R₄, including the regions R₁ and R₂₃, is created. Except for m₄, all the depth maxima included in this region have already been treated and labeled. In this case a compound label L₁₀ equal to L₀ + L₁ (logic OR operation) is affected to the new depth maximum m₄.

In order to characterize a branch, only the deepest maximum with a single label (code with a single bit to 1) is retained. It has been the first one to be labeled in the branch and it is called the feature maximum. For instance, the region R₂₃ and the corresponding vascular branch are characterized by the maximum m₂ with the label L₁. At the end of this process the branching structure describing the scene observed by the virtual sensor is composed of the 3D coordinates of the virtual sensor, the 3D coordinates (M) and the labels (L) of the feature maxima.

3) Detection of the border points

The feature maxima characterize the farthest regions from the virtual sensor likely to correspond to branches. Nevertheless, the depth maxima behave as viewpoint dependent descriptors. The spatial variability of the depth maxima between successive images makes them insufficient to describe and to follow the topological changes in the observed scene. Complementary elements, the border points, that intrinsically feature the frontier between different branches are introduced.

A border point is a depth local minimum separating two branches (represented by feature maxima) in the depth map. Considering the feature maxima in order of decreasing depth, the border point detection process starts with the second maximum (if it exists) and is repeated iteratively for all the following feature maxima. The depth map is binarized according to a first threshold value Th_max which is the same used for the labeling process. For instance, if we consider that the feature maximum m_i has been detected in the scene, the first threshold is defined by: Th_max(m_i) = Th(m_i) with $Th(m_i)=\mathcal{M}(m_i)-b^{\prime }(m_i)$ and b′(m_i) = b(m_i).

According to the definition of the feature maxima, the depth maximum m_i is the single one within a region R_i. By increasing iteratively the value of b′(m_i) (of a step equal to b(m_i)), and consequently decreasing the threshold Th(m_i), the region R_i grows until it merges with another region R (Th(m_i) = Th_min) which contains a feature maximum which has not yet been separated from m_i by a border point. The sought border point is located necessarily at a depth between Th_min and Th_max. The values of Th_min and Th_max (respectively associated to the binary maps F_min and F_max illustrated on Fig. 3(a)–(b)) are then refined by a bisection procedure until δ = Th_max − Th_min becomes lower than a predefined value Δ. The limit value Δ is fixed to 0.1 voxel that corresponds to the resolution step for the detection of the 3D surface points. A search window W is defined adaptively by computing the two mass centers C and C′ respectively associated with the regions R_i and R in F_max. The square window W is placed at the center of [CC′] and its size is defined by: l_w = CC′/2. The border point B is then defined as the local depth minimum of F_min – F_max in the window W (Fig. 3(c)–(d)).

Open in a separate window

Fig. 3

Seeking of the border point: (a) depth map binarized with Th_max (F_max),(b) depth map binarized with Th_min (F_min), (c) subtraction of the depth maps (F_min-F_max), (d) seeking of the depth minimum (border point) in the window W.

B. Structural description

At the local level, i.e. for a given position of the virtual sensor, a data structure can be derived from the feature points (feature maxima, border points) in order to define a structural description of the unknown scene. During the exploration of the anatomical structures, a global description is progressively constructed and updated by the local descriptions.

1) Local description

For each position of the virtual sensor, a structural description is constructed as a tree like structure. A node is given by a border point (featuring the frontier between different branches) and a leaf corresponds to the maximum depth in a branch. Each element of the binary tree is assigned with the 3D feature point (depth maximum or border point), the position and the orientation of the virtual sensor for which the feature point has been obtained, and the address of a possible child element.

The binary tree-like structures are relatively simple to handle and are well appropriate to the description of the anatomical structures that the virtual exploratory navigation is concerned with. A notable exception is encountered when the anatomical structure appears as a single tube-like shape where only one depth maximum is detected and no border point is defined. In this case the local description is reduced to a single element. It is said to be degenerated. The leaf is nonetheless retained within the description.

2) Global description

Different viewpoints along the virtual sensor trajectory lead to different local descriptions. Since border points are implicitly associated with geometrical features of the 3D scene, they are much less sensitive to viewpoint than the depth maxima. Border points are used to perform the matching of the successive local descriptions. Let $B_i^k$ and $B_{i+1}^l$ be the border points given respectively by the local descriptions d_i and d_i+1 at two positions C_i and C_i+1. $B_i^k$ matches $B_{i+1}^l$ if the following condition is verified: $\Vert B_i^k{B}_{i+1}^l\Vert \le d_b$, where d_b is a heuristic parameter fixed at 10 voxels. The associated depth maxima (leaves of the tree) are also matched according to a distance criterion. In order to keep the longest branch in the global representation, a depth maximum is updated only when the new maximum is deeper than the one it is candidate to replace. A border point which is not associated with another one after the matching process highlights a topological variation in the explored structure: either a branch appearance, if the border point belongs to d_i+1, or a branch disappearance if it belongs to d_i.

In addition to the geometrical matching between the local descriptions, some tree editing operations are used to update the global description. Each node is assigned with a weighting factor defined by the number of its occurrences, so that topological events such as a branch appearance or disappearance can be represented with a certain degree of confidence. When a branch disappears, the representation of the corresponding node in the global description is not removed but its occurrence factor is no longer incremented, so that the information is kept (Fig. 4(a)). It might be used to move the virtual sensor backward to an older position in order to explore another branch in the volume. When a new branch appears in the field of view of the virtual sensor, the local description creates a new node in the tree structure (which cannot be matched with a previous node). This node must be inserted in the global description. The operation (Fig. 4(b)) is performed in such a way that the consistency of the description and the heritage relationships are preserved: the new node is inserted between the nodes (in the previous global description) corresponding to its father and its child (in the local description).

Open in a separate window

Fig. 4

Updating of the global description for different topological events (D and Du are respectively the current and the updated global description, d′ and d″ are local descriptions): (a) in the case of the disappearance of a branch, (b) in the case of the appearance of a branch, (c) in the case of a degenerated description.

The tree insertion issue may be a more complex task when the element to be inserted is also a tree structure. Assuming that during its motion, the camera moves from a complex scene (the global description lets appear several nodes) to a simpler one like a single branch and later from this degenerated description to a branching area, the description evolves as described in Fig. 4(c). In this example M′5 corresponds to M5 in the global description, but the branch matching is impossible since d′ is degenerated. d″ does not match with D neither. Therefore, we keep a trace of the sensor displacement regarding D in order to determine an activated link, which in this case is B4 – M5. Then, thanks to the link stored previously, when d″ occurs, the local description is linked to the global one.

IV. NAVIGATION STRATEGY AND PATH PLANNING

The only environment map available to the virtual endoscope is given by the structural description it constructs by itself (scene analysis process). The virtual sensor determines its trajectory by combining information of two kinds (Fig. 5).

Open in a separate window

Fig. 5

Displacement process: the bold contours represent the 3D surface points used for the computation of the displacement direction (C: optical center or virtual sensor position, CM: viewing direction, CG: displacement direction, T: local target point).

One, the viewing direction is determined according to the branch within which the sensor has to navigate. Considering that a 3D target point E has to be reached at the end of the virtual fly-through, a branch characterized by a depth maximum M_s is selected in the local description such as:

The depth maximum M_s is used to set up the viewing direction for the computation of the next virtual image.

Two, the direction of the next displacement within the free space is derived from the mean position G of a reduced set of 3D surface points observed by the virtual sensor. The magnitude of the displacement step MDS is chosen by the user at the beginning of the exploration. Thus, the displacement vector of the optical center C is determined by $\overrightarrow{S}=\overrightarrow{CG}·(\mathit{MDS}/\Vert \overrightarrow{CG}\Vert )$. In order to take into account the nature of the environment (bifurcation, single branch, …), the 3D surface points are selected according to their depth and the proximity or not of a border point. We define l_v as the maximum depth at the periphery of the map. If the depth of the closest border is more than twice the value of l_v, we assume that the virtual sensor is moving within a single branch. G is computed from the observable surface points located at a distance of less than 2l_v from C. It determines a displacement vector noted S_c. If the depth of the closest border point is less than twice the value of l_v, the virtual sensor is moving within a complex area (bifurcation for instance). G is computed from the observable surface points of the branch to follow located at a distance of more than the depth of the border point (region containing M_s). The displacement vector is then noted S_f.

Note that other features could have been chosen to constrain the sensor behavior. For instance, the viewing direction and the sensor displacement direction could have been fused. Our choice is mainly motivated by the fact that when a branch to follow is chosen – more particularly in a branching area –, the path planning has to be defined by the features of this particular branch, whereas the features of the other branch do not have to be considered.

In order to force the virtual sensor to stay inside the cavity, the next location, which is estimated from a current position, is compared to the depth map. To illustrate the collision detection process, we can consider a local 3D target point T and its projection in the depth map p_T. The displacement of the virtual sensor is validated only if $\Vert \overrightarrow{CT}\Vert <\mathcal{M}(p_T)$. When a collision is detected, there is no displacement and only the viewing direction is updated before the computation of the new image. If the collision is still detected, the focal length of the virtual sensor is iteratively reduced by half, in order to enlarge the field of view, and a new image is computed and analyzed at each iteration. If the virtual sensor is jammed, the navigation has to be resumed at a node from which one of the branches has not yet been explored. The restart position is given by the last location where a node with a high weighting factor was observable.

V. RESULTS

Since several parameters intervene in the proposed navigation scheme, we first study their influence both on depth map analysis for different scenes observed in patient datasets, and on virtual endoscope behavior from a phantom model dataset. We then present results of active navigation in the patient CT images.

B. Virtual endoscope behavior

A phantom model representing a vessel with a stenosis and two bifurcations (Fig. 6) was used to illustrate the scene analysis process and to assess the influence of the field of view (FOV), the magnitude of the displacement step (MDS) and the choice of the starting point on the virtual path generated by active navigation. Contrast medium has been injected in the vessel lumen of a maximum of 4 mm in diameter (20 voxels) before the CT acquisition (3 mm thickness slices, reconstructed every 1.5 mm).

Open in a separate window

Fig. 6

Two MIP images of the phantom model (voxel size: 0.2 × 0.2 × 0.2 mm³ after interpolation). The reference trajectories in the main branch (beginning: near points P₁, P₂; end: near point P₄) and in the secondary branch (beginning: near points P₁, P₂; end: near point P₃) are superimposed on the MIP images, a, b, c, d, e, f, g, and h denote eight selected poses along the virtual path in the main branch.

Fig. 7(a) shows the images as well as the local description computed by the virtual sensor inside the main branch. A starting point (P₁) and a target point (P₄) have been defined by the user before launching the virtual exploratory navigation. The virtual sensor passed through the lesion area and automatically reached the target point. At pose a, the vessel bifurcation is not yet detected. At poses b, c, and d, two branches characterized by two depth maxima are detected. At poses e and f, the first bifurcation has disappeared and the second bifurcation is not yet highlighted. At poses g and h, the second bifurcation is detected and represented by the local description. The global description is obtained from the matching of the successive local descriptions (Fig. 7(b)).

Open in a separate window

Fig. 7

(a) Virtual images and local descriptions (gray: 3D path, child node: feature maxima, parent node: border points) constructed by the virtual sensor at 8 selected poses, (b) Global description obtained in the phantom model.

Several virtual paths have been automatically generated for different values of the navigation parameters. For a fixed target point P₃ located in the secondary branch, the variable parameters were: the starting point (with two different positions P₁ and P₂ located in the main branch before the narrowing of the structure), the field of view FOV (120°, 140°, 160°) and the magnitude of the displacement step MDS (1, 2, 4, 8 voxels). The virtual paths have been quantitatively compared to the reference trajectory which is given by the centerline of the phantom CAD model. For each of the 24 combinations of parameters, the mean value and the standard deviation of the euclidean distances between each point of the virtual path and the closest point on the reference trajectory were used to quantify the results obtained by active navigation (Table III). The number of computed points per virtual path varied from 16 (MDS = 8) to 135 (MDS = 1). For both starting points P1 and P2, the distances between the virtual paths and the reference centerlines decrease when FOV increases. Since the surface points observable by the virtual endoscope are determined by the FOV, the best results are obtained when MDS leads to a local target position close to the center of the surface points used for the computation of the displacement direction (low FOV and high MDS, or high FOV and low MDS).

Table III

Euclidean distance (in voxels) between the virtual path and the reference centerline. Mean (Standard Deviation) for different settings of MDS, FOV and starting point.

MDS \ FOV	120°	140°	160°
1 voxel (P1)	2.31 (1.59)	1.9 (1.28)	1.3 (1.04)
2 voxels (P1)	2.58 (1.93)	1.79 (1.32)	1.22 (0.94)
4 voxels (P1)	2.39 (1.88)	1.79 (1.27)	1.25 (0.99)
8 voxels (P1)	2.15 (1.81)	1.7 (1.23)	2.11 (2.02)
1 voxel (P2)	2.74 (1.92)	2.04 (1.14)	1.9 (1.76)
2 voxels (P2)	2.67 (1.76)	2.03 (1.23)	1.32 (0.93)
4 voxels (P2)	2.63 (1.85)	2.09 (1.13)	1.58 (1.28)
8 voxels (P2)	2.32 (1.46)	2.89 (2.22)	2.05 (1.14)

Open in a separate window

As a final evaluation of virtual endoscope behavior, the evolution of the viewing direction of the virtual endoscope along the virtual path is reported in Fig. 8. for three different settings of FOV (120°, 140°, 160°) with MDS = 4 and starting point P1 One can observe that the viewing direction is nearly the same along the virtual path for FOV = 120° and 140°. The peaks between positions 60–80 correspond to the commutations between displacement vectors S_c and S_f which modify the path estimation and therefore the viewing direction. Since, with a FOV of 160°, the virtual endoscope remains closer to the centerline in the single branch area (displacement vector S_c), an important variation of the viewing direction is obtained in the bifurcation area (displacement vector S_f), whereas it becomes more stable in the single branch after the bifurcation.

Open in a separate window

Fig. 8

Evolution of viewing direction (pitch and yaw in degrees) along the virtual paths for different settings of FOV.

C. Virtual fly-through of patient vasculature

Our method was applied to the 20 CT volumes mentioned in Section V.1., in order to automatically navigate between a starting point and different target points defined by the user. The parameters MDS, FOV and SDT were also fixed at the beginning of the exploration. Fig. 9 shows three virtual paths obtained in the carotid arteries with stenosis as well as the detected surface points and border points. It is a typical example of the results obtained in one of the seven datasets of the carotids. Although there was an important variation of the depth maxima positions, the border points remained quite stable. Because of the complex topology in the stenosis area (calcified plaque), secondary branches have been detected. Nevertheless, the main branches were correctly matched. The bifurcation areas were successfully passed through and the virtual paths stayed inside the lumen in the different cases without any collision. Using a PC with a 2 GHz CPU and 1 GB memory, the average computation time was around 0.63s/pose (over a total number of 179 poses) for the virtual path going from P₁ to E₂. It varied locally from 0.6s/pose in a single branch to 1s/pose in a bifurcation area. The computation time at each main step of the algorithm was distributed in the following way: 63–73% for image computation, 4–7% for depth maxima detection, 0–13% for border points detection, and 16–26% for structural description, virtual sensor displacement and other processing tasks.

Open in a separate window

Fig. 9

(a) CT slice at a particular position and the representation (blue) of the viewing direction of the virtual angioscope (CT acquisition parameters: 3 mm thickness slices, reconstructed every 1 mm, voxel size: 0. 24 × 0. 24 × 0.24 mm³ after interpolation). (b) Rendering showing the detected 3D surface points, the border points (red) and three virtual paths (green) computed automatically by virtual active navigation (SDT = 100 HU, FOV = 140°, MDS = 1, Nmax = 5,: starting point:P₁, target points: E₁, E₂, E₃). (c) Selected virtual angioscopy images with the representation of the border points (red).

Fig. 10 illustrates results obtained in one of the six datasets of the aorto-iliac structure. Two different virtual paths going from the supra-renal aorta to the iliac arteries have been computed. Due to the presence of thrombus and many calcified plaques, the number of secondary border points was relatively important. The target point E₂ was successfully reached. Nevertheless, a secondary depth maximum at the level of a calcified plaque caused the virtual angioscope to be jammed and the target point E₁ not to be reached. The backtracking from the last observable node with a high weighting factor (aorto-iliac bifurcation) did not allowed the virtual angioscope to reach the target point. Moreover, we noticed that if we start from the end of the iliac artery towards the previously mentioned starting point, by following only the deepest depth maximum, thus assuming that no bifurcation has to be observed, the virtual path is correctly completed.

Open in a separate window

Fig. 10

Rendering showing the detected 3D surface points, the 2D projection of the border points (red) and two virtual paths computed automatically by virtual active navigation (SDT = 90 HU, FOV = 140°, MDS = 1, N_max = 5, starting point:P₁, target points: E₁, E₂) in the aorta with aneurism and iliac arteries (CT acquisition parameters: 5 mm thickness slices, reconstructed every 4 mm, voxel size: 0.49 × 0.49 × 0.49 mm³ after interpolation).

In the seven datasets of popliteal arteries, the main branches were correctly passed through. Nevertheless, it was impossible to reach some lateral branches with the surface detection threshold fixed at the beginning of the navigation, i.e. at the starting point located in the main branch. Due to low contrast distribution in the bifurcation area, the relative threshold (referring to the difference between the values at the origin and along the ray) fixed in the main branch was either too low or too high to correctly detect the surface points in the small lateral branches (about 8 voxels and 1.7 mm diameter). It was confirmed by modifying the detection algorithm and considering two absolute thresholds (low threshold: detection of soft tissues, high threshold: detection of hard tissues) for which navigation from the main branch to the lateral branches were successfully completed.

VI. DISCUSSION

In the proposed virtual navigation approach, the difficulties of the tasks to be performed depend on the reliability of the detection (performance of the ray-casting algorithm in presence of low contrasted transitions, noise and artifacts), the scene complexity (shape of the surfaces bounding cavities, deviation from pipe-like cavity, hidden parts, multiple tracks, etc …), and the performance of the sensor control (obstacle avoidance, strategy to correct its pose, capability to face complex paths, and loss of control).

The quality of the scene analysis process has been assessed by considering the detection of the border points. Though the maximum probability of detecting a border point when there is an anatomical branching point (sensitivity) is 95 %, the maximum probability that a detected border point is actually an anatomical branching point is 84 %. That means that some secondary nodes do not correspond to anatomical features, notably in the lesion areas. Nonetheless they seem to be more scattered than the nodes which actually correspond to anatomical features, so that most of them might be filtered or discarded according to their weighting factor.

Several definitions of the optimality can be considered when generating virtual paths. They deal with: (i) minimal paths which are mainly concerned with robot motion planning to find the shortest path to the target [44], (ii) centerline paths which always stay inside the tubular structures and which are connected and far from the edges [20], [22], and (iii) reliable paths which ensure that the inner organ surface is fully examined [25]. Although the primary objective of this work was less to define an optimal path than to investigate a path planning approach based on active vision concept, a reference centerline has been considered in the phantom model to quantitatively appreciate the behavior of the virtual active angioscope. The closest paths to the reference centerline were generally obtained for a high field of view.

Although surface detection can be listed as a critical issue (possible failure of the detection process, volume artifacts), the study reported in this paper is not mainly aimed at dealing with this problem. The qualitative results obtained from different patient datasets (popliteal and carotid arteries) show that the interactive settlement of the surface detection threshold in the initialization step provides a simple solution to obtain successful virtual navigation inside vascular structures with lesions. The relatively high computation time compared with other tracking approaches has to be balanced by the fact that the proposed method is automatic. A minimal interaction is required from the user at the initialization step to define the surface detection threshold and to compute a virtual image at the desired starting and target points. The processing algorithms for depth map analysis and other tasks have still to be optimized, but the main time-consuming task remains image computation which may reach 73% of the total computation time. This process may be time-optimized by further considering coherency properties, such as pixel-space, voxel-space and frame-to-frame coherence, thus significantly reducing the number of cast rays and the number of detection operations per ray. Although the structural description may handle different kinds of topological events (branch appearance, disappearance), the existence of cycle (such as bypass) was not taken into account in this study. Moreover some target points were not reached due to the presence of secondary border points. The approach could be improved by reconsidering both the decision process and the surface detection process. The former was mainly concerned in this study with simple decision functions (go on with the navigation towards the target point in case of no collision, enlarge the field of view in case of a collision, select a node with a high weighting factor). It should be able to take the optimal decision to select branching candidates or to correct the virtual sensor trajectory when the target is not directly reached. Concerning the latter point, recent studies [54] allowed us to extend the geometry of the virtual sensor in order to use robust surface detection operators - consequently more time consuming - along rays cast in radial directions. An adaptive threshold could be derived from a reduced set of rays and introduced in the whole image computation process.

VII. CONCLUSION

A new method to guide the path planning in virtual angioscopy has been presented. It is based on the notion of active vision, in which visual perception completely drives the motion of the virtual sensor. The local scene interpretation during the navigation inside vascular structures allows the virtual angioscope to automatically determine its path, without the need of a pre-computed model of the patient data.

Through experiments performed on a phantom model and patient data, the conclusion that can be drawn from this study is the feasibility of the proposed method for guiding autonomous virtual sensor in complex vascular structures with lesions. Further investigations about the selection of branching candidates and the trajectory correction, could help to improve the robustness of the method. Moreover, the computation of quantitative parameters characterizing the anatomical structures from the virtual fly-through should be of main interest in clinical applications, notably for the planning of endovascular interventions. To this end, the difficulty facing us is related to the local orientation of the virtual sensor in the vessel lumen that can affect the results of the quantification. Moreover, further work should also deal with the application of virtual active navigation to other tubular anatomical structures and to spatio-temporal data.

Acknowledgments

References