Rabah Rabah
email: rabah@emn.fr

Grigory Sklyar
email: sklar@univ.szczecin.pl

Pavel Barkhayev
email: barkhayev@ilt.kharkov.ua

Exact null controllability, complete stabilizability and final observability: the case of neutral type systems *

Keywords: Exact null controllability, complete stabilizability, final observability, neutral type system

For abstract linear systems in Hilbert spaces we revisit the problems of exact controllability and complete stabilizability (stabilizability with an arbitrary decay rate), the latter property being related to exact null controllability. We consider also the case when the feedback is not bounded. We obtain a characterization of complete stabilizability for neutral type systems. Conditions of exact null controllability for neutral type systems are discussed. By duality, we obtain a result about continuous final observability. Illustrative examples are given.

Introduction

Consider the controlled neutral type system ż(t) = A -1 ż(t -1) + Lz t (•) + Bu(t),

where

Lz t (•) = 0 -1
[A 2 (θ) ż(t + θ) + A 3 (θ)z(t + θ)] dθ, with z(t) ∈ R n , u(t) ∈ R m , and the matrices A -1 , A 2 , A 3 and B are of appropriate dimensions. The elements of A 2 and A 3 take values in L 2 (-1, 0). System (1) may be represented in a Hilbert space by the equation ẋ(t) = Ax(t) + Bu(t),

where Bu = (Bu, 0) and A is the infinitesimal generator of a C 0 -semigroup e At given in the product space

M 2 (-1, 0; R n ) def = R n × L 2 (-1, 0; R n ),
1 noted shortly M 2 , and defined by

Ax(t) = Lz t (•) dzt(θ) dθ , x(t) = v(t) z t (•) ,
with the domain D(A) given by

D(A) = (v, ϕ) : ϕ(•) ∈ H 1 , v = ϕ(0) -A -1 ϕ(-1)
.

Our purpose is to analyze exact null controllability of delay systems of neutral type [START_REF] Curtain | Well-posedness, stabilizability, and admissibility for Pritchard-Salamon systems[END_REF], to show the relation with the complete stabilizability (exponential stabilizability with an arbitrary decay rate) of the system and, by duality, to give conditions of the exact final observability of such a system with an output y(t) = Cz(t) or y(t) = Cz(t -1), where y(t) takes values in R p .

The problem of exact controllability for systems of neutral type has been widely investigated. References and important results for system (1) can be found in the work of [START_REF] Rabah | The analysis of exact controllability of neutral-type systems by the moment problem approach[END_REF]. A simplification and improvement of some details of the proofs are given in [START_REF] Rabah | On exact controllability of neutral timedelay systems[END_REF]. The duality with exact (continuous) observability is analyzed in [START_REF] Rabah | Exact observability and controllability for linear neutral type systems[END_REF]. For the stabilizability problem, after the first important works [START_REF] Pandolfi | Stabilization of neutral functional differential equations[END_REF][START_REF] Daniel | On stabilization by state feedback for neutral differential equations[END_REF], there were many results on the stabilizability of delay systems (see, for example, [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF][START_REF] Michiels | Stability and stabilization of time-delay systems. An eigenvalue-based approach[END_REF] and references therein) but neutral type systems have been less investigated [START_REF] Pritchard | The linear quadratic control problem for infinite dimensional systems with unbounded input and output operators[END_REF][START_REF] Salamon | Control and observation of neutral systems[END_REF]. In [START_REF] Hale | Strong stabilization of neutral functional differential equations[END_REF] the main scheme of stabilizing neutral type systems and the robustness, with respect to the delays, of the stabilizing feedback were analyzed. The problem of asymptotic nonexponential stabilizability, which appears only for neutral type systems, was treated in [START_REF] Rabah | On strong regular stabilizability for linear neutral type systems[END_REF][START_REF] Rabah | Stability and stabilizability of mixed retarded-neutral type systems[END_REF], such problem occurs for some systems governed by partial differential equations (see for example [START_REF] Grigory | Controlling a non-homogeneous Timoshenko beam with the aid of the torque[END_REF]). This paper is organized as follows. In Section 2 we give results on the relation between exact null controllability and complete stabilizability for abstract systems in Hilbert spaces. In Section 3, we give necessary conditions of exact null controllability and we characterize complete stabilizability for neutral type systems. Then we formulate a conjecture on the equivalence between exact null controllability and complete stabilizability for neutral type systems. Section 4 is concerned with the dual notion of observability: final continuous observability.

Preliminary results

In this section we consider the abstract system

ẋ = Ax + Bu (2)
where the linear operator A, with domain D(A), is the infinitesimal generator of a C 0 -semigroup e At in the Hilbert space X and B is a linear operator, which may be unbounded but admissible (see, for example, [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]), from the Hilbert space U to X.

Bounded input and feedback

Let us first suppose that the operator B is bounded. The solution of the system (2) with the initial condition x 0 and the control u(t) ∈ L loc 2 (R + ; U ) is given by

x(t) = e At x 0 + t 0 e A(t-τ ) Bu(τ )dτ.
The following notions are well known (see for example [START_REF] Curtain | An introduction to infinite-dimensional linear systems theory[END_REF]).

Definition 2.1 System (2) is said to be exactly controllable at time T if for all x 0 , x 1 ∈ X, there is a control u(t) ∈ L 2 (0, T ; U ) such that the corresponding solution of the system verifies x(T ) = x 1 . The system is said to be exactly null controllable if in the preceding definition x 1 = 0.

There are several results about exact (null) controllability. For example, it is well known that if B is compact, particularly if U is finite dimensional, then there is no exact controllability (first proved in [START_REF] Kuperman | On the question of controllability in infinitedimensional spaces[END_REF], see also [START_REF] Curtain | An introduction to infinite-dimensional linear systems theory[END_REF] and references therein). Another condition of exact controllability, in the case of the bounded operator B, is that for all t ≥ 0, the operator e At is onto (surjective) [START_REF] Louis | On exact controllability in Hilbert spaces[END_REF].

In what follows, we need the following criteria of exact (null) controllability [START_REF] Curtain | An introduction to infinite-dimensional linear systems theory[END_REF].

Theorem 2.2 System (2) is exactly null controllable at time T if and only if

∃δ > 0 : ∀x ∈ X, T 0 B * e A * t x 2 dt ≥ δ 2 e A * T x 2 .
For the condition of exact controllability, the operator e A * T must be replaced by the identity I in the right part of the inequality.

The characterization of exact null controllability is due to a result on range inclusion in Hilbert spaces [START_REF] Douglas | On majorization, factorization, and range inclusion of operators on Hilbert space[END_REF]. We also need some notions of stabilizability. 2) is said to be exponentially stabilizable if there is a linear bounded feedback operator F such that the semigroup e (A+BF )t is exponentially stable: there is a ω > 0 such that e (A+BF )t ≤ M ω e -ωt , M ω ≥ 1.

Definition 2.3 System (
(3

)
The system is said to be completely stabilizable (or stabilizable with an arbitrary decay rate) if for all ω > 0 there is a linear bounded feedback F ω such that (3) holds.

The relation between exact controllability and stabilizability is as follows: exact null controllability implies exponential stabilizability. If e At is a group, complete stabilizability implies exact controllability as shown in [START_REF] Zabczyk | Mathematical control theory: an introduction. Systems & Control: Foundations & Applications[END_REF]Theorem 3.4,p. 229]. Note that the original proof was obtained by [START_REF] Zabczyk | Remarks on the algebraic Riccati equation in Hilbert space[END_REF] which extends the result obtained by [START_REF] Megan | On the stabilizability and controllability of linear dissipative systems in Hilbert space[END_REF]. The same result was proved in [START_REF] Rabah | On Exact Controllability and Complete Stabilizability for Linear Systems in Hilbert Spaces[END_REF][START_REF] Zeng | On exact controllability and complete stabilizability for linear systems[END_REF] for the case of a semigroup e At provided that the operators e At are surjective for all t ≥ 0.

We were tempted to extend this latter result to exact null controllability, possibly under some additional conditions. However the situation is not so simple. We have the following implication, but the converse is not true. Theorem 2.4 If system (2) is exactly null controllable, then it is completely stabilizable by a bounded feedback F.

Proof. Suppose that the system is exactly null controllable at time T . Then

∀x 0 ∈ X, ∃u(•) ∈ L 2 (0, T ; U ) : x(T, x 0 , u(•)) = 0,
where x(t) = x(t, x 0 , u(•)) is the solution with the initial condition x 0 and the control u(t):

x(t, x 0 , u(•)) = e At x 0 + t 0 e A(t-τ ) Bu(τ )dτ.
Then for every x 0 ∈ X, there exists u(

•) ∈ L 2 (0, ∞; U ) such that +∞ 0 x(t) 2 + u(t) 2 dt < ∞.
This means that the system is exponentially stabilizable [START_REF] Zabczyk | Mathematical control theory: an introduction. Systems & Control: Foundations & Applications[END_REF]Th. 3.3,p. 227] :

∃F ω 0 ∈ L(U, X) : e (A+BFω 0 )t ≤ M ω 0 e -ω 0 t , ω 0 > 0.
On the other hand, the exact null controllability of system (2) is equivalent to the exact null controllability of the system ẋ = (A + ωI)x + Bu, ω > 0.

This means that for all ω > 0, for some µ ω > 0, there is

F ω ∈ L(U, X) such that e (A+BFω)t ≤ M µω e -(µω+ω)t ≤ M ω e -ωt .
In order to explain the fact that the converse is not true and that the situation is more complicated, we give examples of two systems without control, where the semigroups are exponentially stable with arbitrary decay rate, but where the states may or may not reach the null state in finite time. These examples can be found in [START_REF] Rhandi | Spectral theory for positive semigroups and applications[END_REF] in the spaces of continuous functions.

Example 1.

In the space L 2 (0, +∞), consider the semigroup

S(t)f (x) = e -t 2 2 -xt f (x + t), t ≥ 0, x ≥ 0.
It is not difficult to see that for this semigroup, for all ω > 0, there is a constant M ω ≥ 1 such that S(t) ≤ M ω e -ωt . We have also σ(S(t)) = {0} and then the spectrum of the infinitesimal generator is empty. On the other hand, there are initial conditions f such that S(t)f = 0 for any t ≥ 0.

