Martino Alessandrini
email: martino.alessandrini@creatis.insa-lyon.fr

Olivier Bernard

Adrian Basarab

Hervé Liebgott

Multiscale optical flow computation from the monogenic signal

Introduction

The monogenic signal has been recently introduced by Felsberg and Sommer [START_REF] Felsberg | The monogenic signal[END_REF] as an extension of the analytic signal concept to multiple dimensions. Similarly to the latter, the monogenic signal provides the local amplitude and local phase signal features. It also contains information on the local orientation. These three local features are pointwise orthogonal, which means that they represent independent information: the local amplitude represents the local intensity or dynamics, the local phase describes the local symmetry or grey value transition, and the local orientation describes the direction of the highest signal variance.

Decoupling the local energy from the image structure, accounted for by phase and orientation, has made it possible to derive effective solutions to a number of image-processing problems, in particular when the more traditional pixel intensity cannot be considered as a reliable feature. In this context, we have recently proposed an original optical flow algorithm for the analysis of heart motion based on the monogenic phase [START_REF] Alessandrini | Myocardial motion estimation from medical images using the monogenic signal[END_REF].

Temporal brightness variations are indeed the norm in many medical imaging modalities such as ultrasound and MRI, and this fact negatively affects the performance of intensity based algorithms [START_REF] Alessandrini | Myocardial motion estimation from medical images using the monogenic signal[END_REF].

Thanks to its general formulation, the proposed algorithm is suitable for images from several imaging modalities. In particular, we summarize in this paper some of the results we have obtained on 2D and 3D cardiac ultrasound and on tagged MRI. In all cases, the proposed algorithm was shown to be competitive with state-of-the-art algorithms in the respective fields and extremely effective from a computational point of view.

The paper proceeds as follows. In Section 2.1, the monogenic signal theory is briefly resumed. In Section 2.2, the optical flow estimation algorithm presented in [START_REF] Alessandrini | Myocardial motion estimation from medical images using the monogenic signal[END_REF] is described. In Section 3, the results are presented. Concluding remarks are left to Section 4.

Proposed optical flow estimation algorithm

Monogenic signal computation

Only the most basic aspects of the monogenic signal computation are given here, the reader is otherwise addressed to [START_REF] Felsberg | Optical flow estimation from monogenic phase[END_REF] and the references therein for a theoretically founded derivation.

