A Splitting Algorithm for Image Segmentation on Manifolds Represented by the Grid Based Particle Method

Abstract

We propose a numerical approach to solve variational problems on manifolds represented by the grid based particle method (GBPM) recently developed in Leung et al. (J. Comput. Phys. 230(7):2540–2561, 2011), Leung and Zhao (J. Comput. Phys. 228:7706–7728, 2009a, J. Comput. Phys. 228:2993–3024, 2009b, Commun. Comput. Phys. 8:758–796, 2010). In particular, we propose a splitting algorithm for image segmentation on manifolds represented by unconnected sampling particles. To develop a fast minimization algorithm, we propose a new splitting method by generalizing the augmented Lagrangian method. To efficiently implement the resulting method, we incorporate with the local polynomial approximations of the manifold in the GBPM. The resulting method is flexible for segmentation on various manifolds including closed or open or even surfaces which are not orientable.

Appendices

Appendix 1. Proof of Theorem 1

We first prove the following lemma

Lemma 1

Suppose $\mathcal J _1(\mathbf{v })$ is continuous and convex on a Hilbert space $\mathbb V$, and

$$\begin{aligned} \mathcal J =\mathcal J _1(\mathbf{v })+<\mathbf{p },\mathbf{g }(\mathbf{v }-\mathbf{b })>+\frac{\eta }{2}<\mathbf{v }-\mathbf{b },\mathbf{g }(\mathbf{v }-\mathbf{b })> \end{aligned}$$

where $\eta >0$ and $\mathbf{g }$ is a positive definite and linear symmetric operator with bounded inverse, then the sequence $\{\mathbf{v }^n\}$ produced by the iteration scheme

$$\begin{aligned} \mathbf{v }^{n+1}=&\underset{\mathbf{v }}{\arg \min } ~~\mathcal J (\mathbf{v },\mathbf{p }^n),\end{aligned}$$

(19)

$$\begin{aligned} \mathbf{p }^{n+1}=&\mathbf{p }^{n}+\tau \mathbf{g }(\mathbf{v }^{n+1}-\mathbf{b }), \end{aligned}$$

(20)

converges, i.e. $\mathbf{v }^n\rightarrow \mathbf{v }^{*}$ when $0<\tau <\frac{2\eta }{\Lambda _{\max }}$, where $\mathbf{v }^{*}$ is the saddle point $(\mathbf{v }^*,\mathbf{p }^*)$ of $\mathcal J (\mathbf{v },\mathbf{p })$. Here $\Lambda _{\max }$ is the largest eigenvalue of $\mathbf{g }$.

Proof

Since $(\mathbf{v }^*,\mathbf{p }^*)$ is a saddle of $\mathcal J$, we have $\mathbf{v }^*=\mathbf{b }$ by $\frac{\partial \mathcal J }{\partial \mathbf{p }}|_{(\mathbf{p }^*,\mathbf{v }^*)}=0$. Let $\partial \mathcal J _1(\mathbf{v })$ be the subgradient of $\mathcal J _1$ at $\mathbf{v }$, i.e. $\partial \mathcal J _1(\mathbf{v })=\{\bar{\mathbf{v }}\in \bar{\mathbb{V }}:\mathcal J _1(\mathbf q )-\mathcal J _1(\mathbf v )\geqslant <\bar{\mathbf{v }},\mathbf q -\mathbf v >, \forall \mathbf q \in \mathbb V \}$, where $\bar{\mathbb{V }}$ is the conjugate space of $\mathbb V$. According to the first order optimization conditions of (19), we have

$$\begin{aligned} \begin{array}{rl} \mathbf d ^{n+1}&:=-\mathbf g \mathbf p ^{n}-\eta \mathbf g (\mathbf v ^{n+1}-\mathbf b )\in \partial \mathcal J _1(\mathbf v ^{n+1}),\\ \mathbf d ^{*}&:=-\mathbf g \mathbf p ^{*}-\eta \mathbf g (\mathbf v ^{*}-\mathbf b )\in \partial \mathcal J _1(\mathbf v ^{*}),\\ \end{array} \end{aligned}$$

thus