Example 2.

In the space L 2 (0, 1), consider the semigroup

S(t)f (x) = f (x + t) 0 ≤ t + x ≤ 1, 0 t + x > 1.
It is not difficult to see that for this semigroup, for all ω > 0, there is a constant M ω ≥ 1 such that S(t) ≤ M ω e -ωt . We have also σ(S(t)) = {0} and then the spectrum of the infinitesimal generator is empty. But, for any initial function f ∈ L 2 (0, 1), we have S(t)f (x) = 0 for t > 2. This means that S(t) = 0, for all t > 2. Then, for any control operator B, the corresponding system is exactly null controllable at time T > 2 with the trivial control u = 0.

Unbounded input and feedback operators

For some control systems, the input operator B may not be bounded and it is very restrictive to assume that the feedback operator F is bounded. For a general theory on systems with unbounded control and observation we refer to the paper [START_REF] Salamon | Infinite dimensional linear systems with unbounded control and observation: a functional analytical approach[END_REF]. For the subclass of interest, which includes linear neutral type systems we refer to [START_REF] Pritchard | The linear quadratic control problem for infinite dimensional systems with unbounded input and output operators[END_REF] and [START_REF] Curtain | Well-posedness, stabilizability, and admissibility for Pritchard-Salamon systems[END_REF][START_REF] Guo | Well-posedness and admissible stabilizability for Pritchard-Salamon systems[END_REF]).

As our final goal is to analyze exact null controllability and complete stabilizability for neutral type systems, we will now consider a wider context of systems with unbounded input and output operators. However, the situation is much more complicated, even if some extension may be considered.

Let X 1 be D(A) endowed with the graph norm noted x 1 and X -1 be the completion of the space X with respect to the resolvent norm

x -1 = (λI -A) -1 x X , λ ∈ ρ(A).
We have the following relation

X 1 ⊂ X ⊂ X -1 ,
with continuous dense injections.

Definition 2.5 Let B be a linear operator, bounded from the Hilbert space U to X. We say that B is an admissible input operator for the semigroup e At if there exists t 1 such that

t 1 0 e A(t 1 -τ ) Bu(τ )dτ ∈ X 1 ,
and for some β > 0

t 1 0 e A(t 1 -τ ) Bu(τ )dτ X 1 ≤ β u L 2 (0,t 1 ) .
Definition 2.6 Assume that operator F is a linear operator, bounded from X 1 to the Hilbert space Y . We say that it is an admissible output operator for the semigroup e At if there exists t 1 > 0 such that for some α > 0

Fe A(t 1 -τ ) x L 2 (0,t 1 ) ≤ α x X , x ∈ X 1 .
Admissibility for some t 1 implies admissibility for all t > 0 (see for example [START_REF] Curtain | Well-posedness, stabilizability, and admissibility for Pritchard-Salamon systems[END_REF]). From the general result on the perturbation of semigroup from the Pritchard-Salamon class, we can deduce the following Cauchy formula for the perturbed semigroup e (A+BF )t , for admissible input and output operators B and F:

e (A+BF )t x = e At x + t 0 e A(t-τ ) BFe (A+BF )τ xdτ, (4)
for all x ∈ X 1 . Moreover e (A+BF )t extends to a C 0 -semigroup on X. This means that Definition 2.3 may be reformulated for an admissible input operator and admissible output feedback.

Theorem 2.7 If the system (2) with an admissible operator B is completely stabilizable by an admissible A-bounded feedbacks then it is completely stabilizable by a bounded linear feedback F.

Proof. In [1, Theorem 5.5] (see also [START_REF] Guo | Well-posedness and admissible stabilizability for Pritchard-Salamon systems[END_REF]), in a more general situation, it is shown that system (2) with an admissible operator B is exponentially stabilizable by admissible feedback (in X 1 and X) if and only if it is exponentially stabilizable by a bounded feedback. Hence, we can suppose without loss of generality, that in (4) the operator F is bounded: F ∈ L(X, U ). This means that complete stabilizability by admissible feedbacks holds if and only if there is complete stabilizability by bounded feedbacks.

From this and Theorem 2.4 we can expect to extend the result of Theorem 2.4 for the case of unbounded control and feedback. But, Theorem 2.4 is based on [START_REF] Zabczyk | Mathematical control theory: an introduction. Systems & Control: Foundations & Applications[END_REF]Th. 3.3,227] which needs [START_REF] Zabczyk | Mathematical control theory: an introduction. Systems & Control: Foundations & Applications[END_REF]Th. 4.3,240] based on the assumption of exact null controllability given by Definition 2.1. For the case of unbounded control and feedback, we refer to [START_REF] Pritchard | The linear quadratic control problem for infinite dimensional systems with unbounded input and output operators[END_REF]Theorem 3.3,page 132]. The condition H4 used in this theorem is guaranteed by exact null controllability in X -1 (each initial state from X -1 may be moved to zero by an L 2 control).

Corollary 2.8 If system (2) with admissible operator B is exactly null controllable in X -1 , then it is completely stabilizable by an admissible feedback and then by a bounded feedback F.

A technical Lemma

In the next section we need the following lemma. Lemma 2.9 Let A be a (n × n)-matrix and B a (n × m)-matrix. The following statements are equivalent.

1. For all λ ∈ C, λ = 0, rank λI -A B = n.

The following equality holds

rank B AB • • • A n-1 B = rank B AB • • • A n-1 B A n ,
and this is equivalent to the inclusion:

Im A n ⊂ Im B AB • • • A n-1 B .
Proof. Conditions 1 and 2 of the lemma may be formulated as follows.

1. If there is x = 0 such that A * x = λx and B * x = 0, then λ = 0.

2. If x = 0 is such that B * A * i x = 0, i ∈ N, then A * n x = 0.
Suppose that 1 holds. Let N be the subspace

N = {x : B * A * i x = 0, i ∈ N}.
It is easy to see that N is A * -invariant and contained in Ker B * . The spectrum of the restriction of A * to N is {0} by Condition 1. This means that A * is nilpotent in N . As the dimension of N is k ≤ n, we obtain A * n x = 0 for all x ∈ N . This gives 2.

Let us show the equivalence of these conditions. Suppose now that 2 holds. Let x = 0 be such that A * x = λx and B * x = 0. This implies that B * A * i x = λ i B * x = 0, for all i ∈ N. From Condition 2, we obtain that 0 = A * n x = λ n x.

As x = 0, this implies λ = 0. This gives that statement 1 is verified.

3 The neutral type system: controllability and stabilizability

In this section we analyze exact null controllability and the complete stabilizability (exponential stabilizability with an arbitrary decay rate) of delay system of neutral type [START_REF] Curtain | Well-posedness, stabilizability, and admissibility for Pritchard-Salamon systems[END_REF] and investigate the relation between the two notions. By duality, we give conditions of the exact final observability of such system with outputs

y(t) = Cz(t) or y(t) = Cz(t -1),
where y(t) takes values in R p . The relation between exact controllability and exponential stabilizability for linear neutral type systems may be found in several papers (see, for example [START_REF] Salamon | Control and observation of neutral systems[END_REF][START_REF] Ito | A linear quadratic optimal control for neutral systems[END_REF][START_REF] Daniel | On stabilization by state feedback for neutral differential equations[END_REF][START_REF] Dusser | On exponential stabilizability of linear neutral systems[END_REF] and references therein).

For the analysis of stabilizability, we need the structure of the spectrum of the state operator A of system (1) and the condition of the growth of semigroup e At . Theorem 3.1 [START_REF] Rabah | Stability analysis of neutral type systems in Hilbert space[END_REF] Let ∆ A be the matrix:

∆ A (λ) = λI -λe -λ A -1 - 0 -1 λe λs A 2 (s) + e λs A 3 (s) ds.
The spectrum of A, noted σ(A), consists of eigenvalues only which are the roots of the equation det ∆ A (λ) = 0. The corresponding eigenvectors of A are of the form

v -e -λ A -1 v e λθ v , v ∈ Ker ∆ A (λ).
The spectrum of A contains a non empty set of point of the form

{ln |µ| + i(arg µ + 2πk) + O(1/k), k ∈ Z},
where µ is a non-zero eigenvalue of the matrix A -1 .

The spectrum is countable and the semigroup e At verifies the spectrum growth assumption (see, for example, [START_REF] Hale | Introduction to functional-differential equations[END_REF]): ∀ω > ω 0 = sup Re σ(A), ∃M ω : e At ≤ M ω e ωt . Definition 3.2 System (1) is exactly null controllable if for some T > 0 and for all

x 0 ∈ M 2 , there is a control u(•) ∈ L 2 (0, T ; R m ) such that e At x 0 + T 0 e A(T -τ ) Bu(τ )dτ = 0,
this corresponds to the concept of complete controllability given first by N. N. Krasovskiǐ for retarded systems.

Let R T be the linear operator defined by

R T u(•) = T 0 e A(T -τ ) Bu(τ )dτ, u(•) ∈ L 2 (0, T ; R m ).
The operator R T is bounded from L 2 (0, T ; R m ) to X. Moreover it takes values in D(A) and is bounded from L 2 (0, T ; R m ) to X 1 (see [START_REF] Ito | A linear quadratic optimal control for neutral systems[END_REF]Corollary 2.7] and [START_REF] Rabah | The analysis of exact controllability of neutral-type systems by the moment problem approach[END_REF] for our system).

The exact null controllability may be formulated by the inclusion

Im e AT ⊂ Im R T ,
where Im e AT and Im R T are images of the operators e AT and R T . From the well-known characterization of range inclusion in Hilbert spaces [START_REF] Douglas | On majorization, factorization, and range inclusion of operators on Hilbert space[END_REF] we can obtain the following proposition, which is an extension of Theorem 2.2.

Proposition 3.3 System (1) is exactly null controllable for some T > 0 if and only if there is a constant δ > 0 such that

T 0 B * e A * (T -τ ) x 2 R m dτ ≥ δ 2 e A * T x 2 M 2 , for all x ∈ M 2 .
We can now give the main result of this Section.

Theorem 3.4 If the system (1) is exactly null controllable, then the following two conditions hold