The monogenic signal provides the features local phase φ (x), orientation θ (x) and amplitude A(x) of an image I(x), where x = x, y T is the pixel position. These quantities are computed pointwise from the responses to three 2D spherical quadrature filters (SQF) obtained as follows: one even rotation invariant bandpass filter and two odd filters, computed from the Riesz transform of the even one. Among the existing SQF families, we adopt in this work difference of Poisson (DoP) [START_REF] Felsberg | Optical flow estimation from monogenic phase[END_REF]. Alternatively to the standard pointwise estimate we employ a more robust least-squares technique for the computation of θ [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz laplace wavelet transform[END_REF]. This is obtained by maximizing the directional Hilbert transform H θ I (x) averaged over a local neighborhood, here a Gaussian kernel with variance σ 2 , denoted by υ σ :

θ (x) = arg max θ ′ ∈ [-π,π] R 2 υ σ x ′ -x • H θ ′ I x ′ 2 dx ′ . (1)
As shown in [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz laplace wavelet transform[END_REF], the optimization problem is solved by the eigenvector associated to the maximum eigenvalue of the Riesztransform counterpart of the image structure tensor.

Optical flow computation

Following [START_REF] Felsberg | Optical flow estimation from monogenic phase[END_REF], we compute the displacement d = [d 1 , d 2 ] T between two frames assuming the conservation of the monogenic phase over time φ (x, t + 1) = φ (xd (x) , t). Then, assuming all points translate of the same quantity d 0 within a local window w centered in x 0 , the following linear system of equations is obtained:

J w d 0 = -r t w , J = f nn T (2)
where v w = w (xx 0 ) v (x) dx, r t is the time derivative of the phase vector r = φ • n, n = [cos (φ) , sin φ] T and f is the monogenic frequency, computed as f = (∇φ) T n, with ∇ = ∂ x , ∂ y T . The expressions for the computation of r t and f from SQF filter outputs can be found in [START_REF] Felsberg | Optical flow estimation from monogenic phase[END_REF]. Instead of the simple translation model, we adopt a more general affine model for the local displacement. In addition to translation, this accounts for rotation, expansion, compression and shear. Considering a window centered at the origin, the affine model is expressed as:

d (x) = A (x) u, A = 1 0 x y 0 0 0 1 0 0 x y (3

Multiscale choice of the window size

The choice of the window size is a tedious issue connected with local techniques: the assumed motion model (translational or affine) may not hold when the window is too big, otherwise, the adoption of an excessively small window may result in the well-known aperture problem. To circumvent this issue, in [START_REF] Sühling | Myocardial motion analysis from b-mode echocardiograms[END_REF] Sühling et al. proposed a multiscale strategy for locally choosing the most consistent window size. This is based on the possibility of computing the image moments, i.e., the entries of the system matrix M and the vector b, at multiple scales, by using an efficient B-spline coarse-to-fine strategy. In particular, they are obtained from window functions w that are progressively scaled and subsampled by a factor 2 in each dimension. Among the considered scales considered, the u j producing the smallest residual error Mu jb ℓ2 /|w| ℓ1 is retained as the final displacement estimate. With this strategy, the scale providing the most consistent motion estimate is selected.

Iterative displacement refinement

The hypothesis of small displacements employed in differential techniques may be inadequate whenever the displacement is big or the image intensity profile is non-linear. To overcome this limitation we adopt an incremental coarse-to-fine refinement of the motion estimate: the coarsest scale provides a first estimate of the displacement, this estimate is then used to undo the motion and then the estimator is reapplied to the warped signals at the finer scale to find the residual motion.

Results

We evaluated the proposed algorithm on medical images from different modalities. In order to have a benchmark to compare to, we made use of synthetic sequences in all cases. Endpoint error:

EE = d -d ℓ2 (4)
was used for the evaluation, where d and d denote estimated displacement and benchmark respectively. Unless explicitly mentioned, the displacement has to be intended as the one between two consecutive frames.

For display purposes, we adopted a color encoding of the velocities (Fig. 1), where the color indicates the direction of the displacement and the brightness expresses the magnitude.

Since what is of clinical interest is evaluating the mechanical properties of the cardiac muscle, in all cases we will measure the estimation error within the myocardium only (Fig. 2).

2D cardiac ultrasound

We evaluated our algorithm on custom made realistic synthetic cardiac ultrasound sequences, generated as described in [START_REF] Alessandrini | Simulation of realistic echocardiographic sequences for ground truth validation of motion estimation[END_REF]. Briefly, the simulated sequences are obtained by mimicking the aspect and the motion of a real echocardiographic acquisition taken as a template. For our evaluation, we considered simulated sequences from the two views the most commonly employed in the clinical practice, namely short-axis (SaX) and apical four chambers (A4C). Fig. 2 presents two frames from the simulated SaX sequence along with the benchmark motion.

We have compared our algorithm with the ones of Sühling et al. [START_REF] Sühling | Myocardial motion analysis from b-mode echocardiograms[END_REF] and Felsberg [START_REF] Felsberg | Optical flow estimation from monogenic phase[END_REF]. Sühling's algorithm exploits the standard intensity conservation assumption and, as in our case, employs an affine displacement model and a multiscale window selection. Differently, Felsberg's algorithm makes use of the monogenic phase but employs a simple translation model and a window of fixed size.

Table 1 reports the average errors obtained on the entire simulated sequences (45 frames for the SaX view and 50 frames for the A4C). From Table 1 all the monogenic phase-based algorithms considered perform better than Sühling's algorithm. This confirms that the monogenic phase is a more reliable feature than pixel intensity as far as medical ultrasound is concerned. Also, thanks to the affine model, the multiresolution window selection procedure and the pyramidal refinement, the proposed algorithm outperforms Felsberg's one. For further results on 2D cardiac ultrasound we address the reader to [START_REF] Alessandrini | Myocardial motion estimation from medical images using the monogenic signal[END_REF][START_REF] Alessandrini | Monogenic phase based myocardium motion analysis from cardiac ultrasound with transverse oscillations[END_REF].

3D cardiac ultrasound

Despite our algorithm has been described in this paper for the displacement estimation from 2D images, an extension to 3D volumes can be almost straightforwardly derived [START_REF] Alessandrini | Monogenic phase based optical flow computation for myocardial motion analysis in 3D echocardiography. Statistical Atlases and Computational Models of the Heart Imaging Model Challenges[END_REF]. Also in this case we evaluated our algorithm on synthetic 3D cardiac ultrasound data. In particular, the latter were provided as part of the "Tracking Challenge" promoted by the MICCAI conference [START_REF] Alessandrini | Monogenic phase based optical flow computation for myocardial motion analysis in 3D echocardiography. Statistical Atlases and Computational Models of the Heart Imaging Model Challenges[END_REF].

The heart motion was simulated by displacing a set of point scatterers according to the electromechanical model of the myocardium proposed by Sermesant et al. [START_REF] Sermesant | Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in crt: a preliminary clinical validation[END_REF]. From the time varying scatter map the ultrasound image formation was then simulated with the COLE software developed by Gao et al. [START_REF] Gao | A fast convolution-based methodology to simulate 2-d/3-d cardiac ultrasound images[END_REF]. Each simulated sequence realized a full cardiac cycle from one end-diastole to the next end-diastole. The algorithm evaluation was then performed by measuring the accuracy in tracking the scatterers employed in the simulation. Note that, differently from the previous cases, this implies computing the global displacement from the first frame of the sequence (end-diastole). The global field is therefore computed by accumulating the frameby-frame displacements returned by the proposed algorithm.

In particular, we report here the result on one of the provided sequences representing a healthy heart. On Fig. 3 the error behavior over time is represented as a boxplot. As expected, the largest estimation error corresponds to the end-systolic instant, where the maximum displacement from the rest condition happens. Nevertheless, lets note that the average error remains much smaller with respect to the true maximum displacement that in the end-systolic instant reaches the value of 11.2 mm, also Fig. 4. For a more detailed description of the 3D extension of our algorithm and a more comprehensive performance evaluation we address the reader to [START_REF] Alessandrini | Monogenic phase based optical flow computation for myocardial motion analysis in 3D echocardiography. Statistical Atlases and Computational Models of the Heart Imaging Model Challenges[END_REF].

Tagged MRI

The proposed algorithm is compared with SinMod [START_REF] Arts | Mapping displacement and deformation of the heart with local sine-wave modeling[END_REF], available in the InTag plugin for OsiriX. SinMod was shown to perform better than the state-of-the-art HARP in [START_REF] Arts | Mapping displacement and deformation of the heart with local sine-wave modeling[END_REF]. The evaluation was made on synthetic tMRI sequences, generated with the ASSESS software [START_REF] Clarysse | Simulation based evaluation of cardiac motion estimation methods in tagged-mr image sequences[END_REF]. The synthetic motion is established on the basis of a 2D analytical model taking typical contraction, relaxation, torsion and thickening of the cardiac muscle into account.

The results obtained on nine simulated sequences are summarized in Table 2. These results show that the proposed algorithm systematically returns the estimate with the smallest mean value and variance, which is a proof of precision and reliability.

To better appreciate the difference in performance, it is useful to analyze the local behavior of each algorithm. This is represented on Fig. 5, where the error images obtained on the 4th frame of the two sequences D30F20 and R20F20, implementing a pure contraction/expansion and pure rotation respectively, is displayed. At that instant, the displacement reaches the maximum average value and the greatest spatial variation in both cases: in the first case (first row on the Fig. 5) the angular velocity decreases linearly, passing from the endocardial to the epicardial contour; in the second (second row on the Fig. 5) the radial contraction is null on the epicardium and maximal on the endocardium. As shown by the previous results, SinMod is outperformed by the proposed algorithm.

For further results on tMRI sequences we address the reader to [START_REF] Alessandrini | Myocardial motion estimation from medical images using the monogenic signal[END_REF][START_REF] Alessandrini | Monogenic signal for cardiac motion analysis from tagged magnetic resonance image sequences[END_REF].

Conclusion

We have described a novel algorithm for the analysis of heart motion from medical images. The displacement is estimated from the monogenic phase and is therefore robust to possible variations of the local image energy. A local affine model accounts for the typical contraction, torsion and shear of myocardial fibers. An effective B-spline multiresolution strategy automatically selects the scale returning the most consistent velocity estimate. The multiresolution strategy together with a least-squares estimate of the monogenic orientation make the algorithm robust under image noise. Due to its general formulation, the proposed algorithm is well suited for measuring myocardial motion from images from different modalities. In particular, we have presented an evaluation on simulated 2D and 3D cardiac ultrasound and tagged MRI sequences. The results have shown that the proposed algorithm is a valid alternative to state-of-the-art techniques in the two fields.

The code for the algorithm presented in this paper is freely online available at the url http://www.creatis.insa-lyon.fr/ ustagging/node/13.

  ) where u = d 10 , d 20 , d 1x , d 1y , d 2x , d 2y T is the new unknown vector: d 10 and d 20 correspond to the translation of the window center and d ik = ∂ k d i . Plugging (3) into (2) leads to an underdetermined system of equations, whose least-squares solution is b = Mu, with b = -A T J T r t w and M = A T J T JA w .

Fig. 1 .

 1 Fig. 1. Color encoding of velocities: x and y axes report the horizontal and vertical component of the velocity vector.

Fig. 2 .

 2 Fig. 2. Frames from a synthetic short-axis sequence during ventricular filling (diastole) (a) and ventricular ejection (systole) (b). The benchmark displacement field is superimposed. Note that the benchmark field reflects the physiological expansion and contraction during these two phases of the cardiac cycle.

Fig. 3 .

 3 Fig. 3. Tracking error obtained with the proposed algorithm. Lower and upper limits of the box represent 25th and 75th percentile while the whiskers the 5th and 95th.

Fig. 4 .

 4 Fig.4. Error maps at end-systole, when the displacement is the largest. (a) Real displacement of the scatterers with respect to the first frame of the sequence (end-diastole). This information provides the full-scale of the displacement to be estimated. (b) Estimation error obtained with the proposed algorithm. Note how it is considerably smaller than the maximum displacement reported in (a).

Fig. 5 .

 5 Fig. 5. Error map for the 4-th frame of R20F20 (first row) and D30F20 (second row). (a) and (d) denote the benchmark field.

Table 1

 1 Endpoint error in pixels (µ ± σ).

	ALGORITHM	SEQUENCE
		Apical 4 chambers	Short axis

Sühling 0.395 ± 0.338 0.396 ± 0.346 Felsberg 0.315 ± 0.257 0.364 ± 0.293 Proposed 0.264 ± 0.190 0.313 ± 0.242

Table 2

 2 Endpoint error (µ ± σ) in pixels on nine simulated sequences.

	SEQUENCE	ALGORITHM
		Proposed	SinMod
	D30	0.152 ± 0.121	0.215 ± 0.145
	D30F20	0.082 ± 0.072	0.128 ± 0.112
	D30R10T01P0	0.264 ± 0.149	0.363 ± 0.199
	D30R20T01P0	0.462 ± 0.239	0.970 ± 1.129
	D30R20T01P0F20	0.209 ± 0.139	0.344 ± 0.224
	D30R20T01P3	0.419 ± 0.228	0.911 ± 1.099
	R20F20	0.244 ± 0.164	0.416 ± 0.264
	R10	0.161 ± 0.087	0.220 ± 0.090
	R20	0.104 ± 0.072	0.174 ± 0.122

Acknowledgments

This work was partially supported by US-Tagging grant financed by Agence Nationale de la Recherche (ANR) and by CNRS.

Disclosure of interest

The authors declare that they have no conflicts of interest concerning this article.


I T Lucas
email: ivan.lucas@upmc.fr

S Durand-Vidal

O Bernard

V Dahirel

E Dubois

J.-F Dufrêche

S Gourdin-Bertin

M Jardat

G Meriguet

G Roger

Influence of the volume fraction on the electrokinetic properties of maghemite nanoparticles in suspension

Keywords:

Influence of the volume fraction on the electrokinetic properties of maghemite nanoparticles in suspension

We used several complementary experimental and theoretical tools to characterize the charge properties of well-defined maghemite nanoparticles in solution as a function of the volume fraction. The radius of the nanoparticles is equal to 6 nm. The structural charge was measured from chemical titration and was found high enough to expect some counterions to be electrostatically attracted to the surface, decreasing the apparent charge of the nanoparticle.

DC conductivity measurements were interpreted by an analytical transport theory to deduce the value of this apparent charge, denoted here by "dynamic effective charge". This dynamic effective charge is found to decrease strongly with the volume fraction. On the contrary, the "static" effective charge, defined thanks to the Bjerrum criterion and computed from Monte Carlo simulations turns out to be almost independent of the volume fraction. In the range of Debye screening length and volume fraction investigated here, double layers around nanoparticles actually interact with each other. This strong interaction between nanocolloidal maghemite particles is probably responsible for the experimental dependence of the electrokinetic properties with the volume fraction.

I. Introduction

According to the DLVO theory [START_REF] Verwey | Theory of the Stability of Lyophobic Colloids: The Interactions of Sol Particles Having an Electric Double Layer[END_REF], the long-term stability of charged colloidal suspensions relies on the balance between attractive van der Waals interactions and repulsive electrostatic interactions. The magnitude of these electrostatic interactions is dictated by the charge of the colloidal particle and is strongly affected by the surrounding ionic environment, i.e. the nature and the concentration of the electrolyte. Therefore, understanding the static and dynamical properties of co-ions and counterions around the colloidal particles and their dependence with the medium of dispersion is a cornerstone to predict the long-term stability of charged colloidal suspensions.

For the vast majority of inorganic colloidal systems, the surface charge of the colloidal particle results from surface crystallographic defects (missing valencies at the surface of the particles), the state of the surface depending on the surrounding medium. In the case of oxide particles in contact with an aqueous solution, the surface charge results from proton exchange with the surface oxides and is therefore pH dependent. It can be evaluated from chemical titration of the weak acidities of the colloidal particle surface in contact with strong acidic or alkaline solutions. Depending on the authors, this charge is called the intrinsic, bare, structural, titrated or surface charge. We denote this quantity the structural charge in the following text.

Due to this structural charge, there is a local electric field at the surface of the colloidal particle. If the electrostatic interaction is much stronger than thermal energy, some counterions are strongly attracted to the surface, leading to an accumulation of counterions in the vicinity of the surface of the colloidal particle. This phenomenon is called counterion condensation [START_REF] Belloni | [END_REF]3]. If ions possess a strong chemical affinity to the surface, they can also bind to the surface sites. The structural charge of colloidal particle is partially compensated by these counterions, so that it bears an apparent charge, which we denote by effective charge in what follows. For a given system, the value of the effective charge depends on the method used to probe or to define it. For linear polyelectrolytes, the distinction between condensed ions and free ions has been clearly established by Manning. For colloids, defining a distance below which ions are considered as electrostatically condensed is not straightforward.

Bjerrum and others have made some propositions to define the effective charge based on the radial distribution functions between ions and colloidal particles at equilibrium [START_REF] Belloni | [END_REF]. In the following, we refer to this charge as the "static" effective charge, as it was already proposed in the literature [4,5]. It is also possible to define the effective charge as the charge involved in the linear Poisson-Boltzmann theory, and therefore in the DLVO theory. This DLVO charge ensures that the electric field far from the colloidal particles is equivalent to the field obtained using non-linear Poisson-Boltzmann equation.

The experimental determination of the effective charge is a non-trivial problem, which is still under debates in the community. Various techniques have been used to extract an effective charge. Osmotic pressure measurements yield an estimation of the number of free ions, and therefore are linked to the effective charge as defined by Bjerrum [6,7]. Other techniques enable to characterize the structure of the suspensions, such as X-rays or neutron scattering [8,[START_REF] Lyklema | Solid-Liquid Interfaces[END_REF].This structure depends on the long-range interactions between colloidal particles, and therefore on the DLVO effective charge. Finally, electrokinetic techniques, such as electrophoresis [START_REF] Lyklema | Solid-Liquid Interfaces[END_REF][START_REF] Masliyah | Electrokinetic and colloid transport phenomena[END_REF][START_REF] Delgado | [END_REF][12], measure an electrokinetic charge, which we call a "dynamic" effective charge in what follows, as already proposed in the literature [4,5]. As techniques and related theoretical treatments are different, the obtained results are seldom consistent.

Moreover, some standard results of the electrokinetic theories might not be adapted to the description of colloids of nanometric size. In particular, it is usually assumed that the main contribution to the electrical conductivity of the suspension is due to small ions both in the bulk and in the double layer [START_REF] Masliyah | Electrokinetic and colloid transport phenomena[END_REF]. Conversely, in classical descriptions of electrolyte solutions [START_REF] Robinson | Electrolyte Solutions[END_REF][START_REF] Bernard | [END_REF][15][16] every charge species contribute to the charge transport. For nanoparticles with a radius only 10 to 30 times that of small ions, the ratio between the mobilities at infinite dilution is between 1/30 and 1/10, so that neglecting the contribution of nanoparticles is questionable.

In the present work, we focus on colloidal particles of nanometric size. More precisely, we combine experimental results, analytical theory and numerical simulations to gain more insight into the concepts of static and dynamic effective charges. Other authors have addressed this issue by using numerical simulations: Lobaskin and coworkers used a combination of Lattice Boltzmann simulations and Langevin dynamics to compute effective charges for a model system [4]. They define the dynamic charge as the structural charge decreased by the number of counterions moving with almost the same velocity as the nanoparticle in the presence of the external field. They show that, for nanoparticles of radius equal to 3 nm in a dilute and salt-free solution, strong differences appear between the static and dynamic effective charges. In particular, both charges depend differently on the structural charge of the nanoparticle: While the static charge does almost not vary for increasing values of the structural charge larger than 40, the dynamic charge does not seem to saturate even with strong values of the structural charge (Z = 120).

The dependence of the different effective charges on the volume fraction of colloids is another issue that has been scarcely addressed [5,17]. From the theoretical perspective, it is quite difficult to include the effect of volume fraction beyond the cell model description [START_REF] Masliyah | Electrokinetic and colloid transport phenomena[END_REF],

for which the volume fraction dictates the size of the cells. Nevertheless, in the case where the Debye length is of the order of the size of the colloidal particles, double layers may strongly interact when the volume fraction increases, modifying the electrokinetic properties as compared to the infinite dilution.

The correlations between colloidal particles can be estimated either by using simulation methods, or other theoretical frameworks. We propose in this work to apply an analytical transport theory, which takes into account explicitly electrostatic and excluded volume interactions between every charged species. It previously enabled the accurate interpretation of the direct-current (DC) electrical conductivity of electrolytes or solutions of macroions like micellar species, as a function of their volume fraction. This approach is applied to a homemade and already widely studied system constituted of iron oxide nanoparticles (maghemite) dispersed in aqueous medium [START_REF] Berkovski | Magnetic Fluids and Applications Handbook[END_REF][START_REF] Lucas | [END_REF][20]. The maghemite particles have a mean radius of 6 nanometers. Because of their interesting magnetic properties [START_REF] Berkovski | Magnetic Fluids and Applications Handbook[END_REF][START_REF] Blums | Magnetic Fluids[END_REF][START_REF] Rosensweig | Ferrohydrodynamics[END_REF][START_REF] Odenbach | Ferrofluids : Magnetically Controllable Fluids and Their Applications[END_REF] , such dispersions, known as ferrofluids, are used in an increasing number of applications that strongly rely on their colloidal stability. In this work, we first determined the structural charge of the maghemite particles using chemical titration. Then, we performed dynamic measurements to estimate the dynamic effective charge. These latter were conducted at controlled volume fractions of maghemite nanoparticles and controlled pH. We measured the DC electrical conductivity and deduced the effective charge of nanoparticles by using an analytical transport theory, called the MSA (Mean Spherical Approximation)-transport theory. This analytical theory, which takes into account explicitly interactions between every species, proved indeed to be able to predict not only the electrical conductivity of simple electrolyte solutions [START_REF] Bernard | [END_REF][START_REF] Dufrêche | [END_REF], but also that of micellar systems [25,26]. The electrokinetic potential of maghemite nanoparticles, which is related to the dynamic effective charge, was also measured as a function of the volume fraction from the acoustophoresis technique. Third, the static effective charge was derived from Monte-Carlo (MC) simulations using the structural charge of the maghemite particle deduced from chemical titration as input parameter. The model used accounts for electrostatic and excluded volume interactions and neglects specific chemical binding of ions to the surface. This simulation method is particularly adapted to the determination of equilibrium properties of charged systems [27], in particular for high volume fractions [28].

After a brief recall of the experimental protocols (synthesis of the colloidal suspension, preparation of the samples and electrokinetic measurements) and of the theoretical approaches (MSA-transport theory and MC simulations) in the section II, the results obtained from the different approaches are reported and discussed in section III. As we proceed to show, a dependence of the dynamic effective charge with the volume fraction clearly appears, whereas the static effective charge as computed from MC simulations turns out to be almost independent of the volume fraction.

II. Methods

II.1. Experiments

Synthesis of maghemite colloidal suspensions: Maghemite (g-Fe 2 O 3 ) nanoparticles were synthesized by co-precipitation of iron(II) chloride and iron(III) chloride in alkaline medium followed by acidification with nitric acid and oxidation by iron (III) (ferric nitrate) at 100 ºC [START_REF] Berkovski | Magnetic Fluids and Applications Handbook[END_REF][29][30][31]. The resulting colloidal suspension strongly absorbs in the visible spectrum, appearing red and dark even at low volume fractions. The particle size determined from AFM measurements [START_REF] Lucas | [END_REF] follows a log-normal distribution, with a median diameter d 0 =12 nm (ln d 0 =<ln d>) and a polydispersity of s=0.28. The volume fraction f of the maghemite particles was evaluated from a chemical titration of iron [31] and from atomic absorption spectroscopy.

Determination of the particle's structural charge:

Maghemite nanoparticles through the pH-dependent protonation/deprotonation of their surface oxide groups turn either positively charged, uncharged at the point of zero charge (PZC) or negatively charged, according to the equilibria (1) and ( 2) where M is the metal of the spinel nanoparticles, M-OH 2 + the metaloxide acid site and M-O -the metal-oxide basic site.

(

1) M-OH 2 + = M-OH + H + and (2) M-OH = M-O -+ H +

The number of charged groups on the surface as a function of pH can be determined from titrations. This has been done here using the method already applied in [START_REF] Lucas | [END_REF]: the initial stable suspension at pH around 1.5 is neutralized by addition of a base up to pH = 7, which is the point of zero charge of maghemite [32][33][START_REF] Cornell | The Iron Oxides: Structure, Properties, Reactions, Occurences and Uses[END_REF][START_REF] Kosmulski | Chemical properties of material surfaces[END_REF], where flocculation is complete. The supernatant is removed and the particles are washed with distilled water until the conductivity reaches the one of water. This final state, where the particles are uncharged and the solvent is free of ions, is the starting point of the titrations performed by addition of nitric acid. Note that the particles can be perfectly redispersed in water as soon as their charge is high enough, typically for pH < 3.5. The charge is given by the difference between the number of added H + and the number of H + determined from the pH measurements. Two totally independent experiments are plotted in Figure 1 down to pH=2.3, showing a maximal difference of 10%. Note that, below pH = 2, the accuracy of the determination of the structural charge is poor and should be considered with care since the amount of free H + becomes large compared to the amount of H + on the surface of the particles. At pH=3.1 considered here, the determination is reliable and the structural charge equals Z b = 392. As the radius of the colloid is close to 6 nm, Z b corresponds to a significant surface charge density of 0.87 e nm -2 (0.14 C.m -2 ), consistent with other determinations in similar systems of iron oxides [START_REF] Campos | [END_REF].

Preparation of samples of controlled volume fraction and pH:

The acidic colloidal maghemite suspension (pH range 1.5-1.8) was packed into a dialysis tubing permeable to water molecules and ions (Spectrapore-12000-14000) and immersed into an acidic solution (HNO 3 , pH 2.7) containing 5% in weight of Dextran polymer (500 kg mol -1 Amersham Biosciences).

Even after several days, an asymmetric distribution of protons from both sides of the dialysis tube membrane (pH=3.1 inside vs pH=2.7 outside) is still observed and accounts for the presence of charged nanoparticles inside the tube because of the Donnan effect [37,38].

Successive dilutions of the resulting colloidal suspension at constant pH 3.1 with a nitric acid solution provided a series of samples with volume fractions f between 5.8 % and 0.0001 %.

DC conductivity measurements:

The conductivity of the series of samples prepared in nitric acid HNO 3 at pH = 3.1 was measured without any other additional salt. Conductivity measurements were conducted on each sample by using a Fuoss-Kraus conductivity cell immersed into a water bath with controlled temperature (T = 25±0.1°C). The cell constant, 0.1785 cm -1 , was determined using the method of Lind, Zwolenik, and Fuoss [39] with high purity (99.99 %) KCl aqueous solutions. The impedance of the solution was measured with a high accuracy Wayne-Kerr bridge (model 6425 A) at four frequencies (1, 2, 5 and 10 kHz), and the real part was extrapolated to infinite frequency to separate the contribution of the electrode polarisation and determine the resistance of the solution [40].

Electrokinetic potentials measurements:

The electrokinetic potential of maghemite nanoparticles or zeta-potential in the same series of samples was measured by using the acoustophoresis technique. This technique detects the electrical potential or current generated when a suspension containing charged species is exposed to an ultrasonic wave. These quantities depend both on the density and on the charge of the solutes. We used a DT 1200 (Dispersion Technologies, USA) electro-acoustic spectrometer. This instrument measures the induced current (amplitude and phase) named CVI (Colloidal Vibration Intensity) when an ultrasonic wave is applied, and extracts the dynamic mobility of the colloidal particles.

Calibration of the instrument with a suspension of Ludox ® (silica particles, TM-50 Sigma

Aldrich, w = 10 % in KCl at 0.01mol L -1 and pH=9, z = -38 mV) is required. The electrokinetic potential is deduced from the measured CVI using the extended Helmholtz-Smoluchowski equation [START_REF] Delgado | [END_REF] and several input data for maghemite particles and solvent, i.e.

volume fraction, density, sound speed and dielectric constant. Unfortunately, no access to the raw data is possible from this device, so that the validity of the approximations involved in the treatment of the measures cannot be checked.

II.2. Theoretical approaches

MSA-transport theory:

The DC electrical conductivity of a solution is a non-trivial function of the diffusion coefficients of ions at infinite dilution, of the charge of ions, of their concentration and of the interactions between charged species. The MSA-transport theory [START_REF] Dufrêche | [END_REF][25][26] provides an analytical expression of the conductivity, taking into account the deviations from the ideal behavior, namely the interactions between solutes. In the framework of this transport theory, which describes the solvent as a continuous medium characterized by its dielectric constant and its viscosity, two corrective terms are calculated [41]. The first one accounts for the electrostatic relaxation of the ionic atmosphere around charged species after an external force is applied on the solution (due to the electric field here). The second one accounts for the hydrodynamic couplings mediated by the solvent that tend to equalize the ion drift velocities (electrophoretic phenomena). This electrophoretic effect is generally greater than the first one. The explicit formulation of those corrections is given elsewhere [25,26].

Another interest of the MSA-transport theory is that it allows one to compute the contribution of each species i to the global conductivity.

Several parameters are needed to compute the relaxation and electrophoretic corrections: the self-diffusion coefficient at infinite dilution of charged solutes (denoted by D°), their charge, and the minimal distance of approach between solutes, which can be reduced to an individual radius for each species (denoted by R). For small ions (H + , NO the effective charge of the nanoparticle moving with some of its counterions under the influence of the external electric field.

Monte Carlo simulations: Monte Carlo numerical simulations allow one to compute the average thermodynamical properties of a model system included in a simulation box of length L box with periodic boundary conditions. Here, we were particularly interested in the radial distribution functions between charged solutes in the suspension, with the solvent treated as a continuum. The MC technique is particularly adapted to the simulation of solutions of charged nanoparticles in the presence of small ions [START_REF] Wu | [END_REF][44][45][46]. If one focuses on the pair distribution function g Cj (r) between nanoparticles and ions of type j, the number N Cj of ions j situated at a distance r from the center of a nanoparticle can be easily determined using equation ( 5):

(

) 5
where n j is the bulk concentration of ions j (n j = N j / L box 3 with N j the number of ions j in the simulation box). From this quantity, the static effective charge Z eff-static of the nanoparticle can be deduced at a given distance from the nanoparticle (see below).

The electrostatic potential at a given distance r from a nanoparticle, which is related to the electrokinetic potential, was also computed from these simulations, by solving Poisson's equation within the spherical geometry around the nanoparticle [47]:

1 / 2 2 1 2 2 2 2 1 1 ( ) ( ) C r Lbox
where ε = ε 0 ε r is the permittivity of the medium, R C the radius of the nanoparticle, r e the electrical density of charge around the nanoparticle, and e the elementary charge. In this formula, the electrical density of charge is deduced from the radial distribution functions and includes also other nanoparticles.

Monte Carlo (MC) calculations were carried out at a constant temperature of 25°C to simulate four different volume fractions of nanoparticles surrounded by their co-ions and counterions: dominate over other interactions. The radius of nanoparticles is equal to 6 nm, and that of small ions (nitrate ions and protons) is taken equal to 0.25 nm, as it is usually done to describe simple electrolyte solutions [44]. The proton concentration corresponds to pH=3.1 as established in the experimental protocol. Each nanoparticle has a constant structural charge Z b =392, as found from titration results. The simulation box contains 2 or 3 nanoparticles. One million steps per particle were achieved to ensure the convergence of the radial distribution functions g Cj (r) [44]. The parameters of the simulations are given in Table 1.

f (%) L box (nm) N nanoparticles N NO3- N H+ 0.

III. Results

The DC electrical conductivity of the maghemite suspension as a function of the volume fraction of nanoparticles at constant pH (pH=3.1) is presented on Figure 2 (black squares).

We have used the MSA-transport theory to deduce from these experimental results the dynamic effective charge Z eff-dyn of maghemite nanoparticles (fitting procedure explained in Section II.2). It was not possible to account for the experiments while keeping a constant value of Z eff-dyn as a function of the volume fraction. We have thus repeated the fitting procedure for each volume fraction, so that we obtained the dynamic effective charge as a function of the volume fraction. These theoretical results are also plotted on Figure 2.

Moreover, the contributions of each charged species to the global conductivity were computed and are plotted as a function of the volume fraction in maghemite on Figure 3. The electrokinetic potential of maghemite nanoparticles in the same samples, measured by acoustophoresis for volume fractions larger than 0.5 %, is plotted on Figure 4.

This first set of results calls for several comments. First, at this specific pH, the contribution to the electrical conductivity of protons is of the order of 0.03 S.m -1 , as it can be easily evaluated from the individual electrical conductivity at infinite dilution. Their contribution to conductivity is constant with volume fraction. As shown by the MSA-transport calculations presented in Figure 3, the contributions of maghemite nanoparticles and nitrate ions prevail over that of protons for volume fractions higher than 2 %. Second, the effective charge of nanoparticles as deduced from MSA-transport calculations is much smaller than the structural charge, and is found to decrease with the volume fraction: Z eff-dyn varies between 288 at f=0.4%, and 116 at f=5.6%, whereas we found for the structural charge Z b =392 from chemical titrations. As already stated, the effective dynamic charge as deduced from MSAcalculations accounts for the number of counterions which are sufficiently attracted by a nanoparticle to move with it under the application of an external electric field, decreasing the mobility of the nanoparticle and thus the electrical conductivity of the system. These results are corroborated by the values of the electrokinetic potential, which is found to decrease by a factor two between 0.05 % and 5.8 % (acoustophoresis measurements). The electrokinetic potential, which is the potential difference between the stationary layer of fluid attached to the nanoparticle and the value at long distance, i.e. far from any nanoparticle, depends also on the number of ions situated in the vicinity of the nanoparticle and moving with it. It appears then both from the conductivity and the electrokinetic potential measurements that the number of counterions dynamically coupled to the nanoparticle increases with the volume fraction.

Monte Carlo simulations allow us to probe the number of counterions statically coupled to the nanoparticle according to the Bjerrum criterion. Figure 5 displays the number of counterions at a given distance from the center of a nanoparticle for several volume fractions computed by simulation. We recall that in the model used in simulations, the nanoparticle bears a charge equal to the structural charge of the maghemite particles, and that solute species are treated as charged hard spheres. At a given distance r, one can define an effective charge equal to the structural charge minus the number of counterions N Cj . Several criteria can be used for the determination of the threshold distance where counterions are no longer in strong interaction with the particle [START_REF] Belloni | [END_REF]44]. This optimal threshold distance is often chosen as the point where the first derivative of N Cj (r) is minimum (inflexion point). In the plots displayed in Fig. 5, this corresponds to distances equal to 10.4, 13.4, 15.0 and 16.1 nm from the center of the nanoparticle for volume fractions of f = 5.0, 1.5, 0.5 and 0.1% respectively. At these distances from the center of the nanoparticle the effective charge is found almost constant for the four systems, equal to Z eff-static = 50. Note that this value is consistent with that defined by Belloni in a colloidal suspension [START_REF] Belloni | [END_REF], Z eff-static = 4 r eff / L B with L B the Bjerrum length (we have Z eff-static = 58 for r eff = 10.4 nm). Nevertheless, these results differ from those obtained previously for the dynamic effective charge in two aspects: i) the value of the static effective charge is much lower than those of the dynamic effective charge, ii) the effective charge is found almost independent of the volume fraction.

