Abstract
In this article, we propose a novel image segmentation method called the whole mesh deformation (WMD) model, which aims at addressing the problems of modern medical imaging. Such problems have raised from the combination of several factors: (1) significant growth of medical image volumes sizes due to increasing capabilities of medical acquisition devices; (2) the will to increase the complexity of image processing algorithms in order to explore new functionality; (3) change in processor development and turn towards multi processing units instead of growing bus speeds and the number of operations per second of a single processing unit. Our solution is based on the concept of deformable models and is characterized by a very effective and precise segmentation capability. The proposed WMD model uses a volumetric mesh instead of a contour or a surface to represent the segmented shapes of interest, which allows exploiting more information in the image and obtaining results in shorter times, independently of image contents. The model also offers a good ability for topology changes and allows effective parallelization of workflow, which makes it a very good choice for large datasets. We present a precise model description, followed by experiments on artificial images and real medical data.
Keywords:
Segmentation; Parallelization; Deformable models; Medical imaging1 Introduction
The image segmentation procedure is one of the most important steps involved in the process of image analysis. This refers to the task of partitioning a given image into multiple regions and is typically used to locate and mark objects and boundaries in input scenes. After segmentation, the image represents a set of data far more suitable for further algorithmic processing and decision making. Image segmentation algorithms are a very broad field and they have received significant amount of research interest.
A very interesting family of image segmentation algorithms, that has been gaining a lot of focus for many years, is called deformable models. They are based on the concept of placing a geometrical object in the scene of interest and deforming it until it assumes the shape of objects of interest. This process is usually guided by several forces, which originate in mathematical functions, features of the input images, and other constraints of the deformation process, like object curvature or continuity. A range of very desired features of deformable models include their high capability for customization and specialization for different tasks and also extensibility with various approaches for prior knowledge incorporation. This set of characteristics makes deformable models a very efficient approach, which is capable of delivering results in competitive times and with very good quality of segmentation, robust to noisy, and incomplete data.
Numerous examples of deformable models usage can be found, starting from the earlier years of image processing until the recent research efforts. The first work about that subject has been presented by Kass et al. [1]. These authors have described a method called Snakes, which proposed placing a single contour in the scene of interest and then subjecting it to deformations until it assumes the shape of the objects present in that scene. The deformations have been constrained by the external and internal energies, which described the features of the scene and of the contour itself, respectively. This idea has caused a lot of interest in the field of image segmentation in general and in relation to medicine. Numerous authors started to propose their improvements and changes to the original formulation, including the geometric active contours models [2] (applied to magnetic resonance images of brain) or active contours for segmentation objects without strictly defined edges [3] (tested only on artificial images and pictures) among the most noteworthy ones.
Solutions derived from the original method formulated by Kass et al. [1] usually used an Euler differential equation to determine the solution. A different approach has been presented by Amini et al. [4] in their early publication, where the authors proposed a solution based on dynamic programming. The method allowed the introduction of new type of constraints, called hard constraints, describing rules that could not be violated. It also guaranteed the numerical stability of the solution, thus addressing a serious disadvantage of the Kass method, where the iterations that formed the intermediate steps of execution showed a large level of instability and for the final solution they had to be considered meaningless. The drawback of the Amini solution was a big overhead in terms of memory requirements and execution times. The idea of algorithmic approach was further examined by Williams and Shah [5], who have proposed a greedy algorithm approach and also introduced some advances regarding the energy function estimation. Their solution delivered a significant improvement in terms of execution time and memory needs, as by the definition it considers only local information in each iteration. Therefore, this approach would not guarantee that the found resulting solution would be the globally lowest one. However, those authors argued that the tests of their method have proven its ability to deliver results of a very close quality to those of the dynamic programming version. Tao et al. [6] have extended a very popular formulation of external force, called the Gradient Vector Flow (GVF) [7]. Their formulation, called the fluid vector flow, has improved some features of the GVF that have been not optimal, namely insufficient capture range and poor convergence for concavities. Xie [8] has experimented with gradient vector interaction in order to deal with the initialization dependency problem. An interesting approach has been presented by De Santis and Iacoviello [9]. These authors have experimented with a discrete formulation of the energy function instead of defining the algorithms in the continuum. Their experiments have shown that using the said approach it was possible to obtain more accurate segmentation with significantly lower computational costs.
The ability to change topology of the shape has been a very significant component of deformable models and numerous works have been presented with the aim of classifying and describing different aspects of topology changing [10,11]. McInemey and Terzopoulos [12,13] have considered the incapability of the parametric deformable models for topological transformations without additional mechanisms. They have introduced a model called the Tsnake, which was able to dynamically adapt its topology to that of the target object, flow around objects embedded within the target object, and/or automatically merge with other models interactively introduced by the user. A very important solution has been proposed by Casselles [2], who presented a model based on a curve evolution approach instead of an energy minimization one. It allowed automatic changes in the topology when implemented using the levelsetsbased numerical algorithm [14] and also naturally prevented selfintersecting, which is a costly procedure in parametric deformable models. Their solution has served as the base for numerous following works [1517]. However, the ability to change the topology is not always desired. Especially in the field of medical imaging, we are often dealing with a situation when the topology of the object of interest is known from anatomical knowledge and can be defined upon algorithm execution. In order to provide such functionality a topologypreserving formulation was proposed [18]. It has however achieved the goal by imposing a hard constraint on the number of connected components, which had to be known before the segmentation and could not be modified during it. This was too restrictive for some applications, so a more subtle method was proposed in [19]. The authors have formulated their solution in a way that would allow the components of the segmented shape to merge, split or vanish without changing the genus of the initial deformable model.
Due to their attractive features one of the fields in which deformable models perform very well, and thus are a very popular choice, is the area of medical image analysis. Digital image processing has successfully been applied to this area for more than three decades. The numerous benefits that it offers include, in particular, improvement in the interpretation of examined data, full or nearly full automation of tasks normally performed by a physician, better precision and accuracy of obtained results, and also possibility of exploring new imaging modalities, leading to new anatomical or functional insights.
Despite the large amount of work carried out in this area, deformable models still suffer from a number of drawbacks. Those that have been gaining the most focus are:
• sensitivity to the initial position and shape of the model—without a proper initialization the model might get trapped in local minima and thus not reach the objects of interest or not detect some of their features correctly;
• sensitivity to noise in the input images and to flawed input data;
• problematic topology changing—whenever the scene of interest includes more than one object or when the objects present in the scene contain discontinuities, deformable models need to change the topology of their shape. This is not straightforward in the parametric formulation of the method and requires specific algorithms;
• the need for user supervision over the process.
Some of the above drawbacks have successfully been addressed in an interesting work that has been presented by Barreira and Penedo [20] and further described in [21]. They have formulated a model called The topological active volumes (TAV) which has been introduced as a general model for automatic segmentation of 3D scenes. The TAV model has differed from the first methods based on deformable models, in which the segmentation was performed using solely boundary information and the result included only a contour or a surface describing the shape of the objects of interest. In the case of TAV, the deformed shape was represented with a volumetric mesh, with some nodes being responsible for describing the boundaries of objects and others for modeling their interior structure. The mesh was initialized over the entire image and then converged towards the objects present in the scene. Thanks to the mesh structure, the TAV model was able to describe segmented scenes with more detail and more resemblance to the realworld objects. It also showed a very good potential for the topology changing capabilities and solved the issue of initialization, thanks to the presence of the nodes in the entire image. However, the initialization of the mesh was done in a way, which could be compared to initializing any other Active Contour method over the entire input image. Although it assures covering the entire area of interest, it also introduces disadvantages, like starting the segmentation process from the the most distant location possible. This increases the overall segmentation time and also raises the chances of performing an error during segmentation, because more irrelevant objects and noises are encountered during the process.
Modern medical imaging is focused around highresolution image data, which is capable of delivering increasingly more information to medical practitioners, which in turn leads to improvement in detection rates and increased precision of computeraided surgeries and planning. However, processing such large sets of data presents new challenges, which is partly a result of the current advances in the field of microprocessors. No longer are we witnessing significant growth of possibilities of single processing unit, but rather a trend towards multiunit processing. Numerous attempts have been taken to parallelize the workflow of medical image processing using computer clusters [2225].
In this article, we present our innovative model for 3D image segmentation, called whole mesh deformation (WMD) model. It presents a set of very desired solutions that successfully address the abovementioned disadvantages: it eliminates completely any reliance on the initialization of the process, it allows efficient topology changes and shows low dependence on user interaction. Comparing to the TAV solution it also offers a number of significant advantages, namely much better computational efficiency, high suitability for effective parallelization, and much better solving of the initialization problem. Because of a flexible and parameterized implementation of the energy function the WMD model can also easily be extended with further possibilities, like the ability to incorporate prior knowledge using, e.g., statistical models. The proposed method is designed in a way to be highly suitable for modern applications of image processing, which is a treatment of large datasets composed of 3D, highresolution images and taking advantage of multiprocessor execution environments. The preliminary version of this model has been presented in [26]. In this article, we describe our model with more detail, propose a new mechanism for topology changes and a method for workload parallelization. We also present more experiments using real images as the input.
The remainder of this article is organized as follows: in Section 2, we describe our model for image segmentation, the WMD model, and we describe how can it deal with the most significant problems that can be encountered during 3D image segmentation tasks. In Section 3 we describe the numerical parameters that are used for configuration of the model and we explain the high robustness of the method in case of selection of non optimal values. In Sections 2.5 and 4, we describe the process of shape optimization along with the condition of segmentation finalization and recognition of irrelevant parts of the images. In Section 5, we describe a very important part of the segmentation techniques based on deformable models, namely the topology changes scheme. In Section 6, we describe the parallelization of the model, namely how efficiently it shares the workload when executed on a multiprocessing unit environment. In Section 7, we present our experiments with the WMD model divided into following parts: comparison of the segmentation time with the TAV method using a set of artificial images, comparison of the segmentation time using different levels of precision, performance with real medical images, and the performance gain when executed in a parallel environment. Finally, in Section 8 we present our conclusions.
2 Description of the WMD approach
2.1 General description of the ideas
As mentioned before, the WMD model has been designed using the ideas of deformable models and the TAV method as its foundation. Similar to TAV, the segmentation process is carried out by constructing a threedimensional mesh structure that would cover the entire input image data. Next, the mesh is deformed by moving its nodes in order to detect the objects present in the scene of interest. However, the mesh in the proposed WMD model is defined in a different way and the process of handling its deformations is built upon different assumptions. The comparison of the two approaches is presented in Figure 1.
Figure 1. Visualization of concepts of the TAV (top row) and the WMD models (bottom row).
In the original idea of TAV, the mesh was initialized over the whole volume with all the nodes divided into two groups: internal and external ones. The external nodes, which in fact performed the vital part of the segmentation process, were distributed only over the superficial layer of the volume. The internal nodes would be distributed in the whole area between the external nodes and their function was reduced to model the internal structure of the objects by maintaining an even distribution between the external nodes. We consider this approach to be not optimal for the following reasons:
1. The number of nodes that segment the boundary of the object is usually significantly lower than the number of nodes which model their interior structure. This might not be the case in the exemplary image given on Figure 1, but in realworld scenarios the density of the mesh is significantly higher, resulting in different proportions between the external and internal nodes. Furthermore, because of the mesh structure, every augmentation of the number of external nodes would result in an even bigger increase of the internal nodes number. As a result, the computational power is not efficiently used, because distributing the internal nodes inside of objects of interest is a trivial task and the most precision is required at the boundary detection part.
2. The segmentation process is completed not sooner than when the external nodes place themselves on the borders of objects of interest. This is a very inefficient approach, because the distance that the external nodes need to travel during the segmentation process is relatively large. The mesh is always instanced in a way to cover the entire image (or volume) and the external nodes are instanced on the edges of the mesh. This puts them in the most distant location possible from the objects of interest, which in turn results in longer segmentation times.
3. The shapes present in the scene of interest are approached from image borders with a single layer of external nodes. This makes the segmentation time highly dependent on the contents of the image, because the distance which the nodes need to travel can be varying in different images or even in different parts of a single image. This makes the segmentation time very difficult to predict or to control. That is a much undesired characteristic if we want to perform a parallel implementation of the algorithm or to segment a threedimensional volume composed of a number of slice images. Without the possibility of predicting the segmentation time for each part of the image it will be very hard or impossible to guarantee an even distribution of workload and overall segmentation time in these cases will always be longer than necessary.
In order to take advantage of the great potential of the TAV method but also to cope with the abovedescribed drawbacks, we have proposed the WMD Model, which uses a different organization of the mesh in order to process the scene of interest. An exemplary segmentation using the ideas of our model is presented on the bottom row of Figure 1. The major difference is the lack of division between internal and external nodes of the mesh. All nodes are treated equally and their behavior depends only on the image features in their nearest neighborhood. However, no functionality of the segmentation method is lost. The model maintains the ability to describe segmented scenes with high detail and great resemblance to the realworld objects, offers a very good potential for the topology changing capabilities and independence of initialization. The overall idea is to simulate the behavior of both types of nodes from the TAV model, but without dividing them strictly in two groups. During the segmentation process each node should discover by itself which behavior should it carry out and proceed in the desired way for the rest of the process. As a consequence of the abovedescribed difference, the segmentation process is performed in a different manner: the nodes, which have been instanced near the boundaries of the object, begin to progress towards them, while the remaining nodes change their positions only enough to maintain a stable, regular structure of the mesh. In the end of the process, the nodes that are irrelevant are discarded.
Enforcing the correct behavior on all nodes of the mesh without dividing them strictly in two groups is obtained with an innovative formulation of the energy function of the mesh. The task of defining a single energy function that could incorporate all necessary behaviors and perform a correct segmentation is naturally more complex than in the case of two energy functions, but it opens a possibility to address all of the abovedescribed disadvantages and offer a significant performance improvement.
2.2 Formulation of the energy function
The energy function of the model, which guides the segmentation process, is defined in a way to assume its minimal values when the nodes of the mesh position themselves on the shape of interest. In order to simulate the deforming behavior of the mesh, the TAV model needs two energy functions, one for each type of the mesh nodes. Naturally, the WMD model needs only one energy function, because the mesh is composed of only one type of nodes.
The energy function consists of two groups of forces: internal and external ones, which are responsible, respectively, for preserving the structure of the objects and for applying the features of input images. The division of the energy function into two forces has a purely organizational purpose and it takes its origin from the first publications about deformable models [1].
In order to calculate energy for a given model state, the parameter domain [0,1]×[0,1]×[0,1] is discretized as a regular mesh defined by the internode spacing Gx, Gy, Gz, and each of the contributing forces described in Sections 2.3 and 2.4 is calculated using the appropriate algorithm.
2.3 Internal energy formulation
The internal energy is composed of two forces, namely of continuity and curvature [1], and it is defined as follows:
The parameter k stands for the total number of nodes in the mesh. Parameter represents the average distance between the neighboring nodes (average edge length). The parameters m_{n} and μ_{n} are, respectively, the average distance and angle between the edges incident at a node n and are calculated by dividing the sum of lengths (or angles) of the links in the neighborhood of the node by their number. Only those links that are not longer than the flexibility parameter allows are considered for this step. The flexibility parameter of the links is a feature used for the task of topological changes of the WMD model. This characteristic will be described in more detail in section 5 of this article. The symbols α and β next to the sum expressions represent their weights and serve to balance their impact on the whole equation. As we can see from Equation (1), the continuity force attracts the nodes of the mesh to maintain equal distances between each other, whereas the curvature term attracts the nodes to keep a 90degree angle. This serves to preserve the initial structure of the mesh based on regular, cubic elements.
2.4 External energy formulation
The external energy is composed of three forces and is defined as follows:
The symbol I(v) represents the intensity values taken directly from the input images. The G(v) symbol stands for the GVF [7] values, which are calculated using the GVF algorithm over the input images. Similarly, the Edg(v) symbol represents the values from the edge detector, which are defined using the Canny Edge Detector algorithm [27] over the input images. Similar to Equation (1), the symbols γ,δ, and ε represent the weights of particular components.
2.5 Proposed optimization of the shape of the mesh
The segmentation process is performed by solving an optimization task using the greedy algorithm approach [5]. The energy function is used as the optimized function and the current shape of the mesh is taken as its input. During the procedure, the following is performed for each node N of the mesh: if the coordinates (x_{n},y_{n},z_{n}) describe the position of the node N at the time t, then for time t+1 the coordinates of N would be described with (x_{n+l},y_{n+k},z_{n}) where k,l∈{−1,0,1} and correspond to the lowest possible value of:
which is the energy of the node N calculated with Equations (1) and (2) at these given coordinates. As it can be seen, the node is moved in its nearest neighborhood in X and Y planes and the position with the lowest energy is chosen as the new position of the node. Those steps are repeated for each node of the mesh until the following rule is satisfied:
where μ is a small value near zero. Equation (4) verifies the number of nodes that have changed their position in the last algorithm iteration. Whenever this number is decreased to zero (or very near to zero) the mesh is assumed to be in its stable position and the segmentation is finished.
2.6 Complexity of the WMD model
The proposed WMD model follows a pattern of greedy algorithm optimization for all the nodes of the mesh. The steps performed in that pattern include: (a) calculation of the energy (Equations 1 and 2) in all the neighbor locations of a given node; (b) definition of the location with the lowest energy as a new location for the given node; (c) repetition of steps (a) and (b) for all nodes of the mesh, until the overall energy continues to decrease. From those steps we can see that the execution time of the segmentation task depends linearly on the number of the nodes in the mesh and the necessary number of algorithm repetitions (calculation of new position for each node of the mesh). Furthermore, the number of necessary iterations depends strictly on the distance which the nodes need to travel during the optimization step. In the WMD model, we have focused on decreasing this value as much as possible. As a result, the general complexity of WMD and TAV models is similar but the pessimistic execution time of the WMD model is severely lower than the pessimistic execution time of the TAVbased algorithm.
3 Proposed numerical parametrization of the energy function
The symbols α,β,γ,δ, and ε next to the sum expressions in Equations (1) and (2) represent their weights and serve to balance their impact on the whole equation. In order to achieve the desired behavior of the mesh as shown in Figure 1 we need to choose those values correctly. Our experiments have allowed to establish a relation, which stands correct for every tested execution scenario and therefore highly reduces the risk of choice of incorrect values for the parameters. The said relation is as follows:
Symbols starting from the left represent the weights of: curvature, continuity, GVF, image intensity, and the edge detector terms. The justification for the above assumption is the following:
• The continuity and curvature terms need to be weighted with values just high enough to maintain a high regularity of the mesh. This will not decrease the quality of segmentation as the majority of nodes will have to travel only small distances during the entire process.
• The nodes that are initialized in the proximity of objects of interest need to be attracted to their borders and they should be the only ones attracted—nodes that are initialized farther should maintain their positions and regularity of the mesh. That is why we need to give the GVF force a weight slightly higher than the continuity and curvature ones, but keep the range of GVF on a low level.
• A high weight of the image intensity values serves to impose the correct behavior on the nodes of the mesh that have been initialized inside the objects of interest. The weight needs to be higher than the one of GVF, as the nodes inside the objects should not be attracted by the GVF force. Instead, they should be distributed near their initial locations and maintain a stable structure, thanks to the continuity and curvature forces. This scenario assumes that the objects present in the images are bright and the background is dark. However, this can easily be modified whenever required by the scenario of application.
• The weight of the edge detector energy parameter needs to be set to the highest value, as it will serve to fix permanently the position of the nodes that reach the edges of objects of interest.
During the process of establishing and confirming the correctness of (5), we have also noted that the proposed WMD model proves to be robust and not highly sensitive to changes in the values of the numerical parameters, as long as they follow (5). Each segmentation scenario has a certain choice of optimal values of the parameters, but in most cases within the margin of 2.5% of the optimal values it was not possible to notice any change in the delivered result. The margin of 3.5% introduced small changes but still delivers correct and precise results.
4 Recognition of the unwanted parts of the mesh
As may be seen in Figure 1 some groups of nodes should be discarded at the end of the segmentation process. These would be the parts of the mesh that have been instanced outside of any objects present in the scene of interest. As shown in Figure 1, the behavior that we would want them to maintain would be to reside at their initial positions (or near to them) and gradually disconnect from the parts of the mesh that are attracted to boundaries of object. That would be obtained using the topology changing capabilities of the WMD model. As soon as the optimization step is finished the unwanted parts of the mesh would be recognized and discarded from the result (as shown in Figure 1, bottom right image). To do this, the mesh is divided into submeshes that are isolated from each other, using the information about links that are marked as broken by the topology changing mechanism and the flexibility term. For each submesh, the following values are calculated:
where σ represents the standard deviation of the image intensity values in the submesh and N_{E} stands for the number of nodes finalized on the edges of objects of interest, n is the total number of nodes in the mesh, I(V_{i}) is the intensity value of the input image at the point occupied by the node V_{i} and is the average intensity value for the whole input image. The parts that we want to discard from the final result would show very low values given by those two equations, because most of the nodes in those submeshes would be residing over a background area. Therefore, a simple clustering procedure into two groups is enough to separate the relevant from the irrelevant parts.
5 Proposed dynamic topology changes model
The proposed WMD model requires an efficient and precise scheme for the topology changes. This is due to the fact that some parts of the mesh can become irrelevant during the segmentation process and as a consequence they need to be disconnected from the relevant part and discarded. To cope with those needs we have developed a topology changes mechanism as a component of the WMD model. It allows to change the topology of the mesh during the segmentation process and to create discontinuities in the mesh structure while the optimization step progresses. Whenever a need for mesh reconfiguration would be encountered during the shape optimization, the model should react in a desired way and start creating discontinuities in the mesh structure. On the other hand, the said topology changes feature should be constrained at all times by the mesh energy function, which would guarantee that it will not cause it to progress out of its stable state. This twofold dependency between two processes would assure the correct behavior of the mesh and a stable progression towards its optimal shape. This is performed in contrast to the TAV method, where the task of topology changing is a separate step of the whole segmentation process [20,21].
In principle, the change of topology is carried out by removing connections between nodes of the mesh. During the shape optimization step each node of the mesh can change its position due to the energy minimization process and as a result the lengths of links connecting it with its neighbors are also altered. What follows is the initialization of the topology changing mechanism in order to verify if any discontinuities in the mesh should be created. This is done with the following steps: the lengths of the links between the given node and its neighborhood are checked with the following term:
where L_{n} represents the length of the current link, G_{x} and G_{y} are the lengths of the links in initial distribution of the mesh, and flex is the flexibility parameter, which is defined automatically upon segmentation execution in the range of
taking into consideration the size of discontinuities in the input images. This range has been established experimentally using 30 artificial images as the input. Whenever a certain link fails to satisfy (8) it is marked as broken. We can safely assume that this process is carried out only when desired because the mesh is defined to have a rigid and stable structure and as a consequence the majority of the links of the mesh extend their lengths only by small values during the whole process. When the situation of a link breaking occurs, we know that this is a result of attracting the given node to a border of an object present in the scene of interest and thus the connection breaking should be called.
The broken link will no longer be considered when calculating the continuity term of the energy function in the next iteration. As a result, the node will show a behavior as if this link did not exist at all and it will progress away from its current location more freely. This will usually cause extending and possibly breaking the links of its neighbors, as shown in Figure 2.
Figure 2. Example of a chain reaction during discontinuity detection in the structure of an object.
This sort of chain reaction is much desired as it will trigger the movement of nodes and breaking of connections in a small neighborhood. This will in turn lead to a successful detection of the entire discontinuity in the mesh structure. Such process will be stopped at the right locations, namely where the nodes of the mesh encounter the edges of objects of interest. This is guaranteed by defining the energy of the edge detector to a high value, which will always stop the progression of nodes.
6 Parallel implementation
The WMD model has been designed since the beginning with the aim for execution on parallel architectures. We have implemented a parallel image segmentation algorithm using our model as the foundation. It is able to take advantage of several processing units by dividing the segmentation task into equal parts, which would be carried out by the said units simultaneously. In such a scenario, the workload distribution is performed in the following way: all the nodes of the mesh are divided in as many groups as there are processing units available. Each of the units will have full access to the input data in order to allow full mesh deformation and avoid a situation when the given node is not able to move outside a specific region, because its processing unit does not hold any information about the necessary image data. Naturally this approach imposes the drawback of time required to propagate the input data between all the processing units, but our tests with a computer cluster have shown that this time is not very significant in terms of the whole segmentation time. In a scenario with multicore machine being the execution environment, this time would be even less considerable.
The general idea of the parallel algorithm is presented in Figure 3. Each of the processing units will create a temporary set of variables to use while the segmentation progresses. Those values will be representing the current state of the model on a given processing unit and will not be shared nor synchronized with other participants. This is done to eliminate any dependencies between the processing units and avoid the need for synchronization, which could severely increase the time required for the segmentation. The possibility of performing the optimization in an independent manner is guaranteed by the fact that most of the calculations carried out when minimizing the energy function are performed using only the local data. The only values that are defined as global for the entire model are the average distances and angles between the nodes of the mesh. However, numerous experiments have shown that due to the high rigidity of the WMD mesh structure, their values are modified very slightly during the entire segmentation process and thus there is no need to track their changes very precisely. Therefore, it is safe to assume that the locally calculated value describes well the state of the entire mesh.
Figure 3. General flowchart representation of the parallel algorithm implementation.
The processing units perform the mesh optimization simultaneously using the steps described in Section 2.5 and test Equation 4 in order to decide when to conclude the segmentation.
It is important to note that in the case of TAV the equal distribution of the workload would not be as straightforward as dividing the nodes in equally large groups. This possibility of our model is a result of the approach that we have taken, namely the lack of division between internal and external nodes and by implementing the dynamic topology changes scheme.
Thanks to these feature, as stated before, the time required to perform segmentation with our method depends very little on the contents of the input image. As a consequence, this time is very uniform for different parts of the mesh. On the other hand, in the case of TAV, the distance that the external nodes need to travel during the segmentation task depends highly on the contents of the image, which results in a big variation of the segmentation times for different parts of the mesh. This is the reason why our method is much more suitable for parallelization than TAV or any method based on the concept of contour evolution.
7 Experiments and results
7.1 Performance experiments
We have performed our experiments using two algorithms introduced in the previous section. One of them was an implementation of the TAV model [20] and the other was based on the WMD model. Both models include the ability to perform topological changes. The values of the numerical parameters of energy function that have been used for all the following experiments are the following: α=0.4,β=0.1,γ=0.7,δ=0.5,ε=0.9. As it can be seen, they follow the general rule specified in (5). The algorithm was implemented in C#. The execution environment included a machine equipped with 3GHz processor and 2 GB of RAM. The following tests have been prepared in order to test the most lowlevel features of the model, i.e., the energy formulation, the process of shape optimization and the definition of the mesh, which are the unique features of the WMD model when comparing to the TAV solution or other similar approaches.
7.2 Comparison of segmentation time
As the first input data we have used a set of ten artificial shapes presented in Figure 4. They have been created with the aim of representing features that are often encountered in realworld segmentation tasks, i.e., several objects present in the scene, objects placed inside one another or nonisometric distribution of objects in the image.
Figure 4. A set of artificial images used for the performance experiments.
In each experiment, the 2D image has been used twice, which resulted in a 2slice 3D volume. The results of our experiments are presented in Figure 5. As we can see, the WMD method shows a noteworthy improvement in the segmentation performance. The execution time was always significantly lower than in the case of TAV method, reaching 5.3 times shorter in Scenario 3. The accuracy of segmentation was exactly the same for both approaches. Examples of the segmentation progression, together with the created output, can be seen in Figures 6 and 7. The consecutive images present how the mesh is attracted to the edges of objects and as the result how it deforms and creates discontinuities. The final images present the result with unwanted mesh nodes discarded.
Figure 5. Results of the performance experiments.
Figure 6. Consecutive steps of segmentation algorithm execution using an image with artificial shape as the input. The input image included a simple teapot shape.
This confirms our expectations described in the previous section. What is also very notable is that the WMD model delivered very similar segmentation times in all the scenarios, whereas the TAV model showed a much higher variety in this value. This shows that the segmentation time and the number of iterations are nearly independent from the contents of the images. To explain this let us consider the following.
The size of the mesh in those experiments was 16×16 pixels, which results in 16 pixels of maximum possible distance needed for each node to travel during the segmentation process. This is true if we assume that the segmentation parameters are chosen correctly and every given edge would be approached only by the nodes initialized in their nearest neighborhood. As we can see, the maximum travel distance of a node can be estimated with high confidence without any information about the contents of the input images.
The size of the two input images was 512×512 pixels, which results in mesh of the size 33×33 nodes for each image, giving total number of 1,089 nodes, constructing 1,024 cells of the mesh. The number of mesh cells is relatively high, which lets us assume with high probability that at least in one of the cells a node would have to travel the maximum possible distance. And since the whole segmentation time is defined by the node that has to move for the largest distance of all the nodes, we can assume that this pessimistic scenario is the factual scenario for majority of algorithm executions. Tracing the segmentation process iteration by iteration we have confirmed that the above assumption is correct. Through the initial 16 iterations the nodes are progressing towards the object borders and after that the step of breaking unwanted connections between nodes is taken in order to perform topology changes and separate unwanted parts from the mesh. This part shows to be less predictable and presents some variation, which is mainly responsible for the small differences in required iterations for different scenarios. What can be seen from these experiments is that the WMD model is by definition not highly dependent from the contents of the segmented image and delivers identical results in much shorter times.
7.3 Dependence between the grid density and segmentation speed
We have also compared how the two approaches perform in scenarios of different sizes of the mesh. We have executed the segmentation using Shape 10 from Figure 4 together with 13 different values of internode spacing, namely from 20×20 pixels to 8×8 pixels (in x and yaxis). As it can be seen in Figure 8(left), the TAV method requires exactly the same number of iterations to finish the segmentation, regardless the density of the mesh. In turn, the number of iterations required by the WMD method is decreasing along with the growth of mesh density. The explanation for this behavior lies in the fact that with smaller initial distances between mesh nodes the number of pixels that each node needs to travel is also decreasing, resulting in faster progression towards boundaries of objects. The optimal value of this parameter can be noticed around the initial distance between the nodes of 12×12 pixels. After that threshold the number of necessary iterations starts to depend more on the step of connections breaking and recognition of unwanted parts of the mesh, so further increase in the mesh density will not result in significant decrease of necessary iterations.
Figure 8. Results of segmentation experiments using different values of internode spacing, number of necessary iterations (left) and execution times in seconds (right).
The values of execution times for all those scenarios are depicted in Figure 8(right). The time required to carry out segmentation using the TAV method is growing severely with the growth of mesh complexity, whereas in case of WMD we can notice only a very slight growth of the execution time. As we can see, the smaller number of iterations is balanced with the increased number of nodes to process in the whole mesh. As a result, the execution times are maintained on nearly constant level until the optimal value is reached (12×12 pixels). Then we can see a slightly faster increase of the execution time, as the growing number of nodes in the mesh is no longer balanced by smaller number of iterations.
This low dependance of the segmentation time from mesh density is a much desired feature of the WMD model, as it allows performing precise segmentations of very complex shapes, having almost no drawback of increased execution times. With most of the known segmentation methods, including the TAV, increase of the segmentation accuracy imposes also a significant growth of the segmentation time.
7.4 Performance with realworld medical images
Finally, we wanted to assess the potency of our method when dealing with realworld images. For this purpose, we have used two large 3D volumes presented in part in Figures 9 and 10. The former volume presents a CT scan of human skull. It is composed of 320 slices, 512×512 pixels each. The latter volume presents an XRay scan of a spinal column area and is composed of 380 slices, each of them also of 512×512 dimensions. During the segmentation tasks the mesh density was defined to 8×8 pixels in all the scenarios.
Figure 9. A sample of image slices from the 320slice CT scan volume.
Figure 10. A sample of image slices from the 380slice XRay scan volume.
It is important to notice the large difference in the contents of the consecutive slices in each of the sets. In every case of 3D volume segmentation the whole execution time will highly be dependent from the slice, which requires the biggest amount of time to be processed. In our CT skull scan example, as well as in many similar realworld medical images, we can see that the last slices of the sets contain very little data, located far from the image borders. Performing the volume segmentation using the assumptions of the TAV we would need the nodes of the mesh to travel through nearly the entire image to reach this data. Therefore, the whole segmentation would be very time consuming. We have examined the last slice of the CT volume and found out that the longest distance from the edge of the image to object of interest is 220 pixels. This means that the number of iterations required by the TAV method to complete the segmentation should be around 220. On the other hand, our WMD method is nearly independent from the contents of the image, so the segmentation times should be almost equal for all the slices and the total segmentation time should significantly be lower.
The execution times and the number of iterations for the two test cases are depicted in Figures 11 and 12, respectively. As it can be seen, the WMD model outperforms the TAV model significantly in all the scenarios. The segmentation time of the CT scan volume equalled 689 s for the WMD model and 3320 s for the TAV model. This results in 481% efficiency growth ratio. The number of necessary iterations for WMD and TAV models equalled 45 and 239, respectively. The very large number of iterations required by the TAV model confirms our observations about the required way that the pixels need to travel under the assumptions of this model. The results of the second experiment are of yet higher difference, the execution times equalled 818 s for the WMD model, and 6284 s for the TAV model, which results in 768% efficiency growth ratio. The required iterations equalled respectively 47 and 382. The WMD model maintained a similar number of iterations comparing to the previous experiment and the increased execution time is the result of a larger input volume. The TAV model required a significantly larger number of iterations. This of course resulted in a much longer execution time, which was even further increased by the larger number of slices in the input volume. The higher number of iterations for the TAV model in the second experiment is justified again with the placement of objects of interest in the input volume. In some input slices of the XRay volume (corresponding to the second row of images in Figure 10) we can see that parts of the bone structures are located in the lowerright regions of the image. To close the mesh and to finish the segmentation process, the TAV model needed to move the nodes of the mesh from the topleft corner of the image for the distance of around 370 pixels.
Figure 11. Results of the experiments with real medical images—execution times of the segmentation algorithms.
Figure 12. Results of the experiments with real medical images—iterations required to perform a full segmentation.
These experiments shows how well suited is our model for realworld scenarios, where similar situations can occur often—the slices of a single volume can show a big variety of sizes, which can lead to unreasonably long execution times. Examples of segmentation results can be seen in Figures 13, 14, and 15. The first one presents the result corresponding to the 320slice CT scan volume data. Figures 14 and 15 both present the results corresponding to the 380slice XRay scan volume, the former one presents the results of bone structure segmentation, which was obtained by segmenting the brightest zone of the image, the latter presents the organs, which were obtained by setting the algorithm to segment the medium intensity ranges in the image. It can be seen that results delivered by the WMD model are very precise and reliable, small details are detected successfully. Complex parts of the input images have been segmented correctly and described with surfaces that are smooth and precise. Both WMD and TAV methods present the same level of quality. The change of the shape optimization scheme did not introduce any drawbacks or flaws in the segmentation outcome. WMD model was also much more suitable for working with very high densities of the mesh and very high resolutions of the input images, which results in very precise segmentation results. In fact, there is no reason for which the TAV model could deliver results that are superior to the WMD model, as the segmentation process of these two methods offers exactly the same functionality.
Figure 13. Result of the 320 slice CT scan segmentation.
Figure 14. Results of the 380 slice XRay scan segmentation–bone structure segmentation.
Figure 15. Results of the 380 slice XRay scan segmentation–segmentation of organs.
7.5 Parallel implementation results
The implementation of our parallel algorithm has been performed using the MSMPI [23] message passing interface. The experiments have been executed on a computer cluster constructed from 12 processing units of 3 GHz each, running the Microsoft HPC Server 2008 operating system. As the input image we have used the CT brain scan presented in Figure 9. The times required to carry out the segmentation with different numbers of processing units are presented in Figure 16 left. Furthermore, Figure 16 right presents the speedup obtained with the parallel implementation. This metric is obtained by dividing the execution time of segmentation using single processing unit by the execution time with multiple processing units. The charts show that the parallelization of the method is performed efficiently and the resulting performance improvement is very close to ideal. By this we mean half the execution time for two processing units, quarter for four processing units, and so on. This is a result of the fact that our method does not need any synchronization or exchanging any data between participating processes and also of the fact that the segmentation process takes almost exactly the same time for any part of the input volume.
Figure 16. Execution times and speedup gained with execution of the parallel algorithm on multiprocessing unit environment.
8 Conclusions and future work
In this article, we have presented a novel model for 3D image segmentation, called the WMD Model. This method shows to be highly efficient, giving good segmentation results in short times. More accurately our model shows the following desired features:
• The segmentation procedure is significantly faster than other similar approaches like TAV, thanks to allocation of more nodes into the boundary detection process and to their more optimal distribution.
• The segmentation procedure is predictable and stable. This allows a very efficient parallelization, which is very hard to obtain with other popular methods based on deformable models framework.
• The segmentation time is virtually independent from the density of the mesh. This is a very desirable and also unique feature of our method, as it allows to obtain highquality results without the drawback of increased segmentation time.
• Segmentation time depends also very little from the distribution of the objects in the input image—the only factor that defines the time in a noteworthy manner is the resolution of the input image, but also that has much less impact on the segmentation time than in the case of the TAV model.
• The dependency on user interaction and on initialization of segmentation is significantly lower in the WMD model than in the TAV model.
• Execution on a multiprocessing unit environment is straightforward and offers nearly linear performance gain.
Tests on realworld images have proven that those features allow our solution to perform very well in medical applications. It has been able to detect and segment objects in CT, MRI, and XRay scans of different parts of human body. The WMD model can be implemented on any simple multicore platform, which makes it very easy to apply in any medical office.
Competing interests
The authors declare that they have no competing interests.
Acknowledgements
This study was partially supported by the Portuguese Fundacão para a Ciência e a Tecnologia through the grant contract SFRH/BD/61265/2009, by MSFTMicrosoft Portugal through HIPERCAMBIO project: High Performance Computing Applied to Medical Data and Bioinformatics, by Instituto de Telecomunicacões and by University of Beira Interior.
References

