Skeleton Marching-based Parallel Vascular Geometry Reconstruction Using Implicit Functions

Fast high-precision patient-specific vascular tissue and geometric structure reconstruction is an essential task for vascular tissue engineering and computer-aided minimally invasive vascular disease diagnosis and surgery. In this paper, we present an effective vascular geometry reconstruction technique by representing a highly complicated geometric structure of a vascular system as an implicit function. By implicit geometric modelling, we are able to reduce the complexity and level of difficulty of this geometric reconstruction task and turn it into a parallel process of reconstructing a set of simple short tubular-like vascular sections, thanks to the easy-blending nature of implicit geometries on combining implicitly modelled geometric forms. The basic idea behind our technique is to consider this extremely difficult task as a process of team exploration of an unknown environment like a cave. Based on this idea, we developed a parallel vascular modelling technique, called Skeleton Marching, for fast vascular geometric reconstruction. With the proposed technique, we first extract the vascular skeleton system from a given volumetric medical image. A set of sub-regions of a volumetric image containing a vascular segment is then identified by marching along the extracted skeleton tree. A localised segmentation method is then applied to each of these sub-image blocks to extract a point cloud from the surface of the short simple blood vessel segment contained in the image block. These small point clouds are then fitted with a set of implicit surfaces in a parallel manner. A high-precision geometric vascular tree is then reconstructed by blending together these simple tubular-shaped implicit surfaces using the shape-preserving blending operations. Experimental results show the time required for reconstructing a vascular system can be greatly reduced by the proposed parallel technique.


Introduction
Cardiovascular diseases (CVDs) are the number one cause of death of human beings globally [1,2] . It is anticipated, by 2030, that the death of CVDs will increase to over 23.3 million [3] . Although a large number of CVDs are preventable, they are continually rising due to inadequate measures and the high risk of traditional surgeries [4] . In order to overcome this high risk, computeraided minimally invasive vascular techniques were proposed to assist the diagnosis and surgeries of CVDs [5] . One essential task of computer-aided minimally invasive vascular techniques is the accurate reconstruction of blood vessel geometry out of vascular images. However, blood vessels are of high complexity in their geometric structures, which makes the reconstruction task extremely challenging.
The geometric reconstruction of blood vessels is to re-build the geometric structures of the vascular wall out of medical images. Mathematically, existing techniques for this task can be categorised largely into two groups: explicit modelling and implicit modelling [6] , depending on the form of mathematical representation of the extracted vascular geometric shapes. Explicit modelling typically represents vascular surfaces in the form of triangle meshes, which is a preferred approach when the main objective of generating the vascular shape is for visualisation. However, explicit techniques are very inconvenient for modelling geometric objects with complex branching structures [7] . In contrast, implicit modelling represents a geometric form as a field function, which allows implicitly represented geometric objects to be easily blended to generate more complicated geometries. This easy-blending feature of implicit geometric modelling makes it particularly suitable for vascular modelling [8,9] . Fast patient-specific high-precision vascular geometric modelling based on the scanned volumetric medical image datasets is a difficult task. Since a blood vessel is in general of a tree structure with each of its branches having a thin tubular shape, most vascular modelling techniques are vascular skeleton-based. The skeleton of a blood vessel is an important geometric clue on vascular modelling, which is often represented as the centreline of a vascular tree. However, there are some vascular modelling techniques that are not developed based on vascular skeletons. Whether skeleton extraction is an essential requirement for vascular modelling has been investigated and discussed by some researchers [10][11][12] . In general, existing skeleton-free modelling techniques, such as marching cube [13] and multi-level partition of unity implicits (MPUI) [14] , are difficult to use to create a high-precision geometric reconstruction. In addition, the surface quality created with this type of technique is often poor, and an additional smoothing operation is usually required to further improve the smoothness of the extracted surface, which inevitably makes the accuracy of the extracted surface unguaranteed. Compared with skeleton-free techniques, skeleton-based methods are more intuitive and are straightforward to implement. More importantly, by marching along the skeleton, it is no longer necessary to check those parts of a given medical image that do not contain any vascular vessels. This helps to improve the efficiency and effectiveness of the reconstruction task.
For the skeleton-based vascular reconstruction techniques, sweep surface is one of the representative methods. It intuitively generates a series of cross-sections along the skeleton by sweeping a surface and progressively approximates these cross-sections to be a generalised cylinder. This is a two-step approximation modelling method. A cross-section is first approximated as a generalised circle, and a series of generalised circles are then approximated to be a generalised cylinder.
Inspired by the sweep surface method, this paper proposes a skeleton marching-based vascular geometric reconstruction to produce a high-performance and high-accuracy geometric vascular tree. Instead of using sweep surface, a skeleton marching strategy is followed to build the blood vessel along the skeleton, which not only avoids the two-step approximation of the sweep surface method, but also improves the modelling performance and accuracy. In this method, the skeleton of a blood vessel branch is firstly represented as a curve. This curve is then divided into many shorter segments with small overlaps at the ends of the skeleton segments. For each of these short skeleton segments, a sub-region of the given medical image containing the skeleton segment is identified and a localised image segmentation technique is applied on this sub-image region to extract a small point cloud. A series of connected small point clouds will be extracted in this way along the skeleton of the blood vessel. Since these point clouds are independent of each other, it is possible to approximate them into a set of implicit surfaces in a parallel manner. An unbranched blood vessel will be geometrically reconstructed by blending these implicit surfaces together. In the same way, many unbranched blood vessels can be reconstructed as implicitly represented objects such that a simple implicit blending operation can combine them into an entire geometric vascular tree. This article is an extension based on the conference paper [15] . The proposed technique has been applied to real medical images and very satisfactory results have been achieved. Experimental results show that the proposed method has a much better performance in vascular reconstructions.