The fact that the effective static charge is found to be lower than the dynamic effective charge can be relatively easily understood. Some of the counterions, assumed to be strongly electrostatically coupled to the nanoparticle according to the Bjerrum criterion, actually are not sufficiently "bound" to move with the nanoparticle in the presence of an external field, so that the dynamic effective charge estimated by the MSA-transport theory is larger than the static one. This result was already obtained by Lobaskin et al. [4] from theoretical investigations. The variation of the dynamic effective charge and of the electrokinetic potential with the volume fraction of nanoparticles is not easy to interpret. In standard electrokinetic theories (Helmholtz-Smoluchowki limit) where the Debye layer is very thin, electrokinetic properties of colloids are usually assumed to be independent on the volume fraction. In our case, particles are of nanometric size and in solution with low concentration of added salt (NO 3 -and H + ions in the diluting HNO 3 solution of pH 3.1). To evaluate the Debye length, we can assume either that only the ions of the nitric acid in the bulk participate to the electrostatic screening, or, more realistically, that both counterions of the nanoparticles and ions of the nitric acid in the bulk participate to the screening. In the first case, we find a Debye length k -1 equal to 10.7 nm, whereas in the second case, k -1 varies between 1.6 nm (at f=5%) and 6.3 nm (at f=0.1 %). This means that in every case, and whatever the choice to evaluate the Debye length, we have kR close to one, with R the radius of the maghemite nanoparticle. In other words, we cannot neglect the size of the double layer compared to that of the nanoparticle. Monte Carlo simulation results bring more insights into this issue. We have plotted in Fig. 6 the electrostatic potential created by a nanoparticle and charged species as a function of the distance to the center of the nanoparticle. The electric potential presents in every case a minimum, because several nanoparticles are included in the simulation box. We have translated the computed values so that the minimum of the electric potential corresponds to zero. For the highest volume fraction investigated here, double layers around nanoparticles interact with each other, as attested by the fact that the computed electric potential does not provide any plateau in between nanoparticles. The electric potential is then strongly influenced by the presence of other nanoparticles. This strong interaction between nanocolloidal maghemite particles is probably responsible for the experimental dependence of the electrokinetic phenomena with the volume fraction.