M Kass, A Witkin, D Terzopoulos, Snakes: active contour models. Int. J. Comput. Vis 1(4), 321–331 (1988). Publisher Full Text

V Caselles, Geometric models for active contours. International Conference on Image Processing (NJ, USA: IEEE Press Piscataway, 1995), pp. 9–12

TF Chan, LA Vese, Active contours without edges. IEEE Trans. Image Process 10(2), 266–277 (2001). PubMed Abstract  Publisher Full Text

A Amini, S Tehrani, T Weymouth, Using dynamic programming for minimizing the energy of active contours in the presence of hard constraints. International Conference on Computer Vision (NJ, USA: IEEE Press Piscataway, 5–8 Dec. 1988)

D Williams, M Shah, A fast algorithm for active contours and curvature estimation. CVGIP: Image Understand 55, 14–26 (1992). Publisher Full Text

W Tao, I Cheng, A Basu, Fluid vector flow and applications in brain tumor segmentation. IEEE Trans. Biomed. Eng 56(3), 781–789 (2009). PubMed Abstract  Publisher Full Text

C Xu, J Prince, Snakes, shapes, and gradient vector flow. IEEE Trans. Image Process 7(3), 359–369 (1998). PubMed Abstract  Publisher Full Text

X Xie, Active contouring based on gradient vector interaction and constrained level set diffusion. IEEE Trans. Image Process 19, 154–164 (2010). PubMed Abstract  Publisher Full Text