Implicit geometric reconstruction
In the field of geometric modelling, an implicit geometric model is a mathematical representation of an attribute of a 3D object in the form of a field function. The aim of implicit modelling is to find a scalar function to represent a required geometric shape either as a surface level set or a solid level set of the function. For simple geometric objects such as spheres and cylinders, the corresponding implicit functions can be written out directly. However, finding an implicit function to represent a highly complex real-world shape like the human brain vascular system is obviously difficult.
The activities of geometric modelling can be generally grouped into two types: designing a required shape based on a set of specified constraints and reconstructing a shape from a real-world object. Compared with geometric design, geometric reconstruction is in general a much more difficult task, especially when the object to be reconstructed is of a highly complex geometric structure.
Several implicit geometric reconstruction techniques have been proposed in recent years, such as level set method [16] , moving least square (MLS) method [17,18] , variational implicit surface method [19] , adaptively sampled distance field method [20] , multi-level partition of unity implicits (MPUI) method [14] and radial basis function (RBF) method [21,22] . However, none of these methods can be used directly for reconstructing a complex vascular system.
Fortunately, implicitly modelled geometric objects are easy to blend to generate more complex geometric shapes. Let and be two implicit objects represented by field functions and respectively. The implicit objects corresponding to the intersection , union and subtraction of and can be defined using the maximum function , a special and the simplest intersection blending operator, in the following way: (1) g(x, y) g∩, g∪ g \ For a general blending operation , the blending operators and corresponding to the intersection, union and subtraction of two implicitly represented geo-Q. Qi et al. / Skeleton Marching-based Parallel Vascular Geometry Reconstruction Using Implicit Functions f0 f1 metric shapes with function and are defined in the following way: Though there are various simple ways to define the blending operation , more advanced blending operations are required to blend implicit shapes together smoothly. In this paper, the smooth function is chosen to define . It is a smooth blending operation behaving like the classical function, but can achieve any required continuity at the blending region and the range of the blending region is controllable. Let be the smoothness degree and be the blending rang control parameter, the smooth function takes the following form [9] : and |x| , for n = 0 1 2(n + 1) for n = 1, 2, 3, · · · .
δ > 0 Both sharp and smooth joints can be produced from shape-preserving blending operators. With an additional blending range control parameter , they are an ideal choice for the composition of implicit functions.
f0(x, y) = y − sin(0.5x) f1(x, y) = y− cos(6x) Fig. 1 presents the intersection blending between implicit functions and using the shape-preserving blending operation. Both the two blending results are continuous, but a smaller value gives narrower smooth blending range and better preserves the original shapes of and . The bulge is also smaller at the places where the two implicit shapes intersect. 0 Fig. 2 presents an implicit geometry designed by the level sets of the implicit function where and are sharp blending subtraction and intersection operators, is the -smooth shapepreserving union blending operator with blending range parameter value .
By using blending operations, especially shape-preserving blending operations, the implicit modelling technique can be used to reconstruct very complicated geometric objects by subdividing a complex geometric shape into a set of simple geometric objects, which can be reconstructed independently of each other. Our vascular modelling technique is developed based on this idea.