IV. Discussion and Conclusion

In this study we exposed several complementary experimental and theoretical tools to characterize the ionic condensation on well-defined maghemite nanoparticles in solution.

More precisely, the combination of DC conductivity measurements, acoustophoresis, analytical theory, and numerical simulations allows us to investigate the influence of the volume fraction of nanoparticles on their charge, a key parameter which helps understanding the stability of the suspensions.

Our results show the following trends for such nanocolloidal systems:

(1) The effective charge deduced from dynamic measurements, called here the dynamic effective charge or the electrokinetic charge, is different from the static effective charge deduced from the structure of the ionic cloud at equilibrium. It is indeed found much larger than the static effective charge, which indicates that most of the electrostatically condensed ions (from Bjerrum's definition of condensation) are not fixed relative to the nanoparticle surface and therefore contribute to the electrical current. This means that in the present case, the precise definition of the position of the slip plane where the electrokinetic potential is defined is not straigthforward [START_REF] Delgado | [END_REF].

(2) The dynamic effective charge appears to decrease while the volume fraction of the maghemite nanoparticles increases. This result, determined from conductivity measurements treated within the Mean Spherical Approximation-transport theory, is also confirmed by electroacoustic measurements. It is not usual to observe such a dependence of the properties of the colloidal particles on the volume fraction. This might occur in our systems because it stands at the limits of the colloidal domain. The small size of the colloidal particle together with its important surface charge have indeed several consequences: (i) the concentration of counterions becomes very important even for low nanoparticles volume fractions, and is therefore often larger than the concentration of the added salt, (ii) the average distance between individual particles is much smaller than in the case of micrometer colloidal particles for the same volume fraction, so that at a given Debye screening length, increasing the volume fraction leads to strong interactions between the double layers of neighboring nanoparticles, and (iii) charged nanoparticles participate to the charge current within the solution, and therefore influence the displacements of other nanoparticles. As a consequence, the properties of the solution should be influenced by the volume fraction of colloids. For all these reasons, our choice to resort to a model of the solution that explicitly considers the effect of all charged species (salt ions, charged nanoparticles and their counterions) was necessary.

(3) The static effective charge does not vary with the volume fraction, which shows again that the electrokinetic charges and the static effective charge are not straightforwardly related quantities.

It is generally accepted that the distribution of ions around colloidal particles cannot be simply described by a succession of layers with clear-cut boundaries and well-defined properties [48]. Our results suggest that, for the present system, condensed counterions fall into the following two categories: electrostatically condensed counterions which move with the nanoparticles and do not contribute to the electrical current (those leading to Z eff-dyn ), and the electrostatically condensed ions which participate to the electrical conduction (those leading to Z eff-stat ). While the total number of condensed counterions per particle almost does not vary with the volume fraction of nanoparticles, the populations of each category of ions does.

) 3 -) = 1 . 1 (

 311 , these parameters were taken from the literature[START_REF] Marcus | Ion properties[END_REF]: D°(NO 902 10 -9 m 2 s -1 , D°(H + ) = 9.310 10 -9 m 2 s -deduced from the individual electrical conductivity at infinite dilution of H + ), R(NO 3 -) = 0.19 nm, R(H + ) = 0.095 nm (which corresponds to the size of a water molecule)[START_REF] Robinson | Electrolyte Solutions[END_REF]. We checked that these parameters allowed us to predict the conductivity of dilute HNO 3 solutions with a good precision. The diffusion coefficient at infinite dilution of maghemite particles is calculated thanks to the Stokes-Einstein relation from the radius of the particle (6 nm): D°(maghemite) = k B T / (6phR) with k B the Botzmann constant, T the temperature and h the viscosity of water at ambient temperature. The only unknown parameter is thus the charge of maghemite nanoparticles, taken as an adjustable parameter to model the experimental conductivity. This charge accounts for interactions at long range between moving nanoparticles and is thus an effective dynamic charge, denoted by Z effdyn in what follows, i.e.

  1, 0.5, 1.5 and 5%. Solute particles are modeled as charged hard spheres (primitive model) because long-ranged electrostatic interactions between solutes are assumed to

Figure 1 :

 1 Figure 1: Structural charge Z b of maghemite nanoparticles determined from chemical

Figure 2 :

 2 Figure 2: Experimental electrical conductivity of samples of maghemite nanoparticles

Figure 3 :

 3 Figure 3: Contribution to the electrical conductivity of the different charged species as a

Figure 4 :

 4 Figure 4: Electrokinetic potential or zeta-potential (z) of maghemite nanoparticles prepared

Figure 5 :

 5 Figure 5: Number of NO 3 -counterions as a function of the distance to the center of

Figure 6 :

 6 Figure 6: Electrical potential created by the nanoparticle and small ions as a function of the

Table 1 :

 1

	1	121.86	2	1649	865
	0.5	81.577	3	1436	260
	1.5	56.562	3	1263	87
	5	37.865	3	1202	26

Parameters of Monte Carlo simulations. f (%) volume fraction of nanoparticles, L box the length of the cubic simulation box, N i the number of particles i in the simulation box.


Paolo Robuffo Giordano

Riccardo Spica
email: riccardo.spica@irisa.fr

Franc ¸ois Chaumette
email: francois.chaumette@irisa.fr

F Chaumette

Learning the Shape of Image Moments for Optimal 3D Structure Estimation

The selection of a suitable set of visual features for an optimal performance of closed-loop visual control or Structure from Motion (SfM) schemes is still an open problem in the visual servoing community. For instance, when considering integral region-based features such as image moments, only heuristic, partial, or local results are currently available for guiding the selection of an appropriate moment set. The goal of this paper is to propose a novel learning strategy able to automatically optimize online the shape of a given class of image moments as a function of the observed scene for improving the SfM performance in estimating the scene structure. As case study, the problem of recovering the (unknown) 3D parameters of a planar scene from measured moments and known camera motion is considered. The reported simulation results fully confirm the soundness of the approach and its superior performance over more consolidated solutions in increasing the information gain during the estimation task.

I. INTRODUCTION

The quest for finding a good set of features for visual control and 3D Structure from Motion (SfM) is a classical problem in the visual servoing community, and it has indeed attracted a large body of literature over the last decades. A number of approaches has been proposed over the years for exploiting local geometrical primitives, such as points or lines tracked on the image, as visual features to be measured/controlled. A comprehensive overview of these possibilities can be found in [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. In parallel, another very successful line of research has considered the use of more 'integral' image descriptors able to encode the information contained over a region of interest on the image plane (e.g., enclosed by the contour of a tracked object). This indeed usually results in a (relatively) easier extraction, matching and spatio-temporal tracking across multiple frames of the region of interest, thus generally improving the robustness against image processing errors.

Image moments of binary dense closed regions or of discrete sets of points [START_REF] Tahri | Point-Based and Region-Based Image Moments for Visual Servoing of Planar Objects[END_REF] are a typical (and by now classical) example of integral features exploited for visually controlling the camera pose [START_REF] Hoffmann | Visual servoing with moments of SIFT features[END_REF]- [START_REF] Copot | Predictive control of nonlinear visual servoing systems using image moments[END_REF]. Further extensions of these ideas have dealt with, e.g., the use of particular kernels for evaluating the image moments [START_REF] Kallem | Kernel-based visual servoing[END_REF], or the direct use of pixel intensities [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF] by processing the whole acquired image. Image moments from dense regions [START_REF] Giordano | 3D structure identification from image moments[END_REF] or discrete point clouds [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF], [START_REF] Spica | Experiments of plane estimation by active vision from point features and image moments[END_REF] have also been exploited as visual features for recovering the 3D structure of a planar scene via SfM schemes.

Despite this state-of-the-art, it is however worth noting that the selection of a good set of image moments for 6-dof visual control or SfM is still an open problem. Ideally, one would like to find a unique set of visual features resulting in the 'most linear' control problem with the largest convergence domain, or in maximum observability (i.e., information gain) for a given camera displacement in case of SfM tasks. However, to the best of our knowledge, only local, partial (e.g., depending of the particular shape of the object) or heuristic results are currently available. For instance, [START_REF] Tahri | Point-Based and Region-Based Image Moments for Visual Servoing of Planar Objects[END_REF], [START_REF] Kallem | Kernel-based visual servoing[END_REF], [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF] propose different combinations of image moments able to only guarantee local 6-dof stability of the servoing loop around the desired pose, and with a basin of attraction to be heuristically determined case by case. As for what concerns the SfM case, the choice of which moments to exploit for allowing a converging estimation of the scene structure is also not straightforward. In [START_REF] Giordano | 3D structure identification from image moments[END_REF], [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF] the area a and barycenter coordinates (x g , y g ) of a dense region are successfully fed to a SfM scheme based on the (intuitive) motivation that the same set (a, x g , y g ) is also the typical choice for controlling the camera translational motion in a servoing loop [START_REF] Tahri | Point-Based and Region-Based Image Moments for Visual Servoing of Planar Objects[END_REF]. However, this intuition breaks down when considering moments of a discrete point cloud: in this case, the typical choice for controlling the camera translational motion, that is, the set (x g , y g , µ 20 + µ 02 ) (see [START_REF] Tahri | Point-Based and Region-Based Image Moments for Visual Servoing of Planar Objects[END_REF]), is empirically shown in [START_REF] Spica | Experiments of plane estimation by active vision from point features and image moments[END_REF] to not provide enough information for allowing a converging estimation of the scene structure.

One could argue that the hope of finding a unique set of visual features optimal in all situations might eventually prove to be unrealistic, if not impossible, while it could just be more appropriate (and reasonable) to rely on an automatic and online selection of the best feature set (within a given class) tailored to the particular task at hand. Motivated by these considerations, the goal of this work is to propose an automatic learning strategy able to select online the 'best' set of image moments in order to optimize the SfM performance in recovering the (unknown) 3D parameters of a planar scene.

The rest of the paper is organized as follows: first some preliminary concepts of interests are reviewed in Sect. II, and then the definition of weighted parametric moment, central for the paper developments, is introduced in Sect. III. Section IV discusses a possible strategy for automatically adjusting the shape of the weighted moments so as to maximize an 'observability measure' (or information gain) of the chosen SfM task, and Sect. V reports the results of several simulations where the benefits of the proposed approach can be clearly appreciated. Finally, Sect. VI concludes the paper and draws some future perspectives.

II. PRELIMINARIES

Consider a planar scene with equation P : n T E + d = 0, with n ∈ S 2 being the unit normal vector and d ∈ R the plane distance from the camera center, and let (v, ω) ∈ R 6 represent the linear/angular camera velocity 1 . The (i, j)-th image moment m ij of a collection of N point features p k = (x k , y k , 1) tracked on P is defined as

m ij = N k=1 x i k y j k . (1)
Furthermore, x g = m 10 /N and y g = m 01 /N denote the barycenter coordinates and

µ ij = N k=1 (x k -x g ) i (y k -y g ) j
the (i, j)-th centered moment. By letting χ = -n/d ∈ R 3 , the dynamics of m ij is known to take the following expression linear in

χ ṁij = f ωij (m kl , ω) + f χij (m kl , v)χ, (2)
see [START_REF] Tahri | Point-Based and Region-Based Image Moments for Visual Servoing of Planar Objects[END_REF] for further details. In (2), m kl stands for a generic (k, l)-th moment of order up to i + j + 1 and χ represents the 3D structure of the observed planar scene P (not directly measurable from sole image quantities). A conceptually equivalent derivation can also be obtained for the barycenter, for the centered moments, and for the case of image moments of dense planar regions [START_REF] Tahri | Point-Based and Region-Based Image Moments for Visual Servoing of Planar Objects[END_REF]. Equation ( 2) is at the core of virtually all algorithms for 3D structure estimation with image moments playing the role of measurements. For instance, owing to the linearity of (2) w.r.t. χ, several SfM schemes meant to recover the (unmeasurable) 3D structure χ from the measured s(t) and the known camera motion (v, ω) have been proposed in [START_REF] Giordano | 3D structure identification from image moments[END_REF]- [START_REF] Spica | Experiments of plane estimation by active vision from point features and image moments[END_REF]. However, as discussed in the Introduction, a satisfactory SfM performance requires the preliminary identification of a suitable set of image moments 'rich enough' for ensuring observability of the scene structure. Specifically, given a collection of m measured moments s = (m i1j1 , . . . , m imjm ) and defining (using [START_REF] Tahri | Point-Based and Region-Based Image Moments for Visual Servoing of Planar Objects[END_REF])

Ω = f T χi 1 j 1 (m kl , v) . . . f T χi njn (m kl , v) ∈ R 3×m , (3)
the 'observability matrix' ΩΩ T ∈ R 3×3 must remain full rank during the camera motion [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF].

The weighted parametric moments introduced hereafter are meant to provide an adjustable visual feature set that can be tuned online for automatically coping with this observability requirement.

III. DEFINITION OF WEIGHTED PARAMETRIC

IMAGE MOMENTS Let w = w(x, y, θ) be a smooth function of the coordinates (x, y) on the image plane and of a vector of parameters θ ∈ R p . One can generalize (1) and define a weighted parametric image moment for N observed features p k

m w (θ) = N k=1 w(x k , y k , θ), (4)
with, obviously, m w = m ij for w(x, y, θ) = x i y j . Function w(x, y, θ) can be seen as the class of all the considered image moments (e.g., a quadratic form in x, y) parameterized by vector θ (e.g., the coefficients of the quadratic form). Consider now the following additional definitions

m x wij (θ) = N k=1 x i k y j k ∂w(x, y, θ) ∂x (xk, yk) m y wij (θ) = N k=1 x i k y j k ∂w(x, y, θ) ∂y (xk, yk) m θ w (θ) = N k=1 ∂w(x, y, θ) ∂θ (xk, yk) , (5)
and note that m θ w is a row vector of dimension p. Following the derivations in [START_REF] Tahri | Point-Based and Region-Based Image Moments for Visual Servoing of Planar Objects[END_REF], it is easy to show that the dynamics of m w (θ) takes the expression (reminiscent of (2))

ṁw (θ) =[m A (v, θ) m B (v, θ) m C (v, θ)]χ + [m ωx (θ) m ωy (θ) m ωz (θ)]ω + m θ w (θ) θ (6)
with

                   m A = -m x w10 v x -m y w10 v y + (m x w20 + m y w11 )v z m B = -m x w01 v x -m y w01 v y + (m x w11 + m y w02 )v z m C = -m x w00 v x -m y w00 v y + (m x w10 + m y w01 )v z m ωx =(m x w11 + m y w02 + m y w00 ) m ωy =(-m x w00 -m x w20 -m y w11 ) m ωz =(m x w01 -m y w10 ) , (7)
and χ = (A, B, C) = -n/d. One can then exploit (4-7) for implementing a visual control or estimation algorithm as in the classical case, but with the additional possibility of acting on vector θ (a free parameter) for optimizing any criterium of interest (e.g., the norm of the observability matrix ΩΩ T during an estimation task as it will be discussed in Sect. IV).

A. Design of the weighting function w(x, y, θ)

Clearly, there exist many possibilities for designing the weighting function w(•), i.e., the class of moments spanned by vector θ. A convenient choice, in our opinion, is to take w(•) as some polynomial basis in x and y with θ being the vector of coefficients. Indeed, in this way the weighted moments (4), the expressions in (5) and, eventually, all the terms in (7) will reduce to linear combinations of the unweighted moments m ij in [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. The overall computational complexity will then result equivalent to the classical case [START_REF] Tahri | Point-Based and Region-Based Image Moments for Visual Servoing of Planar Objects[END_REF].

As for which polynomial basis to exploit, many choices are possible depending on the constraints/requirements of the particular application. Within the scope of this work, two possibilities are considered:

1) Polynomial basis of fixed degree: first, one can take w(•) as a polynomial in x and y of a given degree δ ∈ N + , that is,

w(x, y, θ) = δ j=1 j k=0 θ Tj +k x (j-k) y k (8)
with T j = j+1 2 and θ = (θ 1 , . . . . . . θ T δ +δ ) ∈ R T δ +δ . Indeed, this allows (4) to span all the moment linear combinations of order up to δ with coefficients in vector θ. As illustration, by choosing δ = 2 in (8), one obtains the following quadratic polynomial

w(x, y, θ) = θ 1 x + θ 2 y + θ 3 x 2 + θ 4 xy + θ 5 y 2
that, when plugged in (4), yields

m w (θ) = θ 1 m 10 + θ 2 m 01 + θ 3 m 20 + θ 4 m 11 + θ 5 m 02 . (9)
The class ( 9) can then specialize into, e.g., the barycenter coordinate x g for θ = (1/N, 0, 0, 0, 0), the centered moment µ 02 for θ = (0, -y g , 0, 0, 1), and so on. Clearly, the larger the value of the degree δ, the richer the basis representation power in encoding the scene geometry, but at the (well-known) cost of an increasing noise level with the moment order.

2) Constrained polynomial basis: a second possibility is to design a constrained polynomial basis for coping with the possible loss/gain of point features during the camera motion because of the limited camera field of view (fov). Indeed, by imposing that w(•) vanishes (with vanishing derivative) at the image borders, any point feature close to the limits will smoothly enter or leave the image plane and, thus, prevent any discontinuity in the moment dynamics [START_REF] Kallem | Kernel-based visual servoing[END_REF].

Let then x min < x max and y min < y max represent the limits of a rectangular image plane, and consider a weighting function w(•) partitioned as

w(x, y, θ) = w x (x, θ x )w y (y, θ y ), (10)
where w x (x, θ x ) and w x (x, θ y ) are polynomial bases and θ = (θ x , θ y ) ∈ R px+py , p x + p y = p, is the vector of coefficients. Assuming p x ≥ 4 and imposing

     w x (x min , θ x ) = w x (x max , θ x ) = 0 ∂w x (x, θ x ) ∂x xmin = ∂w x (x, θ x ) ∂x xmax = 0 , (11)
one can solve for a set of 4 parameters in vector θ x for shaping w x (x, θ x ) as desired. For instance, by taking

w x (x, θ x ) = θ x 1 x 5 + θ x 2 x 4 + θ x 3 x 3 + θ x 4 x 2 + θ x 5 x + θ x 6 (12
) and by (arbitrarily) choosing the pair (θ x 1 , θ x 2 ) as free parameters in vector θ x , system (11) yields Imposing analogous conditions to function w y (y, θ y ) at y min and y max (with again p y ≥ 4) will then constrain a total of 8 parameters in vector θ, with the remaining p -8 coefficients still free to be exploited for optimization purposes. For the sake of illustration, Fig. 1 shows three examples of weighting functions w(•) smoothly vanishing at the borders of an image plane of size [-1, 1] × [-1, 1] and obtained by picking at random three values for the free parameters in vector θ.

                     θ x 3 = (-3x 2 max -3x 2 min -4x max x min )θ x 1 +(-2x min -2x max )θ x 2 θ x 4 = (2x 3 max + 8x 2 max x min + 8x max x 2 min + 2x 3 min )θ x 1 +(x 2 max + 4x max x min + x 2 min )θ x 2 θ x 5 = (-7x 2 max x 2 min -4x 3 min x max -4x 3 max x min )θ x 1 +(-2x max x 2 min -2x 2 max x min )θ x 2 θ x 6 = (2x 3 max x 2 min + 2x 2 max x 3 min )θ x 1 + x 2 max x 2 min θ x 2 . ( 13
)
We conclude by noting that, compared to the previous case [START_REF] Giordano | 3D structure identification from image moments[END_REF], this latter possibility necessitates of a polynomial basis [START_REF] Spica | Experiments of plane estimation by active vision from point features and image moments[END_REF] with a degree of at least 7. Indeed, as explained, the vanishing conditions at the image border will constrain 4 + 4 coefficients in θ x and θ y , thus forcing both w x (x, θ x ) and w y (y, θ y ) to have a degree of (at least) 3 for ensuring p x ≥ 4 and p y ≥ 4 as required (see [START_REF] Spica | Active Structure from Motion: Application to Point, Sphere and Cylinder[END_REF]). However, for any optimization of the coefficient vector θ to be possible, either p x > 4 or p y > 4 must hold for allowing presence of at least one free coefficient to be optimized besides those already constrained by the vanishing conditions. On the other hand, if either p x > 4 or p y > 4, the final polynomial basis [START_REF] Spica | Experiments of plane estimation by active vision from point features and image moments[END_REF] will necessarily result of at least degree 7.

Therefore, the use of higher-order moments (of at least order 7) is the 'price to pay' for smoothly taking into account the loss/gain of point features during the camera motion 2 . In contrast, the degree of the polynomial basis in [START_REF] Giordano | 3D structure identification from image moments[END_REF] can be chosen at will and thus adjusted, if necessary, for limiting the noise level in the measured moments.

IV. OPTIMIZATION OF THE WEIGHTED PARAMETRIC IMAGE MOMENTS

Let again s ∈ R m represent the set of measured image moments. As explained in Sect. II (and discussed in more detail in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF]), when attempting to estimate the 3D structure χ from the measured s (and the known camera motion (v, ω)), full rankness of the square matrix ΩΩ T from (3) plays the role of a necessary (observability) condition for ensuring a converging estimation. Since in all SfM problems one has Ω = Ω(s, v) (see again [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF]), it is possible to act on the camera linear velocity v in order to maximize some conditioning measure of ΩΩ T for increasing the information gain during the camera motion. This insight has motivated some recent work in the context of active Structure from Motion for planar and 3D scenes [START_REF] Spica | Experiments of plane estimation by active vision from point features and image moments[END_REF]- [START_REF] Spica | Active Structure from Motion: Application to Point, Sphere and Cylinder[END_REF].

