1 Introduction
The incorporation of anatomical knowledge into multiregion medical image segmentation has been the subject of countless articles. Recently, research into the specification and incorporation of anatomyagnostic knowledge structures have been undertaken with varying levels of success. For example, regularization has been used initially to encourage segmentation contiguity or adherence to image edge features, most notably popularized by the graph cuts segmentation approach [1, 2, 3]. This framework represented the image as a finite lattice, incorporating regularization in a globally optimal manner, through the use of edge weighting and minimum cost cuts. Such an approach had direct analogues in Markov Random Field theory [4, 5]
allowing for guided development of cost terms using standard probability theory. Such work has been extended from the discrete domain to a continuous domain without the loss of global optimality
[6] mitigating the effects of differing neighbourhood connectivity and associated metrification artifacts.Regularization was then extended to incorporate or encourage spatial grouping relationships, originally in the form of a full ordering using the discrete Ishikawa model [7] and its continuous counterpart [8]. However, these models enforce that a full ordering be defined a priori for the segmentation problem, and were not applicable to segmentation problems outside of that scope. This scope has been extended recently by Delong et al. [9] for the discrete case and Baxter et al. [10] for the continuous, but maintained constraints on which regularization configurations could be specified.
The motivation behind this work is to extend previous general models to eliminate the constraints on what part/whole or ordering relationships can be defined in the continuous case. Thus, an algorithm for solving continuous maxflow/mincut segmentation problems under a directed acyclic graph overarchitecture is developed, along with a framework for expressing arbitrary superobject regularization using said algorithm.
2 Contributions
As with previous work in extensible maxflow segmentation models, this algorithm displays a high degree of inherent parallelism allowing for acceleration through general purpose graphic processing unit (GPGPU) computation, as well as additional overarching concurrency allowing for additional threading and scheduling to improve performance and multicard use.
3 Directed Acyclic Graphical Model and Previous Work
3.1 Previous Work
Work by Yuan et al. [6] addressed both the continuous binary mincut problem and the convex relaxed Potts model:
This work was further extended by Bae et al. [8] to the continuous Ishikawa model:
using similar variational methods but employed a tiered continuous graph analogous to that used by Ishikawa [7] in the discrete case, that is, with finite capacities on intermediate flows between labels.
Models with limited hierarchical constraints such as that used by Rajchl et al.[11] for myocardial scar segmentation have been posed. They have since been generalized by Baxter et al. [10] in the form:
These techniques both used a continuous maxflow model with augmented Lagrangian multipliers from which efficient solution algorithms could be constructed. Apart from regularization structure, constraints such as starshaped constraints on the various labels[12], as well as volume preserving, and interimage consistency has been incorporated.
3.2 Directed Acyclic Graphical Model
Directed Acyclic Graphical MaxFlow (DAGMF) segmentation relies heavily on the concept of a rooted directed acyclic graph (DAG) with weighted edges for label representation. Rooted directed acyclic graphs extend the notion of part/whole or parent/child relationships used in previous models, specifically hierarchical models[11, 10]. For now, a discussion of interpretation and whole/part relationships will be postponed until Section 6. Figure 1 provides an example of this structure.
Throughout this paper, we will refer to parent and child operators, and respectively. These reflect the direction of flow (always from parent to children) and multiplier accumulation (always from child to parents). In addition, edges in this graph, and corresponding parent/child operators can have multiplicity, that is, multiple edges can exist between two vertices indicating some rational weighting of relative parts. The parentchild operators for the graph in Figure 1 are:
As with the general hierarchical formulation [10], these operators have consistency requirements. That is, they must form a rooted DAG with the properties:

there exists only one label, , has no parent (),

the parent/child relationship is preserved (),

there are no cycles (no label can be a child, grandchild, greatgrandchild, etc… of itself) and that the graph is connected. Note that if this property holds for the child operator, it must also hold for the parent operator due to the above properties, and