Vascular modelling
Skeleton-based implicit vascular modelling usually regards the task of reconstructing the entire tree-structured vascular system as a set of tasks of reconstructing a single branch and reconstructing a single branch can be performed by marching the skeleton of the branch. It is intuitive to sweep a skeleton-orthogonal surface along the skeleton to generate a generalised cylinder [23] . The sweep surface models a blood vessel branch by generating a set of cross-sections along the skeleton of the vessel branch, where each cross-section is approximated by a generalised circle. In order to guarantee the orthogonality between the cross-section and the blood vessel, a local coordinate transformation is introduced as a Frenet frame.
Cross-sections are represented as 2D implicit generalised circles and further extruded along the skeleton tangent as implicit generalised cylinders. These generalised cylinders are treated as an accurate representation of corresponding blood vessels and then blended together to be a vascular tree. Both Hong [11] and Kretschmer [12] adopt this scenario to reconstruct blood vessels. Hong′s method uses 2D piecewise algebraic spline [24] and shape-preserving implicit blending operations [9] . Kretschmer uses Catmull-Rom spline [25] and gradient-based implicit blending [26] . Although the reconstruction accuracy at the cross-sections are preserved, the sweep surface method has several limitations. Firstly, there are two extra operations of the sweep surface method: indispensable coordinate transformation for orthogonality between the blood vessel and the sweep surface, and unavoidable points sorting up on the cross-sections for the approximation of generalised circles. These operations reduced the performance of the modelling. Secondly, the sparsity of cross-sections influences the accuracy of the reconstruction. Data points between two consecutive cross-sections are ignored. The higher the density of the cross-sections, the better the modelling quality will be. Lastly, vascular modelling based on sweep surface uses two-step approximation. Cross-sections are firstly approximated as a series of closed splines, then these splines are approximated as a generalised cylinder. Modelling performance will decrease with this two-step operation.

Skeleton marching
In this paper, an implicit modelling technique, called Skeleton Marching, is proposed analogous to team cave exploration for the geometric reconstruction of blood vessels. With the proposed method, the skeleton of a blood vessel is firstly represented as a spline. This spline is divided into several overlapped shorter segments, which assist a localised image segmentation technique to extract point clouds corresponding these segments. These small point clouds are then fitted and blended into a whole blood vessel represented as an implicit surface. All the locally reconstructed implicit blood vessels are then blended to represent the geometry of the whole vascular system.
Reconstructing a vascular tree with the skeleton marching technique is similar to exploring a cave with multiple branching structures by following the cave skeleton, and thus can be directly implemented in a parallel manner, as each branch and each vascular segment can be reconstructed independently of other branches or vascular segments. The skeleton marching technique consists of subdivision of the skeleton, localised implicit reconstruction, localised implicit object and parallel computing. Fig. 3 illustrates the process of skeleton marching technique.