Consider now a set of m ≥ 3 weighted moments s = (m w (θ 1 ), . . . , m w (θ m )) ∈ R m with θ = (θ 1 , . . . , θ m ) ∈ R p being the stack of all parameters. Plugging the weighted moment dynamics (6-7) in the definition (3), one has

Ω(s, v, θ) =   m A (s, v, θ 1 ) • • • m A (s, v, θ m ) m B (s, v, θ 1 ) • • • m B (s, v, θ m ) m C (s, v, θ 1 ) • • • m C (s, v, θ m )   ∈ R 3×m .
(14) Therefore, when employing the weighted parametric moments (4) instead of the classical moments (1), one gains the additional possibility of also acting on vector θ (i.e., on the 'moment shape') for affecting matrix ΩΩ T .

Different scalar quantities can be taken as a measure of the conditioning of the square (and semi-positive definite) matrix ΩΩ T . For instance, the analysis in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF] shows that its smallest eigenvalue σ 2 1 directly affects the convergence rate of the employed estimator, and thus a reasonable choice is to maximize σ 2 1 over time. However, the evaluation of the derivative/gradient of an eigenvalue is unfortunately not well-defined for repeated eigenvalues [START_REF] Friswell | The derivatives of repeated eigenvalues and their associated eigenvectors[END_REF]. In order to avoid this issue, in this work we chose to take the quantity

