Crouzeix–Raviart Approximation of the Total Variation on Simplicial Meshes

Abstract

We propose an adaptive implementation of a Crouzeix–Raviart-based discretization of the total variation, which has the property of approximating from below the total variation, with metrication errors only depending on the local curvature, rather than on the orientation as is usual for other approaches.

Appendices

A Proof of the Duality Theorem 1

The goal of this section is to provide a proof of the duality theorem, Theorem 1. We prove a more general result for the CR approximation of Sobolev semi-norms.

1.1 A.1 Almost Constant Crouzeix–Raviart Functions

To simplify, we drop the scale parameter h which is useless in this section. We consider a mesh ${\mathcal {T}}$ of d-dimensional simplices in a polygonal domain $\Omega \subset {\mathbb {R}}^d$. We denote $V({\mathcal {T}})$ the non-conforming P1 functions (Crouzeix–Raviart), defined in (10). We then let

$$\begin{aligned} V^0({\mathcal {T}}) = \left\{ u \in CR({\mathcal {T}}): \int _T u(x) \mathrm{d}x=u(c_T)=0 \ \forall T\in {\mathcal {T}}\right\} . \end{aligned}$$

It is the space of the functions which are 0 on average in each simplex (in other words, the average of the middle values of the facets vanishes, equivalently the value at the center of each simplex is zero).

We define $P0({\mathcal {T}})^d\approx ({\mathbb {R}}^d)^{\mathcal {T}}$ as usual as the space of “P0” functions which are constant on each $T\in {\mathcal {T}}$. Endowed with the topology of $L^2(\Omega ;{\mathbb {R}}^d)$, it is a Euclidean space with the weighted scalar product: for $p,q\in P0({\mathcal {T}})$,

$$\begin{aligned} \int _\Omega p(x)q(x) \mathrm{d}x = \sum _{T\in {\mathcal {T}}} |T|p_T\cdot q_T. \end{aligned}$$

Then, we consider the gradients

$$\begin{aligned} GV^0({\mathcal {T}}) = \left\{ Du : u\in V^0({\mathcal {T}}) \right\} \subset P0({\mathcal {T}}). \end{aligned}$$

We want to characterize this space and its orthogonal. In order to do this, we consider the space $RT0({\mathcal {T}})$ of the first-order Raviart–Thomas vector fields subject to the mesh ${\mathcal {T}}$ (cf Sect. 2.3, these are defined by their fluxes through the edges of the elements $T\in {\mathcal {T}}$). As before we also let $RT0_0({\mathcal {T}})\subset RT0({\mathcal {T}})$ the RT0 fields with zero flux through $\partial \Omega$.

We know that cf Lemma 2.4:

$$\begin{aligned}&\left\{ Du : u\in V({\mathcal {T}})\right\} = \left\{ p\in P0({\mathcal {T}})^d: \int _\Omega p(x)\cdot \sigma (x)\mathrm{d}x=0\right. \\&\quad \left. \ \forall \sigma \in RT0_0({\mathcal {T}}) , \hbox {div}\,\sigma =0\right\} . \end{aligned}$$

More generally, if $u\in V({\mathcal {T}})$ and $\sigma \in RT0_0({\mathcal {T}})$, one has thanks to (11):

$$\begin{aligned} \int _\Omega \sigma (x)\cdot Du(x) \mathrm{d}x = - \int _\Omega \hbox {div}\,\sigma (x) u(x) \mathrm{d}x. \end{aligned}$$

As Du and $\hbox {div}\,\sigma$ are constant on each triangle, this can also be written:

$$\begin{aligned} \sum _{T\in {\mathcal {T}}} |T|\sigma (c_T) \cdot Du(T) = -\sum _{T\in {\mathcal {T}}} |T|\hbox {div}\,\sigma (T) u(c_T) \end{aligned}$$

(41)

where $c_T$ is the center of mass of the simplex T (so that given any affine function a(x), $\int _T a(x)\mathrm{d}x = |T|a(c_T)$).