1. rank ∆ A (λ) B = n for all λ ∈ C, 2. rank µI -A -1 B = n for all µ ∈ C, µ = 0.
Proof. Suppose that system (1) is exactly null controllable. The necessity of condition 1 is trivial. Let us show that condition 2 is verified. We follow a method used in [START_REF] Metel | Criteria for the constructive identifiability and complete controllability of linear time-independent systems of neutral type[END_REF] (see also [START_REF] Khartovskiȋ | Complete controllability and controllability for linear autonomous systems of neutral type[END_REF]). Then, for some T , for all initial conditions, in particular for all ϕ ∈ H 1 (-1, 0; R n ), there is a control u(•) ∈ L 2 (0, T ; R), u(t) = 0 for t > T , such that z(t) = 0, t > T . We can suppose that T > n.

The function z(t) is absolutely continuous and then almost everywhere differentiable. Then we have ż(t) = A -1 ż(t -1) + Lz t + Bu(t).

Replacing ż(t -1) in this equation, we obtain

ż(t) = A -1 (A -1 ż(t -2) + Lz t-1 + Bu(t -1)) + Bu(t).
Without loss of generality one can suppose that the time t is such that the function u is well defined at these points. Repeating this procedure, we obtain

ż(t) = A N -1 ż(t -N )+ N -1 k=0 A k -1 (Lz t-k + Bu(t -k)) .
Putting t = N ≥ T , and using the continuity of z(t) we obtain

0 = A N -1 ( ż(+0) -ż(-0)) + N -1 k=0 A k -1 (Bu(N -k + 0) -Bu(N -k -0)) . (5)
As z(t) for t > 0 is the solution of equation ( 1) we have ż(+0) = A -1 ż(-1) + Lz +0 (•) + Bu(+0).

Then, replacing this expression in [START_REF] Guo | Well-posedness and admissible stabilizability for Pritchard-Salamon systems[END_REF], and putting the initial condition z 0 (θ) = ϕ(θ), we obtain

A N -1 (A -1 φ(-1) + Lϕ(θ) -φ(-0)) + A N -1 Bu(+0)+ N -1 k=0 A k -1 (Bu(N -k + 0) -Bu(N -k -0)) = 0.
As φ(-0) ∈ R n may be chosen arbitrarily, we obtain

Im A N -1 ⊂ Im B A -1 B • • • A N -1 -1 B .
This may be written as

rank B A -1 B • • • A N -1 -1 B = rank B A -1 B • • • A N -1 -1 B A N -1 .
By the Cayley-Hamilton theorem, this gives

rank B A -1 B • • • A n-1 -1 B = rank B A -1 B • • • A n-1 -1 B A n -1 .
Now, using Lemma 2.9, we obtain Condition 2.

The necessary conditions of exact null controllability characterize in fact the complete stabilizability property.

Theorem 3.5 System (1) is completely stabilizable by a feedback law of the form

u(t) = F -1 ż(t -1) + F z t (•), (6)
where

F z t (•) = 0 -1 [F 2 (θ) ż(t + θ) + F 3 (θ)z(t + θ)] dθ if and only if 1. rank ∆ A (λ) B = n for all λ ∈ C, 2. rank µI -A -1 B = n for all µ ∈ C, µ = 0.
Proof. We give a short and direct proof of the necessity even if it may be obtained from Corollary 5.1.3 of [START_REF] Salamon | Control and observation of neutral systems[END_REF]. If Condition 1 is not verified, then there is an eigenvalue λ 0 of the operator A which cannot be modified by the control operator B. This implies the lack of complete stabilizability.

If Condition 2 is not verified, then there is a nonzero eigenvalue µ 0 of the matrix A -1 which cannot be modified. Then the spectral set

{ln |µ 0 | + i(arg µ 0 + 2πk) + O(1/k), k ∈ Z} ⊂ σ(A),
which belongs to a vertical strip, cannot also be modified. This means that complete stabilizability is not possible.

Let us show now that the two conditions are sufficient for complete stabilizability feedback laws of the form [START_REF] Hale | Introduction to functional-differential equations[END_REF].

Suppose that Condition 2 is satisfied. Let us fix an arbitrary ω > 0. As all the non-zero poles of the matrix A -1 are controllable by Condition 2, then a matrix F -1 can be found such that the spectrum σ(A -1 + BF -1 ) verifies ∀µ ∈ σ(A -1 + BF -1 ), µ = 0, ln |µ| < -ω.

Consider now the neutral type system

ż(t) = (A -1 + BF -1 ) ż(t -1) + Lz t + Bu. ( 7
)
Let A 1 be the generator of system [START_REF] Hale | Strong stabilization of neutral functional differential equations[END_REF]. From the structure of the spectrum of neutral type systems like (1), we have only a finite number of eigenvalues λ ∈ σ(A 1 ) such that Re λ ≥ -ω. Now, using Condition 1 of the theorem, a feedback u(t) = F 1 z t (•), where

F 1 z t (•) = 0 -1 [F 2 (θ) ż(t + θ) + F 3 (θ)z(t + θ)] dθ,
can be found (see for example [START_REF] Pandolfi | Stabilization of neutral functional differential equations[END_REF][START_REF] Pritchard | The linear quadratic control problem for infinite dimensional systems with unbounded input and output operators[END_REF][START_REF] Rabah | On strong regular stabilizability for linear neutral type systems[END_REF][START_REF] Rabah | Stability and stabilizability of mixed retarded-neutral type systems[END_REF]) such that all the eigenvalues λ of the system

ż(t) = (A -1 + BF -1 ) ż(t -1) + (L + BF 1 )z t verify Re λ < -ω. If we denote by F the global feedback u(t) = F -1 ż(t -1) + F 1 z t (•)),
then we obtain e (A+BF )t ≤ M e -ωt , M ≥ 1.

Since ω has been arbitrarily taken, this means that the system is completely stabilizable by a feedback of the form (6).

A similar result was obtained for modal controllability (assignment of characteristic quasipolynomial) for neutral system with multiple discrete delays in [START_REF] Metel | Criteria for modal controllability of linear systems of neutral type[END_REF].

In view of Corollary 2.8, Theorem 3.4 and Theorem 3.5, one can formulate the following natural conjecture. Conjecture. System (1) is exactly null controllable if the following two conditions hold

1. rank ∆ A (λ) B = n for all λ ∈ C, 2. rank µI -A -1 B = n for all µ ∈ C, µ = 0.
This means that exact null controllability is equivalent to complete stabilizability for neutral type systems.

It is well known that the Conjecture is verified for some class of neutral type systems with discrete delays [START_REF] Metel | Criteria for the constructive identifiability and complete controllability of linear time-independent systems of neutral type[END_REF][START_REF] Khartovskiȋ | Complete controllability and controllability for linear autonomous systems of neutral type[END_REF] and in the case of retarded systems [START_REF] Pandolfi | Null controllability of a class of functional-differential systems[END_REF]. It seems to us, that one can use the conditions of complete stabilizability to show the result of the conjecture. But at this moment, we have not a satisfactory formal proof.

Final exact observability

The dual notion of exact null controllability in Hilbert space is the notion of final continuous observability. Sometimes the term continuous is replaced (by analogy) by the term exact. In [START_REF] Rabah | Exact observability and controllability for linear neutral type systems[END_REF], the duality between exact controllability and exact observability was analyzed. In the present section we give the result for null exact controllability and the corresponding notion of observability.

We consider the finite dimensional observation

y(t) = Cx(t), (8)
where C is a linear operator and y(t) ∈ R p is a finite dimensional output. There are several ways to design the output operator C [START_REF] Salamon | Neutral functional differential equations and semigroups of operators[END_REF][START_REF] Salamon | Control and observation of neutral systems[END_REF][START_REF] Metel | Criteria for the constructive identifiability and complete controllability of linear time-independent systems of neutral type[END_REF]. One of our goals in this paper is to investigate how to design a minimal output operator like

Cx(t) = Cz(t) or Cx(t) = Cz(t -1), (9)
where C is a p × n matrix. More general outputs, for example with several and/or distributed delays are not considered here. We want to use some results on exact controllability in order to analyze, by duality, the exact observability property in the infinite dimensional setting like, for example, in [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF].

The operator C defined in ( 9) is linear but not bounded in M 2 . However, in both cases it is admissible in the following sense:

T 0 Ce At x 0 2 R p dt ≤ κ 2 x 0 2 M 2 , ∀x 0 ∈ D(A),
because e At x 0 ∈ D(A), t ≥ 0 (see for example [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]).

Definition 4.1 Let K be the output operator

K : M 2 -→ L 2 (0, T ; R p ), x 0 -→ Kx 0 = Ce At x 0 .
System (1) is said to be exactly finally observable or continuously finally observable [START_REF] Salamon | Control and observation of neutral systems[END_REF] if

Kx 0 2 L 2 = T 0 Ce At x 0 2 R p dt ≥ γ 2 e AT x 0 2 M 2 , (10)
for some constant γ > 0, and for all x 0 ∈ D(A). We say that the system is exactly (or continuously) observable if in (4.1) in the second term of the inequality e AT x 0 is replaced by x 0 .

Exact observability means that one continuously determinates the initial state z 0 (•) from the observation on [0, T ], final exact observability that we can continuously determinate the final state z T (•).

The exact (final) observability depends essentially on the topology of the space. We can expect that, the given neutral type system is not exactly observable if we consider x 0 ∈ D(A), with the norm of the graph and no longer in the topology of M 2 . In fact, we obtain the final observability in the initial norm but we need some delay in the observation in the general case.

In order to use the duality between observability and controllability, we need the expression of the adjoint operator K * in the duality with respect to the pivot space M 2 in the embedding

X 1 ⊂ X = M 2 ⊂ X -1 ,
where X 1 = D(A) with the graph norm noted x 1 and X -1 the completion of the space M 2 with respect to the resolvent norm x -1 = (λI -A) -1 x M 2 . The duality relation is

Kx 0 , u(•) L 2 (0,T ;R p ) = x 0 , K * u(•) X 1 ,X d -1
, where X d -1 is constructed as X -1 with A * instead of A (see [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF] for example): X d -1 is the completion of the space M 2 with the resolvent norm corresponding to the operator A * .