ρ = det(ΩΩ T ) (15)
as a conditioning measure for matrix ΩΩ T . Indeed, from classical linear algebra [START_REF] Bernstein | Matrix Mathematics: Theory, Facts, and Formulas[END_REF] the following relationship holds for a square matrix A

d det(A) dt = tr adj(A) d(A) dt (16)
with tr(•) and adj(•) being the trace and adjugate operators, respectively. Contrarily to the derivative of an eigenvalue, the relationship ( 16) is always well-defined with, in particular, no possible ill-conditioning due to repeated eigenvalues. By then applying (16) to matrix ΩΩ T and expanding the various terms, one obtains

ρ = i tr adj ΩΩ T ∂(ΩΩ T ) ∂v i vi + i tr adj ΩΩ T ∂(ΩΩ T ) ∂θ i θi + i tr adj ΩΩ T ∂(ΩΩ T ) ∂s i ṡi = J v v + J θ θ + J s ṡ
(17) where the Jacobian matrixes

J v = . . . tr adj ΩΩ T ∂(ΩΩ T ) ∂v i . . . ∈ R 1×3 J θ = . . . tr adj ΩΩ T ∂(ΩΩ T ) ∂θ i . . . ∈ R 1×p J s = . . . tr adj ΩΩ T ∂(ΩΩ T ) ∂s i . . . ∈ R 1×m (18)
are function of (s, v, θ) (all available quantities). We stress that all the terms in (18) can be computed in closed-form.