Subdivision of skeleton
The basic idea of skeleton marching is localisation. It not only simplifies the modelling complexity but also accelerates the modelling speed. Since a vascular image contains highly complicated blood vessels together with a large amount of non-vascular information, the vascular modelling in global scope cannot, in general, produce satisfactory results for both large and small blood vessels. A localised region marching on the skeleton of the blood vessel can be used to subdivide a vascular modelling into small modellings for qualified reconstruction results. A localised region on the skeleton contains a short segment of a blood vessel. This segment is going to be reconstructed as a localised implicit object, which is a small generalised cylinder whose position and size are determined by the length and curvature of the skeleton enclosed in it.
In order to subdivide a blood vessel into shorter vascular segments with a reasonable shape, the skeleton of  this blood vessel is regarded as a parametric curve such that its curvature and length can be obtained conveniently for the subdivision.
Suppose a 3D curve is represented as a parametric function , then the skeleton of a blood vessel can be represented as a parametric curve [27] γ The curvature of the curve is given by γ(t0) The length of the curve started from point is given by With the curvature and length, a long and curvy skeleton can be subdivided into shorter and less curvy segments such that the corresponded vascular segments are in simple shapes. The neighboured segments share an overlapped part which is used for maintaining the fitting accuracy during the blending process after these segments are reconstructed as localised implicit objects.
This subdivision is delineated in Algorithm 1, which generates a set of knot pairs . One knot pair marks the two ends of a skeleton segment. The shape of this segment is controlled by and such that the curve of the segment is neither too long nor too curvy. Consequently the localised region surrounding the skeleton segment is simple enough. The data points in this small region can be easily extracted and reconstructed to be a localised implicit object.
are two empirical values. They are dependent on the choice of the surface fitting discussed in the next section. This paper uses the RBF surface fitting with ellipsoid constraint, which will give redundancies if the objective shape is long and curvy. However, thanks to the shape-preserving blending operations in (6), the reconstruction accuracy will not be influenced. The performance is not affected either since parallel computing is insensitive on the fluctuation of the dataset when it is big enough.
Algorithm 1. Division of an unbranched skeleton L Parameter: : threshold of length In order to smoothly blend the neighbouring localised implicit objects together, Algorithm 1 is designed in such a manner that the start endpoint of knot pair is in the middle of knot pair . This overlapped region is long enough to avoid unwanted bulges when blending. Fig. 4 gives an example of this subdivision. A medical volume containing a blood vessel is rendered in Fig. 4(a). This blood vessel is very curvy and is expected to be subdivided into shorter and less curvy segments. In Fig. 4(b), a skeleton is shown as a curve inside the blood vessel. Nine markup knots are located with Algorithm 1, each markup knot is paired with such that two consecutive segments share a common part. These nine markup points produce seven segments shown in Fig. 5.