$$\begin{aligned} \mathbf d ^{n+1}-\mathbf d ^{*}=-\mathbf g (\mathbf p ^{n}-\mathbf p ^{*})-\eta \mathbf g (\mathbf v ^{n+1}-\mathbf v ^{*}). \end{aligned}$$

Taking the inner product with $\mathbf v ^{n+1}-\mathbf v ^*$ for both sides of the above equation, it becomes

$$\begin{aligned} <\mathbf p ^{n}\!-\!\mathbf p ^{*},\mathbf g (\mathbf v ^{n\!+\!1}\!-\!\mathbf v ^*)>\!=\!-<\mathbf d ^{n\!+\!1}\!-\!\mathbf d ^{*},\mathbf v ^{n\!+\!1}\!-\!\mathbf v ^*>\!-\!\eta <\mathbf v ^{n\!+\!1}\!-\!\mathbf v ^{*},\mathbf g (\mathbf v ^{n\!+\!1}\!-\!\mathbf v ^{*})>.\nonumber \\ \end{aligned}$$

(21)

By the iteration Eq. (20) and the fact that $\mathbf v ^*-\mathbf b =0$, we have

$$\begin{aligned} \mathbf p ^{n+1}-\mathbf p ^{*} =\mathbf p ^{n}-\mathbf p ^{*}+\tau \mathbf g (\mathbf v ^{n+1}-\mathbf v ^{*}). \end{aligned}$$

Taking the norm for the both sides of the above equation, we get

$$\begin{aligned} ||\mathbf p ^{n+1}-\mathbf p ^{*}||^2 =||\mathbf p ^{n}-\mathbf p ^{*}||^2+\tau ^2||\mathbf g (\mathbf v ^{n+1}-\mathbf v ^{*})||^2+2\tau <\mathbf p ^{n}-\mathbf p ^{*},\mathbf g (\mathbf v ^{n+1}-\mathbf v ^{*})>.\nonumber \\ \end{aligned}$$

(22)

Substituting (21) into (22), we have

$$\begin{aligned}&||\mathbf p ^{n+1}-\mathbf p ^{*}||^2 -||\mathbf p ^{n}-\mathbf p ^{*}||^2\nonumber \\&\quad \quad =-2\tau <\mathbf d ^{n+1}-\mathbf d ^{*},\mathbf v ^{n+1}-\mathbf v ^*>-<\mathbf v ^{n+1}-\mathbf v ^{*},\left(-\tau ^2\mathbf g ^2+2\tau \eta \mathbf g \right)(\mathbf v ^{n+1}-\mathbf v ^{*})>.\nonumber \\ \end{aligned}$$

(23)

Since $\mathcal J _1$ is convex, $\mathbf d ^{n+1}\in \partial \mathcal J _1(\mathbf v ^{n+1})$ and $\mathbf d ^{*}\in \partial \mathcal J _1(\mathbf v ^{*})$, thus

$$\begin{aligned} <\mathbf d ^{n+1}-\mathbf d ^{*},\mathbf v ^{n+1}-\mathbf v ^*>\geqslant 0. \end{aligned}$$

(24)

With the condition $0<\tau <\frac{2\eta }{\Lambda _{\max }}$, we conclude that the operator $-\tau ^2\mathbf g ^2+2\tau \eta \mathbf g$ is positive definite.

Now, using both (23) and (24), we have $||\mathbf p ^{n+1}-\mathbf p ^{*}||^2 -||\mathbf p ^{n}-\mathbf p ^{*}||^2<0$, which implies that the sequence $||\mathbf p ^{n}-\mathbf p ^{*}||^2$ is monotonic decreasing with a lower bound $0$ and so it must be convergent.

Finally from (23), we conclude that $\mathbf v ^{n}\rightarrow \mathbf v ^{*}$. $\square$

Now, we can use this Lemma to prove Theorem 1. It is easy to check that $\lambda \int _{\mathcal{M }}\sqrt{\mathbf{v}^\mathrm{T }\mathbf{gv}} \,\mathrm d M$ is convex when $\mathbf g$ is positive definite. For any fixed $u$ and $\mathbf c$, let $\mathbf b =\nabla _s u$ and follow the lemma, one can show that $\{\mathbf{v}^n\}$ produced by Eqs. (8) and (10) converges to the saddle point $(\cdot ,\mathbf v ^*,\mathbf{p}^*,\cdot )$ of $L(\cdot ,\mathbf v ,\mathbf p ,\cdot )$, i.e. $\mathbf v ^n\rightarrow \mathbf v ^*$. Since $(\cdot ,\mathbf v ^*,\mathbf p ^*,\cdot )$ is the saddle point, we have $\mathbf v ^{*}=\nabla _s u$, which completes the proof.