Exact null controllability is dual with exact final observability in the corresponding spaces and with the corresponding topologies. It is expected that the operator K * corresponds to a control operator for some adjoint system. However, the situation is not so simple, as it was pointed out in the paper [START_REF] Rabah | Exact observability and controllability for linear neutral type systems[END_REF], from which we take our main considerations on duality. Proposition 4.2 [START_REF] Rabah | On strong regular stabilizability for linear neutral type systems[END_REF][START_REF] Rabah | Exact observability and controllability for linear neutral type systems[END_REF] The adjoint operator A * is given by

A * w ψ(•) = (A * 2 (0)w + ψ(0) - d[ψ(θ)+A * 2 (θ)w] dθ + A * 3 (θ)w
, with the domain D(A * ) consisting of (w, ψ(•)) ∈ M 2 such that:

ψ(θ) + A * 2 (θ)w ∈ H 1 , A * -1 (A * 2 (0)w + ψ(0)) = ψ(-1) + A * 2 (-1)w.
Let x be a solution of the abstract equation

ẋ = A * x, x(t) = w(t) ψ t (θ) . (11)
Then the function w(t) is the solution of the neutral type equation

ẇ(t + 1) = A * -1 ẇ(t)+ 0 -1 [A * 2 (τ ) ẇ(t + 1 + τ ) + A * 3 (τ )w(t + 1 + τ )] dτ. ( 12
)
This means that the form of the adjoint system is not a simple transposition of the initial one [START_REF] Curtain | Well-posedness, stabilizability, and admissibility for Pritchard-Salamon systems[END_REF]. Let us now specify the relation between the solutions of the neutral type equation [START_REF] Megan | On the stabilizability and controllability of linear dissipative systems in Hilbert space[END_REF] related to the adjoint system [START_REF] Louis | On exact controllability in Hilbert spaces[END_REF] and the transposed neutral type equation

ż(t) = A * -1 ż(t -1)+ 0 -1 [A * 2 (τ ) ż(t + τ ), A * 3 (τ )z(t + τ )] dτ, (13)
with initial z 0 (θ). Let A † be the infinitesimal generator of the semigroup corresponding to equation [START_REF] Metel | Criteria for modal controllability of linear systems of neutral type[END_REF].

Let us put w(t) ψ t (θ) = e A * t ξ 0 = e A * t w(0) ψ 0 (θ) , and the conditions of

v(t) z t (θ) = w(t + 1) -A * -1 w(t) w(t + 1 + θ) = e A † t v(0) z 0 (θ) ,
where z 0 (θ) = w(θ + 1) and v(0) = z 0 (0) -A -1 z 0 (-1). We can give the explicit relation between the initial conditions ξ 0 and x 0 :

ξ 0 = w(0) ψ 0 (θ) , x 0 = v(0) z 0 (θ) .
The formal relation between these vectors is

ξ 0 = w(0) ψ 0 (θ) = Φx 0 = Φ w(1) -A -1 w(0) w(θ + 1) ,
and we have the following result.

Theorem 4.3 [START_REF] Rabah | Exact observability and controllability for linear neutral type systems[END_REF] The operator Φ representing the relation between initial conditions x 0 and ξ 0 corresponding to neutral type systems ( 11) -( 12) and ( 13) is linear bounded and bounded invertible from X d 1 to M 2 , where X d 1 is D(A * ) with the graph norm.

Let us now consider the reachability operator of the transposed controlled system:

ẋ(t) = A † x(t) + C † u(t),
where C † = C * 0 . This operator is given by

R † T u(•) = T 0 e A † (T -τ ) C † u(τ )dτ.
The operator K may be written using R † T and the semigroup e A † of system (13) as follows (see [START_REF] Rabah | Exact observability and controllability for linear neutral type systems[END_REF]):

Kx 0 = R † * T Φx 0 if Cx(t) = Cz(t -1), R † * T e A † * Φx 0 if Cx(t) = Cz(t). (14)
We can now formulate the main result of this section.

Theorem 4.4 System (1) with the output y = Cz(t -1) is exactly (continuously) finally observable if and only if system (13) is exactly null controllable. A necessary condition of exact final observability is given by two conditions:

1. Ker ∆ A (λ) C = {0} for all λ ∈ C, 2. Ker λI -A -1 C = {0} for all λ ∈ C, λ = 0.
Proof. According to relation [START_REF] Metel | Criteria for the constructive identifiability and complete controllability of linear time-independent systems of neutral type[END_REF] we have

Kx 0 L 2 = T 0 C * e A † * (T -τ ) Φx 0 2 dτ 1 2
.

As system (13) is exactly null controllable, we obtain

Kx 0 L 2 ≥ δ e A † * T Φx 0 , for all x 0 ∈ D(A). It is easy to see from [21] that e A † * T Φx 0 = Φe A † * (T -τ ) x 0 = e AT x 0 .
This gives Kx 0 L 2 ≥ δ e AT x 0 , which means that exact final observability holds.

For the case of the output y = Cz(t) we cannot say anything if det(A -1 ) = 0. If A -1 is not singular, then e At is a group and exact final observability coincides with exact observability [START_REF] Rabah | Exact observability and controllability for linear neutral type systems[END_REF].

Examples

To illustrate our results and hypothesis we give here 3 examples. The first one shows that for continuous observability a delay in the output is needed if the semigroup is not a group. The second one is taken from [START_REF] Metel | Criteria for the constructive identifiability and complete controllability of linear time-independent systems of neutral type[END_REF] and it is shown that in fact we have exact controllability (not only exact null controllability). The last example illustrates our Conjecture on equivalence between exact controllability and complete stabilizability. All examples are given in the form of a system with one discrete delay: ż(t) = A -1 z(t -1) + A 0 z(t) + A 1 z(t -1) + Bu(t).

Example 3.

System [START_REF] Michiels | Stability and stabilization of time-delay systems. An eigenvalue-based approach[END_REF] with

A 0 = 0 1 0 0 , A 1 = 0, A -1 = 0 1 0 0 , B = 0 1 .
It is easy to see that, for all λ ∈ C,

rank ∆ A (λ) B = λ -λe -λ -1 0 0 λ 1 = 2.
Moreover, for all λ ∈ C, rank λI -A -1 B = n, then the system is exactly controllable (not only to zero). The transposed system ż1 (t) = 0 ż2 (t) = ż1 (t -1) + z 1 (t) is continuously observable with the output y = z 2 (t -1) but not with y(t) = z 2 (t).

Example 4.

The following system was given for exact null controllability and continuous final observability in [START_REF] Metel | Criteria for the constructive identifiability and complete controllability of linear time-independent systems of neutral type[END_REF]:

A 0 = 0, A 1 = 0 1 0 0 , A -1 = 0 -1 0 1 .
In fact, for this system the initial condition is exactly observable by the output y = Cz(t -1), C = 1 0 , and the transposed system is exactly controllable because, for all λ ∈ C,

rank λI -A * -1 C * = rank λ 0 1 1 λ -1 0 = 2.
However, the initial system is not exactly observable by the output y = Cz(t) = z 1 (t), because the initial function z 0 (θ), θ ∈ [0, 1[ cannot be determined.

Example 5.

System [START_REF] Michiels | Stability and stabilization of time-delay systems. An eigenvalue-based approach[END_REF] with

A 0 = 0 0 1 0 , A 1 = 0, A -1 = 1 0 0 0 , B = 1 0 .
We have for all λ ∈ C,

rank ∆ A (λ) B = λ -λe -λ 0 1 -1 λ 0 = 2,
and for all complex λ = 0, rank λI -A -1 B = λ -1 0 1 0 λ 0 = 2.

The system is exactly null controllable by Lemma 2.9 and result in [START_REF] Metel | Criteria for the constructive identifiability and complete controllability of linear time-independent systems of neutral type[END_REF]. It is completely stabilizable by Theorem 6. Consider now the transposed system ż1 (t) = z 2 (t), ż2 (t) = ż2 (t -1).

This system is continuously finally observable by the feedback y = z 1 (t -1) by Theorem 4.4.

Conclusion

We gave some relations between exact null controllability and complete stabilizability of abstract systems in Hilbert spaces. A characterization of complete stabilizability has been given for a large class of linear neutral type systems. Necessary conditions of exact null controllability are given, which conjectured to be also sufficient for neutral type systems, even if they are not in the general case. This also enables the final continuous observability of such systems to be characterized. The following step is to prove the conjecture and to extend such results to the problem of detectability, which is dual with stabilizability.

Acknowledgment

An earlier version of this work was presented at the International Conference "System Analysis: Modelling an Control", organized in the memory of Academician Arkady Kryazhimskiy at the Krasovskii Institute of Mathematics, Ural Branch of RAS, Ekaterinburg, Russia. The first author thanks Professor Vyacheslav Maksimov and all the Organizing Committee of the Conference.


Chafik Samir

Sebastian Kurtek

Anuj Srivastava

Michel Canis

Elastic Shape Analysis of Cylindrical Surfaces for 3D/2D Registration in Endometrial Tissue Characterization

Keywords: Elastic geodesics, 3D-2D registration, elastic deformation, generalized cylinders, re-parametrization, shape analysis, endometriosis, multi-modal registration

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Endometriosis is a complex gynecological medical disease in which endometrial cells (glands and stroma) appear outside their usual locations in the uterine cavity, as illustrated using Fig. I(a). The main symptoms of endometriosis may depend on the site of active endometriosis and are influenced by hormonal changes. This disease affects approximately 10% of women in the reproductive age group and may cause chronic pelvic pain, severe dysmenorrhea, and several digestive problems including diarrhea, occasional constipation, rectal bleeding and infertility. The endometriosis can be found in the pelvic cavity region, specifically in the pelvic peritoneum and pelvic organs. Currently there is no ethologic cure for endometriosis, but it can be treated in a variety of ways, including using pain medication, hormonal treatments, and laparoscopic surgery in severe cases.

As shown in Fig. I(a), the endometriosis tissue frequently penetrates into other tissues and organs. An accurate diagnosis must be made in order to obtain important medical information, such as the depth of infiltration and exact location. Several methods have been used for the purpose of assessing whether a patient has endometriosis. These methods use multiple imaging modalities, such as MRI, transvaginal ultrasound (TVUS), pelvic ultrasound, and laparoscopy to determine the stage of endometriosis in a patient. If MRI and ultrasound techniques fail to detect endometrial lesions, then laparoscopy becomes necessary. Laparoscopy involves minor surgery under general anesthesia and, for many years, laparoscopy was the only method to diagnose endometriosis [START_REF] Chamie | Findings of pelvic endometriosis at Transvaginal US, MR imaging, and Laproscopy[END_REF]. However, there are some limitations in this approach. First, it does not provide any Manuscript received xxx; revised xxx. information about the depth of infiltration. Second, it might not detect the endometrial lesions which are located under other organs. Finally, the place of insertion of the laparoscope into the patient's body depends on the prediction of the location of endometriosis by gynecologists. To avoid these drawbacks, an attractive solution is to combine information collected from TVUS and MR imaging and reach a comprehensive inference on endometrial lesions.

The MRI provides a global map of a nodule's location and gives a good contrast between normal and malignant tissues, as shown in Fig. I(c). However, MR imaging does not include information about the elasticity of endometrial lesions [START_REF] Chamie | Findings of pelvic endometriosis at Transvaginal US, MR imaging, and Laproscopy[END_REF]. On the other hand, TVUS imaging can provide information about the depth of infiltration and, importantly, about the tissue deformation resulting from the transducer's pressure. An example of a TVUS image is shown in Fig. I(b). The resulting nonlinear deformation is highly dependent on the transducer's position, area, shape, and pressure. This scheme of imaging is depicted in Fig. 2 where a surface (shown in the left) is observed using a transducer that changes the surface shape locally (in the middle) and captures a planar image shown in the right. The limitations of TVUS imaging include limited field of view and low signal to noise ratio [START_REF] Kennedy | Transvaginal sonography in reproductive endocrinology and infertility[END_REF]. Due to their different operating principles, TVUS and MRI can provide complementary information for diagnosing endometriosis. In order to counter the individual limitations of MR and TVUS imaging modalities, the registration and fusion of MR/TVUS data must be performed. The main advantage of MR/TVUS fusion is the combination of accurate localization (from MR) and estimation of the depth of infiltration (from TVUS), which can help the surgeon accurately remove pathological tissues. However, registration of this imaging data across modalities is challenging due to differences in the tissue's intensity, dimensionality, patient positioning (translation, rotation, and scaling) and the organ's elasticity (re-parametrization) introduced by the transducer's pressure. Furthermore, registration and fusion of these two images is a challenging task due to the disparity in gray-level intensities, the existence of artifacts including deformation and noise, and the large non-overlapping fields of view. Therefore, shape analysis of structures becomes an important tool in this context. Before we describe our shapebased approach to registration of the two modalities, we summarize related research from past literature.

A. Related Work

The general problem of registering 2D or/and 3D images is especially hard because: (1) often data are derived from different modalities, (2) they represent deformable organs with different shapes in different images, and (3) there is no initial correspondence between shapes contained therein. The problem of 3D-2D registration is a special case of the general image registration problem, with a high interest to clinicians. In a radiological context, the precise knowledge of position and elasticity of some specific organs is very important. With the help of 3D-2D registration methods the surgeon can obtain such knowledge, and benefit in time efficiency and safety during the laparoscopic intervention. In this section, we provide a short summary of some relevant works in 3D-2D registration in the context of diagnosis, treatment, and guided-intervention. A large majority of existing ideas are concerned with registering pre-operative data such as MRI, CT, MR Angiography (MRA), Digital Subtracted Angiography (DSA), etc., to the inter-operative data, i.e. data available during the surgery. Although our problem context is different, some of these methods are relevant as they have the same goal, i.e. register a 3D object to a 2D one. Based on the representations used for registration, one can roughly categorize the most common approaches into feature-based and intensity-based.

Intensity-based methods:

Intensity-based registration methods are quite common in medical image analysis as they can use pixel/voxel values to register single or multiple modalities. Most of the intensity-based methods consist of minimizing a cost function with two terms: (1) a data term that consists of one of the following measures: mutual information, entropy, cross correlation or others, and (2) a regularization term that represents a regularization of the deformation of the first or second order [START_REF] Modersitzki | Regularization[END_REF]. A common approach is to integrate similarity metrics to match image representations. There are two classes of such methods: first, images are represented by their intensity distribution on the whole images [START_REF] Learned-Miller | Joint mri bias removal using entropy minimization across images[END_REF] or on some characteristic subregions [START_REF] Kadir | Saliency, scale and image description[END_REF]. This is the case when a surface is deformed using a parametric model to simulate mammographic compression [START_REF] Hopp | Automatic multimodal 2d/3d breast image registration using biomechanical fem models and intensity-based optimization[END_REF], or when a probabilistic model is defined to estimate a cutting plane as a registration between 3D MRAs and 2D X-ray DSAs [START_REF] Hipwell | Intensity-based 2d-3d registration of cerebral angiograms[END_REF]. The other types of methods represent both images in an artificial modality using morphological tools [START_REF] Maintz | Registration of 3d medical images using simple morphological tools[END_REF], color variances [START_REF] Andronache | Nonrigid registration of multi-modal images using both mutual information and cross-correlation[END_REF], gradient directions [START_REF] Roche | Rigid registration of 3d ultrasound with mr images: A new approach combining intensity and gradient information[END_REF], etc. The main drawback of methods of the first type is that they are highly dependent on the images' intensity correlation. A big difference between intensity values, such as in TVUS and MRI of endometriosis, can make those approaches prone to large errors. Methods of the second type are even more challenging for TVUS images, which vary significantly in their appearance due to multiple echoes and repetition times.

Feature-based methods: Feature-based methods are those that make use of landmarks or other extracted geometric features to register images. Such methods rely on lowdimensional transformations, such as rigid or affine, for registration. For example, Cosse et al. [START_REF] Cosse | Diffeomorphic surface-based registration for mr-us fusion in prostate brachytherapy[END_REF] proposed an ICP-based algorithm to register MR/US prostate and rectum surfaces. Vermandel et al. [START_REF] Vermandel | A 2d/3d matching based on a hybrid approach: improvement to the imaging flow for avm radiosurgery[END_REF] used a hybrid approach by combining both intensity and skeleton-based methods to register 3D MRAs and 2D DSAs. The registration process consists of estimating the six (rigid-body) motion parameters between the 2D image and the projected views of the 3D structure. More recently, Biesdorf et al. [START_REF] Biesdorf | Segmentation and quantification of the aortic arch using joint 3d model-based segmentation and elastic image registration[END_REF] used a combination of 3D modelbased segmentation and elastic image registration to quantify the aortic arch. Their method consists of minimizing an energy functional corresponding to a segmented 3D cylindrical tubular model as a first term and a 2D registered circular cross-section as a second term. However, they used a parametric spline model to ensure smoothness of the resulting displacement field, which makes the method lose its generality in function spaces. Although this paper has a similar setting to the current paper, the focus is somewhat different. In that paper, one assumed the knowledge of projection from 3D to 2D and used this information to segment a 3D volume. In our case, the cutting plane for TVUS is unknown and the projection is not available. Furthermore, we have an unknown deformation between the 3D and 2D objects. We assume the availability of segmentation to estimate the cutting plane, the registration, and the deformation between objects. Some of the current approaches consider registration using locally-affine transformations. Here, one splits the deformation vector field into small regions. For example, Zikic et al. [START_REF] Zikic | Linear intensity-based image registration by markov random fields and discrete optimization[END_REF] used Markov random fields to estimate affine transformations that align 3D to 2D data by comparing the 2D projected views of the volume and the available 2D images. Also, this way the deformations can be computed locally and then interpolated. Buerger et al. [START_REF] Buerger | Hierarchical adaptive local affine registration for fast and robust respiratory motion estimation[END_REF] exploited the same idea by using an adaptive region-based method to estimate wholethorax deformations during a respiratory motion. They split the 3D volume into hierarchical structures including adaptive sub-blocks. A local affine transformation is then estimated on each block and a multilevel B-spline interpolation is used to estimate a full displacement (a deformation field on the whole volume). To cope with the limitations of intensitybased approaches, Zhang et al. [START_REF] Zhang | Supir: Surface uncertaintypenalized, non-rigid image registration for pelvic ct imaging[END_REF] combine intensity and uncertainty weighted surfaces to register 3D CBCT and 3D FBCT prostate images. Their cost function includes three terms: the similarity between images, the similarity between spherical surfaces including uncertainty, and the regularization term.

Among all intensity or hybrid based techniques, the methods relying on mutual information have been shown to robustly align images from different imaging modalities. However, matching two image distributions assumes that they are sufficiently correlated which is not the case when images are produced from very distinct modalities, resulting in different data dimensions.

B. Our Approach

Since the main goal of this paper is to study endometriosis, we focus on this anatomical structure directly. The purpose of this work is to characterize endometrial tissues in MR images and, since the nodules can be deformed in the TVUS images, there is no automatic correspondence between features in both modalities. Thus, we propose a new method based on shape analysis, since a shape feature is invariant to translation, scale, rotation, and re-parametrization. We assume that we can extract, manually or automatically, the boundaries of endometrial tissue in MRI and TVUS images. In case of 3D MRI images, we obtain a 2D surface as the boundary, while for TVUS images we get a simple closed curve. Although the boundary of an endometrium forms a closed surface, the MRI image is collected at such a low resolution that in the absence of the closing slices on both ends, the resulting surfaces are topologically cylindrical. So, our modified goal is now to focus on the shapes of these two objects -a cylindrically-parameterized surface from MRI and a planar closed curve from TVUS -and perform registration of these extracted structures. When comparing two objects in different dimensions, as is the case here, there are two possibilities. One is to project the higher-dimensional object to the lower space and then compare with the other object. However, the parameters of projection (or the plane of intersection) are unknown and can make a big difference in the final result. The other possibility is to extrapolate the lower-dimensional object to the higher space and perform a more exhaustive comparison. Although this approach can be computationally expensive, it is more stable with respect to the choice of parameters. So, we will take this approach by extending the TVUS curve into a generalized cylinder and comparing the two shapes. In this setup, we pose the following questions:

1) Which closed curve on the MRI surface optimally corresponds to the TVUS curve? 2) What is the optimal deformation between the two curves and how to quantify the elasticity of endometriosis using a measure of this deformation? 3) Can we extend this deformation from the optimal curve to the full MRI surface to help a physician in locating and assessing endometrium shape variability? To answer these and related questions, we propose a surfacebased elastic registration and deformation method for this task. This approach is a modification of shape analysis of certain parametric surfaces developed by Kurtek et al. [START_REF] Kurtek | Elastic geodesic paths in shape space of parametrized surfaces[END_REF], [START_REF] Kurtek | Parameterization-invariant shape comparisons of anatomical surfaces[END_REF] in previous papers. While previous papers dealt with spherical and quadrilateral surfaces, the current paper will extend that mathematical framework to cylindrically-parameterized surfaces. The main strength of this framework is that it solves the problems of registration and optimal deformation in a joint framework.