The relation (17) can then be exploited for affecting ρ(t) over time by acting on v (the camera linear acceleration) and/or θ (the parameter vector). Among the many possibilities, we considered here the following update rules

           v = k v I - vv T v T v J T v θ = k θ I - θθ T θ T θ J T θ , k v > 0, k θ > 0 (19)
which are meant to maximize ρ(t) by following its gradient w.r.t. (v, θ) projected on the null-spaces of the constantnorm constraints v(t) = const and θ(t) = const. As explained in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF], the constraint v(t) = const is meant to prevent a better conditioning of matrix ΩΩ T only due to a faster camera motion while observing the scene. Indeed, since, roughly speaking, ΩΩ T is monotonically increasing with v 2 , the faster the camera motion the faster the SfM convergence regardless of any other optimization action. The second constraint θ(t) = const is motivated by similar arguments: an increasing θ(t) would (artificially) magnify ρ(t) at the cost of an increased noise level (all the terms in ( 7) would just result amplified).

The optimization action (19) will then maximize the observability measure det(ΩΩ T ) for the SfM task at hand by (i) adjusting the direction of the camera linear velocity v and, at the same time, (ii) by adapting the shape of the m weighted moments s = (m w (θ 1 ), . . . , m w (θ m )) as only a function of the perceived scene and camera motion. We also remark that (19) (or any other equivalent strategy) assumes the possibility of acting at will on the direction of the linear camera velocity v. There could be cases where this is not (fully) possible, and v (or components of it) are given (for example during a combined estimation/servoing loop as in [START_REF] Spica | Coupling Visual Servoing with Active Structure from Motion[END_REF]). In all these cases, it is obviously still possible to just keep on optimizing θ(t) during the (given/known) camera motion in order to adapt, as best as possible, the moment shape. Finally, since Ω(t) (and thus ρ(t)) does not depend on the camera angular velocity, one can freely choose ω to fulfil any additional goal of interest.

V. SIMULATION RESULTS

All the following simulations consider a free-flying camera observing a planar scene P consisting of N = 30 points, and with plane parameters n = (0, 0, -1) and d = 1.5 [m] in 0 F C (the initial camera frame at t = t 0 ). The initial estimations of the plane normal and distance are always taken as n(t 0 ) = (-0.87, 0, -0. As for the SfM algorithm providing an estimation χ(t) = -n(t)/ d(t) of the unknown 3D structure χ(t) = -n(t)/d(t), we made use of the general scheme recently discussed in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF] in the context of active SfM and already exploited in, e.g., [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF], [START_REF] Spica | Experiments of plane estimation by active vision from point features and image moments[END_REF], [START_REF] Spica | Active Structure from Motion: Application to Point, Sphere and Cylinder[END_REF] in a number of different applications. The reader is referred to these works for full details on the inner machinery of the algorithm: for our goals it is just worth mentioning that, for a given choice of the estimation gains, the convergence rate of the structure estimation error defined as

z(t) = χ(t) -χ(t)
is directly determined by the norm of matrix ΩΩ T . Thus, the larger the value of ρ(t) from ( 15) during the camera motion, the faster the expected convergence of z(t) → 0.

A. Unconstrained polynomial basis

We start with the results obtained by making use of the unconstrained polynomial basis (8) of fixed degree δ introduced in Sect. III-A.1. In particular, we tested our method by considering a set of m = 3 weighted moments s = (m w (θ 1 ), m w (θ 2 ), m w (θ 3 )) ∈ R 3 with degree δ = 2 defined as in [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF], with then θ i ∈ R 5 , i = 1 . . . 3, and θ = (θ 1 , θ 2 , θ 3 ) ∈ R p , p = 15. This choice was meant to provide a direct comparison against the use of:

1) the more 'classical set' (x g , y g , µ 20 + µ 02 ) that, as explained, is known to be an optimal choice for controlling the camera translational motion but also to yield poor results when employed for SfM purposes; 2) the set of five moments (x g , y g , µ 20 , µ 11 , µ 02 ) which, as reported in [START_REF] Spica | Experiments of plane estimation by active vision from point features and image moments[END_REF], does allow for a converging estimation but at cost of an increased complexity (need of propagating five image moments). The goal of the comparison is to prove that estimation of vector χ is, instead, fully possible when a suitable combination of just three moments of order up to 2 is selected.

