The incorporation of anatomical knowledge into multi-region medical image segmentation has been the subject of countless articles. Recently, research into the specification and incorporation of anatomy-agnostic 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 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  and its continuous counterpart . 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.  for the discrete case and Baxter et al.  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 max-flow/min-cut segmentation problems under a directed acyclic graph over-architecture is developed, along with a framework for expressing arbitrary super-object regularization using said algorithm.
As with previous work in extensible max-flow 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 over-arching concurrency allowing for additional threading and scheduling to improve performance and multi-card use.
3 Directed Acyclic Graphical Model and Previous Work
3.1 Previous Work
Work by Yuan et al.  addressed both the continuous binary min-cut problem and the convex relaxed Potts model:
This work was further extended by Bae et al.  to the continuous Ishikawa model:
using similar variational methods but employed a tiered continuous graph analogous to that used by Ishikawa  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. for myocardial scar segmentation have been posed. They have since been generalized by Baxter et al.  in the form:
These techniques both used a continuous max-flow model with augmented Lagrangian multipliers from which efficient solution algorithms could be constructed. Apart from regularization structure, constraints such as star-shaped constraints on the various labels, as well as volume preserving, and inter-image consistency has been incorporated.
3.2 Directed Acyclic Graphical Model
Directed Acyclic Graphical Max-Flow (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 parent-child operators for the graph in Figure 1 are:
As with the general hierarchical formulation , 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, great-grandchild, 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 non-negative and normalized for each child, that is .
As with the Potts model  we would like to define the set of end-labels, , as a partition of the image, that is:
Unlike those in the Ishikawa  model and generalized hierarchical max-flow model , the intermediate labels do not immediately lend themselves to a set-theoretic interpretation since they are not constrained to be the union of end-labels. (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]
for labels in . In the case where the DAG forms a hierarchy, the same constraints apply as in the generalized hierarchical model. Otherwise, the non-unit multiplicands render the union operator meaningless. In terms of labeling function, the intermediate nodes are defined as:
These yield the convex relaxed generalized hierarchical model:
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 end-labels, as demonstrated by the NP-hardness of the Potts model, even over finite lattices.
As with previous approaches, we only consider the case where the end-labels have non-zero data terms. To prove that this does not limit the applicability of the method, a similar linear-time data term pushdown proof can be constructed along the lines presented by Baxter at al.
4 Continuous Max-Flow 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:
subject to the constraints
This is equivalent to a multi-flow 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. considering hierarchies to be a specific class of rooted DAG.
4.2 Primal-Dual Formulation
The primal model can be converted to a primal-dual model through the use of Lagrangian multipliers over the flow conservation constraint yielding the Lagrangian:
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.  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. 
4.3 Dual Formulation
As implied in the previous section, we can find the saddle point through the optimization of the sink-flows, , working bottom-up 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:
when . (If , the function can be arbitrarily maximized by .) Working through , every label, , where and can be isolated in (8) as:
at the saddle point defined by . Lastly, the source flow, , can be isolated in a similar manner, that is:
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 well-studied form  as:
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 Primal-Dual Formulation
To address the optimization problem, we can find this saddle point by augmenting the Lagrangian function :
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 2-5, 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 Max-Flow Algorithm
To improve the convergence rate, we perform an initialization step that ensures optimality for the zero-smoothness 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 bottom-up approach used by Baxter et al.  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 Super-Objects
As demonstrated in the introduction, continuous max-flow 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:
in which could refer to either an end-label in or a subset with . We refer to this problem as having arbitrary super-object 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, Ishikawa, and General Hierarchical models. This section aims to illustrate this as a subclass of that represented by DAGMF.
To show how arbitrary super-object 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 end-label in . To do so, let us consider to be the set of super-objects represented in the problem, not including end-labels. 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 end-labels in . Lastly, we must ensure that all end-label vertices have the same number of parents, so we can add edges between and each vertex (allowing multiplicity) until they do. (Note that each end-label 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:
which is equivalent to the following by multiplying the labelling functions for by
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).
In this paper, we present an algorithm for addressing max-flow segmentation problems where the underlying architecture is a directed acyclic graph. Using this architecture, we can achieve arbitrary super-object regularization, making it more general version of Generalized Hierarchical Max-Flow 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.
-  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.
-  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.
-  Y. Boykov and V. Kolmogorov, “An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision,” Pattern Analysis and Machine Intelligence, IEEE Transactions on 26(9), pp. 1124––1137, 2004.
Y. Boykov, O. Veksler, and R. Zabih, “Markov random fields with efficient
Computer vision and pattern recognition, 1998. Proceedings. 1998 IEEE computer society conference on, pp. 648––655, 1998.
-  S. Z. Li, Markov Random Field Modeling in Image Analysis, Springer, Jan. 2009.
-  J. Yuan, E. Bae, and X.-C. Tai, “A study on continuous max-flow and min-cut approaches,” in Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on, pp. 2217––2224, 2010.
-  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.
-  E. Bae, J. Yuan, X.-C. Tai, and Y. Boykov, “A fast continuous max-flow approach to non-convex multi–labeling problems,” 2011.
-  A. Delong and Y. Boykov, “Globally optimal segmentation of multi-region objects,” in Computer Vision, 2009 IEEE 12th International Conference on, pp. 285–292, IEEE, 2009.
-  J. S. Baxter, M. Rajchl, J. Yuan, and T. M. Peters, “A continuous max-flow approach to general hierarchical multi-labelling problems,” arXiv preprint arXiv:1404.0336 , 2014.
-  M. Rajchl, J. Yuan, J. White, E. Ukwatta, J. Stirrat, C. Nambakhsh, F. Li, and T. Peters, “Interactive hierarchical max-flow segmentation of scar tissue from late-enhancement cardiac MR images,” IEEE Transactions on Medical Imaging , 2014.
-  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.
-  R. B. Potts, “Some generalized order-disorder transformations,” in Proceedings of the Cambridge Philosophical Society, 48, pp. 106––109, 1952.
-  J. Yuan, E. Bae, X.-C. Tai, and Y. Boykov, “A continuous max-flow 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.
-  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.
-  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.
-  I. Ekeland and R. Temam, “Convex analysis and variational problems,” 1976.
-  E. Giusti, Minimal surfaces and functions of bounded variation, vol. 80, Birkhauser, 1984.
-  D. P. Bertsekas, “Nonlinear programming,” 1999.
-  A. Chambolle, “An algorithm for total variation minimization and applications,” Journal of Mathematical imaging and vision 20(1-2), pp. 89––97, 2004.