We summarize the goal of this paper in precise mathematical form. As stated above, the variability between TVUS and MRI images results from two transformations. Firstly, there is a change in patient position due to the fact that MR imaging is performed in supine position, whilst TVUS imaging is performed in lithotomy position. As shown in Fig. 2, the first transformation, denoted by T , introduces a change in translation and rotation of the 3D MRI surface f MRI (depicted as a cylinder). Secondly, there is a nonlinear deformation Φ due to the ultrasound transducer's pressure, resulting in the deformed surface Φ(T (f MRI )). Thus, a planar curve on the MRI surface may be transformed to a non-planar one on Φ(T (f MRI )), and vice versa. Let β US be a planar curve observed using the intersection of a cutting plane P with Φ(T (f MRI )), i.e. β US = P ∩ Φ(T (f MRI )). The problem of interest is: Given β US and f MRI , estimate a closed curve βMRI ∈ f MRI that optimally corresponds to β US under elastic deformations. Note that T , Φ, and P are unknown nuisance variables and one needs to estimate them, at least partially, in the estimation process.

The rest of the paper is organized as follows. We describe a conceptual overview of our framework in Section 2 and the mathematical framework for registration of cylindricallyparametrized surfaces in Section 3. We use that framework to compute optimal deformations (geodesics) between surfaces in Section 4. The application of this framework to analyze MRI surfaces and TVUS curves is presented in Section 5, with some validation analysis in Section 6. We close this paper with a brief conclusion in Section 7.