Hence in particular, for any ${\mathbf {p}}\in GV^0({\mathcal {T}})$ and $\sigma \in RT0_0({\mathcal {T}})$, $\int _\Omega \sigma \cdot {\mathbf {p}}\mathrm{d}x = 0$. A natural question is whether this is an if and only if; that is, if the orthogonal of $\{ (\sigma (c_T))_T\in P0({\mathcal {T}}): \sigma \in RT0_0({\mathcal {T}})\}$ is $GV^0({\mathcal {T}})$.

Assume ${\mathbf {p}}\in P0({\mathcal {T}})^d$ is such that $\int _\Omega \sigma \cdot {\mathbf {p}}\mathrm{d}x = 0$ for all $\sigma \in RT0_0({\mathcal {T}})$. Then in particular it is orthogonal to fields with zero divergence and there exists $u\in V({\mathcal {T}})$ with $Du={\mathbf {p}}$, thanks to Lemma 2.4. In particular, because of (41) we have

$$\begin{aligned} 0 = -\sum _{T\in {\mathcal {T}}} |T|\hbox {div}\,\sigma (T) u(c_T) \end{aligned}$$

for all Raviart–Thomas field $\sigma$ (vanishing on $\partial \Omega$). Now, given any inner facet $S\subset \partial T\cap \partial T'$, $T,T'\in {\mathcal {T}}$, $T\ne T'$, we can introduce the field $\sigma$ with flux 1 through S from T to $T'$ and zero through all other facets. Then the formula becomes

