Invariance properties of generalized polarization tensors and design of shape descriptors in three dimensions

Accepted Manuscript Invariance properties of generalized polarization tensors and design of shape descriptors in three dimensions Habib Ammari, Daewon Chung, Hyeonbae Kang, Han Wang

PII: DOI: Reference:

S1063-5203(14)00082-7 10.1016/j.acha.2014.05.004 YACHA 982

To appear in:

Applied and Computational Harmonic Analysis

Received date: 16 December 2012 Revised date: 20 March 2014 Accepted date: 30 May 2014

Please cite this article in press as: H. Ammari et al., Invariance properties of generalized polarization tensors and design of shape descriptors in three dimensions, Appl. Comput. Harmon. Anal. (2014), http://dx.doi.org/10.1016/j.acha.2014.05.004

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Invariance Properties of Generalized Polarization Tensors and Design of Shape Descriptors in Three Dimensions∗ Habib Ammari†

Daewon Chung‡

Hyeonbae Kang‡

Han Wang†

Abstract We derive transformation formulas for the generalized polarization tensors under rigid motions and scaling in three dimensions, and use them to construct an infinite number of invariants under those transformations. These invariants can be used as shape descriptors for dictionary matching. AMS subject classifications (2010). 35R30, 35B30 Key words. Generalized polarization tensors, rigid and scaling transformations, invariant, shape descriptor.

1

Introduction

The shape of a domain can be represented in terms of various physical and geometric quantities such as eigenvalues, capacity and moments. The generalized polarization tensor (GPT) is one of them. GPTs are an (infinite) sequence of tensors associated with inclusions (domains) and they appear naturally in the far field expansion of the perturbation of the electrical field in the presence of the inclusion. They are geometric quantities in the sense that the full set of GPTs completely determines the shape of the inclusion. Moreover, one can use a first few terms of GPTs to recover a good approximation of the actual shape of an inclusion [3, 4]. When the domain is transformed by a rigid motion and a scaling, the corresponding GPTs change according to certain rules and it is possible to construct as combinations of GPTs invariants under these transformations. This property makes GPTs suitable for the dictionary matching problem. The dictionary matching problem is to identify the object in the dictionary when the target object is identical to one of the objects in the dictionary up to shifting, rotation, and scaling. The standard method of dictionary matching is to construct invariants, called shape descriptors, under rigid motions and scaling, and to compare those invariants. In [1], invariants are constructed using GPTs in two dimensions. It is the purpose of this letter to extend results of [1] to construct invariants using GPTs in three dimensions. In fact, using contracted GPTs (CGPT), which are harmonic combinations of GPTs, we are able to construct an infinite number of shape descriptors. Using transformation formulas of spherical harmonics under rigid motions [6], we are able to derive transformation formulas obeyed by the CGPTs. We then use these formulas and the method for registration to construct invariants which can be used for target shape description. ∗ This work was supported by the ERC Advanced Grant Project MULTIMOD–267184 and NRF grants No. 2010-0017532. † Department of Mathematics and Applications, Ecole Normale Sup´ erieure, 45 Rue d’Ulm, 75005 Paris, France ([email protected], [email protected]). ‡ Department of Mathematics, Inha University, Incheon 402-751, Korea ([email protected], [email protected]).

1

2

The CGPT block matrix

∗ For D  R3 with Lipschitz boundary, the Neumann-Poincar´e operator KD is defined by  x − y, ν(x) 1 ∗ p.v. [φ](x) = φ(y)dσ(y), x ∈ ∂D . KD 4π |x − y|3 ∂D

Here ·, · denotes the scalar product and ν(x) the unit outward normal vector x ∈ ∂D. Let λ be a real number such that |λ| > 1/2. The generalized polarization tensor (GPT) Mαβ for multi-indices α = (α1 , α2 , α3 ) and β = (β1 , β2 , β3 ) associated with λ and D is defined by  ∗ −1 y β (λI − KD ) [ν · ∇y α ]dσ. Mαβ (λ, D) = ∂D