II. MATHEMATICAL FRAMEWORK

Let β US : S 1 → R 2 represent a closed, planar curve coming from the TVUS image. Also, let f MRI : S 1 × [0, 1] → R 3 be the surface obtained from the MRI image. We will develop a novel Riemannian framework for shape analysis of cylindrically-parametrized surfaces. In particular, the developed tools will give us a principled approach for registra-tion and deformation of surfaces. Furthermore, the resulting distance (that satisfies symmetry and triangle inequality) will be based on an elastic metric, which explicitly captures the amount of "stretching" and "bending" needed to deform one surface into another. In order to use this framework, we form an ultrasound surface using f US :

S 1 ×[0, 1] → R 3 by replicat- ing β US in the axial direction, i.e. f US (u, v) = β US (u), ∀ v ∈ [0, 1].
This surface represents all of the information provided by β US coming from the ultrasound image. Our approach will be to compare shapes of f US and f MRI in an elastic framework. This requires the definition of the elastic metric and a set of tools to register and compute geodesics (and geodesic distances) between surfaces. We present these details next.

A. Representation of Surfaces

Let F be the space of all smooth embeddings of a cylinder in R 3 and let Γ be the set of all boundary-preserving diffeomorphisms from S 1 × [0, 1] to itself. It acts naturally on F by composition: (f, γ) = f • γ, where • represents the composition of f and γ. In other words, it is a reparametrization of the surface f using γ. Since F is a vector space, the tangent space at

f ∈ F, T f (F ), is F itself. Using two tangent vectors v 1 , v 2 ∈ T f (F ) we can define the standard inner product on F as follows ⟨v 1 , v 2 ⟩ = S 1 ×[0,1] ⟨v 1 (s), v 2 (s)⟩ ds,
where ds is the standard Lebesgue measure on S 1 × [0, 1] and the L 2 norm between any two

surfaces f 1 , f 2 ∈ F is S 1 ×[0,1] |f 1 (s) -f 2 (s)| 2 ds 1 2 .
As was previously shown in many works ( [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF], [START_REF] Kurtek | A novel Riemannian framework for shape analysis of 3D objects[END_REF], [START_REF] Kurtek | Parameterization-invariant shape comparisons of anatomical surfaces[END_REF], [START_REF] Kurtek | Elastic geodesic paths in shape space of parametrized surfaces[END_REF]), such a framework is not appropriate for analyzing shapes of surfaces because the action of Γ does not preserve distances, i.e. ∥f 1 -

f 2 ∥ ̸ = ∥f 1 • γ -f 2 • γ∥.
One solution is to restrict the analysis to area-preserving diffeomorphisms [START_REF] Gu | Matching 3D shapes using 2D conformal representations[END_REF], but this is very limiting in practice. Another approach, taken by [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF], [START_REF] Kurtek | A novel Riemannian framework for shape analysis of 3D objects[END_REF], [START_REF] Kurtek | Parameterization-invariant shape comparisons of anatomical surfaces[END_REF], [START_REF] Kurtek | Elastic geodesic paths in shape space of parametrized surfaces[END_REF], is to utilize a different representation of surfaces such that the action of Γ preserves L 2 distances. Such representations are motivated by Riemannian metrics that: (1) make the action of Γ on F by isometries, and (2) help quantify deformations needed to transform surfaces into each other. Jermyn et al. [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF] defined a specific Riemannian metric and the corresponding mathematical representation to perform registration of spherical and quadrilateral surfaces; in this framework, the registration of points across surfaces is accomplished by re-parameterization of surfaces. The use of this metric is very intuitive in our application as one of our goals is to define a measure of elasticity of the ultrasound surface. We extend their work by considering cylindricallyparametrized surfaces, and in addition, we define a numerical technique termed path-straightening to compute geodesic paths and distances in F similar to [START_REF] Kurtek | Elastic geodesic paths in shape space of parametrized surfaces[END_REF]. Let n(s) = ∂f ∂u (s)× ∂f ∂v (s) ∈ R 3 denote the normal vector to the surface at the point s = (u, v) ∈ S 1 × [0, 1]. Using n(s), we define a mathematical representation of surfaces termed square-root normal fields (SRNFs) [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF]. Definition 1. For a surface f : S 1 × [0, 1] → R 3 define a mapping Q : F → L 2 and the square root normal field (SRNF)

q : S 1 × [0, 1] → R 3 as Q(f )(s) = q(s) = n(s) √ |n(s)|
.

The L 2 -norm of q equals the surface area of f and thus the space of SRNFs is a subset of L 2 (S 1 ×[0, 1], R 3 ) from now on referred to as L 2 . This representation of surfaces has a number of important properties:

1) As shown in [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF], a complicated elastic metric on F simplifies to the standard L 2 metric on the space of SRNFs. This greatly simplifies the development of certain shape analysis tools such as registration and deformation. It also provides explicit interpretations of the types of deformations measured by the distance function introduced later.

2) The SRNF q does not depend on a translation of f , since it only depends on spatial derivatives of f . 3) SO( 3) is the special orthogonal group, which contains all 3 × 3 rotation matrices. For O ∈ SO(3) and f ∈ F, the SRNF of a rotated surface Of is Oq. Γ is the reparameterization group. For a γ ∈ Γ and f ∈ F, the SRNF of a re-parameterized surface f • γ is (q, γ) = (q •γ) J γ , where J γ is the determinant of the Jacobian of γ. Furthermore, the actions of SO(3) and Γ on L 2 commute. 4) The action of the product group SO(3) × Γ on L 2 is by isometries. In other words, if we rotate and reparametrize any two surfaces in the same way, the L 2norm between their SRNFs remains unchanged:

∥q 1 - q 2 ∥ = ∥(Oq 1 , γ) -(Oq 2 , γ)∥ [19].
Since our notion of shape is invariant to translation, scale, rotation and re-parameterization, we need to ensure these invariances. The SRNFs are already invariant to translations and we can achieve scale invariance by re-scaling all surfaces to have unit surface area. The re-sclaling is problematic in our application because we do not know the relative scale of β US and the third coordinate in f US . Thus, we take a different approach. We re-scale all of the surface coordinates of f US individually to lie in [0, 1]. In order to define a rotation and parameterization invariant distance on L 2 we first define an equivalence class of an SRNF as all possible rotations and re-parameterizations of this function: [q] = {(Oq, γ)|O ∈ SO(3), γ ∈ Γ}. The distance between shapes of surfaces is then given by:

d([q 1 ], [q 2 ]) = inf (O,γ)∈SO(3)×Γ ∥q 1 -(Oq 2 , γ)∥. (1)
Note that while this distance is important in quantifying differences between shapes of surfaces it defines an extrinsic distance. It does not provide a recipe for computing geodesic paths between surfaces due to difficulties in computing the inverse mapping of Q.

It is also important to point out that the distance in Eqn. 1 is fundamentally different from the distance proposed in [START_REF] Kurtek | Parameterization-invariant shape comparisons of anatomical surfaces[END_REF]. The distance in [START_REF] Kurtek | Parameterization-invariant shape comparisons of anatomical surfaces[END_REF] is based on a different representation of surfaces termed the q-map (defined as q(s) = |n(s)|f (s)). The advantages of the SRNF representation compared to the q-map representation are described in more detail in [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF].

B. Registration of Cylindrical Surfaces

The computation of d in Eqn. 1 requires solving the joint optimization problem on SO(3) × Γ. This procedure, albeit with different mathematical representations of surfaces, has been previously provided in [START_REF] Kurtek | A novel Riemannian framework for shape analysis of 3D objects[END_REF], [START_REF] Kurtek | Parameterization-invariant shape comparisons of anatomical surfaces[END_REF], [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF] for different parameter domains (unit sphere, unit square, etc.) and we give a brief summary next for convenience. The novelty of this paper lies in extending these methods to the case of the cylindrical domain

S 1 × [0, 1].
For a fixed γ ∈ Γ, the minimization over SO(3) can be easily performed directly using Procrustes analysis. Let q 2 denote (q 2 , γ); then, the optimal value of O is obtained as follows. Compute the 3 × 3 matrix A = S 1 ×[0,1] q 1 (s)q 2 (s) T ds. Using singular value decomposition A = U ΣV T , we can define the optimal rotation as O * = U V T (if the determinant of A is negative, the last column of V changes sign).

