Markov random fields, specifically in the form of the Ishikawa model
, have previously been used to address the issue of reconstructing an image from noisy data using probabilistic priors and maximum a posteriori probability optimization. The benefit of such frameworks is their submodularity and therefore global optimizability via graph cuts. The primary improvement of Ishikawa models over previous graph cut approaches in reconstructing an image or depth field is that the topology of the intensity space, an ordered linear topology, is fundamentally built into the functional.
Continuous max-flow approaches have traditionally been seen as continuous analogs to earlier discrete models. Along these lines, a continuous Ishikawa model  was proposed which maintains the same linear topology over the intensity space but using a continuous, rather than discrete domain. This is especially crucial in image reconstruction as the mitigation of discrete metrification artifacts achieved by the continuous formulation  has an immediate visual effect on the reconstruction.
Certain types of fields, however, do not have ordered linear intensities. For example, hue in HSL and HSV colour spaces or phase in a complex-valued image represents a cyclic, rather than linear, topology. These require a specialized reconstruction optimizer which, as the Ishikawa model was to linear intensity models, sensitive to the intensity’s cyclic topology. Although this could in theory be accomplished through a directed acyclic graph continous max-flow (DAGMF) model as described by Baxter et al., the quadratic growth of said models is undesirable, especially for medical applications such as phase processing in susceptibility weighted magnetic resonance imaging in which high intensity resolution is necessary.
The main contribution of this work is the development of a continuous max-flow model describing cyclic field image reconstruction, that is, reconstruction of an image in which the image intensity is defined over a cyclic topology rather than a linear topology as explored by Bae et al.. This reconstruction algorithm is better suited for the reconstruction of cyclic fields such as colour/hue in natural images or phase in complex-valued magnetic resonance images.
3 Previous Work
Previous work by Yuan et al.  has addressed the continuous binary min-cut problem:
as well as the convex relaxed continuous Potts Model:
Both of these techniques both used a continuous max-flow model with augmented Lagrangian multipliers. In the case of the convex-relaxed continuous Potts model, the source flow had infinite capacity and the costs in the functional corresponded with constraints on the sink flows.
Bae et al.  extended the work on the continuous binary min-cut problem to the continuous Ishikawa model:
using similar variational methods but a tiered continuous graph analogous to that used by Ishikawa  in the discrete case, that is, with finite capacities on intermediate flows between labels.
Both continuous Potts and Ishikawa models are important to consider as they were the first max-flow models to be used in image reconstruction. However, they are not topologically suitable for cyclic fields because the Potts model disallows any notion of label topology, and the Ishikawa model is defined only over linear orderings. Additionally, cyclic orderings are non-hierarchical, making them inaccessible to generalized hierarchical max-flow solvers and directed acyclic graph max-flow models of these orderings require quadratic time to handle linear increases in the intensity resolution.
4 Cyclic Topology in Continuous Max-Flow Segmentation
4.1 Continuous Max-Flow Model
The desired optimization model for cyclic field reconstruction would be:
in which is defined over a cyclic manifold to prevent modulation errors.
One can re-express the reconstructed intensity with respect to a series of fuzzy indicator functions:
Using these indicator functions, the data term can be re-written as:
and the smoothness term can be re-written as:
yielding the alternative minimization formulation:
The benefit of this new formulation via indicator functions is that it resembles the previously explored continuous Potts model with the exception that the gradient magnitude operator is defined over the Cartesian product of the spatial domain with a cyclic intensity domain, . Thus, the spatial domain augmented with the direction forms a high dimensional cylinder as shown in Figure 1.
4.2 Primal Formulation
As in the previous work, the primal formulation will be expressed as a max-flow optimization, specifically:
subject to the capacity constraints:
and the flow conservation constraint:
, the spatial flow now has an additional dimension, the direction. The direction is equipped with a cyclic topology and therefore is not susceptible to modulation artifacts. Similar to , makes use of the spatial domain augmented with the cyclic domain of as shown in Figure 1.
4.3 Primal-Dual Formulation
The primal model can be converted to a primal-dual model through the use of Lagrangian multipliers on the flow conservation constraint . This yields the equation:
To ensure that this function meets the criteria of the minimax theorem, we must ensure that it is convex with respect to , considering to be fixed, and concave with respect to with fixed.  For the first, it is sufficient to note that if are fixed, then is fixed as well, meaning that (4) is linear and therefore convex with respect to . It should also be noted that is a linear function of , implying the linearity (and concavity) of (4) with respect to . Thus, a saddle point must exist and the formulations are equivalent. 
4.4 Dual Formulation
To show the equivalence of the primal-dual formulation (Eq. (4)) and the desired dual formulation (Eq. (1)), we must show that the original formula can be reconstructed from the definition of the saddle-point. First, consider the primal-dual formulation with the sink flows, isolated:
which implies that as if , the minimization is unbounded as . This reconstructs the data term portion of the original formulation using the fuzzy labeling function.
The source flow, in Eq. 4 can be isolated as:
at the saddle point defined by . This reconstructs the labeling constraint, that is, that only one intensity value is assigned to each voxel.
Considering coupling the spatial domain with the cyclic domain of via a simple Cartesian product with members taking the form . In this coupled space, the maximization of the spatial flow functions can be expressed  as:
This yields the smoothness component in Eq. (1) thus showing the equivalence of the two models.
4.5 Solution to Primal-Dual Formulation
The optimization problem being addressed is the augmented form of Eq. (4):
in which is a non-negative constant which has no effect on Eq. (4) when the constraint holds, but serves as an additional penalty. Optimizing Eq. (8) can be achieved by performing the following steps iteratively, optimizing Eq. (8) with respect to each individual variable:
The specific augmented Lagrangian based algorithm used is:
4.6 Pseudo-Flow Formulation
The approach taken here for the use of proximal Bregman or pseudo-flow methods is derived from Baxter et al.. The corresponding non-smooth pseudo-flow model can be written as:
which can be derived from Eq (4) by:
Now that a pseudo-flow representation is developed, we can take advantage of Bregman proximal projections to optimize this formula. Consider the distance between labeling functions and as:
which can be verified to be a Bregman distance (when ) using the entropy function:
If we consider a feasible labeling, , we can find another proximal labeling, , which has a lower energy by addressing the optimization:
where is a positive constant. Using a Lagrangian multiplier on the constraint , we can solve for analytically as:
noting that this answer fulfills the constraint provided the same holds for . By letting the distance weighting parameter approach 0, if . Using that fact,
illustrating that the Bregman proximal method is a smoothed version of the non-smooth equation Eq (9).
4.7 Solution to Pseudo-Flow Formulation
By taking the gradient of equation (11) with respect to (with substituted by equation (12)) one can derive the appropriate Chambolle iteration scheme for maximizing (11) with respect to the spatial flow variables :
where is a positive gradient descent parameter. Thus, the proximal Bregman, or pseudo-flow, algorithm relies on the alternation between the analytic optimization of the label functions, through Eq. (12) and the specified Chambolle iteration.
The proximal Bregman based algorithm proposed is therefore:
5 Discussion and Conclusions
In this paper, we present two algorithms for addressing continuous max-flow image reconstruction in which the intensity being reconstructed is fundamentally organized along a cyclic, as opposed to linear, topology. Such a reconstruction method is more efficient than that developed from previous continuous max-flow methods capable of incorporating said topology, specifically DAGMF, by displaying linear rather than quadratic growth in complexity in terms of the intensity resolution. This solver has been implemented using GPGPU acceleration due to its inherent parallelizability. This solver is available open-source for both two- and three-dimensional images at www.advancedsegmentationtools.org.
Future work in continuous max-flow based image reconstruction includes the robust automatic detection of gradient direction which may allow for anisotropic smoothness terms to be employed. These terms could be used to estimate level sets in the image, maintaining them while allowing for faster or less-constrained gradients in the direction orthogonal to said level set.
Acknowledgements.The authors would like to acknowledge Zahra Hosseini and Dr. Martin Rajchl for their invaluable discussion and editing.
-  H. Ishikawa, “Exact optimization for markov random fields with convex priors,” Pattern Analysis and Machine Intelligence, IEEE Transactions on 25(10), pp. 1333–1336, 2003.
-  V. Kolmogorov and R. Zabin, “What energy functions can be minimized via graph cuts?,” Pattern Analysis and Machine Intelligence, IEEE Transactions on 26(2), pp. 147–159, 2004.
-  Y. Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” Pattern Analysis and Machine Intelligence, IEEE Transactions on 23(11), pp. 1222–1239, 2001.
-  V. Kolmogorov and R. Zabih, “Multi-camera scene reconstruction via graph cuts,” in Computer Vision ECCV 2002, pp. 82–96, Springer, 2002.
-  E. Bae, J. Yuan, X.-C. Tai, and Y. Boykov, “A fast continuous max-flow approach to non-convex multi-labeling problems,” in Efficient Algorithms for Global Optimization Methods in Computer Vision, pp. 134–154, Springer, 2014.
J. Yuan, E. Bae, and X.-C. Tai, “A study on continuous max-flow and min-cut
Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on, pp. 2217–2224, IEEE, 2010.
-  J. S. Baxter, M. Rajchl, J. Yuan, and T. M. Peters, “A continuous max-flow approach to multi-labeling problems under arbitrary region regularization,” arXiv preprint arXiv:1405.0892 , 2014.
-  E. M. Haacke, Y. Xu, Y.-C. N. Cheng, and J. R. Reichenbach, “Susceptibility weighted imaging (swi),” Magnetic Resonance in Medicine 52(3), pp. 612–618, 2004.
-  J. Yuan, E. Bae, X.-C. Tai, and Y. Boykov, “A continuous max-flow approach to potts model,” in Computer Vision–ECCV 2010, pp. 379–392, Springer, 2010.
-  J. S. Baxter, M. Rajchl, J. Yuan, and T. M. Peters, “A continuous max-flow approach to general hierarchical multi-labeling problems,” arXiv preprint arXiv:1404.0336 , 2014.
-  I. Ekeland and R. Temam, “Convex analysis and 9 variational problems,” 1976.
-  E. Giusti, Minimal surfaces and functions of bounded variation, no. 80, Springer Science & Business Media, 1984.
-  A. Chambolle, “An algorithm for total variation minimization and applications,” Journal of Mathematical imaging and vision 20(1-2), pp. 89–97, 2004.
-  J. S. Baxter, M. Rajchl, J. Yuan, and T. M. Peters, “A proximal bregman projection approach to continuous max-flow problems using entropic distances,” arXiv preprint arXiv:1501.07844 , 2015.
-  L. M. Bregman, “The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming,” USSR computational mathematics and mathematical physics 7(3), pp. 200–217, 1967.