$$\begin{aligned} 0 = u(c_{T'})-u(c_T). \end{aligned}$$

This shows that u must take the same value in all the centers of the vertices. As u is defined up to a constant (in each connected component of $\Omega$), we can assume this value is zero, and we have shown that in $P0({\mathcal {T}})^d$,

$$\begin{aligned} \left\{ (\sigma (c_T))_{T\in {\mathcal {T}}} : \sigma \in RT0_0({\mathcal {T}})\right\} ^\perp = GV^0({\mathcal {T}}). \end{aligned}$$

(42)

In particular, as a consequence, also

$$\begin{aligned} GV^0({\mathcal {T}})^\perp = \left\{ (\sigma (c_T))_{T\in {\mathcal {T}}} : \sigma \in RT0_0({\mathcal {T}})\right\} . \end{aligned}$$

(43)

1.2 A.2 Duality for Discrete Sobolev Semi-norms

Let us introduce, for $u\in V({\mathcal {T}})$, $p\in [1,\infty )$,

$$\begin{aligned} J_p(u) = \frac{1}{p}\int _\Omega |Du|^p \mathrm{d}x \end{aligned}$$

(44)

and, for ${\bar{u}}\in P0({\mathcal {T}})$

$$\begin{aligned} J_p^0({\bar{u}}) = \min \left\{ J_p(u): u\in V({\mathcal {T}}), u(c_T)={\bar{u}}_T \ \forall T\in {\mathcal {T}}\right\} .\nonumber \\ \end{aligned}$$

(45)

By a slight abuse of notation, we also let $J_p^0(u) =J_p^0((u(c_T)_T)$ for $u\in V({\mathcal {T}})$. Then:

Theorem 2

For all $p\in ( 1,\infty )$ and any $u\in V({\mathcal {T}})$,

$$\begin{aligned} J_p^0(u)= & {} \sup \Big \{ \int _\Omega \sigma (x)\cdot Du(x)\mathrm{d}x \nonumber \\&-\, \frac{1}{p'}\sum _{T\in {\mathcal {T}}} |T| |\sigma (c_T)|^{p'}: \sigma \in RT0_0({\mathcal {T}}) \Big \}, \end{aligned}$$

(46)

where $p'=p/(p-1)$, and for $p=1$,

$$\begin{aligned} J_p^0(u)= & {} \sup \Big \{ \int _\Omega \sigma (x)\cdot Du(x)\mathrm{d}x : \sigma \in RT0_0({\mathcal {T}}) , |\sigma (c_T)|\nonumber \\&\le 1 \Big \}. \end{aligned}$$

(47)

In particular, Theorem 1 corresponds to the particular case $p=1$.

Proof

We consider the case $p>1$. The case $p=1$ is then recovered as a limit problem. The “$\ge$” inequality is obvious. To show the reverse, we assume that u is a minimizer in (45). Then, for any $v\in V^0({\mathcal {T}})$ (cf Appendix A.1), $u+v$ is admissible in the problem and one has $J_p(u+v)\ge J_p(u)$. Hence, taking the derivative $\lim _{t\downarrow 0} (J_p(u+tv)-J_p(v))/t$, it follows that

$$\begin{aligned} \int _\Omega |Du|^{p-2}Du\cdot Dv \mathrm{d}x = 0 \end{aligned}$$

for all $v\in V^0({\mathcal {T}})$. That is, the field $(|Du(T)|^{p-2}Du(T))_{T\in {\mathcal {T}}}\in P0({\mathcal {T}})$ is orthogonal to $GV^0({\mathcal {T}})$, and hence thanks to (43), there exists $\sigma \in RT0_0({\mathcal {T}})$ such that $\sigma (c_T) = |Du(T)|^{p-2}Du(T)$ for all $T\in {\mathcal {T}}$. Clearly, $|\sigma (c_T)|^{p'} = |Du(T)|^p$ and the conclusion easily follows. The case $p=1$ can be recovered as follows: one builds, for $p>1$, a field $\sigma _p\in RT0_0({\mathcal {T}})$ optimal in (46), and letting then $p\rightarrow 1$, one checks that it converges to a field which is optimal in (47). $\square$

Remark A.1

It is quite easy to derive, as a more general result, that given $f:{\mathbb {R}}^d\rightarrow {\mathbb {R}}$ a convex, lower semi-continuous function and any ${\bar{u}}\in P0({\mathcal {T}})$ one has

$$\begin{aligned}&\inf \left\{ \int _\Omega f(\nabla u)\mathrm{d}x: u\in V({\mathcal {T}}), u(c_T)={\bar{u}}_T \ \forall T\in {\mathcal {T}}\right\} \\&\quad = \sup \Big \{ -\int _\Omega {\bar{u}}\hbox {div}\,\sigma \,\mathrm{d}x -\sum _{T\in {\mathcal {T}}}|T| f^*(\sigma (c_T)) : \sigma \in RT0_0({\mathcal {T}}) \Big \}. \end{aligned}$$

Fig. 13

Left: the variant (49), right: the ACR result of Fig. 8

B A Variant with One Node Per Pixel

For imaging application, one drawback of our approach could be the need to introduce more nodes in the representation than the number of pixels.

Given a (gray-level) $n\times m$ image $(u_{i,j})_{i=1,\ldots ,n}^{j=1,\ldots ,m}$ (to simplify, we assume that the scale $h=1$ in this section), even if one rotates the grid by $45^\circ$ and considers the pixels (i, j) as centers of edges of larger squares (for instance, (1, 1), (1, 2), (2, 1), (2, 2) would be the centers of the edges of the square $[(3/2,1/2),(5/2,3/2),(3/2,5/2),(1/2,3/2)]$), one still needs to introduce an additional node in the center of each square (in the above example, at (3/2, 3/2)) and introduce fictitious values $u_{i+1/2,j+1/2}$ (i, j both even or both odd) at these nodes. On average, this increases the dimension of the problems by roughly 50%.

Unfortunately, it seems there is no simple strategy to eliminate this additional node. To illustrate this issue, let us first concentrate on one square. We consider the four vertices $\{0,1\}^2$ as the middle points of the edges of the square (of area 2)

$$\begin{aligned} C= \left[ \left( \tfrac{1}{2},-\tfrac{1}{2}\right) , \left( \tfrac{3}{2},\tfrac{1}{2}\right) , \left( \tfrac{1}{2},\tfrac{3}{2}\right) , \left( -\tfrac{1}{2},\tfrac{1}{2}\right) \right] \end{aligned}$$

and a fifth vertex in (1/2, 1/2) in the middle, which is the middle of both the vertical and horizontal edges cutting the square into two halves. Then, given the values $u_{\alpha ,\beta }:=u(\alpha ,\beta )$, $(\alpha ,\beta )\in \{0,1\}^2$, and c the value at the center, the Crouzeix–Raviart total variation inside the square is

$$\begin{aligned}&\max \Big \{ \sqrt{2}\sqrt{(u_{00}-c)^2+(u_{10}-c)^2}\\&\quad +\,\sqrt{2}\sqrt{(u_{11}-c)^2+(u_{01}-c)^2 }, \\&\qquad \sqrt{2}\sqrt{(u_{00}-c)^2+(u_{01}-c)^2}\\&\quad +\,\sqrt{2}\sqrt{(u_{11}-c)^2+(u_{10}-c)^2 } \Big \}. \end{aligned}$$

(Each gradient norm is multiplied by the area 1 of the corresponding triangle, and we have used that the distance between a vertex of the cube and the center is $1/\sqrt{2}$.) A possibility to eliminate c is to minimize this quantity with respect to c. In the “inpainting” problems of Fig. 8, this would give the same results (since this is precisely what is done automatically by the minimization in this case). Unfortunately, we have, at this point, no idea of how to solve this problem explicitly. It means that to compute the “proximity” operator of the corresponding energy on a whole image, we need to solve subproblems which keep the additional central variable.

Fig. 14

Left: original “clock” image, middle: with a Gaussian blur of SD 1.5 and a $1\%$ Gaussian noise, right: TV-regularized deblurred image with (49)

A simpler possibility is to first optimize with respect to the value c and then, afterward, pick the best orientation. In that case, one needs to solve

$$\begin{aligned} \sqrt{2} \max \Big \{&\min _c \sqrt{(u_{00}-c)^2+(u_{10}-c)^2} \\&+\sqrt{(u_{11}-c)^2+(u_{01}-c)^2}, \\&\quad \min _c \sqrt{(u_{00}-c)^2+(u_{01}-c)^2} \\&\qquad +\sqrt{(u_{11}-c)^2+(u_{10}-c)^2}, \Big \}. \end{aligned}$$

A careful analysis shows that this value is given by the function

$$\begin{aligned}&J^4((u_{00},u_{10},u_{01},u_{11})) \nonumber \\&\quad \,:= \sqrt{2} \max \Big \{ \sqrt{ (u_{11}-u_{00})^2+(u_{10}-u_{01})^2}, \nonumber \\&\qquad \qquad \qquad \quad \qquad \sqrt{(u_{01}-u_{00})^2+(u_{11}-u_{10})^2}, \nonumber \\&\qquad \qquad \qquad \qquad \qquad \sqrt{(u_{10}-u_{00})^2+(u_{11}-u_{01})^2} \Big \}. \end{aligned}$$

(48)

One can use (48) to define, given u defined by its pixel values $(u_{i,j})_{i=1,\ldots ,n}^{j=1,\ldots ,m}$, a discrete total variation as

$$\begin{aligned} J(u)&:= \sum _{(i,j)\hbox { even}} J^4((u_{i,j},u_{i+1,j},u_{i,j+1},u_{i+1,j+1})) \nonumber \\&\quad +\,\sum _{(i,j)\hbox { odd}} J^4((u_{i,j},u_{i+1,j},u_{i,j+1},u_{i+1,j+1})). \end{aligned}$$

(49)

We remark this is a variant of the energy defined in [24] (see also [26] for a theoretical study), which can be optimized by an efficient alternating descent method as soon as one knows how to solve explicitly the subproblems given by the proximity operator of $J^4$, on each square.

Unfortunately, our implementation shows that it does not perform as well as the ACR technique introduced in this paper. Figure 13 compares this to the ACR result in Fig. 8: we obtain a very diffusive solution, with practically no improvement over a non-optimized Crouzeix Raviart implementation.

On the other hand, as is expected, this approximation (which in any case is still based on a hidden, underlying Crouzeix–Raviart discretization), yields to a quite precise approximation of the energy and is a reasonable regularizer for standard inverse problems, cf. Fig. 14.

C The Proximity Operator of (22)

We describe in this section how to implement the proximity operator of the function f in (22). The problem we need to solve is as follows, given ${\bar{\xi }}=({\bar{\xi }}_{mn})_{m=1,\ldots ,4}^{n=1,2}\in {\mathbb {R}}^{4\times 2}$ and $\tau >0$:

$$\begin{aligned} \min _{\xi =(\xi _{mn})_{m=1,\ldots ,4}^{n=1,2}\in {\mathbb {R}}^{4\times 2}} f(\xi ) + \frac{1}{2\tau } \Vert \xi -{\bar{\xi }}\Vert ^2. \end{aligned}$$

(50)

We call $\hbox {prox}\,_{\tau f}({\bar{\xi }})$ the solution of (50). We recall that the prox of the convex conjugate $f^*$ is also easily recovered, once (50) is solved, using Moreau’s identity:

$$\begin{aligned} x = \hbox {prox}\,_{\tau f}(x) + \tau \hbox {prox}\,_{\frac{1}{\tau }f^*}(\tfrac{x}{\tau }). \end{aligned}$$

To solve (50), we first make the following obvious observation: denoting

$$\begin{aligned} x_1 = \sqrt{\xi _{1,1}^2+\xi _{2,1}^2}, \quad&(\xi _{1,1},\xi _{2,1})^T = x_1\eta _1, \\ x_2 = \sqrt{\xi _{3,1}^2+\xi _{4,1}^2}, \quad&(\xi _{3,1},\xi _{4,1})^T = x_2\eta _2, \\ x_3 = \sqrt{\xi _{1,2}^2+\xi _{2,2}^2}, \quad&(\xi _{1,2},\xi _{2,2})^T = x_3\eta _3, \\ x_4 = \sqrt{\xi _{3,2}^2+\xi _{4,2}^2}, \quad&(\xi _{3,2},\xi _{4,2})^T = x_4\eta _4, \end{aligned}$$

(and the same for ${\bar{\xi }}$), it is equivalent to solve:

$$\begin{aligned}&\min _{(x_i)\ge 0,(\eta _i)} \max \{ |x_1|+|x_2|,|x_3|+|x_4|\} \\&\quad +\, \frac{1}{2\tau }\sum _{i=1}^4 |x_i\eta _i - {\bar{x}}_i{\bar{\eta }}_i|^2. \end{aligned}$$

We obtain at the minimum that $\eta _i={\bar{\eta }}_i$, $i=1,\ldots ,4$ and the problem boils down to

$$\begin{aligned} \min _{x=(x_i)_{i=1}^4\in {\mathbb {R}}_+^4} \max \{ |x_1|+|x_2|,|x_3|+|x_4|\} + \frac{1}{2\tau } |x-{\bar{x}}|^2 \end{aligned}$$

where $|x-{\bar{x}}|^2=\sum _{i=1}^4 |x_i-{\bar{x}}_i|^2$. Remark that here, ${\bar{x}}_i\ge 0$ and it is equivalent to look for $x\in {\mathbb {R}}_+^4$ or in ${\mathbb {R}}^4$.

We now explain how to solve this 4-dimensional convex problem. We can rewrite it as

$$\begin{aligned}&\min _{x} \max _{\begin{array}{c}\mu _{12}+\mu _{34}=1\\ \mu _{12}\ge 0,\mu _{34}\ge 0 \end{array}} \mu _{12}(|x_1|+|x_2|) + \mu _{34}(|x_3|+|x_4|) \\&\quad +\, \frac{1}{2\tau } |x-{\bar{x}}|^2 \end{aligned}$$

and then we exchange $\min$ and $\max$. We obtain 4 problems of the form

$$\begin{aligned} \min _{x_1} \mu _{12}|x_1|+\frac{1}{2\tau }|x_1-{\bar{x}}_1|^2. \end{aligned}$$

This is well known to be solved by $x_1=({\bar{x}}_1-\tau \mu _{12})^+$ and with value

$$\begin{aligned} \mu _{12}({\bar{x}}_1-\tau \mu _{12})^+ + \frac{1}{2\tau }{\left\{ \begin{array}{ll} |{\bar{x}}_1|^2 &{}\quad \hbox { if } {\bar{x}}_1\le \tau \mu _{12} \\ |\tau \mu _{12}|^2 &{}\quad \hbox { else.} \end{array}\right. } \end{aligned}$$

When ${\bar{x}}_1\le \tau \mu _{12}$, this is $|{\bar{x}}_1|^2/(2\tau )$, otherwise

$$\begin{aligned} \mu _{12}{\bar{x}}_1 - \frac{\tau }{2}|\mu _{12}|^2 = \frac{1}{2\tau }|{\bar{x}}_1|^2 - \frac{1}{2\tau }|{\bar{x}}_1-\tau \mu _{12}|^2. \end{aligned}$$

We end up with the dual problem

$$\begin{aligned}&\max _{\begin{array}{c}\mu _{12}+\mu _{34}=1\\ \mu _{12}\ge 0,\mu _{34}\ge 0 \end{array}} \frac{1}{2\tau } \Big (\sum _{i=1}^4 |{\bar{x}}_i|^2 - \big (|({\bar{x}}_1-\tau \mu _{12})^+|^2 \\&\quad +\, |({\bar{x}}_2-\tau \mu _{12})^+|^2 +\, |({\bar{x}}_3-\tau \mu _{34})^+|^2 + |({\bar{x}}_4-\tau \mu _{34})^+|^2 \big )\Big ), \end{aligned}$$

whose optimality reads, if $0<\mu _{12}<1$,

$$\begin{aligned}&({\bar{x}}_1-\tau \mu _{12})^+ +({\bar{x}}_2-\tau \mu _{12})^+ \\&\quad = ({\bar{x}}_3-\tau \mu _{34})^+ +({\bar{x}}_4-\tau \mu _{34})^+ \end{aligned}$$

with $\mu _{34}=1-\mu _{12}$.

Without loss of generality, assume that ${\bar{x}}_2\ge {\bar{x}}_1$ and ${\bar{x}}_4\ge {\bar{x}}_3$. We recast the problem as

$$\begin{aligned}&\min _{0\le \mu \le 1} |({\bar{x}}_1-\tau \mu )^+|^2 + |({\bar{x}}_2-\tau \mu )^+|^2 +|({\bar{x}}_3-\tau + \tau \mu )^+|^2 \\&\quad +\, |({\bar{x}}_4-\tau + \tau \mu )^+|^2 \end{aligned}$$

by letting $\mu :=\mu _{12}$ and $\mu _{34}=1-\mu$.

By convexity of the objective, $\mu \in [0,1]$ is optimal if and only if:

$$\begin{aligned} {\left\{ \begin{array}{ll} ({\bar{x}}_1-\tau \mu )^+ + ({\bar{x}}_2-\tau \mu )^+ - ({\bar{x}}_3-\tau +\tau \mu )^+ - ({\bar{x}}_4-\tau +\tau \mu )^+ \le 0 &{} \hbox { if } \mu <1\,; \\ ({\bar{x}}_1-\tau \mu )^+ + ({\bar{x}}_2-\tau \mu )^+ - ({\bar{x}}_3-\tau +\tau \mu )^+ - ({\bar{x}}_4-\tau +\tau \mu )^+ \ge 0 &{} \hbox { if } \mu >0. \end{array}\right. } \end{aligned}$$

(51)

Hence, one sees that if one knows which terms are positive in the above sums, $\mu$ is found by solving the above equations with “$=0$” instead of “$\ge /\le 0$” and then projecting the value onto the interval [0, 1]. For instance, if all values are positive,

$$\begin{aligned} \mu = \left( 0\vee \frac{{\bar{x}}_1+{\bar{x}}_2-{\bar{x}}_3-{\bar{x}}_4+2\tau }{4\tau }\right) \wedge 1. \end{aligned}$$

(52)

Whenever $\mu \in (0,1)$, of course, (51) reads

$$\begin{aligned} ({\bar{x}}_1-\tau \mu )^+&+ ({\bar{x}}_2-\tau \mu )^+ = ({\bar{x}}_3-\tau +\tau \mu )^+ \nonumber \\&+ ({\bar{x}}_4-\tau +\tau \mu )^+ . \end{aligned}$$

(53)

Hence, the problem is solved by exhaustion of the following cases:

• 1. if ${\bar{x}}_2+{\bar{x}}_4\le \tau$, then clearly one can find $\mu \in [0,1]$ such that all terms of the sums in (51) are zero, hence the solution is $x_1=x_2=x_3=x_4=0$.

• 2. if ${\bar{x}}_2+{\bar{x}}_4>\tau$ then:

  • (a) either both${\bar{x}}_2-\tau \mu >0$ and ${\bar{x}}_4-\tau +\tau \mu >0$,

  • (b) or one side of (53) is zero so that one must be in a case of strict inequality in (51), and $\mu \in \{0,1\}$.

  The second case 2b can be first easily eliminated by checking whether $\mu =0$ or $\mu =1$ is a solution of the optimality condition: one has

  $$\begin{aligned} {\bar{x}}_1 + {\bar{x}}_2 \le ({\bar{x}}_3-\tau )^+ + ({\bar{x}}_4-\tau )^+&\ \Leftrightarrow \ \mu =0, \\ ({\bar{x}}_1-\tau )^+ + ({\bar{x}}_2-\tau )^+ \ge {\bar{x}}_3 + {\bar{x}}_4&\ \Leftrightarrow \ \mu =1. \end{aligned}$$

• 3. Otherwise, we must be in the first case 2a, where ${\bar{x}}_2-\tau \mu >0$ and ${\bar{x}}_4-\tau +\tau \mu >0$, equality (53) holds, and which is then split into four possible cases:

  • (a) $\mu$ given by (52), and ${\bar{x}}_1\ge \tau \mu$, ${\bar{x}}_3 \ge \tau (1-\mu )$, then $x_1={\bar{x}}_1-\tau \mu$, $x_2={\bar{x}}_2-\tau \mu$, $x_3={\bar{x}}_3-\tau (1-\mu )$, $x_4={\bar{x}}_4-\tau (1-\mu )$;

  • (b) $\mu =\frac{{\bar{x}}_1+{\bar{x}}_2-{\bar{x}}_4+\tau }{3\tau }$ and ${\bar{x}}_1\ge \tau \mu$, ${\bar{x}}_3\le \tau (1-\mu )$, ${\bar{x}}_3 \ge \tau (1-\mu )$, then $x_1={\bar{x}}_1-\tau \mu$, $x_2={\bar{x}}_2-\tau \mu$, $x_3=0$, $x_4={\bar{x}}_4-\tau (1-\mu )$;

  • (c) $\mu =\frac{{\bar{x}}_2-{\bar{x}}_3-{\bar{x}}_4+2\tau }{3\tau }$ and ${\bar{x}}_1\le \tau \mu$, ${\bar{x}}_2\ge \tau \mu$, ${\bar{x}}_3 \ge \tau (1-\mu )$, then $x_1=0$, $x_2={\bar{x}}_2-\tau \mu$, $x_3={\bar{x}}_3-\tau (1-\mu )$, $x_4={\bar{x}}_4-\tau (1-\mu )$;

  • (d) $\mu =\frac{{\bar{x}}_2-{\bar{x}}_4+\tau }{2\tau }$ if all the previous cases fail to hold, and then $x_1=x_3=0$, $x_2={\bar{x}}_2-\tau \mu$, $x_4={\bar{x}}_4-\tau (1-\mu )$.