In order to solve the optimization problem over Γ (for a fixed rotation), we will use a gradient approach. To specify the gradient for iterative updates, we define the cost function at any iteration to be:

E : Γ → R, E[γ] = ∥q 1 -(q 2 , γ)∥ 2 = ∥q 1 -φ(γ)∥ 2
, where φ : Γ → [q 2 ] is defined to be φ(γ) = (q 2 , γ). Let b be a unit vector in T γ id (Γ) for γ id (s) = s. The tangent space T γ id (Γ) is the set of all smooth vector fields tangent to the cylinder with appropriate boundary conditions. Using an orthonormal basis for T γ id (Γ) we can compute the full gradient of E at γ id and use it to update the current estimate of the optimal parameterization of q 2 . In principle, this tangent space is infinite dimensional, but we restrict to a finite number of basis elements. This basis is formed by considering the coordinates on S 1 and [0, 1] separately. Tangent vector fields on S 1 can be written as {sin(k

1 u), 1 - cos(k 1 u), u, 1 -u | k 1 = 1, . . . , n 1 , u ∈ S 1 }. Tangent vector fields on [0, 1] can be written as {sin(2πk 2 v), 1 - cos(2πk 2 v) | k 2 = 1, . . . , n 2 , v ∈ [0, 1]}.
We take all possible products of the elements of these two sets and orthonormalize them using Gram-Schmidt. Similar basis construction for the sphere and unit square domains is described in [START_REF] Kurtek | Parameterization-invariant shape comparisons of anatomical surfaces[END_REF], [START_REF] Kurtek | A novel Riemannian framework for shape analysis of 3D objects[END_REF]. Some basis elements (vector fields on the cylinder) are shown in Fig. 3. With a basis of T γ id (Γ) , the directional derivative of E at γ id , in the direction of b, is given by ⟨q 1 -φ(γ id ), φ * (b)⟩b, where φ * : T γ id (Γ) → T q2 ([q 2 ]) is the differential of φ. (An expression for φ * can be found in [START_REF] Kurtek | Parameterization-invariant shape comparisons of anatomical surfaces[END_REF]. ) Furthermore, the full gradient can be approximated using b ⟨q 1φ(γ id ), φ * (b)⟩b.

To solve the joint minimization problem in Eqn. 1, we alternate between minimizing over SO(3) and Γ and stop when we converge to a solution. Fig. 4 shows an example of this optimization for two surfaces -these surfaces have two bumps each but the bumps are placed at different points along the axial coordinate. The final registration provides a better matching of features, as depicted using identical color maps, across the two surfaces. The importance of registration is further emphasized later in the bottom row of Fig. 5 where we show geodesic paths between the same two surfaces, before and after the registration.

III. GEODESICS BETWEEN CYLINDRICAL SURFACES

As mentioned earlier, it is not a simple task to compute the inverse mapping of Q. This prevents us from computing the geodesic between surfaces by simply mapping the straight line path from the space of SRNFs to F . For this reason, we develop a numerical technique, termed path-straightening, to compute geodesic paths and geodesic distances (these will now be intrinsic distances between surfaces) directly in F . We will do this using the pullback metric from L 2 .

In order to define the pullback metric on F we begin by defining the differential of Q, Q * :

T f (F ) → T Q(f ) (L 2 ) (for w ∈ T f (F )): Q * ,f (w)(s) = d dϵ Q(f + ϵw)(s)| ϵ=0 = n w (s) |n(s)| - (n(s) • n w s)n(s) 2|n(s)| 5/2 , (2)
where n w (s) = ∂f ∂u (s) × ∂w ∂v (s) + ∂w ∂u (s) × ∂f ∂v (s). Note that w represents a vector field on the surface f . Using this expression, we can define a Riemannian metric on F for w 1 , w 2 ∈ T f (F ) as follows:

⟨⟨w1, w2⟩⟩ f = ⟨Q * ,f (w1), Q * ,f (w2)⟩ = (nw 1 (s) • nw 2 (s)) |n(s)| - 3 4 (n(s) • nw 1 (s))(n(s) • nw 2 (s)) |n(s)| 3 (3)
In order to compute a geodesic path between two given surfaces f 1 , f 2 ∈ F, we begin by defining the length of a path given by F : [0, 1] → F, F(0) = f 1 , F(1) = f 2 using the pullback metric (all arguments have been suppressed for brevity):

L(F ) = 1 0 ⟨⟨Ft, Ft⟩⟩ F dt = 1 0 ⟨Q * ,F (Ft), Q * ,F (Ft)⟩dt = 1 0 [0,1]×S 1 -3(N • Nt) 2 4|N | 3 + |Nt| 2 |N | dsdt. (4)
In the above expression, we use N (t) to denote ∂F (t) ∂u (s) ×

∂F (t)
∂v (s). N (t) is a normal vector field on the surface f = F (t). Also, we denote partial derivatives using subscripts (i.e. N t ). It is a well known result that a critical point of this path energy provides a geodesic between f 1 and f 2 in F . Equivalently, one can consider the following energy, which is computationally easier to work with:

Ẽ(F ) = 1 0 ⟨⟨F t , F t ⟩⟩ F dt = 1 0 ⟨Q * ,F (F t ), Q * ,F (F t )⟩ dt = 1 0 [0,1]×S 1 -3(N • N t ) 2 4|N | 3 + |N t | 2 |N | dsdt.
(

) 5
Thus, a geodesic path between two surfaces f 1 and f 2 is the solution to the following minimization problem:

F * = argmin F :[0,1]→F , F(0)=f1, F(1)=f2 Ẽ(F ). ( 6
)
We will compute the solution to this problem using a gradient descent approach. We will approximate ∇ Ẽ using directional derivatives ∇ ẼF (G) where G ∈ G and G is an orthonormal basis that spans all possible perturbations of a path F . We begin by computing ∇ ẼF (G):

∇ ẼF (G) = d dϵ Ẽ(F + ϵG)|ϵ=0 = [ -3(N • δGNt + Nt • δGN )(N • Nt) -2|Nt| 2 (N • δGN ) 2|N | 3 + 9(N • δGN )(N • Nt) 2 4|N | 5 + 2(Nt • δGNt) |N | ]dsdt, (7)
where

δ G N = G u × F v + F u × G v and δ G N t = G tu × F v + F tu × G v + G u × F tv + F u × G tv (
subscripts are used for time and spatial derivatives here). Then, in order to approximate the gradient we will use a finite number of basis elements G i , i = 1, ..., B: ∇ Ẽ ≈ B i=1 ∇ ẼF (G i )G i . We will update the current estimate of the geodesic path according to F = F -η∇ Ẽ, where η > 0 is a small step size. Once the algorithm has converged, it results in the geodesic path (F * ) as well as the unique shooting vector field (F t (0)), which can be used to optimally deform f 1 to f 2 under our metric.

The basis set {G i } is created as follows. Each element G : (S 1 ×[0, 1])×[0, 1] → R 3 has three arguments u, v, t where (u, v) are the spatial coordinates on S 1 × [0, 1] and t is the time index along the path. We first form bases for each of the u, v, t components separately:

P 1 (u) = {sin(iu), cos(iu) - 1, u, 1 -u | i = 1, . . . , n 1 , u ∈ S 1 }, P 2 (v) = {sin(2πjv), cos(2πjv), 1 | j = 1, ..., n 2 , v ∈ [0, 1]}, P 3 (t) = {sin(2πkt), cos(2πkt) -1 | k = 1, . . . , n 3 , t ∈ [0, 1]}.
In order to define the basis P : (S 1 ×[0, 1])×[0, 1] → R we utilize all possible products of the three component bases P 1 (u), P 2 (v) and P 3 (t). The final step in forming the full basis is to utilize a copy of P for each of the Euclidean coordinates (in R 3 ) and to orthonormalize this basis using the Gram-Schmidt procedure.

Fig. 5 shows examples of minimization of Ẽ using path straightening between surfaces before (top) and after (bottom) registration. Fig. 6 shows a geodesic between two cylindrically-parametrized surfaces f 1 and f 2 at the end points and uniformly spaced samples along the optimal path F * . IV. APPLICATION TO 2D-3D REGISTRATION Now that we have laid out tools for registering and deforming cylindrically-parameterized surfaces, we return to the four major goals outlined earlier and show how they are achieved using this framework. Starting with a surface f MRI and a curve β US , we first form a generalized cylinder f US from β US . Then, we perform a registration of f US and f MRI using the framework outlined in Sections II-A and II-B. Given this optimal registration, we compute the geodesic distance and geodesic path between these surfaces using tools developed in Section III. This provides an optimal deformation between f US and f MRI under the elastic metric and the shooting vector field on f US (we always use f US as the source surface and f MRI as the target surface). We then define a function

H : [0, 1] → R + as H(v) = u∈S 1 |F * t (0)(u, v)| 2 du
, where Ḟ * t (0) is the initial velocity vector of the geodesic path from f US to the re-parameterized f MRI . We remind the reader that the shooting vector is sufficient to characterize the geodesic path or the deformation from one shape to the other. Let v * = argmin v∈[0,1] H(v). The quantity H(v * ) represents our estimate of a quantification of the elasticity of the ultrasound surface. We also use v * to find the curve on f MRI that most corresponds to β US , as follows. Select the curve on f MRI that registers with f US (•, v * ) under the optimal registration obtained above and call it β MRI . It is important to note that β MRI can be non-planar. This is a reasonable possibility because the corresponding transformation to optimally register the ultrasound and MRI images is often nonlinear.

Figure 7 provides a depiction of the entire pipeline of our framework and presents results for one of the subjects in our study. In this case, we found that the measure of elasticity of the endometrial tissue was given by H(v * ) = 13.54 (it corresponds to the red point in Figure 7(a)). Fig. 7. Pipeline of our framework. The curve β US leads to a generalized cylinder f US that is registered and deformed into f MRI using elastic shape geodesic. The geodesic provides (a) H(v), the magnitude of the shooting vector field at each level of the US surface, with its minimum value H(v * ) marked using a red dot, (b) a visualization of the magnitude of this vector field on f US (hot colors correspond to higher magnitude), (c) the deformation vector field on the selected β US and (d) its magnitude, (e) the optimal curve β MRI on the MRI surface, and (f) the deformation between the β US and β MRI curves.