Figures 2(a-b) start showing the results obtained by employing the set (x g , y g , µ 20 + µ 02 ) for estimating vector χ while, at the same time, optimizing the camera linear velocity v by implementing the first row of (19) with k v = 1. b-d)). For a comparison, in all plots dashed lines correspond to the use of the five moments (xg, yg, µ20, µ11, µ02): it is worth noting how, despite the increased measurement set (five moments vs. only three), the estimation convergence results still slower than in the weighted moment case camera motion. As expected, and even despite the velocity optimization, the value of ρ(t) keeps (numerically) very close to 0 with a maximum of ≈ 1.2 • 10 -9 (Fig. 2(a)). Thus, the chosen set (x g , y g , µ 20 + µ 02 ) is not able to provide enough information for allowing convergence of the SfM scheme, and indeed the estimation error z(t) = χ(t) -χ(t) even starts diverging (Fig. 2(b)).

The linear velocity was initially set to

v(t 0 ) = [0 0.1 0] T [m/s] with then v(t) = v(t 0 ) = 0.1 [m/s] during the
On the other hand, exploiting the three weighted moments of degree 2 yields a much more satisfactory estimation performance: Figs. 3(a-d) report in solid lines the results obtained by implementing (19) with k v = 1 and k θ = 3, and by taking again v(t 0 ) = [0 0.1 0] T [m] as in the previous case. The parameter vector θ was instead chosen at random under the constraint θ i (t 0 ) = 1, i = 1 . . . 3.

Looking at Fig. 3(a) one can then verify how now ρ(t) attains an overall quite larger value compared to Fig. 2(a), with a maximum of ≈ 1.8 • 10 -2 (thus, more than 10 7 times larger than in the previous case). As a result, the estimation error z(t) is able to quickly converge towards 0 in about 4 seconds (Fig. 3(b)). For a better appreciation of the estimation performance, Figs. 3(c-d) also report the behavior of e d (t) = d(t) -d(t) (the error in estimating the plane distance d) and e n = arccos(n T (t) n(t)) (the angular error in estimating the direction of the plane normal n) with d = 1/ χ and n = -χ/ χ . Finally, Figs. 3(ad) superimpose in dashed lines the behavior of ρ(t) and of the estimation errors when instead relying on the set of m = 5 moments (x g , y g , µ 20 , µ 11 , µ 02 ) for estimating χ: in this case, the estimation error does actually converge (as expected), but nevertheless at a slower rate compared to the weighted moment case (indeed, the maximum value of ρ(t) is now 'only' ≈ 5.9 • 10 -4 ). We then believe these results clearly show the advantages of the proposed approach: the SfM scheme has its best performance when relying on the optimization (19) for automatically selecting (online) the best combination of three moments of order up to 2.

As an additional evaluation, Figs. 4(a-d) show the results obtained when only optimizing the parameter vector θ while keeping a constant linear velocity v(t) = v(t 0 ) during the whole motion (thus, by setting k v = 0 in (19)). This case is meant to assess the optimization performance in a situation in which the camera velocity cannot be arbitrarily adjusted but it must be considered as 'given' by an external source. Thus, the only possibility for improving the conditioning of the observability matrix ΩΩ T is to act on vector θ, i.e., on the moment shape. Nevertheless also in this situation ρ(t) still reaches a range of values comparable with the previous case, with indeed max ρ(t) ≈ 1.05•10 -2 against the previous 1.6•10 -2 (thus, still ≈ 10 7 times larger than when employing the classical set (x g , y g , µ 20 + µ 02 )). As a result, vector z(t) keeps converging to 0 in about 4 [s] (Fig. 4(b)) even if slightly more slowly w.r.t. the previous case of Fig. 3(b) (as one could expect because of the smaller value of ρ(t)). In any case, we believe it is worth noting how the sole optimization of the moment shape (via vector θ) is still able to yield a very satisfactory SfM performance even for a non-optimal camera motion.

We finally remark that in all simulations the camera angular velocity ω was exploited for keeping the centroid of the observed point features p i at the center of the image b-d): behavior of the error quantities z(t), e d (t) and en(t). Note how all quantities keep behaving smoothly despite the frequent loss of tracked points.

In dashed lines, the behavior that all quantities would have had in case no optimization of θ had been performed plane (we recall that matrix Ω and, thus, ρ(t) do not depend on ω that can then be freely chosen without affecting the estimation performance).

B. Constrained polynomial basis

We now address the case of the constrained polynomial basis (10-11) described in Sect. III-A.2 and meant to smoothly take into account the loss/gain of point features because of the camera limited fov. We consider again a set of m = 3 weighted moments s = (m w (θ 1 ), m w (θ 2 ), m w (θ 3 )) ∈ R 3 with both functions w x (•) and w y (•) taken as the fifthorder polynomials given in [START_REF] Spica | Active Structure from Motion: Application to Point, Sphere and Cylinder[END_REF] with, therefore, a total of p x + p y -8 = 4 parameters to be optimized. The initial camera velocity v(t 0 ) was set as in the previous cases, and the optimization action (19) was again implemented with k v = 1 and k θ = 3. The camera angular velocity ω was Figure 5(a) shows the camera trajectory during the estimation task, and Fig. 5(b) reports the number N p (t) of tracked features over time: after about 5 [s] some points start being lost, dropping from a total of 30 to a minimum of 6 at t = 10 [s]. Nevertheless, thanks to the adopted constrained weighted moments, the scene structure is still correctly estimated without suffering from discontinuities or numerical instabilities because of the lost features. Figure 6(a) reports again the behavior of ρ(t) (blue solid line) that reaches a maximum of ≈ 2.9•10 -4 before starting to decrease at t ≈ 5 [s] because of the fewer tracked points. As a comparison, Fig. 6(a) also reports the superimposed behavior of ρ(t) in case no optimization of vector θ had been performed (the almost horizontal red dashed line). In this case, the maximum attained value for ρ(t) would have been ≈ 1.2 • 10 -6 (100 times smaller), thus proving again the importance of properly optimizing the shape of the chosen weighting function w(•). Figures 6(b-d) then show (in solid lines) the behavior of the estimation error z(t) and of the corresponding quantities e d (t) and e n (t) that smoothly reach convergence in about 10 [s] of motion despite the loss of point features. Again, for a comparison, Figs. 6(b-d) also report the superimposed behavior (in dashed lines) of the estimation errors in case of no optimization of vector θ (all quantities have a slower convergence rate as expected).

Finally, Fig. 7 depicts four snapshots of the shape of function w(x, y, θ 1 ) used to compute the first constrained weighted moment. The reader can then appreciate how the function shape evolves over time and, in particular, automatically polarizes its peaks around the location of the tracked point features. A video showing an animation of the reported simulations is also attached to the paper.

VI. CONCLUSIONS

In this paper we addressed the problem of automatically selecting an optimal set of image moments for improving the performance in estimating the 3D structure of a planar scene. This is achieved by replacing the classical definition of unweighted moments with the novel concept of parametric weighted moments that can span a whole class of image moments as a function of a parameter vector. By then optimizing online this parameter vector, the moment shape can be automatically adjusted as a function of the perceived scene. A number of simulation results, involving the estimation of the 3D parameters of a planar scene from measured moments (and known camera motion), fully confirms the effectiveness of the proposed approach. We are currently working towards an experimental validation of these results, as well as the application of these ideas in the context of image-based visual servoing.

Fig. 1 :

 1 Fig. 1: Three examples of the constrained polynomial basis w(x, y, θ) in (10-13). Note how w(x, y, θ) smoothly vanishes at the image borders of size [-1, 1] × [-1, 1]

  49) and d(t 0 ) = 1 [m], thus representing an initial incertitude of ≈ 60 [deg] on the real normal direction and of 0.5 [m] on the real distance to the plane. Finally, the point features p k , k = 1 . . . N , are sampled at 60 Hz and then corrupted component-wise by a uniformly distributed random noise of magnitude 2 pixels before being processed for evaluating the image moments. The camera motion (and the optimization (19)) is instead updated at 100 Hz.

Fig. 2 :

 2 Fig. 2: Results obtained by employing the classical moment set (xg, yg, µ20 + µ02) and optimizing for the camera linear velocity v. In this case ρ(t) keeps very close to zero (the scale of Fig. (a) is 10 -9 ) and as a consequence the estimation error z(t) does not converge (Fig. (b))

Fig. 3 :

 3 Fig.3: Solid lines: results obtained by employing three weighted moments of degree δ = 2 defined as in[START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF], and optimizing for both vector θ and the camera linear velocity v. Fig.(a): the value of ρ is now ≈ 10 7 times larger than in the previous case of Fig.2(a). This then allows for a quick convergence of the error quantities z(t), e d (t) and en(t) (Figs.(b-d)). For a comparison, in all plots dashed lines correspond to the use of the five moments (xg, yg, µ20, µ11, µ02): it is worth noting how, despite the increased measurement set (five moments vs. only three), the estimation convergence results still slower than in the weighted moment case

Fig. 4 :

 4 Fig. 4: Results obtained by employing three weighted moments of degree δ = 2 defined as in (9) and only optimizing for vector θ. Fig.(a): note how ρ(t) still reaches a range of values comparable with the previous case of Fig. 3(a) despite the lack of any optimization of the camera velocity v. Figs. (b-d): behavior of the error quantities z(t), e d (t) and en(t)

Fig. 5 :Fig. 6 :

 56 Fig. 5: Results obtained by employing the constrained weighted moments (10-11), a camera with limited fov, and optimizing for both vector θ and the linear velocity v. Fig. (a): camera trajectory and direction of the optical axis during the estimation task. Fig. (b): number Np of tracked point features over time

Fig. 7 :

 7 Fig. 7: Four snapshots of the weighting function w(x, y, θ1) taken during the camera motion. It is interesting to visualize how function w(•) automatically adjusts its shape as a function of the observed scene (e.g., it tends to peak around clusters of points)

All quantities are expressed in the camera frame, and the camera is assumed calibrated.

Of course, the use of different functional bases, also non-polynomial, could be possible.