Localised implicit reconstruction
Subdivision of the vascular skeleton simplifies the reconstruction of a blood vessel by breaking down a big modelling task into many small sub-tasks. Due to the simplicity of each local vascular shape, the localised reconstruction of the region surrounded a skeleton segment is simple, accurate and efficient.
A skeleton segment is a short curve inside a localised blood vessel. The shape of the blood vessel segment is small and simple such that the point data of the vascular segment can be collected with many different image segmentation algorithms. In this paper, the implicit deformable model segmentation algorithm [28] is chosen to extract the data points out of a short blood vessel. This segmentation method evolves a level set inside a blood ves- sel from its start point to the endpoint and collects the surface points of the blood vessel on the go. For a small blood vessel segment, this algorithm gives efficient and effective segmentation results. Fig. 6 presents localised segmentation results of seven vascular segments from the curvy blood vessel in Fig. 4. Each segmentation result is a scatter point cloud. From Figs. 6(a) to 6(g), meshes are rendered with the data points for better observation. All the data points shown with different colours are given in Fig. 6(h).
The localised segmentation collects the surface point clouds of the small blood vessels corresponding to the skeleton segments. Each point cloud is a small dataset with a capsule-like shape and is expected to be reconstructed as a surface. Because of the simplicity of the underlying geometry represented by the small dataset, an uncomplicated surface fitting algorithm is preferred.
As shown in Fig. 6, the localised segmentation result is always a closed point cloud without holes and self-intersections. Although some of them have curvy parts, the general shape is always a capsule or deformed ellipsoid. Under such conditions, the direct RBF surface fitting with ellipsoid constraints is selected to reconstruct these localised blood vessel segments, which is a high-accuracy surface fitting method designed for small datasets [22] .
The direct RBF surface fitting with ellipsoid constraints is based on an assumption that the RBF-based surface fitting problem can be regarded as a blending of two implicit surfaces: a surface fitted by radial basis functions and a surface modelled by a low degree polynomial. Since the implicit surface of a radial basis function is always a blobby model, a shape-closed polynomial constraint, such as an ellipsoid, is an ideal choice to deform the blob to be the desired shape.
Given a set of points from a surface, the direct RBF fitting with ellipsoid constraint is to fit points set as an implicit function : where is a radial basis function and is a polynomial always representing an ellipsoid, is a general point. Assume , then can be expressed in the following matrix form: where The solution to (12) can be obtained by solving an eigensytem subject to the condition that always represents an ellipsoid. A full explanation of (11) and (12) can be found in [22].
The direct RBF surface fitting with ellipsoid constraint is especially suitable for the reconstruction of small datasets like the localised segmentation results shown in Fig. 6. As a one-step fitting algorithm based on the solution of the eigensystem in (12), this direct fitting does not require additional information such as surface normals and extra fitting layers. The fitting is fast when the dataset is small. In addition, because of the ellipsoid constraint, this method always gives bounded fitting results which is the expected shape of a vascular segment. As the fitting results are expressed as implicit functions, it is easier to combine the reconstructed vascular segments together using implicit blending operations. Furthermore, since each closed vascular segment has a simple shape, the corresponding implicit shape has a simple form, whose computational cost can be dramatically reduced. This will be discussed in the next section. Last but not least, this method is a high-accuracy surface fitting method when the shape of the datasets are simple. The fitting errors are very small and negligible [22] . Fig. 7 presents localised modelling results of the curvy blood vessel corresponding to the segmentation results of Fig. 6. Both the reconstructed surfaces and the data points are rendered. The whole curvy blood vessel is given in Fig. 7(h). It is the blended result of the implicit objects shown in other seven subfigures. Although the blending smoothness can be any degree by using the shape-preserving blending operators in (6), a smoothness is sufficient for the purpose of geometric modelling. In this paper the smoothness degree is set to with a smoothness span controller , where is from an empirical validation. Table 1 gives the distance errors of the localised modelling results of the seven subfigures from Figs. 7(a) to 7(g). The distance is measured from each data point to the reconstructed surface. The distance error is evaluated by the standard deviation of all the distances of a localised reconstruction result. It is shown that the errors are so small that the data points can be regarded as on the surface.

Oi
The localised reconstruction process of vascular segments generates a group of bounded small implicit objects , which are defined by a set of implicit functions in an image space . Although each individual implicit object takes only a  is crossing the whole image space . Its evaluation is expected to be limited inside as the implicit object is only a locally defined implicit object.
Definition 1. (Localised implicit object) An implicit object defined by an implicit function is said to be a localised implicit object if there exists a subset of , such that it is a bounded object defined in , and that , where .
and is expected to be blended with another localised implicit object defined in . However, these two implicit objects cannot be blended together directly because the domains of the corresponding implicit functions are different. In order to achieve this blending, their domains have to be extended to the origin domain . The truncated implicit function is used to do this extension.f