Appendix 2. Derivation of (11)

The directional derivative

$$\begin{aligned} \begin{array}{rl} \frac{\mathrm{d }\tilde{L}(u+\tau w)}{\mathrm{d } \tau }|_{\tau =0}=&\frac{\mathrm{d }}{\mathrm{d } \tau }\left\{ \int _{\mathcal{M }}(f-c_1)^2(u+\tau w) \,\mathrm d M+ \int _{\mathcal{M }}(f-c_2)^2[1-(u+\tau w)] \,\mathrm d M\right.\\&\left.+\frac{\eta }{2}\int _{\mathcal{M }}\left(\mathbf v -\nabla _{\mathbf{s }}(u+\tau w)+\frac{1}{\eta }\mathbf{p }\right)^{\mathrm{T }} \mathbf{g }\left(\mathbf{v }-\nabla _{\mathbf{s }}(u+\tau w)+\frac{1}{\eta }\mathbf{p }\right)\,\mathrm d M\right\} |_{\tau =0}\\ =&\frac{\eta }{2}\frac{\mathrm{d }}{\mathrm{d }\tau }\left\{ {\int _{\smallint }}\sqrt{EG-F^2}\left(\mathbf{v }-\nabla _\mathbf{s }(u+\tau w)+\frac{1}{\eta }\mathbf{p }\right)^{\mathrm{T }}\mathbf{g }\left(\mathbf{v }-\nabla _\mathbf{s }(u+\tau w)+\frac{1}{\eta }\mathbf{p }\right)\,\mathrm d \mathbf{s }\right\} |_{\tau =0}\\&+ \int _{\mathcal{M }}\left[(f-c_1)^2-(f-c_2)^2\right]w\,\mathrm d M\\ =&-\eta {\int _{\smallint }}\sqrt{EG-F^2}<\mathbf{g }\mathbf{v }-\mathbf{g }\nabla _{\mathbf{s }}u+\frac{1}{\eta }\mathbf{g }\mathbf{p },\nabla _\mathbf{s }w>\,\mathrm d \mathbf s + \int _{\mathcal{M }}\left[(f-c_1)^2-(f-c_2)^2\right]\\&\quad w\,\mathrm d M\\ =&\eta {\int _{\smallint }}\nabla _\mathbf{s }\cdot \left(\sqrt{EG-F^2}(\mathbf g \mathbf v -\mathbf g \nabla _\mathbf{s }u+\frac{1}{\eta }\mathbf g \mathbf p )\right)w\,\mathrm d \mathbf s + \int _{\mathcal{M }}\left[(f-c_1)^2-(f-c_2)^2\right]w\,\mathrm d M\\ =&\eta \int _{\mathcal{M }}\frac{1}{\sqrt{EG-F^2}}\nabla _\mathbf{s }\cdot \left(\sqrt{EG-F^2}(\mathbf g \mathbf v -\mathbf g \nabla _\mathbf{s }u+\frac{1}{\eta }\mathbf g \mathbf p )\right)w\,\mathrm d M \\&+\int _{\mathcal{M }}\left[(f-c_1)^2-(f-c_2)^2\right]w\,\mathrm d M.\\ \end{array} \end{aligned}$$

By the variational formulation $\displaystyle \biggl .\frac{\mathrm{d }\tilde{L}(u+\tau w)}{\mathrm{d } \tau }|_{\tau =0}=<\frac{\delta \tilde{L}}{\delta u},w>$, one gets

$$\begin{aligned} \frac{\delta \tilde{L}}{\delta u}&= -\frac{\eta }{\sqrt{EG-F^2}}\nabla _\mathbf{s }\cdot \left(\sqrt{EG-F^2}\mathbf{g }\nabla _{\mathbf{s }}u\right) +\frac{\eta }{\sqrt{EG-F^2}}\nabla _\mathbf{s }\cdot \\&\left(\sqrt{EG-F^2} \mathbf{g }\left(\mathbf{v }+\frac{1}{\eta }\mathbf{p }\right)\right) +(f-c_1)^2-(f-c_2)^2 \, , \end{aligned}$$