AD Santis, D Iacoviello, A discrete level set approach to image segmentation. Signal, Image Video Process 1(4), 303–320 (2007). Publisher Full Text

H Delingette, J Montagnat, Shape and topology constraints on parametric active contours. Comput. Vis. Image Understand 83(2), 140–171 (2000)

T Mcinerney, D Terzopoulos, Deformable models in medical image analysis: a survey. Med. Image Anal 1, 91–108 (1996). PubMed Abstract  Publisher Full Text

T Mcinerney, D Terzopoulos, Topologically adaptable snakes. Proceedings of the International Conference on Computer Vision (ICCV’95) Cambridge, MA (NJ, USA: IEEE Press Piscataway, 1995), pp. 840–845

T Mcinerney, D Terzopoulos, Topology adaptive deformable surfaces for medical image volume segmentation. IEEE Trans. Med. Imag 18(10), 840–850 (1999). Publisher Full Text

S Osher, J Sethian, Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. J. Comput. Phys 79, 12–49 (1988). Publisher Full Text

K Fritscher, R Schubert, 3D image segmentation by using statistical deformation models and level sets. Int. J. Comput. Assist. Radiol. Surg 1(3), 123–135 (2006). Publisher Full Text

M Leventon, W Eric, L Grimson, O Faugeras, Statistical shape influence in geodesic active contours. Proceeding of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (NJ, USA: IEEE Press Piscataway, 2000), pp. 316–323