Definition 2. (Truncated implicit function) A field function
is said to be a truncated implicit function (TIF) of an implicit function if there exists a sub-domain , such that the implicit solid object is a localised implicit object.
For a normal implicit function , let be the domain on which it is truncated, then the truncated implicit function on can be expressed in the following form:f c where is a sufficient large constant.
With the truncated implicit function, a group of localised implicit objects can be quickly evaluated inside their localised defining domains and blended together in a global domain. It dramatically improves the performance of the geometric reconstruction process when the number of the localised implicit objects is large and their sizes are small. Fig. 8 demonstrates how the truncated implicit function works. In Fig. 8(a), isocontours of implicit function are shown and the level set is marked. Since an implicit object is an interior implicit solid object , the evaluation outside can be neglected. In Fig. 8(b), the implicit function has been truncated to be inside a sub-domain enclosing such that the unwanted evaluation outside this sub-domain is eliminated.
Ui Oi Ui However, since an implicit function in general cannot explicitly represent an implicit surface, there is no direct way to generate points or surfaces of an implicit function [7] . In order to locate the sub-domain of the localised implicit object , the location of the localised segmentation will be used as a substitution. Consider the small size of the localised segmentation, a minimum cubic bounding box is chosen to be the sub-domain . Fig. 9 shows the differences between a globalised implicit object and a localised implicit object with their minimum bounding box. The branched blood vessel has been divided into 41 simple vascular subsections and their data points on the vascular walls, which are rendered as blue points, have been extracted by the implicit deformable model segmentation method [28] . The red capsule-like objects in the two subfigures are reconstructed small vascular sections. Although the two objects are geometrically identical, they are generated from different domains. The fitting of Fig. 9(a) uses the global domain and costs   Table 2 compares the fitting speeds of the blood vessel segments in Fig. 7. The globalised fitting costs much more time than the localised fitting since the localised fitting uses truncated implicit function such that only the sub-domain of the localised implicit object is evaluated and the fitting is much faster. Fig. 10 compares the fitting speed of each blood vessel segment. The speed remains at a low level of the globalised fitting using the normal implicit function. On average nine points are fitted per second. In contrast, the fitting speed of using the localised implicit function is about ten times faster, and the speed grows when the number of points increases.
By using the truncated implicit function, the fitting speed of the localised implicit object is dramatically improved at the algorithmic level. Since the fitting of one localised implicit object is independent of the others, the fitting efficiency can be further improved by using parallel computing techniques.

Parallel computing
The subdivision of the vascular skeleton turns a big modelling task into many small subtasks, each of them is regarded as a localised reconstruction and generates a localised implicit object. A notable feature of the localised implicit objects is that they can be built independently of the others, thus, localised implicit reconstructions work well with parallel computing techniques.
Parallel computing is a technique that improves computational efficiency by using multiple processing elements simultaneously to solve a problem. [29] . There are several different types of parallel computing, but only data parallelism is used in this research.
Data parallelism distributes the data across multiple processors for parallel computing. It has received much attention with the fast development of general-purpose computing on graphics processing units (GPGPU). The massive parallel processors and memories on a GPGPU are able to process complicated tasks in a parallel manner and significantly improve the computational efficiency [30] .
Data parallelism has been used on the RBF fitting with ellipsoid constraint to improve the performance of the vascular reconstruction. For a point cloud from the localised segmentation, the implicit function in (12) is expressed and distributed over the GPGPU as a kernel of parallel computing. The solution is transferred back to the host programme on CPU for further processing. Table 3 compares the fitting time of the blood vessel segments based on Table 2. In this table, Columns 3 and 4 show the globalised fitting time on CPU and GPU. The fitting speed of data points on GPU is significantly improved. Columns 5 and 6 give the localised fitting on CPU and GPU, but as can be observed, the fitting performance on GPU is not improved, instead, it is even slower. This is because the localised implicit objects of the blood vessel segments are too small to be accelerated. Data transformation between CPU and GPU slows down the acceleration. The fitting speed comparison is given in Fig. 11. The experiments are remotely running on an high performance computer (HPC) with an Intel(R) Xeon(R)   Besides the data parallelism of the modelling on one implicit object, multiple localised modellings can be processed concurrently. This is also regarded as task parallelism. Localised implicits and localised modelling are designed for this purpose. Because the calculation of one localised implicit object is independent of the others, different localised implicits can be configured into different streams for concurrent computing. Futhermore, higher parallelism can be achieved when multiple GPU devices are available.

