Reconstruction of domains with algebraic boundaries from generalized polarization tensors

TL;DR

This paper demonstrates the stability of reconstructing smooth planar domains with algebraic boundaries using a finite set of generalized polarization tensors, extending previous work to less regular domains.

Contribution

It introduces a stable recovery method for algebraic boundary domains from polarization tensors without regularity assumptions, building on prior theoretical identification of boundary-defining polynomials.

Findings

Recovery procedure is stable and effective for domains with algebraic boundaries.

No regularity assumptions are needed for the domain in the reconstruction process.

Performance and limitations of the method are demonstrated through examples.

Abstract

This paper aims at showing the stability of the recovery of a smooth planar domain with a real algebraic boundary from a finite number of its generalized polarization tensors. It is a follow-up of the work [H. Ammari et al., Math. Annalen, 2018], where it is proved that the minimal polynomial with real coefficients vanishing on the boundary can be identified as the generator of a one dimensional kernel of a matrix whose entries are obtained from a finite number of generalized polarization tensors. The recovery procedure is implemented without any assumption on the regularity of the domain to be reconstructed and its performance and limitations are illustrated.

Equations

$${\rm Y}=k\chi (D)+\chi ({\mathbf{R}}^2\setminus \overline{D}),$$

$$\displaystyle\left\{\begin{array}[]{ll}\nabla\cdot\Upsilon\nabla u(x,z)=\delta_{z}(x)&\mbox{in }{\mathbf{R}}^{2},\\ u(x,z)-u_{0}(x,z)=O(|x|^{-1})&\mbox{as }|x|\to\infty,\end{array}\right.$$

$$u(x,z)-u_0(x,z)=\sum _{|\alpha |,|\beta |=1}^{\infty }\frac{(-1{)}^{|\alpha |+|\beta |}}{\alpha !\beta !}{\partial }^{\alpha }{u}_0(x){\mathbf{M}}_{\alpha \beta }{\partial }^{\beta }{u}_0(z)\text{as }|x|\to +\infty ,$$

$$x^{\alpha }=x_1^{{\alpha }_1}{x}_2^{{\alpha }_2},\alpha =({\alpha }_1,{\alpha }_2)\in {\mathbf{N}}^2,\text{ and}|\alpha |={\alpha }_1+{\alpha }_2.$$

$$\partial G\subset V:=\{x\in {\mathbf{R}}^2:g_1(x)=\cdots =g_m(x)=0\}.$$

$$g(x)=\sum _{|\alpha |\le 2}{g}_{\alpha }{x}^{\alpha }$$

$${\mathcal{G}}^{\ast }:=\left\{G\in \mathcal{G}:G=\text{int}\overline{G}\right\}.$$

$$I(X_j)=(P_j),j\in J;$$

$$p(x)=\sum _{|\alpha |\le n}{p}_{\alpha }{x}^{\alpha },$$

$$\mathbf{p}=(p_{\alpha })\in {\mathbf{R}}_{r_d};\sum _{|\beta |\le d}{\mathbf{M}}_{\alpha \beta }(\lambda ,G)p_{\beta }=0\text{for }|\alpha |\le 2d;p_{{\alpha }^{\ast }}=1,{\alpha }^{\ast }=\underset{p_{\alpha }\ne 0}{\max }\alpha .$$

$$({\mathbf{M}}_{\alpha \beta }(G,\lambda ){)}_{|\alpha |\le 2d,0<|\beta |\le d}=({\mathbf{M}}_{\alpha \beta }(\tilde{G},\tilde{\lambda }){)}_{|\alpha |\le 2d,0<|\beta |\le d}\text{iff}G=\tilde{G}\text{mbox}and\lambda =\tilde{\lambda }.$$

$${\mathcal{G}}_0^{\ast }:=\left\{G\in {\mathcal{G}}_0^{\ast }:G\subset B_R(0),I(\partial G)=(g),g(0)=0,\text{deg}(g)=d,\Vert \mathbf{g}\Vert \le M_0,\underset{\partial G}{\min }\Vert \nabla g\Vert \ge \kappa \right\},$$

$$|g|,\Vert \nabla g\Vert ,\Vert H(g)\Vert \le M\text{on}{B}_R(0)$$

$${\mathbf{d}}_{\mathbf{H}}(K_1,K_2)=\max \left\{\underset{x\in K_1}{\sup }\mathbf{d}(x,K_2),\underset{x\in K_2}{\sup }\mathbf{d}(x,K_1)\right\},$$

$$\sum _{|\alpha |\le 2d,0<|\beta |\le d}{\Vert {\mathbf{M}}_{\alpha \beta }({\lambda }^{\ast },G)-{\mathbf{M}}_{\alpha \beta }({\lambda }^{\ast },\tilde{G})\Vert }^2={\epsilon }^2<1,$$

$${\mathbf{d}}_{\mathbf{H}}(\partial G,\partial \tilde{G})\le C{\epsilon }^{\eta }.$$

$${\mathcal{K}}_G^{\ast }[\varphi ](x)=\frac{1}{2\pi }\text{mbox}p.v.{\int }_{\partial G}\frac{⟨x-y,{\nu }_G(x)⟩}{\Vert x-y{\Vert }^2}\varphi (y)d\sigma (y),x\in \partial G,$$

$${\varphi }_{\alpha }(y):=(\lambda I-{\mathcal{K}}_G^{\ast }{)}^{-1}\left[{\nu }_G(x)\cdot \nabla x^{\alpha }\right](y),y\in \partial G.$$

$${\mathbf{M}}_{\alpha \beta }(\lambda ,G):={\int }_{\partial G}{y}^{\beta }{\varphi }_{\alpha }(y)d\sigma (y).$$

$${\mathbb{M}}_{\alpha \beta }(\mu ):={\int }_{\partial G}(I-\mu {\mathcal{K}}_G^{\ast }{)}^{-1}\left[{\nu }_G(x)\cdot \nabla x^{\alpha }\right]y^{\beta }d\sigma (y),$$

$${\tilde{\mathbb{M}}}_{\alpha \beta }(\mu ):={\int }_{\partial \tilde{G}}(I-\mu {\mathcal{K}}_{\tilde{G}}^{\ast }{)}^{-1}\left[{\nu }_{\tilde{G}}(x)\cdot \nabla x^{\alpha }\right]y^{\beta }d\sigma (y).$$

$$F(\mu ):=\sum _{|\alpha |\le 2d,0<|\beta |\le d}{\Vert {\mathbb{M}}_{\alpha \beta }(\mu )-{\tilde{\mathbb{M}}}_{\alpha \beta }(\mu )\Vert }^2.$$

$${{\mathbf{d}}_{\mathbf{H}}}^2(\partial G,\partial \tilde{G})\le C{F}^{1/2}(0).$$

$$\Vert g-\tilde{g}{\Vert }_{L^2(\partial G)}^2+\Vert g-\tilde{g}{\Vert }_{L^2(\partial \tilde{G})}^2\le 2{\kappa }^{-1}{M}^2{F}^{1/2}(0).$$

$${\mathbf{q}}^t\left(\mathbb{M}(0)-\tilde{\mathbb{M}}(0)\right)\mathbf{p}={\int }_{\partial G}{\nu }_G(y)\cdot \nabla q(y)p(y)d\sigma (y)-{\int }_{\partial \tilde{G}}{\nu }_{\tilde{G}}(y)\cdot \nabla q(y)p(y)d\sigma (y),\forall p\in {\mathbf{R}}_d[x],q\in {\mathbf{R}}_{2d}[x],$$

$${\mathbf{q}}^t\left(\mathbb{M}(0)-\tilde{\mathbb{M}}(0)\right)\mathbf{p}={\int }_{\partial G}\frac{\nabla g}{\Vert \nabla g\Vert }\cdot \nabla q(y)p(y)d\sigma (y)-{\int }_{\partial \tilde{G}}\frac{\nabla \tilde{g}}{\Vert \nabla \tilde{g}\Vert }\cdot \nabla q(y)p(y)d\sigma (y),\forall p\in {\mathbf{R}}_d[x],q\in {\mathbf{R}}_{2d}[x].$$

$$(\tilde{\mathbf{g}}\mathbf{g}{)}^t\left(\mathbb{M}(0)-\tilde{\mathbb{M}}(0)\right)(\mathbf{g}+\tilde{\mathbf{g}})={\int }_{\partial G}\Vert \nabla g\Vert {\left(g-\tilde{g}\right)}^{2}d\sigma +{\int }_{\partial \tilde{G}}\Vert \nabla \tilde{g}\Vert {\left(g-\tilde{g}\right)}^{2}d\sigma ,$$

$${\int }_{\partial G}\Vert \nabla g\Vert {\left(g-\tilde{g}\right)}^{2}d\sigma +{\int }_{\partial \tilde{G}}\Vert \nabla \tilde{g}\Vert {\left(g-\tilde{g}\right)}^{2}d\sigma \le 2\left(\Vert \mathbf{g}{\Vert }^2+\Vert \tilde{\mathbf{g}}{\Vert }^2\right)F^{1/2}(0).$$

$${\mathcal{O}}_r:=\left\{y+s{\nu }_G(y);y\in \partial G,s\in (-r,r)\right\}.$$

$$|g(x)|\ge \frac{\kappa }{2}r,\forall x\in \partial {\mathcal{O}}_r.$$

$$g(y\pm r{\nu }_G(y))=\pm \nabla g(y)\cdot {\nu }_G(y)r+\frac{r^2}{2}{\nu }_G^t(y)H(g)(y\pm s_0{\nu }_G(y)){\nu }_G(y),$$

Full text

This work was supported in part by the grant ANR-17-CE40-0029 of the French National Research Agency ANR (project MultiOnde), the LabEx PERSYVAL-Lab (ANR-11-LABX- 0025-01), and the Swiss National Science Foundation grant number 200021–172483.

1. 1 Introduction

2. 2 Real algebraic domains

3. 3 Uniqueness and stability estimates

4. 3.1 Uniqueness

5. 3.2 Stability estimates

6. 4 Algorithm description and numerical examples

7. 4.1 Algorithm

8. 4.2 Examples

1. Introduction

Let $D$ be a bounded connected Lipschitz domain in ${\mathbf{R}}^{2}$, and assume that its boundary $\partial D$ contains the origin. Let $\Upsilon$ be the conductivity distribution in ${\mathbf{R}}^{2}$ given by

[TABLE]

where $\chi$ denotes the indicator function, and $k$ is a fixed constant in $(0,1)\cup(1,+\infty)$. Let $u_{0}(x)=\frac{1}{2\pi}\ln|x|$, be the fundamental solution to the Laplacian in ${\mathbf{R}}^{2}$. For a given position $z$ in ${\mathbf{R}}^{2}$, we consider the following conductivity equation

[TABLE]

where $\delta_{z}$ is the Dirac function at $z$ and $u_{0}(x,z):=u_{0}(x-z)$. The system (1) has a unique solution $u$ which is the total voltage potential generated by the point source placed at $z$ [7]. The function $-\nabla u_{0}(x,z)$ represents the background electric field while $u(x,z)-u_{0}(x,z)$ is the perturbation of the voltage potential due to presence of the inclusion $D$. Then, the far-field perturbation of the voltage potential due to the presence of $D$ is given by [7]

[TABLE]

where and throughout this paper, we use the conventional notation:

[TABLE]

We also use the graded lexicographic order: $\alpha,\,\beta\in\mathbf{N}^{2}$ verifies $\alpha\leq\beta$ if $|\alpha|<|\beta|$, or, if $|\alpha|=|\beta|$, then $\alpha_{1}\leq\beta_{1}$ or $\alpha_{1}=\beta_{1}$ and $\alpha_{2}\leq\beta_{2}.$

The quantities $\mathbf{M}_{\alpha\beta}$ that appear naturally in the multi-polar asymptotic expansion  (2), are called Generalized Polarization Tensors (GPTs). We emphasize that GPTs are not dependent on the positions $x$ and $z$. In fact they only depend on the inclusion $D$ and the conductivity ratio $1/k$ or conductivity contrast $\lambda:=\frac{k+1}{2(k-1)}$. For a fixed contrast $\lambda$, the GPTs are indeed geometric quantities associated with the shape of the domain $D$ such as eigenvalues, capacities, and moments. The notion of GPTs has been used in diverse fields of academic research as well as of engineering applications such as the theories of composites, inverse problems, bio-medical imaging, bio-sensing, nano-sensing, and electro-sensing [5, 9, 11, 12, 13, 10, 15, 26].

From the asymptotic expansion (2), we deduce that the knowledge of all the GPTs is equivalent to knowing the far-field responses of the inclusion for all harmonic excitations. It is well known that in that case the inverse problem of recovering $(\lambda,D)$ has a unique solution [6], and a number of algorithms have been proposed for its numerical treatment [7, 3, 4, 8]. However, in applications, the GPTs are usually only measured with finite accuracy and only a finite number of them can be determined from noisy data. Hence, studying the well-posedness of the inverse problem when only a finite number of GPTs are available is of importance.

The purpose of this paper is to evaluate how much information one can get from the knowledge of a finite number of these GPTs. Precisely, assuming that the domain has an algebraic boundary, we are interested in the inverse problem of recovering its position, its shape and the contrast for given a finite number of its GPTs. Recently the uniqueness to this inverse problem was established by the same authors [1]. Our goal in the present paper is twofold: (i) to quantify the stability of the inversion and (ii) to implement the inversion procedure and apply it to much more general cases than those discussed in [1]. In particular, we show here how to recover the true domain (with possibly nonsmooth boundary) from the recovered polynomial level set even in the case where several candidate domains have the same polynomial level set. In doing so, we resolve key numerical issues which include handling of bifurcation points, segmentation points, and arc sets. It is worth emphasizing that the stability estimates proved in this paper holds for algebraic domains with smooth boundaries. Their generalization to the nonsmooth case is technically quite challenging.

The paper is organized as follows. In Section 2, we introduce the class under consideration of domains with algebraic boundaries. Stability issues are studied in Section 3. The main stability estimates are given in Theorem 3.2. Section 4 is devoted to the presentation of our new numerical algorithm which is designed to recover algebraic domains from finite numbers of their associated GPTs. It is worth mentioning that based on the density with respect to Hausdorff distance of algebraic domains among all bounded domains, the proposed algorithm can be extended via approximation beyond its natural context. This observation has already turned algebraic curves into an efficient tool for describing shapes and reconstructing them from their associated moments [19, 20, 21, 22, 23, 25].

2. Real algebraic domains

In this section, we introduce the class of bounded open subsets in ${\bf R}^{2}$ with real algebraic boundaries. We recall the following definition.

Definition 2.1**.**

An open set $G$ in ${\mathbf{R}}^{2}$ is called real algebraic (or simply algebraic) if there exists a finite number of real coefficient polynomials $g_{i}(x),i=1,\cdots,m$, such that

[TABLE]

The ellipse is a simple example of an algebraic domain, since its general boundary coincides with the zero set of the quadratic polynomial function

[TABLE]

for given real coefficients $(g_{\alpha})_{{|\alpha|\leq 2}}$ and proper signs in the top degree part.

We further denote by $\mathcal{G}$ the collection of bounded algebraic domains. It is well-known that the differential structure of the boundary $\partial G$ consists of algebraic arcs joining finitely many singular points, see for instance [16].

As mentioned in [1], since the connectedness of the respective sets is not accessible by the linear algebra tools we developed for reconstructing an algebraic domain from a finite number of its generalized polarization tensors, we drop such a constraint here. Nevertheless, we call "domains" all elements $G\in\mathcal{G}$.

Following [23] we consider a particular class of algebraic domains which are better adapted to the uniqueness and stability of our inverse shape problem. Let

[TABLE]

An element of $\mathcal{G}^{\ast}$ is called an admissible domain, although it may not be connected.

The assumption that $G=\textrm{int}\,\overline{G}$ implies that $G$ contains no slits or $\partial G$ does not have isolated points. If $G\in\mathcal{G}^{*}$, the algebraic dimension of $\partial G$ is one, and the ideal associated to it is principal. To be more precise, $\partial G$ is contained in a finite union of irreducible algebraic sets $X_{j},\ j\in J,$ of dimension one each. The reduced ideal associated to every $X_{j}$ is principal:

[TABLE]

see, for instance, [14, Theorem 4.5.1]. We assume that each $P_{j}$ is indefinite, i.e., it changes sign when crossing $X_{j}$. Therefore, one can consider the polynomial $g=\prod_{j\in J}P_{j}$, vanishing of the first-order on $\partial G$, that is $|\nabla g|\not=0$ on the regular locus of $\partial G$. According to the real version of Study’s lemma (cf. [16, Theorem 12]) every polynomial vanishing on $\partial G$ is a multiple of $g$, that is $I(\partial G)=(g)$. We define the degree of $\partial G$ as the degree of the generator $g$ of the ideal $I(\partial G)$. For a thorough discussion of the reduced ideal of a real algebraic surface in ${\mathbf{R}}^{d}$, we refer the reader to [17].

Throughout this paper, we denote by $g(x)$ the single polynomial vanishing on $\partial G$ which is the generator of $I(\partial G)$ and satisfying the following normalization condition $g_{\alpha^{*}}=1$, where $\alpha^{*}=\max_{g_{\alpha}\not=0}\alpha$. We further assume that $G\in\mathcal{G}^{*}$.

3. Uniqueness and stability estimates

In this section, we first recall the uniqueness result obtained in [1] and then derive stability estimates for the inversion procedure for smooth algebraic domains.

3.1. Uniqueness

Let ${\mathbf{R}}[x]$ be the ring of polynomials in the variables $x=(x_{1},x_{2})$ and let ${\mathbf{R}}_{n}[x]$ be the vector space of polynomials of degree at most $n$ (whose dimension is $r_{n}=(n+1)(n+2)/2$). Any polynomial function $p(x)\in{\mathbf{R}}_{n}[x]$ has a unique expansion in the canonical basis $x^{\alpha},|\alpha|\leq n$ of ${\mathbf{R}}_{n}[x]$, that is,

[TABLE]

for some vector coefficients ${\bf p}=(p_{\alpha})\in{\mathbf{R}}_{r_{n}}$. The following results are established in [1].

Theorem 3.1**.**

Let $G\in\mathcal{G}^{*}$ with $\partial G$ Lipschitz of degree $d$, and let $g(x)=\sum_{|\alpha|\leq d}g_{\alpha}x^{\alpha}$, be a polynomial function that vanishes of the first-order on $\partial G$, satisfying $I(\partial G)=(g),$ $g_{\alpha^{*}}=1$, and $g(0)=0$, where $\alpha^{*}=\max_{g_{\alpha}\not=0}\alpha$. Then, there exists a discrete set $\Sigma\subset\mathbf{C}_{0}:=\mathbf{C}\setminus[-1/2,1,2]$, such that for any fixed $\lambda\in\mathbf{C}_{0}\setminus\Sigma$, ${\bf g}=(g_{\alpha})\in{\mathbf{R}}_{r_{d}}$ is the unique solution to the following normalized linear system:

[TABLE]

Corollary 3.1**.**

Let $G,\widetilde{G}\in\mathcal{G}^{*}$ be Lipschitz of degree $d$. Let $g$ and $\widetilde{g}$ be two polynomials that vanish respectively of the first order on $\partial G$ and on $\partial\widetilde{G}$ satisfying $I(\partial G)=(g)$   and  $I(\partial\widetilde{G})=(\widetilde{g})$. Assume that $g(0)=\widetilde{g}(0)=0$ and $\|\nabla g\|,\,\|\nabla\widetilde{g}\|>0$ on respectively $\partial G$ and $\partial\widetilde{G}$. Moreover, assume that $G$ is the unique element of $\mathcal{G}^{*}$ containing [math] such that $\partial G\subset\{g=0\}\cup B_{r}(0)$, where $B_{r}(0)$ is the disk of center [math] and radius $r$ large enough. Let $\lambda$ and $\widetilde{\lambda}$ be fixed in $\mathbb{C}_{0}$ such that $\lambda\notin\Sigma,\,\widetilde{\lambda}\notin\widetilde{\Sigma}$, where the sets $\Sigma(\partial G)$ and $\widetilde{\Sigma}=\Sigma(\partial\widetilde{G})$ are as defined in Theorem 3.1. Then, the following uniqueness result holds:

[TABLE]

Proof.

The result is a direct consequence of Theorem 3.1. Since the generalized polarization tensors coincide, and $\lambda\notin\Sigma$, we can deduce from Theorem 3.1 that $g=\widetilde{g}$. The fact that $g(0)=\widetilde{g}(0)=0$ and $\|\nabla g\|,\,\|\nabla\widetilde{g}\|>0$ on respectively $\partial G$ and $\partial\widetilde{G}$ implies that $G=\widetilde{G}$. A straightforward calculation shows then that $\lambda=\widetilde{\lambda}$, which finishes the proof. ∎

3.2. Stability estimates

In this section we derive, under some regularity assumption, stability estimates for the considered inverse problem. For fixed integer $d>0$, and constants $R>0$, $M_{0}>0$, $\kappa>0$, define a reduced set of algebraic domains $\mathcal{G}_{0}^{*}$ by

[TABLE]

where deg denotes the degree. It is not difficult to show that there exists a constant $M>M_{0}$, that only depends on $\mathcal{G}^{*}_{0}$, such that

[TABLE]

for all $g$ satisfying $I(\partial G)=(g)$, where $G\in\mathcal{G}_{0}^{*}$ and $H(g)$ is the Hessian matrix of $g$.

Let $K_{1}$ and $K_{2}$ be two compact sets in ${\mathbf{R}}^{2}$. Recall that the Hausdorff distance between $K_{1}$ and $K_{2}$ is defined by

[TABLE]

where ${\bf d}(x,K_{i})=\inf_{y\in K_{i}}\|x-y\|,\,i=1,2$. Let $\|\;\|$ denote the Euclidean norm of tensors.

Theorem 3.2**.**

Let $G\in\mathcal{G}^{*}_{0},\widetilde{G}\in\mathcal{G}^{*}_{0}$ with respectively $\partial G$ and $\partial\widetilde{G}$. Let $\delta>0$ be a fixed constant and $\lambda_{0}\in{\mathbf{R}}$ satisfying $B_{\delta}(\lambda_{0})\Subset\mathbf{C}\cap\left\{|\lambda|>\frac{1}{2}\right\}$. Then there exists $\lambda^{*}\in(\lambda_{0}-\delta,\lambda_{0}+\delta)$, constants $\eta=\eta(\lambda_{0},\delta,\mathcal{G}^{*}_{0})\in(0,1)$, and $C=C(\lambda_{0},\delta,\mathcal{G}^{*}_{0})>0$, such that if

[TABLE]

then the following stability result holds:

[TABLE]

In order to prove Theorem 3.2, we need to show several intermediate results. Let $g(x)=\sum_{|\alpha|\leq d}g_{\alpha}x^{\alpha}$ and $\widetilde{g}(x)=\sum_{|\alpha|\leq d}\widetilde{g}_{\alpha}x^{\alpha}$ be respectively polynomial functions that vanish respectively of the first-order on $\partial G$ and $\partial\widetilde{G}$ satisfying $I(\partial G)=(g),\;g_{\alpha^{*}}=1$,  $g(0)=0$,  and  $I(\partial\widetilde{G})=(\widetilde{g}),\;\widetilde{g}_{\alpha^{*}}=1,$   $\widetilde{g}(0)=0$.

Further, we shall use standard notation concerning Sobolev spaces. For a density $\phi\in H^{-1/2}(\partial G)$, define the Neumann-Poincaré operator: $\mathcal{K}^{*}_{G}:\,H^{-1/2}(\partial G)\rightarrow H^{-1/2}(\partial G),$ by

[TABLE]

where p.v. denotes the principal value, $\nu_{G}(x)$ is the outward unit normal to $\partial G$ at $x\in\partial G,$ $\langle\,,\,\rangle$ denotes the scalar product in ${\mathbf{R}}^{2}$, and $\|\;\|$ denotes the Euclidean norm in ${\mathbf{R}}^{2}$.

The following lemma characterizes the resolvent set $\rho({\mathcal{K}}^{*}_{G})$ of the operator ${\mathcal{K}}^{*}_{G}$, see, for instance,  [7] and [18].

Lemma 3.1**.**

We have $\mathbb{C}\setminus(-1/2,1/2]\subset\rho({\mathcal{K}}^{*}_{G})$. Moreover, if $|\lambda|\geq 1/2$ , then $(\lambda I-\mathcal{K}^{*}_{G})$ is invertible on $H_{0}^{-1/2}(\partial G):=\{f\in H^{-1/2}(\partial G):\langle f,1\rangle_{-1/2,1/2}=0\}$. Here, $\langle\;,\;\rangle_{-1/2,1/2}$ denotes the duality pairing between $H^{-1/2}(\partial G)$ and $H^{1/2}(\partial G)$.

For $|\lambda|>1/2$ and a multi-index $\alpha=(\alpha_{1},\alpha_{2})\in\mathbf{N}^{2}$, define $\phi_{\alpha}$ by

[TABLE]

The GPTs $\mathbf{M}_{\alpha\beta}$ for $\alpha,\beta\in\mathbf{N}^{2}\;(|\alpha|,|\beta|\geq 1)$, associated with the contrast $\lambda$ and the domain $G$ can be rewritten as [7]

[TABLE]

Denote by $\mathbf{C}_{\star}:=\mathbf{C}\setminus(-\infty,-2]\cup[2,+\infty),$ and let $\mu=\lambda^{-1}\in\mathbb{C}_{\star}$. Define respectively $\mathbb{M}(\mu)$ and $\widetilde{\mathbb{M}}(\mu)$ to be the rectangular matrices with coefficients:

[TABLE]

Note that $\mathbf{M}_{\alpha\beta}(\lambda,G)=\lambda\mathbb{M}_{\alpha\beta}({1}/{\lambda})$ and $\mathbf{M}_{\alpha\beta}(\lambda,\widetilde{G})=\lambda\widetilde{\mathbb{M}}_{\alpha\beta}({1}/{\lambda})$.

Recall the following result from [1].

Lemma 3.2**.**

The functions $\mu\rightarrow\mathbb{M}(\mu),\widetilde{\mathbb{M}}(\mu)\in\mathcal{L}\left({\mathbf{R}}^{r_{d}},{\mathbf{R}}^{r_{2d}}\right)$ are holomorphic matrix-valued on $\mathbf{C}_{\star}$. In addition, $\textrm{ker}(\mathbb{M}(0))=\left\{c{\bf g};\;\;c\in{\mathbf{R}}\right\}$ and $\textrm{ker}(\widetilde{\mathbb{M}}(0))=\left\{c\widetilde{{\bf g}};\;\;c\in{\mathbf{R}}\right\}$.

The proof of Theorem 3.2 has two main steps. In the first step, using the normalized linear system (4), we estimate ${\bf g}-{\bf\widetilde{g}}$ in terms of $\mathbb{M}(0)-\widetilde{\mathbb{M}}(0)$. The second step consists in applying the unique continuation of holomorphic functions on $\mathbb{M}(\mu)-\widetilde{\mathbb{M}}(\mu)$ to "propagate the information" from [math] to $\mu=\lambda^{-1}$.

Let

[TABLE]

We remark that $F(\mu)$ is a real positive function on $\mathbf{C}_{\star}\cap{\mathbf{R}}$. We deduce from Lemma 3.2 that $F(\mu)$ is holomorphic on $\mathbf{C}_{\star}$ and that $F(0)=0$ implies $\mathbf{g}=\widetilde{\mathbf{g}}$. We next estimate how much $\mathbf{g}$ is close to $\widetilde{\mathbf{g}}$ when $F(0)$ is very small.

Proposition 3.1**.**

Let the constants $\kappa$ and $M$ be defined by (6) and (7), respectively. Let $\varepsilon_{0}=\frac{\kappa^{5}}{65M^{4}}$ and $C=64\frac{M^{4}}{\kappa^{5}}$. Assume that $F(0)\leq\varepsilon_{0}^{2}$. Then the following inequality holds:

[TABLE]

In order to prove Proposition 3.1 we need the following three lemmas.

Lemma 3.3**.**

We have

[TABLE]

Proof.

From the definition of the matrices $\mathbb{M}(0)$ and $\widetilde{\mathbb{M}}(0)$, we have

[TABLE]

where the superscript $t$ denotes the transpose.

Since $g$ and $\widetilde{g}$ are in $\mathcal{G}^{*}$ defined by (3) and they respectively generate the ideals associated to $\partial G$ and $\partial\widetilde{G}$, we have $\nu_{G}(x)=\frac{\nabla g(x)}{\|\nabla g(x)\|},\,x\in\partial G$ and $\nu_{\widetilde{G}}(x)=\frac{\nabla\widetilde{g}(x)}{\|\nabla\widetilde{g}(x)\|},\,x\in\partial\widetilde{G}$. Then (14) becomes

[TABLE]

By taking $q(x)=\widetilde{g}(x)g(x)$, $p(x)=g(x)+\widetilde{g}(x)$, and considering the fact that $g(x)$ and $\widetilde{g}(x)$ respectively vanish on $\partial G$ and on $\partial\widetilde{G}$, one finds that

[TABLE]

which in turn implies that

[TABLE]

Hence, (13) holds. ∎

For $r>0$ small, let $\mathscr{O}_{r}\subset{\mathbf{R}}$ being the tubular domain along $\partial G$, defined by

[TABLE]

Lemma 3.4**.**

Assume that $0<r\leq\frac{\kappa}{M}$. Then

[TABLE]

Proof.

Let $x=y+\pm r\nu_{G}(y)\in\mathscr{O}_{r}$, for some $y\in\partial G$ be fixed. From the regularity of $g$, it follows that the function $s\to g(y\pm s\nu_{G}(y))$ is $C^{2}$ and satisfies the following Taylor expansion of order two at zero:

[TABLE]

where $H(g)(y)$ is the Hessian matrix of $g$ at $y$, and $s_{0}$ is some constant in between [math] and $\pm r$. Recalling that $\nu_{G}(x)=\frac{\nabla g(x)}{\|\nabla g(x)\|}$, we therefore obtain that

[TABLE]

which finishes the proof.

∎

The proof of Lemma 3.4 shows that if the zero level set of $g$ is isolated, that is, $\|\nabla g\|\not=0$ on $\partial G$, then the polynomial $g$ behaves as a weighted signed distance function to the boundary $\partial G$ in the small tubular neighborhood domain $\mathscr{O}_{r}$.

Lemma 3.5**.**

Let $r^{*}=\frac{\kappa}{M}$ and $\varepsilon_{0}=\frac{\kappa^{5}}{65M^{4}}$. Assume that $F(0)\leq\varepsilon_{0}^{2}$. Then

[TABLE]

Proof.

Let $\widetilde{x}(t)$ be the parametric representation of the boundary $\partial\widetilde{G}$ ($\partial{\widetilde{G}}=\left\{\widetilde{x}(t),\;t\in{\mathbf{R}}_{+}\right\}$) satisfying

[TABLE]

where $J$ is the counter-clockwise rotation matrix by $\pi/2$. Since $\widetilde{g}$ is smooth, $\widetilde{x}(t)$ is the unique solution to the system (17), which is in addition of class $C^{1}$ and is periodic on ${\mathbf{R}}_{+}$.

Now we shall prove that $\widetilde{x}(t)$ lies indeed in $\mathscr{O}_{r^{*}}$, for all $t\in{\mathbf{R}}_{+}$. Assume that $\partial\widetilde{G}$ is not entirely included in $\mathscr{O}_{r^{*}}$, and define

[TABLE]

Since $0\in\partial G$, $t_{0}>0$ is well defined, is finite, and verifies $\widetilde{x}(t_{0})\in\partial\mathscr{O}_{r^{*}}$. Lemma 3.4 then implies that

[TABLE]

In view of the regularity of $g$ and since $\widetilde{x}$ verifies (17), we have

[TABLE]

Combining inequalities (18) and (19), we obtain that

[TABLE]

for all $t$ satisfying $|t-t_{0}|\leq\frac{\kappa}{4M^{2}}$. Whence

[TABLE]

This together with (13) entail

[TABLE]

which is in contradiction with the fact that $F(0)\leq\varepsilon_{0}^{2}$. Then the inclusion (16) is satisfied.

∎

Proof of Proposition 3.1.

Now, we are ready to prove Proposition 3.1. We further assume that $F(0)\leq\varepsilon_{0}^{2}$. Let $\widetilde{x}(t)$ be defined by (17). Since $\partial\widetilde{G}\subset\mathscr{O}_{r^{*}}$, for each $t>0$, there exists $r(t)\in(0,r^{*})$ and $y(t)\in\partial G$, such that $x(t)=y(t)\pm r(t)\nu_{G}(y(t))$. Noting that $x(t)\in\mathscr{O}_{r(t)}$, we get from Lemma 3.4 the following estimate:

[TABLE]

Following the same arguments as those in the proof of Lemma 3.5, we get

[TABLE]

for all $s$ satisfying $|t-s|\leq\frac{\kappa}{4M^{2}}$. Whence

[TABLE]

for all $t\in{\mathbf{R}}_{+}$. Then

[TABLE]

for all $t\in{\mathbf{R}}_{+}$, which implies

[TABLE]

Repeating the same steps by interchanging $G$ and $\widetilde{G}$, we also get

[TABLE]

Finally, combining inequalities (20), (21), and (13), we obtain the final result of Proposition 3.1. ∎

The second step in proving Theorem 3.2 consists in showing the following proposition.

Proposition 3.2**.**

Let $\delta>0$ be a fixed constant and $\lambda_{0}\in{\mathbf{R}}$ satisfying $B_{\delta}(\lambda_{0})\Subset\mathbf{C}\cap\left\{|\lambda|>\frac{1}{2}\right\}$. Then, there exist constants $\theta=\theta(\lambda_{0},\delta)>0$ and $C=C(\lambda_{0},\delta)>0$ such that

[TABLE]

Proof.

Let $\omega\Subset B_{2}(0)$ be the image of $(\lambda_{0}-\delta,\lambda_{0}+\delta)$ by the complex function $\lambda\mapsto 1/\lambda$. Then there exists a constant $r_{0}\in(0,2)$ such that $\omega\Subset B_{r_{0}}(0)$. Denote by $M_{1}=\|\mu\mapsto F(\mu)\|_{L^{\infty}(B_{r_{0}}(0))}$, and let $w$ be the harmonic measure satisfying

[TABLE]

Since $\mu\mapsto F(\mu)$ is holomorphic on $B_{r_{0}}(0)$, the function $\mu\mapsto\log|F(\mu)|$ is subharmonic, and we can deduce from the Two constants Theorem [24] the following inequality:

[TABLE]

Then by taking $\theta=w(0)$, and $C=M_{1}^{1-w(0)}$, we obtain the result. ∎

Proof of Theorem 3.2.

Finally, we are now in a position to prove Theorem 3.2. Let $\lambda^{*}\in(\lambda_{0}-\delta,\lambda_{0}+\delta)$. By combining estimates (12) and (22) together with the fact that $\mathbf{M}_{\alpha\beta}=\lambda\mathbb{M}_{\alpha\beta}$, we finally obtain the desired stability result stated in Theorem 3.2. ∎

4. Algorithm description and numerical examples

4.1. Algorithm

Before we can dive into the algorithm for recovering algebraic domains from finitely many of their GPTs we must first define a processed form of the GPTs that will form our starting point. In [1, Algorithm 6.2] the GPTs $(\textbf{M}_{\alpha\beta})_{|\alpha|\leq 2d,|\beta|\leq d}$ are flattened out into a linear system. We define one such system explicitly here. For doing so, we use the notation $\textbf{M}_{\alpha\beta}=\textbf{M}_{[\alpha_{1},\alpha_{2}],[\beta_{1},\beta_{2}]}$, where $\alpha=(\alpha_{1},\alpha_{2})$ and $\beta=(\beta_{1},\beta_{2})$.

Definition 4.1**.**

The GPT tessera of order (m,n) is given by

[TABLE]

Definition 4.2**.**

The Tesselated GPT (TGPT) of order (d) is given by

[TABLE]

Our algorithm has in total nine steps. The detail of each step is given algorithmically below with an accompanying description and diagrams. The main steps consist in first recovering the polynomial level set from the given GPTs, then then reconstructing the domain candidates and finally selecting one of the domain candidates in order to minimise the discrepancy between its GPTs and those of the true domain. Our algorithm goes far beyond the stability estimates established in the previous section. Here there is no need to assume that the curve to be recovered is smooth. Nevertheless, in order to reconstruct the domain candidates, several issues need to be carefully resolved. These include bifurcation points, segmentation points, and arc sets.

There are also tuning parameters scattered throughout the various processes and for the most part they are fixed. These tuning parameters should not distract from the otherwise straightforward process.

4.2. Examples

In this section, we apply the algorithm described in the previous subsection to a few examples. We demonstrate its performance by means of a well chosen examples. We also show where the algorithm fails.

In the first example, Figure 4.7 present the possible seven domain candidates corresponding to a disk with a sector missing shown in Figure 4.6. The true domain is recovered by Algorithm 4.9. Here, it corresponds to the one with relative error $0.021$.

In the second example, we consider domains with the same recovered level set. These are discerned from each other by using some boundary information and matching the associated TGPTs. Figures 4.9, 4.10, and 4.11 show three of the six distinct domains. We call these domains "conjoined circles", "crescent" and "intersection of circles" respectively to indicate the shape. All of these shape have the same level set namely two overlapping circles as seen in Figure 4.8. Among the candidates of the conjoined circles the best candidate was found to have relative error $0.01$, see Figure 4.9. Among the candidates of the crescent the best candidate was found to have relative error $0.053$, see Figure 4.10. And among the candidates of the intersection of circles shape the best candidate was found to have relative error $0.044$, see Figures 4.11.

In the third example, we present in Figure 4.12 a square with sinusoidal sides and its recovery from a single domain candidate.

Finally, we show in Figure 4.13 that sometimes the recovered polynomial simply does not give the right domain. The true domain is in blue while the level set of the reconstructed polynomial from the GPTs is in red. This failure to recover the level set could stem from several reasons. The first reason is that higher degree domains are more unstable due to the higher powers taken in computing their GPTs. The second reason is that the proximity of the origin to a bifurcation point could cause instability. This however is still under investigation. We invite the reader to play around with the algorithm which is open source and available at https://github.com/JAndriesJ/ASPT.