which leads to Eq. (11).

Appendix 3. Derivation of (12)

Denote $\mathbf{q }^n=\nabla _\mathbf{s }u^{n+1}-\frac{\mathbf{p }^n}{\eta }$, we have

$$\begin{aligned} \frac{\delta \tilde{L}}{\delta \mathbf{v }}=\frac{\lambda }{\sqrt{\mathbf{v }^ \mathrm{T }\mathbf{g }\mathbf{v }}}\mathbf{g }\mathbf{v }+\eta \mathbf{g }(\mathbf{v }-\mathbf{q }^n). \end{aligned}$$

Solving $\frac{\delta \tilde{L}}{\delta \mathbf{v }}=0$ for $\mathbf{v }$, we get

$$\begin{aligned} \left(\frac{\lambda }{\sqrt{(\mathbf v ^{n+1})^ \mathrm T \mathbf g \mathbf v ^{n+1}}}+\eta \right)\mathbf v ^{n+1}=\eta \mathbf q ^n. \end{aligned}$$

(25)

Taking modulus $||\cdot ||_\mathbf{g }$ for the two sides of (25), it becomes

$$\begin{aligned} ||\mathbf v ^{n+1}||_\mathbf{g }=||\mathbf q ^{n}||_\mathbf{g }-\frac{\lambda }{\eta } \end{aligned}$$

if $||\mathbf q ^{n}||_\mathbf{g }\geqslant \frac{\lambda }{\eta }$. Plugging it back to the above (25), we get

$$\begin{aligned} \mathbf{v }^{n+1}=\frac{\mathbf{q }^n}{||\mathbf{q }^n||_\mathbf{g }}\left(||\mathbf{q }^n||_\mathbf{g }-\frac{\lambda }{\eta }\right). \end{aligned}$$

(26)

On the other hand, if $||\mathbf q ^n||_\mathbf{g }<\frac{\lambda }{\eta }$, then

$$\begin{aligned} \tilde{L}(\mathbf{v },\cdot )&= \biggl .\lambda \int _{\mathcal{M }}\sqrt{\mathbf{v }^{\mathrm{T }}{\mathbf{g }}\mathbf{v }}\,\mathrm d M +\frac{\eta }{2}\int _{\mathcal{M }}(\mathbf{v }-\mathbf{q })^{\mathrm{T }}{\mathbf{g }}(\mathbf{v }-{\mathbf{q }}^n)\,\mathrm d M\\&= \biggl .\lambda \int _{\mathcal{M }} \sqrt{\mathbf{v }^{\mathrm{T }}{\mathbf{g }}{\mathbf{v }}}\,\mathrm d M +\frac{\eta }{2}\int _{\mathcal{M }}\left(|| \mathbf{v } ||_{\mathbf{g }}^2+||\mathbf{q }^n||_{\mathbf{g }}^2\right)\,\mathrm d M -\eta \int _{\mathcal{M }}\mathbf{v }^{\mathrm{T }}\mathbf{g }{\mathbf{q }}^n\,\mathrm d M\\&\geqslant \biggl .\lambda \int _{\mathcal{M }} \sqrt{\mathbf{v }^{\mathrm{T }}{\mathbf{g }}{\mathbf{v }}}\,\mathrm d M +\frac{\eta }{2}\int _{\mathcal{M }}\left(||\mathbf v ||_{\mathbf{g }}^2+||\mathbf q ^n||_{\mathbf{g }}^2\right)\,\mathrm d M -\eta \int _{\mathcal{M }}||\mathbf v ||_{\mathbf{g }}\cdot ||\mathbf q ^n||_{\mathbf{g }}\,\mathrm d M\\&> \biggl .\frac{\eta }{2}\int _{\mathcal{M }}\left(||\mathbf v ||_{\mathbf{g }}^2+||\mathbf{q }^n||_{\mathbf{g }}^2\right)\,\mathrm d M.\\ \end{aligned}$$

This means that $\mathbf v ^{n+1}=0$ must be the minimizer of $\tilde{L}$ with respect to $\mathbf v$. Summarizing these two results, we get the $\mathbf g$-shrinkage operator (12).