α

y1α1 y2α2 y3α3 .

Here y = Throughout this letter λ is fixed, so we use Mαβ (D) for Mαβ (λ, D). The contracted GPTs (CGPT) are harmonic combinations of GPTs. Let Ynm , −n ≤ m ≤ n, be the (complex) spherical harmonic of homogeneous degree n and order m. We define the CGPT Mnmkl by      l k ∗ −1 ∂ n m ry Yn (θy , ϕy ) (y)dσ(y) (2.1) ry Yl (θy , ϕy )(λI − KD ) Mnmlk (D) = ∂ν ∂D ∂D with y = ry (cos ϕy sin θy , sin ϕy sin θy , cos θy ). We introduce the matrix Mln of dimension (2l + 1) × (2n + 1) by (Mln )km := Mnmlk , For an integer K, we define

⎡

−l ≤ k ≤ l, −n ≤ m ≤ n.

M11 ⎢ M21 ⎢ M := ⎢ . ⎣ ..

M12 M22 .. .

··· ··· .. .

⎤ M1K M2K ⎥ ⎥ .. ⎥ . . ⎦

MK1

MK2

···

MKK

(2.2)

We call Mln and M the CGPT matrix and the CGPT block matrix of order K, respectively. The following propsoition holds [2]. Proposition 2.1 The CGPTs block matrix M is Hermitian. Furthermore, the matrices Mnn are invertible for all n ≥ 1.

3

Transformation formulas

In this section, we derive transformation formulas of the CGPT matrix Mln defined in (2.2) and the CGPT block matrix under rigid motions and scaling. By a change of variables in (2.1), these formulas play an essential role in finding the invariant under these transformations in the sequel. We first consider the scaling of Mln . The following lemma holds [2]. Lemma 3.1 (Scaling) For any positive integers l, n and scaling parameter s > 0, the following holds: Mln (sD) = sl+n+1 Mln (D) . (3.1) To deal with a shifting of Mln , we need the shifting of the regular spherical harmonics rn Ynm (θ, ϕ). For y = (r, θ, ϕ), z = (rz , θz , ϕz ), and yz = y + z = (r , θ , ϕ ) , we have the following expression of the shifting of the regular spherical harmonic [6], (n,m)

rn Ynm (θ , ϕ ) =



m−μ Cνμnm rzn−ν Yn−ν (θz , ϕz )rν Yνμ (θ, ϕ) ,

(ν,μ)

2

(3.2)

where  Cνμnm =

4π(2n + 1)(n − m)!(n + m)! (2n − 2ν + 1)(2ν + 1)(n − ν − m + μ)!(n − ν + m − μ)!(ν − μ)!(ν + μ)!

1/2 .

Here we use special summation notation: (n,m)



(ν,μ)

=

n μ=min(ν,−ν+n+m)

.

ν=0 μ=max(−ν,ν−n+m)

From (3.2), we obtain (l,k) (n,m)

Mnmlk (Dz ) =



k−j m−μ Cijlk rzl−i Yl−i (θz , ϕz )Mijνμ (D)Cνμnm rzn−ν Yn−ν (θz , ϕz ) .

(3.3)

(i,j) (ν,μ)

For a pair of integers i and l, define (2l + 1) × (2i + 1) matrix Gli = Gli (z) by

k−j Cijlk rzl−i Yl−i (θz , ϕz ) if max(−i, i − l + k) ≤ j ≤ min(i, −i + l + k), (Gli )kj = Gijlk := 0 otherwise (3.4) for −l ≤ k ≤ l. Then, from (3.3) and taking the fact that Mln = 0 if n = 0 or l = 0 into account, we immediately have the following property. Lemma 3.2 (Shifting) For any positive integer l, n and the shifting parameter z, the following holds: n l Gli (z)Miν (D)Gnν (z)t . (3.5) Mln (Dz ) = i=1 ν=1