V. EXPERIMENTAL RESULTS

Now we apply our framework on real data to estimate the quantities of interestβ MRI , elasticity measure H(v * ), etc.and to achieve the stated goals. However, since an important part of the experimentation is the evaluation of results, we are interested in quantifying our estimation performance. The difficulty in this evaluation comes from the fact that the ground truth is unknown. It is even difficult to visually evaluate the registration due to the complexity introduced by different dimensions (2D versus 3D) and different modalities (TVUS versus MRI). Therefore, it is important to generate some semisynthetic data where the ground truth is known by construction and one can easily compare the estimated quantities with the ground truth. In terms of evaluating inferences on real data, we will rely on annotations made by the experts, for β MRI , and compare them with the estimated β MRI . These experimental results and their evaluations are described next.

A. Validation on Semi-Synthetic Data

As stated earlier, the non-rigid deformation Φ is due to the transducer's pressure and, thus, the deformation location and magnitude depend on: 1) the position of the transducer, especially the surface of contact between the transducer and the organ; 2) the magnitude and the direction of its pressure. In order to imitate the real-world process, we choose an elastic deformation model that follows these physical principles and generate a controlled deformation as follows:

1) Pick up a random point on the surface around which we define a contact area between the transducer and the organ. 2) Apply different strengths of deformations, ranging from weak to strong, controlled by the deformation model. 3) Choose random directions to simulate the direction of pressure. Once the original surface f MRI is deformed, using one of the deformation fields Φ described above, we randomly choose a point on the deformed surface Φ(f MRI ) and a normal vector to define an intersecting plane P . The intersection P ∩ Φ(f MRI ) is a planar curve β US , which plays the role of a manually extracted curve from the TVUS image. Fig. 8 illustrates the experimental design for the construction of synthetic deformations, for which the correspondence between surfaces is known. In the process, we generate an estimate βMRI of β MRI on the original surface f MRI . An example of geodesic paths between the original MRI-surfaces and the reconstructed TVUS-surfaces is given in Fig. 9. Since we have the original curve β MRI and its estimate βMRI , we can quantify the estimation performance quite naturally. In the top panel of Fig. 10, we present an evaluation of the differences between these curves for ten simulated cases under four different metrics -root-mean squared (RMS or STD), supremum norm, and fiducial registration error (FRE). As we can see, especially using the supremum norm, the error between the curves is very small in most of the cases and our framework is able to estimate β MRI faithfully. For comparison purposes, we have applied the same experimental protocol to register f MRI and f US using an ICP-based approach and present these results in the bottom panel of Fig. 10. Our approach significantly outperforms this method.

B. Qualitative Validation

We start by displaying results of geodesic deformations between two surfaces in Fig. 11. In each example, we display the geodesic path, and the deformation between the ultrasound curve and the optimal MRI curve chosen using our method. The correspondence of points across surfaces is mapped onto the whole geodesic where points corresponding to each other are shaded with the same colors. In addition, we plot the optimal curves along the geodesic path in black color.

As a main qualitative approach to validation of our estimation, we investigate results on real data from ten patients. As described earlier, we obtained f MRI and β US for each of these cases using expert segmentation and used them to compute the estimate βMRI . In this case, the ground truth (β MRI ) is not available and we asked two experts to evaluate the results. The experts were shown the estimated curve βMRI , the deformation between β US and the estimated curve, and a plane that contains as much of βMRI as possible. The evaluation is given in Fig. 12, and examples of H(v) (with H(v * ) highlighted in red) and the estimated magnitude of the vector field on β U S are shown in Fig. 13. To evaluate the performance of the proposed method, we setup a registration quality score ranging from 1 (bad) to 5 (excellent). The experts rated eight out of ten as 4, one as 3, and one as 5 -a very high evaluation of the estimation performance. We note that the experts' evaluation is based on the anatomical features, tissue structures, and the surrounding organs, like the uterus. We observe that the deformation magnitude provides a nice display of elasticity along the curve. Such information is very useful when coupled with the localization of the tissue (surrounding organs) as it will help a surgeon evaluate the phase of the disease and the tissue's mobility. Since our approach provides the optimal curve βMRI on the MRI surface that matches the original TVUS curve, we can locate βMRI in the MRI volume. In Fig. 14, we show this result on one patient. Given the segmented MRI surface f MRI (Fig. 14(a)) and the TVUS image (Fig. 14(b)) we describe the clinical utility of our approach as following: first, we have constructed f US and then applied our approach to register them and to estimate βMRI on f MRI ; second, we have computed its barycentre and its normal to define the cutting plane P ; third, we have interpolated the corresponding plane P MRI between the MRI volume and P ; finally, we have projected βMRI on P MRI with small changes to show its closest position on P MRI (Fig. 14(d)). To show the performance of the registration, Fig. 14(c) displays the overlap of the original curve β US and the estimated curve βMRI .

As stated in the introduction, our clinical goals in this work were to combine complementary information from MRI and TVUS (transvaginal ultrasound) images of rectovaginal endometriotic nodules, to locate the US intersecting plane in the MRI volume as shown in Fig. 14(d), and to quantify the deformability of endometrial tissues and surrounding structures as shown in Fig. 13. These goals were reached using two models: (1) a clinical model from TVUS and MRI images (10 patients) and (2) a semi-synthetic model (original data were deformed using a biomechanical model to simulate physical deformations). A segmentation (delimitation) of the nodule was performed in the two modalities (MRI and TVUS) and a recent shape analysis approach was applied, with appropriate modifications. The study resulted in an accurate quantification and localization of the deformable endometriotic tissues.

In short, using both semi-synthetic simulated data and the original real data from ten patients, we were able to visually improve the preoperative diagnosis by facilitating interpretation of ultrasound data. Specifically, combining images from MRI and TVUS allows us to compensate for differences in patient positioning during MRI and TVUS observations, and the potential deformation of endometrial tissue during the TVUS measurement process (transducer's pressure). An immediate impact of these results will be in improving the planning of surgical procedures. Furthermore, we believe that combining complementary information from both modalities, by locating the US plane inside the MRI volume, and quantifying deformability of endometriosis under a given level of applied stress can open new avenues for diagnosing and monitoring complex cases.

Nevertheless, there are some limitations in the performed study. One is the use of boundaries that were manually extracted by the experts to represent the true shape of the endometrial tissue. As a consequence, the deformation vector field depends on the quality of the segmentation. This can introduce a bias when the quality of the imaging modalities decreases or when the endometrial tissues are present in very small areas. Another application of this methodology (with sufficient adaptation to the new problem) can be the registration of cylindrical representations of the colon in different modalities to prevent the development of advanced colorectal cancer. Such study will allow practitioners to quantify the uncontrolled cell growth in the colon or rectum that may cause colorectal cancer.

VII. SUMMARY

In this work, we have proposed an elastic shape analysis framework for 2D-3D TVUS-MR registration to localize and characterize endometrial tissues. The information is available in the form of parametrized surfaces in the 3D MRI data and segmented closed curves in the 2D TVUS images. Our approach extrapolates the US curve into a generalized cylinder and applies an elastic shape analysis technique to obtain registration and deformation between the two surfaces. This way, we obtain the optimal curve on the MRI surface that corresponds to the original TVUS curve. We also obtain deformations and quantification of the deformations as measures of elasticity at different regions on the tissue.

Using both semi-synthetic and real data, we are able to validate the estimation results and demonstrate the appropriateness of elastic shape analysis in this context. In the semi-synthetic data, the closeness of estimated quantities (curves and deformations) with the ground truth underscore the success of this framework despite a loss of information inherent in the 2D nature of US imaging. A qualitative evaluation by two experts, that involves more syntactical information, generally rated the estimation results from very good to excellent.

Fig. 1 .

 1 Fig. 1. Example of endometrial tissue localization outside the pelvic cavity shown in laparoscopic (a), TVUS (b), and MR (c) images.

Fig. 2 .

 2 Fig. 2. An illustration of the mapping from an MRI surface f MRI to the TVUS curve β US . Left: an open cylinder depicting an MRI surface f MRI . Middle: Φ(T (f MRI )) after a nonlinear deformation Φ and its intersection with a cutting plane. Right: Resulting in β US .

Fig. 3 .

 3 Fig. 3. Examples of basis elements of Tγ id (Γ) that are used in incremental optimization over Γ. For each example we display the full vector field (top) and a zoom-in on a specific part (bottom).

Fig. 4 .

 4 Fig. 4. Using identical colors to denote registration, the figure shows improvement in registration of similar features across the two surfaces during optimization over Γ.

Fig. 5 .

 5 Fig. 5. Geodesic computation using path straightening between surfaces before (top) and after (bottom) registration.

Fig. 6 .

 6 Fig. 6. Geodesic between two cylindrically-parametrized surfaces.

Fig. 8 .

 8 Fig. 8. An illustration of semi-synthetic data. (a) f MRI , (b) a deformation vector field generating Φ, (c) the deformed surface Φ(f MRI ), and (d) a cutting plane P resulting in β US .

Fig. 9 .

 9 Fig. 9. The geodesic between f US and f MRI surfaces, with the estimated curve βMRI marked in black on f MRI .

Fig. 12 .

 12 Fig.[START_REF] Vermandel | A 2d/3d matching based on a hybrid approach: improvement to the imaging flow for avm radiosurgery[END_REF]. Qualitative evaluation by two experts. For each patient, the performance is evaluated by a score from 1 (bad) to 5 (excellent).

Fig. 13 .

 13 Fig. 13. Estimation of βMRI : for each patient, we display (a) H(v) and (b) the corresponding magnitude of the deformation vector field on β US .

Fig. 14 .

 14 Fig. 14. (a) The MRI surface, (b) the TVUS image, (c) the overlap between β US from (b) and the estimated curve βMRI on (a), and (d) βMRI on the estimated intersecting plane P .

ACKNOWLEDGEMENTS

Anuj Srivastava was supported in part by the NSF DMS 1208959 and Chafik Samir was supported in part by Auvergne Regional Council grants.