Appendix 4. Minimizer of (14)

Let $\mathbf H =\sum _{j=1}^m \left(\mathbf x ^j-\mathbf x ^i\right)\left(\mathbf x ^j-\mathbf x ^i\right)^\mathrm{T }=\mathbf U \begin{pmatrix} \lambda _1&\,&\\&\lambda _2&\\&\,&\lambda _3\\ \end{pmatrix} \mathbf U ^\mathrm{T }$, where $\lambda _1\geqslant \lambda _2\geqslant \lambda _3$ and $\mathbf U =\begin{pmatrix} U_1&U_2&U_3 \\ \end{pmatrix}$ is an orthogonal matrix. Then

$$\begin{aligned} \begin{array}{rl} E=&\sum _{j=1}^m\left(<\mathbf x ^j,\tilde{\mathbf{n }}>-<\mathbf x ^i,\tilde{\mathbf{n }}>\right)^2=\tilde{\mathbf{n }}^\mathrm{T }\mathbf H \tilde{\mathbf{n }} =\left(\mathbf U ^\mathrm{T }\tilde{\mathbf{n }}\right)^\mathrm{T } \begin{pmatrix} \lambda _1&\,&\\&\lambda _2&\\&\,&\lambda _3\\ \end{pmatrix} \left(\mathbf U ^\mathrm{T }\tilde{\mathbf{n }}\right).\\ \end{array} \end{aligned}$$

Since $\left(\mathbf U ^\mathrm{T }\tilde{\mathbf{n }}\right)^\mathrm{T } \left(\mathbf U ^\mathrm{T }\tilde{\mathbf{n }}\right)=1$, thus we have $E\geqslant \lambda _3$. If $\tilde{\mathbf{n }}=U_3$, $E=U_3^\mathrm T \mathbf H U_3=\lambda _3$. Thus we conclude that $\mathbf n =U_3$ is a minimizer of $E$ such that $|\mathbf n |=1$, which completes the proof.

Appendix 5. Explicit formula for the gradient descent update $\Delta u^i$

In this Appendix, we state explicitly the formula to compute the gradient descent update of $u^i$ as a function of the coefficients of the second order polynomial approximation in the GBPM representation. For convenience, let us denote $h_j=v_j^n+\frac{1}{\eta }p_j^n$ and write the coefficients of the approximated second-degree polynomial at $\varvec{x}^i$ for functions $h_1$ and $h_2$ as $\gamma _{\tau _1\tau _2}^i, \delta _{\tau _1\tau _2}^i (0\leqslant \tau _1+\tau _2\leqslant 2,\tau _1,\tau _2\in \mathbb N )$, respectively.

The PDE (11) in the explicit form is given by

$$\begin{aligned} -\frac{\eta }{EG-F^2}\left( S_1 + S_2 S_3 + S_4 S_5 \right) +\frac{\eta }{EG-F^2}\left[ S_6 +S_7 +S_8 + \frac{1}{EG-F^2} S_9 \right]\\ +(f-c_1)^2-(f-c_2)^2&= 0 \, , \end{aligned}$$

where