the edge weights are nonnegative and normalized for each child, that is .
As with the Potts model [13] we would like to define the set of endlabels, , as a partition of the image, that is:
(1) 
Unlike those in the Ishikawa [7] model and generalized hierarchical maxflow model [10], the intermediate labels do not immediately lend themselves to a settheoretic interpretation since they are not constrained to be the union of endlabels. (Although, this case will be explicitly explored in Section 6 with corresponding DAG structure.) Nevertheless, the labeling function, has the following properties:[6, 14, 15]
(2) 
(3) 
for labels in . In the case where the DAG forms a hierarchy, the same constraints apply as in the generalized hierarchical model. Otherwise, the nonunit multiplicands render the union operator meaningless. In terms of labeling function, the intermediate nodes are defined as:
(4) 
These yield the convex relaxed generalized hierarchical model:
(5)  
subject to  
These formulations can be solved with global optimality for probabilistic labels, which will be demonstrated by this paper. However, it is only an approximation algorithm for models under integrality constraints on the endlabels, as demonstrated by the NPhardness of the Potts model, even over finite lattices[16].
As with previous approaches, we only consider the case where the endlabels have nonzero data terms. To prove that this does not limit the applicability of the method, a similar lineartime data term pushdown proof can be constructed along the lines presented by Baxter at al.[10]
4 Continuous MaxFlow Model
4.1 Primal Formulation
The modeling approach is derived from those presented by Yuan et al.[6, 14]] and follows the same format, using duality through an augmented Langrangian formulation. The primal model represents network flow maximization through a large graph with only the sink flows constrained. The dual of this formulation is the DAGMF equation (5) as we shall prove in this section. We can write the primal model as:
(6) 
subject to the constraints
(7)  
This is equivalent to a multiflow problem over a large graph constructed from the image dimensions and the provided directed acyclic graph as overall architecture. Other than constraints put on the magnitude of the spatial flows, and the capacity of the sink flows ( where ), the system is assumed to have infinite capacity. This is a strict generalization of the hierarchical formed explored by Baxter et al.[10] considering hierarchies to be a specific class of rooted DAG.
4.2 PrimalDual Formulation
The primal model can be converted to a primaldual model through the use of Lagrangian multipliers over the flow conservation constraint yielding the Lagrangian:
(8)  
First, we must ensure that equation (8) is convex with respect to , considering to be fixed, and concave with respect to with fixed, as to meet the requirements of the minimax theorem. [17] Considering as fixed, is obviously fixed as well, implying that equation (8) is linear over and therefore convex. It should also be noted that is a linear function of , meaning that (8) is again linear and therefore concave with respect to This implies the existence of a saddle point and the equivalence of the formulation regardless of the order of the prefix max and min operators under the minimax theorem. [17]
4.3 Dual Formulation
As implied in the previous section, we can find the saddle point through the optimization of the sinkflows, , working bottomup and the spatial flows within each label. For the sake of simplicity, we refer to an ordering generated by topologically sorting the graph . We proceed through the graph in that order. Starting with any label, , such that we can isolate in (8) giving:
(9)  
when . (If , the function can be arbitrarily maximized by .) Working through , every label, , where and can be isolated in (8) as:
(10) 
at the saddle point defined by . Lastly, the source flow, , can be isolated in a similar manner, that is:
(11) 
at the saddle point defined by . These constraints combined yield the labeling constraints in the original formulation. The maximization of the spatial flow functions can be expressed in a wellstudied form [18] as:
(12)  
The above implies that we can express the saddle point of equation (8) as the original energy functional, (5), and therefore, finding the saddle point of (8) is equivalent to solving the DAGMF segmentation problem.
5 Solution to PrimalDual Formulation
To address the optimization problem, we can find this saddle point by augmenting the Lagrangian function [19]:
(13)  
where is a positive penalty parameter encouraging faster convergence to solutions that fulfill the optimization constraints. Using this formula, we can iteratively maximize each component. The solution steps are:
Note that within each step there exists a large amount of inherent parallelism, that is, each voxel can be accounted for completely independently of all other voxels in steps 25, and with only a dependence on a local neighbourhood in step 1. This inherent parallelism allows for GPGPU acceleration of each step with relative ease.
5.1 Directed Acyclic Graphical MaxFlow Algorithm
To improve the convergence rate, we perform an initialization step that ensures optimality for the zerosmoothness condition. This is achieved by initializing the system with optimal flows and multipliers under the assumption that all spatial flows are zero. To ensure faster convergence, we order the tasks using a topological sort over the graph. In this ordering, each child label occurs only after all of its parents. The inverse ordering is the opposite. This is equivalent to the bottomup approach used by Baxter et al. [10] over hierarchies. In addition, each label is equipped with two ‘working’ buffers, and , for the purposes of accumulation.
which makes use of the following function definitions:
For the sake of conciseness, the ‘for do ’ loops surrounding each assignment operation have been replaced with the prefix in both the algorithm and the function definitions.
6 Regularization of Arbitrary SuperObjects
As demonstrated in the introduction, continuous maxflow segmentation models have been progressing towards more and more general regularization structures. In terms of a discrete analogue, the most general structure possible would be:
(14) 
in which could refer to either an endlabel in or a subset with . We refer to this problem as having arbitrary superobject regularization in that it incorporates all possible regularization while maintaining the discrete analogue for intermediate labels. It is easy to see it as a generalization of Potts[13], Ishikawa[7], and General Hierarchical[10] models. This section aims to illustrate this as a subclass of that represented by DAGMF.
To show how arbitrary superobject regularization can be implemented with DAGMF, we must consider the construction of a DAG with associated transformations on smoothness parameters. First, let us add a vertex to the graph for the source node, , and one for each endlabel in . To do so, let us consider to be the set of superobjects represented in the problem, not including endlabels. We will associate a vertex in the DAG to each element, , in this set. Each of these vertices has as their sole parent, and their only children will be vertices corresponding to the endlabels in . Lastly, we must ensure that all endlabel vertices have the same number of parents, so we can add edges between and each vertex (allowing multiplicity) until they do. (Note that each endlabel will have parents.) An example of this is given in Figure 2.
Before the algorithm can be used, we must determine edge weights and eliminate multiplicities in the edges. We do this by associated the unnormalized edge weight with the multiplicity of the edge, followed by a normalization step to ensure the weights are valid. The corresponding example is given in Figure 3.
Returning to the general case, if we consider the equation associated with this graph, we will find it is:
(15)  
subject to  
which is equivalent to the following by multiplying the labelling functions for by
(16)  
subject to  
This indicates that the intended regularization field for must be multiplied by for correct scaling, and for the union operator to be a discrete analogue for addition (since the multiplicand is removed in the constraint equations).
7 Conclusions
In this paper, we present an algorithm for addressing maxflow segmentation problems where the underlying architecture is a directed acyclic graph. Using this architecture, we can achieve arbitrary superobject regularization, making it more general version of Generalized Hierarchical MaxFlow[10] and thus both Potts and Ishikawa models. Such a regularization mechanism allows for the definition of more flexible part/whole relationships, where any individual part could belong to more than one whole, a capability not present in any earlier extendable model. These also represent the most general form of ordering relationships possible while maintaining an analogous purely discrete problem.
This solver has been implemented using the NVIDIA Compute Unified Device Architecture (CUDA) taking advantage of the inherent parallelism within each optimization step. Additional concurrency between steps has been exploited to improve computational speed and allow for multiple graphics cards to be used simultaneously on a single segmentation problem.
Acknowledgements.
The authors would like to acknowledge Dr. Elvis Chen, Kamyar Abhari, and Jonathan McLeod for their invaluable discussion, editing, and technical support.References
 [1] Y. Y. Boykov and M.P. Jolly, “Interactive graph cuts for optimal boundary & region segmentation of objects in ND images,” in Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on, 1, pp. 105––112, 2001.
 [2] Y. Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” IEEE Transactions on Pattern Analysis and Machine Intelligence 23(11), pp. 1222–1239, 2001.
 [3] Y. Boykov and V. Kolmogorov, “An experimental comparison of mincut/maxflow algorithms for energy minimization in vision,” Pattern Analysis and Machine Intelligence, IEEE Transactions on 26(9), pp. 1124––1137, 2004.