Rotations in three dimensions may be described in many different ways. Among them we use the Euler angles which can be conveniently used to represent the rotation formula for spherical harmonics. The rotation R is given by ⎤ ⎡ ⎤⎡ ⎤⎡ cos α − sin α 0 cos γ − sin γ 0 cos β 0 − sin β cos γ 0⎦ ⎣ 0 1 0 ⎦ ⎣ sin α cos α 0⎦ . R = ⎣ sin γ 0 0 1 sin β 0 cos β 0 0 1 That is, we rotate by an angle α about z-axis, and by an angle β about the new y axis, and finally by an angle γ about the new z-axis. Let R be a rotation matrix. Since the homogeneous polynomials and Laplace operator are invariant under rotation, Ynm (Rξ) is also spherical harmonic of degree n, and Ynm (Rξ) can be written as follows: Ynm (Rξ) =

n





,m m ρm Yn (ξ) , n

(3.6)

m =−n

where







,m ,m = eim γ dm (β)eimα . ρm n n

Here, 

,m (β) = [(n + m )!(n − m )!(n + m)!(n − m)!]1/2 dm n (−1)m +m+s (cos β )2(n−s)+m−m (sin β )2s−m+m 2 2 , ×  − m + s)!(n − m − s)! (n + m − s)!s!(m s

3

(3.7)

where the sum is over values s such that the factorials are nonnegative: max(0, m − m ) ≤ s ≤ min(n − m , n + m) . We want to find the transformation formula for Mnmlk under rotation, i.e., Mnmlk (DR ). We have t Mnmlk (DR ) = qkl Mln (qm n) ,

where

(3.8)

−n,m , . . . , ρn,m qm n = (ρn n ).

Let for each positive integer n ⎡

ρ−n,−n n ⎢ρ−n,−n+1 ⎢ n Qn = Qn (R) := ⎢ ⎣ ··· ρ−n,n n

ρ−n+1,−n n ρ−n+1,−n+1 n ···

ρ−n+1,n n

··· ··· .. . ···

⎤ ρn,−n n ⎥ ρn,−n+1 n ⎥ ⎥. ··· ⎦ ρn,n n

(3.9)

The rotation matrix Qn is called Wigner D-matrix and known to be unitary [6]. Lemma 3.3 (Rotation) For a unitary matrix R the following relation holds: Mln (DR ) = Ql (R)Mln (D)Q(R)t ,

(3.10)

where t denotes the transpose. In addition to the CGPT block matrix, we define two block matrices: Let K be the truncation order, s be a scaling parameter, z shifting factor, and R rotation, and define ⎡ ⎤ ⎡ ⎤ sQ1 0 ··· 0 G11 0 ··· 0 2 ⎢ G21 G22 · · · ⎢ 0 0 ⎥ s Q2 · · · 0 ⎥ ⎢ ⎥ ⎢ ⎥ G(R) := ⎢ . , Q(s, R) := ⎢ . (3.11) ⎥ . .. ⎥ , . . . . .. .. .. .. .. ⎦ ⎣ .. ⎣ .. . ⎦ GK1 GK2 · · · GKK 0 0 · · · sK QK where Gln = Gln (z) and Qn = Qn (R) are defined by (3.4) and (3.9), respectively. Theorem 3.4 Let Tz be the shift by z, T s scaling by s, and R a unitary matrix. Then the CGPT block matrix M satisfies M(Tz T s R(D)) = sG(z)Q(s, R)M(D)Q(s, R)t G(z)t .

(3.12)

Proof. We have from (3.5) that M(Tz T s R(D)) = GM(T s R(D))Gt . We also have from (3.1) and (3.10) that ⎡ 2 s Q1 M11 (D)Qt1 ⎢ .. M(T s R(D)) = s ⎣ . sK+1 QK MK1 (D)Qt1

··· .. . ···

⎤ s1+K Q1 M1K (D)QtK ⎥ .. ⎦. . s2K QK MKK (D)QtK 

So we obtain (3.12).

4

4

Shape descriptors

In this section, we construct the invariants using CGPT matrices under shifting, scaling, and rotation which can be used as shape descriptors. Let B be a reference domain and D be the one obtained by rotating B by R, scaling by s, −1 and shifting by z, i.e., D = Tz T s R(B). Since Qn is unitary, we have Qn = Qtn . One can easily see that M11 and Q1 commute, i.e., M11 Q1 = Q1 M11 . Since Gnn = I2n+1 ((2n + 1) × (2n + 1) identity matrix), we have M11 (D) = s3 Q1 M11 (B)Qt1 = s3 M11 (B) , M21 (D) = G21 M11 (D) + s4 Q2 M21 (B)Qt1 . Define

UD = M21 (D)M11 (D)−1

Then we have

4.1

and

UB = M21 (B)M11 (B)−1 .

UD = G21 + sQ2 UB Qt1 .

(4.1)

Invariance by registration

We first recall that the matrices Gln and Qn are determined by the shift factor z and rotation R, respectively. Moreover, we may write G21 in terms of rectangular coordinates as ⎡ ⎤ −(z1 + z2 i) 0  0 ⎢ ⎥ z3 − 12 (z1 + z2 i) 0 ⎢ ⎥ ⎢ ⎥    √ ⎢ ⎥ 1 4 1 G21 (z) = 5 ⎢ 6 (z1 − z2 i) (4.2) ⎥. z − (z + z i) 3 1 2 3 6 ⎢ ⎥  ⎢ ⎥ 1 ⎣ ⎦ z3 0 2 (z1 − z2 i) 0 0 z1 − z2 i We present here a method for registration to construct invariants. This method is based on a nontrivial linear mapping u : M5×3 (C) → C3 (M5×3 (C) is the collection of 5 × 3 complex matrices) which satisfies (4.3) u(G21 (z)) = z and u(Q2 (R)UQ1 (R)t ) = Ru(U) for any U ∈ M5×3 (C) and any rotation R. Such a linear mapping does exist and an example is √ √ ⎡ 3 ⎤ 3 U53 − 10 U11 + 10√32 U31 − 103√2 U22 + 103√2 U42 − 10√32 U33 + 10 √ √ 1 ⎢ 3 ⎥ 3 u(U) = √ ⎣i( 10 U11 + 10√32 U31 + 103√2 U22 + 103√2 U42 + 10√32 U33 + 10 U53 )⎦ . √ 5 3 3 3 10 U21 + 5 U32 + 10 U43

(4.4)

(4.5)

One can easily check that this linear transformation satisfies (4.3) using (4.2). We check that it satisfies (4.4) through symbolic computations using Matlab. Since the computation is lengthy, we omit the detail here. In the following we use the shorthand uD = u(UD ). Applying u to (4.1), we have uD = z + sRuB ,

(4.6)

so uD can be used as a registration point. The shape descriptors are constructed after shifting by the registration point. We emphasize that since UD = M21 (D)M11 (D)−1 , uD can be computed using CGPTs of D. 5

The first invariant we introduce is Jln (D) := Mln (T−uD D) = Gl (−uD )Mln (D)Gl (−uD )t ,

l, n = 1, 2, . . . .

Proposition 4.1 For all indices l, n the quantity Jln is invariant by shifting: Jln (Tz D) = Jln (D) for any z ∈ R3 .

(4.7)

Proof. Let Dz = Tz D. By definition of Jln , and using formula (3.5), we have Jln (Tz D) = Mln (T−uDz Dz ) = Mln (Tz−uDz D) = Gl (z − uDz )Mln (D)Gn (z − uDz )t . The relation (4.6) says that uDz = z + uD , and hence we have Jln (Tz D) = Gl (−uD )M(D)Gn (−uD )t = Jln (D). 

This completes the proof. We also have the following lemma.

 = T−u D. For any indices l, n and scaling parameter s > 0 and rotation R, Lemma 4.2 Let D D we have  = Mln (sR(D))  = sl+n+1 Mln (R(D)).  (4.8) Jln (sR(D)) = Jln (sR(D)) In particular, we have and

Jln (sD) = sl+n+1 Jln (D),

(4.9)

Jln (R(D)) = Ql (R)Jln (D)Qn (R)t .

(4.10)

Proof. Since Jln is shifting invariant, we have for any z ∈ R3 : Jln (sR(D)) = Jln (TsRz sR(D)) = Jln (sR(Tz D)). Then by taking z = −uD , we obtain the first identity in (4.8). Moreover, we have from (4.6) uD = −uD + uD = 0.  = sRu  = 0 , which implies the second identity in (4.8). The third one is So we have u(sR(D)) D  (3.1). By taking R = I, we have (4.9). (4.10) follows from the definition of Jln and (3.10). In particular, it can be seen from (4.8) that all Jnn are invertible. So we define the second invariants: Sln (D) := Jnn (D)−1 Jnl (D)Jll (D)−1 Jln (D).

(4.11)

It is worth emphasizing that Sln (D) is a square matrix of dimension (2n + 1). Proposition 4.3 For any indices l, n, the quantity Sln (D) is shifting and scaling invariant: Sln (Tz sD) = Sln (D) for any z ∈ R3 , and s > 0.

(4.12)

Moreover, for any rotation R: Sln (R(D)) = Qn (R)Sln (D)Qn (R)t .

6

(4.13)

Proof. Since Jln ’s are shifting invariant, so are Sln ’s. Scaling invariance of Sln follows from (4.9). Formula (4.13) can be seen using (4.10) as follows: Sln (R(D)) = Jnn (R(D))−1 Jnl (R(D))Jll (R(D))−1 Jln (R(D)) = Qn Jnn (D)−1 Qtn Qn Jnl (D)Qtl Ql Jll (D)−1 Qtl Ql Jln (D)Qtn = Qn Jnn (D)−1 Jnl (D)Jll (D)−1 Jln (D)Qtn = Qn Sln (D)Qtn . This completes the proof.  Since the Frobenius norm of a matrix remains unchanged after the multiplication by a unitary matrix, so finally we define the shape descriptors Iln (D) = Sln (D)F ,

l, n = 1, 2, . . . ,

(4.14)

which is clearly invariant by any rotation, scaling and shifting.

References [1] H. Ammari, T. Boulier, J. Garnier, W. Jing, H. Kang, and H. Wang, Target identification using dictionary matching of generalized polarization tensors, Found. Comp. Math., 14 (2014), 27–62. [2] H. Ammari, J. Garnier, W. Jing, H. Kang, M. Lim, K. Sølna, and H. Wang, Mathematical and Statistical Methods for Multistatic Imaging, Lecture Notes in Mathematics, Vol. 2098, Springer-Verlag, Heidelberg, 2013. [3] H. Ammari, J. Garnier, H. Kang, M. Lim, and S. Yu, Generalized polarization tensors for shape description, Numer. Math., 126 (2014), 199–224. [4] H. Ammari and H. Kang, Polarization and Moment Tensors: with Applications to Inverse Problems and Effective Medium Theory, Vol. 162, Springer-Verlag, New York, 2007. [5] K. Atkinson and W. Han, Spherical Harmonics and Approximations on the Unit Sphere: An Introduction, Lecture Notes in Mathematics, Vol. 2044, Springer-Verlag, Heidelberg, 2012. [6] E.O. Steinborn and K. Ruedenberg, Rotation and translation of regular and irregular solid spherical harmonics, Adv. Quantum Chem., 7 (1973), 1–81.

7