R Malladi, J Sethian, B Vemuri, Shape modeling with front propagation: a level set approach. IEEE Trans. Pattern Anal. Mach. Intell 17(2), 158–175 (1995). Publisher Full Text

X Han, C Xu, J Prince, A topology preserving level set method for geometric deformable models. IEEE Trans. Pattern Anal. Mach. Intell 25(6), 755–768 (2003). Publisher Full Text

F Segonne, J Pacheco, B Fischl, Active contours under topology control genus preserving level sets. Computer Vision for Biomedical Image Applications (New York, NY, USA: Springer US, 2005), pp. 135–145

N Barreira, MG Penedo, Topological active volumes. EURASIP J. Appl. Signal Process 2005, 1939–1947 (2005). Publisher Full Text

N Barreira, M Penedo, C Alonso, C Marino, F Ansia, Handling topological changes in the topological active volumes model. Proceedings of the 5th International Conference on Image Analysis and Recognition (Berlin, Germany: Springer Berlin Heidelberg, 2008), pp. 121–131

R Natarajan, D Krishnaswamy, A case study in parallel scientific computing: the boundary element method on a distributedmemory multicomputer. Proceedings of the IEEE/ACM Supercomputing Conference, San Diego, CA, USA (NY, USA: ACM New York, 1995), pp. 33–33

RI Taniguchi, Y Makiyama, N Tsuruta, S Yonemoto, D Arita, Software platform for parallel image processing and computer vision. in SPIE Proceedings Parallel and Distributed Methods for Image Processing, ed. by Shi H, Coffield PC (UK: SPIE, Cardiff, 2001), pp. 2–10

C Nicolescu, P Jonker, A data and task parallel image processing environment. Parallel Comput 28(7–8), 945–965 (2002)

Y Ataseven, Z AkalinAcar, C Acar, N Gencer, Parallel implementation of the accelerated BEM approach for EMSI of the human brain. Med. Biol. Eng. Comput 46(7), 671–679 (2008). PubMed Abstract  Publisher Full Text

P Lenkiewicz, M Pereira, M Freire, J Fernandes, The whole mesh deformation model for 2D and 3D image segmentation. International Conference on Image Processing (NJ, USA: IEEE Press Piscataway, 2009), pp. 4045–4048

J Canny, A computational approach to edge detection. IEEE Trans. Pattern Anal. Mach. Intell 8(6), 679–698 (http://portal, 1986), . acm.org/citation.cfm?id=11275 webcite PubMed Abstract