[4]
Y. Boykov, O. Veksler, and R. Zabih, “Markov random fields with efficient
approximations,” in
Computer vision and pattern recognition, 1998. Proceedings. 1998 IEEE computer society conference on
, pp. 648––655, 1998.  [5] S. Z. Li, Markov Random Field Modeling in Image Analysis, Springer, Jan. 2009.
 [6] J. Yuan, E. Bae, and X.C. Tai, “A study on continuous maxflow and mincut approaches,” in Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on, pp. 2217––2224, 2010.
 [7] H. Ishikawa, “Exact optimization for markov random fields with convex priors,” IEEE Transactions on Pattern Analysis and Machine Intelligence 25(10), pp. 1333–1336, 2003.
 [8] E. Bae, J. Yuan, X.C. Tai, and Y. Boykov, “A fast continuous maxflow approach to nonconvex multi–labeling problems,” 2011.
 [9] A. Delong and Y. Boykov, “Globally optimal segmentation of multiregion objects,” in Computer Vision, 2009 IEEE 12th International Conference on, pp. 285–292, IEEE, 2009.
 [10] J. S. Baxter, M. Rajchl, J. Yuan, and T. M. Peters, “A continuous maxflow approach to general hierarchical multilabelling problems,” arXiv preprint arXiv:1404.0336 , 2014.
 [11] M. Rajchl, J. Yuan, J. White, E. Ukwatta, J. Stirrat, C. Nambakhsh, F. Li, and T. Peters, “Interactive hierarchical maxflow segmentation of scar tissue from lateenhancement cardiac MR images,” IEEE Transactions on Medical Imaging , 2014.
 [12] J. Yuan, W. Qiu, E. Ukwatta, M. Rajchl, Y. Sun, and A. Fenster, “An efficient convex optimization approach to 3D prostate MRI segmentation with generic star shape prior,” Prostate MR Image Segmentation Challenge, MICCAI , 2012.
 [13] R. B. Potts, “Some generalized orderdisorder transformations,” in Proceedings of the Cambridge Philosophical Society, 48, pp. 106––109, 1952.
 [14] J. Yuan, E. Bae, X.C. Tai, and Y. Boykov, “A continuous maxflow approach to potts model,” in Computer Vision – ECCV 2010, K. Daniilidis, P. Maragos, and N. Paragios, eds., Lecture Notes in Computer Science, pp. 379–392, Springer Berlin Heidelberg, Jan. 2010.
 [15] E. Bae, J. Yuan, and X.C. Tai, “Global minimization for continuous multiphase partitioning problems using a dual approach,” International journal of computer vision 92(1), pp. 112––129, 2011.
 [16] V. Kolmogorov and R. Zabih, “Computing visual correspondence with occlusions using graph cuts,” in Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on, 2, pp. 508––515, 2001.
 [17] I. Ekeland and R. Temam, “Convex analysis and variational problems,” 1976.
 [18] E. Giusti, Minimal surfaces and functions of bounded variation, vol. 80, Birkhauser, 1984.
 [19] D. P. Bertsekas, “Nonlinear programming,” 1999.
 [20] A. Chambolle, “An algorithm for total variation minimization and applications,” Journal of Mathematical imaging and vision 20(12), pp. 89––97, 2004.
Comments
There are no comments yet.