$$\begin{aligned} S_1&= G\frac{\partial ^2 u}{\partial s_1^2}-2F\frac{\partial ^2 u}{\partial s_1 \partial s_2}+E\frac{\partial ^2 u}{\partial s_2^2} \, , \\ S_2&= <\frac{\partial ^2 \varvec{{x}}}{\partial s_1 \partial s_2}, \frac{\partial \varvec{{x}}}{\partial s_2}>-<\frac{\partial ^2 \varvec{{x}}}{\partial s_2^2},\frac{\partial \varvec{{x}}}{\partial s_1}>-\frac{1}{EG-F^2}(T_1G-T_2F) \, , \\ S_3&= \frac{\partial u}{\partial s_1} \, , \\ S_4&= <\frac{\partial ^2 \varvec{{x}}}{\partial s_1 \partial s_2}, \frac{\partial \varvec{{x}}}{\partial s_1}>-<\frac{\partial ^2 \varvec{{x}}}{\partial s_1^2},\frac{\partial \varvec{{x}}}{\partial s_2}>-\frac{1}{EG-F^2}(T_2E-T_1F) \, , \\ S_5&= \frac{\partial u}{\partial s_2} \, , \\ S_6&= h_1\left(<\frac{\partial ^2 \varvec{x}}{\partial s_1\partial s_2},\frac{\partial \varvec{x}}{\partial s_2}>-<\frac{\partial ^2 \varvec{x}}{\partial s_2^2},\frac{\partial \varvec{x}}{\partial s_1}>\right) \, , \\ S_7&= h_2\left(<\frac{\partial ^2 \varvec{x}}{\partial s_1\partial s_2},\frac{\partial \varvec{x}}{\partial s_1}>-<\frac{\partial ^2 \varvec{x}}{\partial s_1^2},\frac{\partial \varvec{x}}{\partial s_2}>\right) \, , \\ S_8&= G\frac{\partial h_1}{\partial s_1}-F\left(\frac{\partial h_2}{\partial s_1}+\frac{\partial h_1}{\partial s_2}\right) +E\frac{\partial h_2}{\partial s_2} \, , \\ S_9&= \left(Gh_1-Fh_2\right)T_1+\left(-Fh_1+Eh_2\right)T_2 \end{aligned}$$

and

$$\begin{aligned} T_1&= <\frac{\partial ^2 \varvec{x}}{\partial s_1^2},\frac{\partial \varvec{x}}{\partial s_1}>G+<\frac{\partial ^2 \varvec{x}}{\partial s_1\partial s_2},\frac{\partial \varvec{x}}{\partial s_2}>E-\left(<\frac{\partial ^2 \varvec{x}}{\partial s_1^2},\frac{\partial \varvec{x}}{\partial s_2}>+<\frac{\partial ^2 \varvec{x}}{\partial s_1\partial s_2},\frac{\partial \varvec{x}}{\partial s_1}>\right)F \, ,\\ T_2&= <\frac{\partial ^2 \varvec{x}}{\partial s_1\partial s_2},\frac{\partial \varvec{x}}{\partial s_1}>G+<\frac{\partial ^2 \varvec{x}}{\partial s_2^2},\frac{\partial \varvec{x}}{\partial s_2}>E-\left(<\frac{\partial ^2 \varvec{x}}{\partial s_1\partial s_2},\frac{\partial \varvec{x}}{\partial s_2}>+<\frac{\partial ^2 \varvec{x}}{\partial s_2^2},\frac{\partial \varvec{x}}{\partial s_1}>\right)F \, .\\ \end{aligned}$$

Now, replacing all local geometry by the local polynomial least square approximation, we have

$$\begin{aligned} \Delta u^i&= \frac{\eta }{E^iG^i-(F^{i})^2} \left( S_1^i + S_2^i S_3^i + S_4^i S_5^i\right) -\frac{\eta }{E^iG^i-(F^i)^2}\\&\left[ S_6^i + S_7^i + S_8^i + \frac{1}{E^iG^i-(F^i)^2} S_9^i \right] -(f-c_1)^2+(f-c_2)^2\, , \end{aligned}$$

where