Experimental results and discussions
A skeleton marching algorithm is proposed in this paper using the localised implicit objects for the reconstruction of blood vessels. The segmentation, surface fitting and parallel computing of a single localised implicit object have been discussed to illustrate the proposed method. However, it is worth stressing that this method is designed for the reconstruction of vascular trees rather than simple blood vessel segments. This section discusses the advantages and limitations of the proposed method with further experiment results.

Unbranched blood vessel
A short curvy blood vessel has been reconstructed us- ing the proposed method when discussing skeleton marching. This section presents the reconstruction of another unbranched blood vessel in Fig. 12. The volume rendering of this blood vessel is shown in Fig. 12(a). Its skeleton is represented as a cubic Hermite spline and ordered knots have been positioned on it using Algorithm 1. Spline segments are divided in such a way that knot is paired with . Each segment overlaps its neighbours, this is to make sure that the adjacent locally fitted implicits are blended smoothly. Fig. 12(b) gives the localised segmentation of each spline segment. Random colours distinguish different point clouds where small and simple shapes can be observed. Fig. 12(c) shows the localised reconstruction result of one point cloud. Fig. 12(d) presents the blending result of all modelling results with point clouds rendering. The shape-preserving blending operation is used [9] . Fig. 12(e) shows the reconstructed blood vessel without point clouds. Fig. 13 shows a branched blood vessel reconstructed with the proposed method. Figs. 13(a) to 13(c) give the point clouds, point clouds on the blood vessel wall and the blood vessel wall as implicit surface, respectively.

Branched blood vessel
This branched blood vessel consists of 59 segments from 7 smaller single blood vessels. Each segment corresponds to a small point cloud extracted with localised segmentation. Small and simple-shaped implicit objects are fitted from these point clouds and then blended together to create the final result. The segmentation and surface fitting of one point cloud has no relation to the others and therefore is parallel-friendly. The underlying implicit function of the surface fitting is represented as the truncated implicit function in (13) which makes the fitting very fast.

Blood vessel tree
The reconstruction of a blood vessel tree is similar to the reconstruction of branched blood vessels. Fig. 14 presents a reconstructed blood vessel tree with 294 segments from 58 unbranched blood vessels. The original image is rendered as a volume in Fig. 14(a). In Fig. 14(b), the point clouds of the 294 segments are rendered with random colours. In Fig. 14(c), both the points and the reconstructed surface are rendered. Fig. 14(d) gives the reconstructed vascular tree without the points.
The flowchart in Fig. 15 illustrates this reconstruction. Localised segmentation is applied on short blood vessel segments guided by skeleton segments, then localised modelling fits the segmentation results with small localised implicits. Implicit blending operation combines them together to get the final vascular tree. Compared with Fig. 3, this chart emphasises the parallel computing used in the Skeleton Marching technique.
Parallel computing plays an important role in the proposed reconstruction method. For a single localised modelling on one point cloud, data parallelism aims at the acceleration of the RBF fitting. Multiple localised model-  Table 4 compares the performances of reconstructing blood vessel(s) on CPU and GPU with globalised fitting and localised fitting. Datasets A, B and C in the first column correspond to the reconstructed blood vessel(s) in Figs. 12(c), 12(e) and 13(c), respectively, each of them contains 316 data points, 4 196 data points and 17 503 data points. Fig. 16 gives a comparison of the modelling speed of the three examples in Table 4. The first bar shows the globalised fitting speed on the CPU remaining at a stable level. The second bar reflects the rapid growth of globalised fitting speed on the GPU when the number of data points increases. In contrast, the modelling speed differences between the localised fitting over the CPU and GPU are not that significant. The performance of the localised fitting on GPU is even worse than the CPU counterpart when the dataset is small, but its advantage is emerging with the increasing of the size of the dataset. This is because the localised fitting is less computationally intensive as the dataset is small. In general, parallel computing gives limited improvement for non-computational intensive tasks.