$$\begin{aligned} S_1^i&= 2\beta _{20}^iG^i-2\beta _{11}^iF^i +2\beta _{02}^iE^i \, , \\ S_2^i&= \alpha _{11}^i\left(\alpha _{01}^i+\alpha _{11}^i\bar{x}_1^i+2\alpha _{02}^i\bar{x}_2^i\right) -2\alpha _{02}^i\left(\alpha _{10}^i+\alpha _{11}^i\bar{x}_2^i+2\alpha _{20}^i\bar{x}_1^i\right)\\&\quad - \frac{1}{E^iG^i-(F^i)^2}\left(T_1^iG^i-T_2^iF^i\right) \, , \\ S_3^i&= \beta _{10}^i+\beta _{11}^i\bar{x}_2^i+2\beta _{20}^i\bar{x}_1^i \, , \\ S_4^i&= \alpha _{11}^i \left(\alpha _{10}^i+\alpha _{11}^i\bar{x}_2^i+2\alpha _{20}^i\bar{x}_1^i\right) -2\alpha _{20}^i\left(\alpha _{01}^i+\alpha _{11}^i\bar{x}_1^i+2\alpha _{02}^i\bar{x}_2^i\right)\\&\quad - \frac{1}{E^iG^i-(F^i)^2}\left(T_2^iE^i-T_1^iF^i\right) \quad , \\ S_5^i&= \beta _{01}^i+\beta _{11}^i\bar{x}_1^i+2\beta _{02}^i\bar{x}_2^i \quad , \\ S_6^i&= h_1^i\left[\alpha _{11}^i\left(\alpha _{01}^i+\alpha _{11}^i\bar{x}_1^i+2\alpha _{02}^i\bar{x}_2^i\right)-2 \alpha _{02}^i\left(\alpha _{10}^i+\alpha _{11}^i\bar{x}_2^i+2\alpha _{20}^i\bar{x}_1^i\right) \right] \quad , \\ S_7^i&= h_2^i\left[\alpha _{11}^i\left(\alpha _{10}^i+\alpha _{11}^i\bar{x}_2^i+2\alpha _{20}^i\bar{x}_1^i\right)-2 \alpha _{20}^i\left(\alpha _{01}^i+\alpha _{11}^i\bar{x}_1^i+2\alpha _{02}^i\bar{x}_2^i\right) \right] \quad , \\ S_8^i&= G^i\left(\gamma _{10}^i+\gamma _{11}^i\bar{x}_2^i+2\gamma _{20}^i\bar{x}_1^i\right) +E^i\left(\delta _{01}^i+\delta _{11}^i\bar{x}_1^i+2\delta _{02}^i\bar{x}_2^i\right) \\&-F^i\left(\gamma _{01}^i+\gamma _{11}^i\bar{x}_1^i+2\gamma _{02}^i\bar{x}_2^i+\delta _{10}^i+\delta _{11}^i\bar{x}_2^i+2\delta _{20}^i\bar{x}_1^i\right) \quad , \\ S_9^i&= \left(G^ih_1^i-F^ih_2^i\right)T_1^i+\left(-F^ih_1^i+E^ih_2^i\right)T_2^i \end{aligned}$$

and

$$\begin{aligned} T_1^i&= 2\alpha _{20}^i\left(\alpha _{10}^i+\alpha _{11}^i\bar{x}_2^i+2\alpha _{20}^i\bar{x}_1^i\right)G^i+ \alpha _{11}^i\left(\alpha _{01}^i+\alpha _{11}^i\bar{x}_1^i+2\alpha _{02}^i\bar{x}_2^i\right)E^i \\&-\left[2\alpha _{20}^i\left(\alpha _{01}^i+\alpha _{11}^i\bar{x}_1^i+\displaystyle \biggl .2\alpha _{02}^i\bar{x}_2^i\right)+ \alpha _{11}^i\left(\alpha _{10}^i+\alpha _{11}^i\bar{x}_2^i+2\alpha _{20}^i\bar{x}_1^i\right)\right]F^i \, ,\\ T_2^i&= \alpha _{11}^i\left(\alpha _{10}^i+\alpha _{11}^i\bar{x}_2^i+2\alpha _{20}^i\bar{x}_1^i\right)G^i+ 2\alpha _{02}^i\left(\alpha _{01}^i+\alpha _{11}^i\bar{x}_1^i+2\alpha _{02}^i\bar{x}_2^i\right)E^i \\&-\left[2\alpha _{02}^i\left(\alpha _{10}^i+\alpha _{11}^i\bar{x}_2^i+\displaystyle \biggl .2\alpha _{20}^i\bar{x}_1^i\right)+ \alpha _{11}^i\left(\alpha _{01}^i+\alpha _{11}^i\bar{x}_1^i+2\alpha _{02}^i\bar{x}_2^i\right)\right]F^i,\\ h_1^i&= \sum _{\tau _1=0}^2\sum _{0\leqslant \tau _1+\tau _2\leqslant 2}\gamma _{\tau _1\tau _2}^i\left(\bar{x}_1^i\right)^{\tau _1}\left(\bar{x}_2^i\right)^{\tau _2} \, ,\\ h_2^i&= \sum _{\tau _1=0}^2\sum _{0\leqslant \tau _1+\tau _2\leqslant 2}\delta _{\tau _1\tau _2}^i\left(\bar{x}_1^i\right)^{\tau _1}\left(\bar{x}_2^i\right)^{\tau _2}. \end{aligned}$$