Discussions
The skeleton marching-based vascular reconstruction is a high-accuracy and high-performance vascular reconstruction method. Compared with the sweep surface method, it has several distinctive advantages in blood vessel reconstruction: 1) No coordinate transformation is required. Each vascular segment is reconstructed to be a localised implicit object inside the image space. This reconstruction does not need the Frenet coordinate.
2) The dataset of each localised reconstruction does not need to be sorted out. The reconstruction of a vascular segment is a surface fitting on an unorganised dataset. The coordinates of the data points are the only prerequis-ite of this fitting.
3) All data points are considered by the proposed method. The reconstruction of a blood vessel segment is based on a 3D point cloud rather than a fraction of surface points. No data point is neglected.
4) The fitting of a localised implicit object is a onestep fitting. The data points of a vascular segment are directly reconstructed as an implicit object. No post-fitting operation is required.
5) The proposed method is a high-performance vascular reconstruction technique. On the one hand, the localised implicit object limits the evaluation of each vascular segment reconstruction inside a small region such that the fitting efficiency is dramatically improved because of the significantly recused fitting size. On the other hand, the proposed method is parallel computing-friendly. All localised implicit objects can be reconstructed in a parallel manner to save more computational cost.
6) This vascular reconstruction technique is of high accuracy. First, the direct RBF fitting with ellipsoid constraint is a high-accuracy surface fitting method. Various shapes can be accurately represented. In addition, the blended shape of the localised implicit objects can preserve the original shapes as much as possible by using the shape-preserving implicit blending operations. Both the degree of smoothness and the blending range can be flexibly controlled to adjust the blending results.
The main drawback of the proposed method is the overuse of data points. In order to blend the neighbouring vascular segments together, each segment has overlapped regions with its neighbours and the data points in these regions will be processed twice. This overuse of data points guarantees the smoothness and accuracy of blending results but slightly decreases the modelling efficiency. However, the size of the overlapped region is controlled by the length and curvature of the vascular segments such that this negative issue can be suppressed as much as possible. This overuse can be observed in Fig. 6(h). The size of the overlapped regions is to be optimised in the future work.

Conclusions
The skeleton-based vascular modelling is the mainstream approach in accurate reconstruction of blood vessels out of medical images. Among the existing skeletonbased vascular modelling techniques, the sweep surface method allows high-accuracy modelling results at the cross-sections but it is far less accurate at reconstructing the parts between two consecutive sweep surfaces. Furthermore, it is a two-step approximation technique, whose performance can drop significantly with the increase of the number of cross-sections considered. In order to avoid this problem and improve the performance of geometric vascular modelling, this paper proposes the skeleton marching-based vascular reconstruction to give a  With the proposed technique, a vascular tree is firstly divided along the vascular skeleton into short simple blood vessel sections and a small simple point cloud is extracted from each of these vessel sections using a segmentation technique. Then a set of simple implicit shapes are reconstructed in a parallel manner using these point clouds. A high-accuracy vascular tree is finally generated by blending these localised implicit objects together using shape-preserving blending operations. Compared with the CPU-intensive implementations, much less time is required using the proposed parallel implicit fitting technique. The proposed method can be used to reconstruct scene from scanned "3D information" similar to the "P-RM method" proposed in [31].