F Vialla

Y Chassagneux
email: yannick.chassagneux@lpa.ens.fr

R Ferreira

C Roquelet

C Diederichs

G Cassabois

Ph Roussignol

J S Lauret

C Voisin

Unifying the low-temperature photoluminescence spectra of carbon nanotubes: the role of acoustic phonon confinement

At low temperature the photoluminescence of single-wall carbon nanotubes show a large variety of spectral profiles ranging from ultra narrow lines in suspended nanotubes to broad and asymmetrical line-shapes that puzzle the current interpretation in terms of exciton-phonon coupling. Here, we present a complete set of photoluminescence profiles in matrix embedded nanotubes including unprecedented narrow emission lines. We demonstrate that the diversity of the low-temperature luminescence profiles in nanotubes originates in tiny modifications of their low-energy acoustic phonon modes. When low energy modes are locally suppressed, a sharp photoluminescence line as narrow as 0.7 meV is restored. Furthermore, multi-peak luminescence profiles with specific temperature dependence show the presence of confined phonon modes.

Single-Wall Carbon Nanotubes (SWNTs) are fascinating nanostructures with unique properties that make them promising for a number of opto-electronic applications. The understanding of their basic photo-physical properties has been the focus of an intense research effort with major leaps such as the measurement of an exceptional exciton binding energy [1], the detection of trions and bi-excitons [2,3] or the observation of photon anti-bunching due to the localization of excitons [4]. In contrast, the understanding of the low temperature photoluminescence (PL) spectra of carbon nanotubes remains a pending issue, with apparently inconsistent results reported in the literature. For instance, the low temperature luminescence profiles of SWNTs range from broad (more than 10 meV) asymmetrical lines to narrow (less than 40 µeV) and symmetrical lines [5][6][7][8][9].

The general framework of the coupling of a localized exciton to an acoustic phonon bath was developed by Krummheuer et al. [10]. In the case of nanotubes, Galland et al. [7] showed that this model predicts intrinsically broad and asymmetrical line-shapes. This so-called Ohmic model accounts for a limited set of luminescence profiles but fails to explain the variety of profiles reported in the literature. In this model, the red shifted (resp. blue shifted) wing of the PL line corresponds to the radiative recombination of the exciton assisted by the emission (resp. absorption) of an acoustic phonon. In contrast to the case of epitaxial quantum dots embedded in a 3 dimensional (3D) phonon bath where a sharp purely electronic Zero Phonon Line (ZPL) accompanied by weak phonon side-bands is observed [11,12], a strong prediction of this model in the case of a 1D phonon bath is the complete merging of the ZPL into the phonon wings [13]. Thus, for a 1D phonon bath, the dominant optical process near the electronic resonance is the simultaneous emission of a photon and a phonon. As a consequence, the photoluminescence line-shape is driven by this non perturbative coupling to long wavelength acoustic phonons rather than by the eigen properties of the electronic two-level system. This peculiar exciton-phonon coupling provides a new means to probe the coupling between exciton and low-energy acoustic phonons with an unprecedented sensitivity.

In this Letter, we present a complete set of PL spectra spanning a large variety of profiles and line-widths. We show that, due to phonon confinement, the coupling of excitons to acoustic phonons can deviate strongly from the Ohmic regime especially for phonon energies lower than a few meV. This has drastic effects on the photoluminescence spectra recovering a sharp ZPL. Phonon confinement allows us to explain all the observed spectra both for matrix embedded and suspended SWNTs, provided that the much longer dephasing time of suspended nanotubes is properly taken into account.

The PL spectra were recorded on individual SWNTs using a homemade confocal microscope. The surfactant (Sodium Cholate) embedded CoMoCat SWNTs were deposited onto the flat surface of a high refractive index solid immersion lens (SIL) to improve the excitation and collection efficiencies. Prior to SWNTs deposition, the SIL was functionalized with a poly-L-lysine layer to increase the linkage of the micelles to the substrate. The sample was placed on the cold finger of an optical helium cryostat.

Optical excitation was performed with a cw Ti:Sapphire laser with intensities of 1kW/cm 2 or below. Exposure times of typically one minute were used.

Fig. 1 (a)-(d) show the typical low-temperature PL profiles observed in our samples. The first one (a) consists of an asymmetrical profile and is similar to those reported by Galland et al. [7]. The second one (b) exhibits a sharp peak onto an asymmetrical pedestal. Note that this kind of profile is closer to the ones observed for suspended SWNTs [5,6]. We also found a continuous set of spectra (such as in Fig. 1 (c)) bridging the narrow and the broad ones. Finally, the last example (d) shows a multi-peak feature onto an asymmetrical pedestal.

In fact, we note that the pedestal (Fig. 1 (b)-(d)) has -in all spectra-a shape and a width (a few meV) comparable to that of the asymmetrical spectrum (Fig. 1 (a)). The temperature dependence of the PL line shape of tubes (c) and (d) is given in Fig. 2 and shows striking similarities : when raising the temperature, the blue wing tends to increase leading to a symmetrized profile. Importantly, we found that the relative height of the red and blue side peaks in Fig. 1 (d) fits to a Boltzmann factor with an activation energy of 2.5 meV equal to the splitting between the peaks [13]. This rules out that these peaks would arise from spectral diffusion. Actually, spectral diffusion would not lead to such temperature dependence since there is no simple relationship between the trapping energy of a charge and the magnitude of the related Stark shift [13]. More generally, we managed to fit the relative weight of the blue and red wings to a Boltzmann factor for all spectra, showing that these wings are related to phonon assisted processes [13]. Thus, we speculate that all the PL profiles originate from the coupling of excitons to phonons.

In order to extract the coupling between the exciton and acoustic phonons, we used the calculations developed for semiconductor quantum dots in [10] and adapted for SWNTs in [7]. We assume that excitons couple only to the 1D acoustic phonons of the nanotube due to the strong sound velocity mismatch between the nanotube (v ≃ 2×10 4 m/s) and the substrate (≃ 3 × 10 3 m/s). In carbon nanotubes, three acoustic phonon branches are coupled to the exciton: the stretching, twisting and radial breathing (RBM) modes [14,15]. The acoustic phonons couple to excitons through the difference of the off-diagonal terms of the deformation potential for electrons and holes D = D e -D h . Due to the complementary dependences with the chiral angle, the different branches can be merged into an effective acoustic mode [16]. The exciton-phonon coupling matrix element is given by [10]:

hg(ω) = hω 2ρLv 2 F (ω) (1)
L is the nanotube length, ρ is the linear mass density, v is the sound velocity, hω is the energy of the phonon mode and F (ω) is a form factor. In the case of a perfectly 1D acoustic mode, one has F = F 0 (Ohmic coupling) with :

F 0 (ω) = D dz|ψ(z)| 2 e iqz (2)
where ψ(z) is the exciton center-of-mass envelope wave function and where ω = v|q|. Note that in this Ohmic model one has F 0 (0) = D regardless of the exciton wave function, which leads to the merging of the ZPL into the phonon wings. In the case of a Gaussian exciton envelope, ψ(z) =

1 π 1/4 σ 1/2 exp(-z 2 2σ 2 )
, where σ is the exciton localization length, the form factor reads :

F 0 (ω) = D exp[-(ωσ/2v) 2 ] (3)
This implies a high energy cut-off of the exciton-acoustic phonon coupling at hω σ ≈ 2hv/σ that simply results from momentum conservation within the width of the exciton envelope in the k-space.

The emission spectra can be computed from the coupling matrix element by summing up the contributions of all phonons modes weighted by their occupation numbers [7,10]. Conversely, the form factor F (ω) can be traced back from the experimental spectra with the procedure described in [13]. The form factor extracted from the broad and asymmetric profile is presented in Fig. 1 (e). This form factor fits well to a Gaussian such as F 0 , as expected in the Ohmic model [7].

Even if the experimental spectra show a large variety of profiles (Fig. 1 (a-d)), the extracted form factors share common features. Importantly, we note that they do not depend on the temperature (as expected from the model). In addition, all the form factors display a Gaussian baseline, from which we can extract a deformation potential D between 12 and 15 eV, and an exciton localization length σ of around 3 nm, similar to the ones found in [7]. The main deviation from the Ohmic coupling occurs for low phonon energies, typically below 3 meV. For instance the narrow spectrum of Fig. 1 (b) corresponds to a complete suppression of the coupling (F (ω) = 0) for low energy phonons (Fig. 1 (f)). Again, we stress that this reduced coupling cannot be accounted for by a modification of the exciton wave function in the Ohmic model, since the form factor goes to D for vanishing ω regardless of the envelope shape (Eq. 2). In terms of emission spectrum, this means that changing the exciton envelope function (shape or localization length) cannot account simultaneously for the observed narrow ZPL (requiring a delocalized exciton) and a broader pedestal (requiring a tightly localized exciton). Therefore, the lowenergy modification of the form factor must arise from the phonon modes themselves (as discussed at the end of the paper).

The key to the understanding of the strong variations of the PL line shapes actually lies in the low-energy exciton-phonon coupling. When F takes sizable values for vanishing phonon energies, the ZPL is strongly suppressed to the profit of the phonon wing. In contrast, if the form factor is reduced for low phonon energies (ω < ω c ) a sharp ZPL is restored (Fig. 1 (b)). Practically, this effect is observable if the cut-off phonon energy hω c is greater than the natural line-width of the ZPL : ω c > 2/T 2 , where T 2 is the ZPL dephasing time. Experimentally, we found a width of ≈ 0.65 meV for the ZPL and hω c ≈ 1.2 meV in Fig. 1 (a), which explains why the phonon sideband and the ZPL are not completely split off.

More generally, the effect of the dephasing time on the spectral line shape is presented in Fig. 3. In the case of the Ohmic coupling, the line shape is roughly insensitive to the value of the electronic dephasing time. This is a peculiarity of the 1D geometry, where the ZPL is completely merged into the phonon wings. In the case of a low-energy gap in the form factor (in Fig. 3 we used a 1 meV gap), a sharp ZPL is restored. The amplitude of the ZPL is directly proportional to the dephasing time T 2 since the area of the line must be constant, reflecting the oscillator strength of the transition. In contrast the phonon wings hardly depend on T 2 . Longer dephasing times such as those observed in suspended nanotubes therefore lead to an apparent magnification of the ZPL. For example in the Fig. 1 of [5], a sharp ZPL and a weak phonon sideband are observable [17]. The ZPL / sideband amplitude ratio is of the order of 100 for a ZPL width of 40 µeV. Consistently, the spectrum displayed in Fig. 1 (b) shows a ratio of the order of 5, for a ZPL width of 600 µeV.

We note that some suspended NT spectra have no apparent phonon wings. We believe that this is due to experimental limitations. In fact, reaching a given signal to noise ratio requires much shorter integration times for sharp ZPLs, which may loose the weak sidebands in the noise background. [18] Let us now examine possible origins of the reduction of the coupling to the low-energy phonon modes. Equivalently, this low-energy cut-off can be seen as a reduction of the low energy phonon density of states probed by the exciton. Generally speaking, such local modifications of the phonon modes can be accounted for by introducing acoustic barriers, with height hω c ≃ 2 meV [13]. Actually, since 1D structures are much more sensitive to perturbations than their higher dimension coun- terparts, several ordinary effects may cause important phonon modifications. In the case of suspended SWNTs, amorphous carbon deposited by places or structural defects can break the 1D translational symmetry, leading to low energy phonon modifications [19]. In the case of embedded SWNTs, the non-covalent contact to the substrate or matrix is known to occur on nanometer scale zones [20] and leads to the opening of few-meV-phonon gaps [21] due to a local hardening of the nanotube lattice.

The different types of PL profiles then arise from different geometries of the barrier. We make use of the general expression for the form factor F :

F (ω i ) = D √ L v/ω i dz|ψ(z)| 2 W ′ i (z) (4)
where W i (z) are the orthonormal acoustic modes, solutions of :

-ω 2 i W i (z) = v 2 W ′′ i (z) -ω 2 b (z)W i (z) (5)
For the sake of simplicity, we assume here that the acoustic barrier amplitude hω b (z) is either 0 or hω c along the nanotube (fig. 1(i-l) . Without barrier, free propagating modes are obtained and F = F 0 (Fig. 1 (a) and (e)).

If the barrier width is much larger than the exciton localization length, a true gap is found in the form factor (Fig. 1 (b) and (f)). For an intermediate barrier width, typically a few tens of nanometers, the coupling is only partially canceled, owing to the finite attenuation length of the evanescent low energy phonon modes inside the barrier (Fig. 1 (c) and (g)) [13].

An additional case of interest corresponds to the situation where the exciton is located between two acoustic barriers separated by a few tens of nanometers. Our model predicts that in this case side-peaks should appear in the phonon wing because the exciton is coupled to a resonant vibration mode between the barriers. We actually found experimental evidences for such resonant coupling (Fig. 1 (d) and (h)). The observed splitting (2.5 meV) corresponds to a distance between the barriers of the order of 20 nm, consistent with the scale of variations of the environment of the nanotube. This original observation corresponds to the interesting situation of disorder-induced confinement of vibrations at the nanoscale.

In conclusion, the low temperature PL profiles of carbon nanotubes (ranging from sharp sub-meV lines to broad asymmetrical lines) are well accounted for within a single model describing the coupling of a localized exciton and a one-dimensional phonon bath. The key to the diversity of the profiles lies in tiny modifications of the low-energy phonons. The observation of a narrow lineshape indicates the presence of a gap in the excitonphonon coupling. Furthermore, we pointed out the key role of the ZPL dephasing rate in the overall profile of the PL, explaining the apparent discrepancy between suspended and embedded nanotubes. This reduction of the coupling to acoustic phonons is a first important step towards the control of the emission line-width of matrix embedded carbon nanotubes. Finally, we showed that phonons may be confined by disorder with specific fingerprints in the PL spectrum.

Supplemental Material for : "Unifying the low-temperature photoluminescence spectra of carbon nanotubes: the role of acoustic phonon confinement".

SINGLE EMITTER

The sample was characterized by Atomic Force Microscopy (AFM). We observed surfactant aggregates along with a few well individualized tubes. Their length ranges from 100 nm to 1 µm. Since the density of SWNTs can strongly vary from one location to another at the surface of the sample, optical studies were performed only in zones with a rather low density of emitting tubes (≃ 10 -1 µm -2 ).

At cryogenic temperatures, the PL signals exhibit blinking, spectral jumps and spectral diffusion which are typical of single emitters. Extensive optical studies were performed on the most stable ones. In the following we present observations made on the SWNTs described in Fig. 1 (c) and Fig. 1 (d) of the paper. Similar behaviors were observed for several dozens of SWNTs. All the spectra were recorded at 10K.

Long time traces with short exposure time steps show a rather stable emission energy, with a spectral diffusion smaller than the linewidth, and a stable emission intensity (Fig. S-1 and Fig. S-2). In the particular case of the multi-peak spectrum, we observe that the magnitude of the spectral diffusion (2 pixels = 0.6 meV) remains much weaker than the typical splitting between the components of the spectrum (of the order of 8 pixels) and that this splitting is constant up to 1 pixel (0.3meV) whatever the absolute position of the lines.

In addition, spatial PL maps show clear localized emission spots (Fig.

S-3 and Fig. S-4). A 2-dimensional

Gaussian fit gives a full width at half-maximum of 410 ± 30 nm, close to the diffraction limited resolution of the setup. However, one cannot discriminate between a nanometer scale localized exciton in a long tube and a more delocalized exciton in a non resolved very short nanotube.

Moreover, we always observed strongly anisotropic excitation and emission polarization diagrams (Fig. S-5), which is typical of nanotubes for which the oscillator strength is much larger for a polarization along the tube axis.

All these observations support the fact that the spectra arise from single nanotubes rather than from several emitters.

Special care has been taken to come up to this conclusion in the case of the multi-peak spectrum. Indeed, spatial maps of the 3 different peaks show a perfect overlap within the Gaussian fit error bars of 30 nm (Fig. S-4). In addition, we can note the correlated polarization be- Importantly, we note that the characteristic asymmetrical and multi-peak profiles reported in the paper clearly show up for any integration time. Together with the singular temperature dependence of the spectra, we can safely rule out that these feature arise from spectral diffusion.

BRANCHING RATIO BETWEEN ZPL AND PHONON WING

The Hamiltonian of a two-level excitonic system coupled to an acoustic phonon mode reads:

H = hΩ|e e| + hω i b † i b i + hG i |e e| ⊗ (b i + b † i ) (6)
where |g and |e are the ground and excited electronic state, b i (resp. b † i ) is the annihilation (resp. creation) acoustic phonon operator, hΩ is the energy difference between the ground and excited electronic state and hω i is the phonon energy. The bulk electron -acoustic phonon coupling G i is given by:

G i = D ωi 2ρV v 2 h
, where D is the deformation potential, ρ the mass density, V a normalisation volume and v the sound speed (for an unidimensional system, ρ is the linear mass density and V should be replace by L the length of the system). H can be diagonalized into:

H = hω i |g g|⊗b † i b i + |e e|⊗(hΩ+hω i D † ξ b † i b i D ξ )-hω i ξ 2 (7)
where ξ = G i /ω i and D ξ is the displacement operator given by D ξ = exp(ξb † i -ξb i ). The eigenstates of this hamiltonian are given by |g, n and D † ξ |e, n . The probability of optical transition is proportional to

g, n ′ | d |D † ξ |e, n 2
, where d is the dipolar operator acting only on the electronic part of the wavefunction. The branching ratio (br) between the emission of a photonphonon pair (i.e. in the wings) and the emission of a photon alone (i.e. in the ZPL) (cf Fig. S-6) is given by:

br i = 1|D † ξ |0 0|D † ξ |0 2 = ξ 2 exp(|ξ| 2 ) exp(|ξ| 2 ) = G 2 i ω 2 i ( 8
)
When considering the total coupling to all phonon modes, the previous result is multiplied by the density of acoustic phonon states, given by: DOS ∝ ω Dim-1 , where Dim is the dimensionality of the phonon bath.

From the dependency of the coupling G i with ω i , one obtains that the branching ratio goes as :

br ∝ ω Dim-2 (9)
This ratio becomes diverging at low energy for a onedimensional phonon bath, leading to a complete merging of the ZPL into the phonon wings.

TEMPERATURE DEPENDENCE OF THE PHONON WINGS

For each spectrum, we evaluate the amplitude ratio between the blue and the red wings. We introduce the ratio I +ǫ /I -ǫ where I +ǫ (resp. I -ǫ ) is the PL intensity at an energy E = E max + ǫ (resp. E = E max -ǫ), E max being the energy at maximum PL intensity (Fig. S-7). In a simple single phonon coupling interpretation, the blue wing involves the creation of an acoustic phonon proportional to n ǫ the occupation number of a phonon mode with energy ǫ. Conversely, the red wing involves the annihilation of a phonon and is proportional to n ǫ +1. We obtain I +ǫ /I -ǫ = n ǫ /(n ǫ + 1) = exp(-ǫ/k B T ). An effective temperature T ef f can be deduced from this relation. We note that T ef f is accurately evaluated only for temperatures T lower than a few ǫ/k B (∼ 30 K) since a I +ǫ /I -ǫ ratio close to 1 (within the PL noise) gives diverging uncertainties for T ef f . We find a good agreement between T ef f and the experimentally measured temperatures (Fig. S-8). The slightly higher values of T ef f , especially at very low temperature (T ≤ 10 K) may come from local heating due to the excitation laser beam. We note that the temperature dependence is better reproduced by the full model (as shown in the main text) because it includes the ZPL dephasing time, which tends to smooth the spectra and also because it includes the contribution of multi-phonons emission (or absorption), which generates deviations from the simple Boltzmann factor.

FORM FACTOR: NUMERICAL ESTIMATION

Following the methods presented in [7,10], the spectrum can be computed from the form factor. The relationship between the spectrum and the form factor is a mathematically reversible operation, and thus it ensures that a given spectrum corresponds to an unique form factor. We note that the noise on the experimental spectrum leads to numerical instabilities in the inversion process. Therefore, we choosed to use a random optimization procedure to extract the form factor. We initialized the procedure with a flat form factor. Then, by means of an iterative procedure, we searched for random modifications that make the computed spectrum closer to the experimental one. This procedure gives reproducible form factors, at the cost of some high frequency noise.

ACOUSTIC BARRIERS

The electron -acoustic phonon coupling Hamiltonian is given by : H e-AP = D div(u) (10) where the displacement u reads :

u = i h 2ρω i (b i + b † i )W i (z) (11)
The W i (z) form an orthonormal basis of the acoustic modes. For an exciton, the coupling Hamiltonian H X-AP can be simplified in :

Ψ X |D div(u)|Ψ X = i hg ωi (b i + b † i ) (12)
with the matix element g ωi given in the main text.

hg ωi = hω i 2ρLv 2 F (ω i ) (13)
where F (ω i ) reads :

F (ω i ) = D √ L v/ω i Ψ X dW i (z) dz Ψ X (14)
To find the acoustic mode W i one has to solve the following equation :

-ω 2 i W i (z) = v 2 d 2 W i (z) dz 2 - 0 , no barrier ω 2 c W i (z) , barrier (15)
The free propagating mode corresponds to the upper line, whereas the lower one describes the acoustic barrier due to the local coupling to the environment (for instance a local increase of the spring constant). The discretization of the eigen frequencies ω i and eigen modes W i come from the finite length of the tube L = 1 µm in the calculation. For the sake of simplicity, we assume that the acoustic barrier amplitude is either 0 or hω c along the length of the nanotube.

Once the acoustic modes W i are normalized, the form factor F can be computed. Simplified analytical expressions can be obtained for simple symmetrical configurations with the exciton at the centre (position z = 0 in Fig. 1 (i-l)). The free propagating case (no acoustic barrier F = F 0 ) and the case of an infinitely long acoustic barrier (F = F ∞ ) are represented in Fig. 1 (i) and (j) respectively. Approximated expressions can be obtained when the barrier length a is finite but much larger than the exciton size σ (L >> a >> σ, F a , Fig. 1 (k)). Similarly, an approximated analytical expression can be derived for two identical barriers (of width a) spaced by a distance l (F l,a , Fig. 1 (l)).

F 0 (ω) = D exp - σ 2 4v 2 ω 2 (16)
F ∞ (ω) = F 0 (ω) Re[ω] ω (17)
F a (ω) ≈ F 0 (ω) 1 + ω 2 c sin 2 (aω/2v) ω2 Where ω is the complex number ω = ω 2 -ω 2 c , and α and β are given by α = aω/v and β = lω/2v.

For each case, a fit to the form factor extrated from experimental spectra presented in Fig. 1 in the main text is performed. The best fits were obtained with the set of parameters listed in Table I. The deformation potential D and exciton enveloppe size σ are similar to the one evaluated with the Gaussian baseline. The barrier height hω c is found to be between 1 and 2 meV. The spatial structures (barrier widths and distance between barriers) are of the order of 10 to 35 nanometers. These fits to the form factor were further used to simulate the emission spectra with an excellent accuracy, as shown in Fig. 1 (a)-(d) of the main text. I. Best fit parameters for the data displayed in Fig. 1. The acoustic barrier is defined with its height hωc and width a. In the case of the multipeak spectrum (Fig. 1 (d)), the gap width between the two identical barriers is l.

Finally, we can deduce from the statistical occurrence of the different cases presented in Figure 1, that there is no preferential localization sites of the exciton with respect of the mechanical contact zones. Among a few tens of specimens, we found that most spectra fall in the category of Figure 1 (c), corresponding to an intermediate coupling scheme between the limit cases (a) and (b). The observation of confined phonon modes such as in Figure 1 (d) represents less than 10% of the cases.

FIG. 1 .

 1 FIG. 1. (a) to (d) PL spectra of individual SWNTs recorded at 10K (open dots). Chiral species are indicated in parenthesis.The red lines are simulated spectra using the fit of the form factor of panels (e) -(h). (e) to (h) Form factors (see text) extracted from the corresponding spectra. The continuous lines are fits to the acoustic barrier model[13]. (i)-(l) Schematic representation of two exemplary local phonon modes (hω = 1 meV (red line) and hω = 2.5 meV (blue line)) resulting from the presence of acoustic barriers (depicted in green). For the spectrum simulation, the exciton is supposed to be located at z = 0.

FIG. 2 .FIG. 3 .

 23 FIG. 2. Evolution of the normalized PL spectra (vertically shifted for clarity) as a function of temperature (black open dots). The red lines are the calculated spectra from the extracted form factor at the lowest temperature (Fig.1(e) to (h)). The measured base plate temperature is indicated together with the effective temperature used for the fitting.

  FIG. S-1. (a) Time trace of the PL spectrum of the nanotube presented in Fig. 1 (c), recorded every second. (b) PL spectra recorded every 10 s with a 1 s exposure time. The color stripe shows the integration energy range (1.5 meV) used for Fig. S-3. .

  FIG. S-2. (a) Time trace of the PL spectrum of the nanotube presented in Fig. 1 (d), recorded every 2 seconds. (b) PL spectra recorded every 20 s with a 2 s exposure time. The color stripes show the integration energy range (1.5 meV) used for Fig. S-4-5. .

3 .

 3 FIG. S-3. Confocal spatial PL map of the nanotube in Fig. 1 (c). For each point of the map, a PL spectrum was recorded with a 1 s exposure time and the integrated intensity over the energy range given in Fig. S-1(b) is then reported according to the color scale. Spatial steps are 100 nm long .

FIG. S- 4 .

 4 FIG. S-4. Confocal PL maps of the 3 different peaks observed in the spectrum of the nanotube in Fig. 1 (d). For each point of the map, a single PL spectrum was recorded with a 2 s exposure time and the integrated intensity over each energy range given in Fig. S-2(b) is then reported according to the color scale. (a) (b) and (c) correspond respectively to the peak at 1.416 eV, 1.418 eV and 1.42 eV. Spatial steps are 100 nm long .

FIG. S- 5 .

 5 FIG. S-5. Excitation polarization diagram of the integrated PL of the nanotube of Fig. 1 (d). The 3 different peaks are represented according to their respective color attributed in Fig. S-2(b) .

FIG. S- 7 .FIG. S- 8 .

 78 FIG. S-7. Evolution of 3 different spectra presented in the main article as a function of the temperature. Spectra are normalized and shifted vertically for clarity. The PL intensities I-ǫ and I+ǫ are represented by a vertical bar for an energy shift ǫ of 2.3 meV.

Figure

  Figure 1 (e) (f) (g) (h) D [eV] 14 12 12 15 σ [nm] 2.4 2.4 2.6 3.0 hωc [meV] -1.1 1.3 1.9 a [nm] --34 8 l [nm] ---20

Acknowledgement

This work was supported by the GDR-I GNT, the grant "C'Nano IdF TENAPO". CV and GC are members of "Institut Universitaire de France".


Jan Heller

Didier Henrion

Tomáš Pajdla

Stable radial distortion calibration by polynomial matrix inequalities programming

Polynomial and rational functions are the number one choice when it comes to modeling of radial distortion of lenses. However, several extrapolation and numerical issues may arise while using these functions that have not been covered by the literature much so far. In this paper, we identify these problems and show how to deal with them by enforcing nonnegativity of certain polynomials. Further, we show how to model these nonnegativities using polynomial matrix inequalities (PMI) and how to estimate the radial distortion parameters subject to PMI constraints using semidefinite programming (SDP). Finally, we suggest several approaches on how to incorporate the proposed method into the overall camera calibration procedure.

Introduction

Radial distortion modeling is the most important non-linear part of the camera calibration process [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. The first works on the topic came from the photogrammetric community [START_REF] Brown | Decentering distortion of lenses[END_REF][START_REF] Brown | Close-range camera calibration[END_REF][START_REF] Chester | Manual of photogrammetry[END_REF]. Since then, a plethora of models has been suggested in the literature [START_REF] Sturm | Camera models and fundamental concepts used in geometric computer vision[END_REF]. Among the proposed models, the ones based on polynomial and rational functions are the most popular. This popularity undoubtedly stems from the fact that these function are easily manipulated and yet provide sufficient fitting power for wide range or distortions. Unfortunately, the extrapolation qualities of polynomials can be quite unpredictable in situations where little or no data is available. However, even if data points are missing, the overall shape of the distortion is known a priori in many calibration scenarios, e.g., the lens introduces barrel or pincushion distortions. Based on such a priori information, the shape of the polynomial and rational distortion functions can be controlled by enforcing nonnegativity of certain polynomials. For example, in the case of pincushion distortion we can accomplish the desired shape by enforcing nonnegativity of the first and the second derivatives of the distortion function on the whole field of view of the camera.

In this paper, we propose a radial distortion calibration procedure where a polynomial cost function, e.g., reprojection error, is minimized subject to such shape constraints. This shape optimization procedure is designed to stabilize the shape of the distortion function. It is based on polynomial matrix inequalities (PMI) programming and can be easily incorporated into an existing camera calibration procedure.

In Section 2, we formally introduce the radial distortion function and present several extrapolation issues arising while using polynomial and rational distortion models. Next, in Section 3 we provide a minimal theoretical background needed for our shape stabilization approach. In Section 4, we demonstrate the proposed method on three types of radial distortion shapes and models and show how to incorporate the method into an overall camera calibration procedure. Finally, in Section 5 we experimentally validate our approach and show that the method guarantees the correct shape of a distortion function without compromising the quality of the overall camera calibration as measured by the reprojection error.

Camera Radial Distortion

Let us suppose that a set of scene points X i ∈ R 3 , i = 1, . . . , n is observed by a camera. If R ∈ SO(3), t ∈ R 3 are the camera extrinsic parameters, a scene point X i gets projected into an image point (x i , y i , 1) ⊤ :

λ i (x i , y i , 1) ⊤ = RX i + t, λ i ∈ R.
In reality, some amount of radial distortion is always present and the camera observes a point (x i , ŷi , 1) ⊤ which does not coincide with the ideal (and unobservable) point (x i , y i , 1) ⊤ . In pixel coordinates, the camera observes a point K(x i , ŷi , 1) ⊤ , where K ∈ R 3×3 is the matrix of intrinsic camera parameters, the so-called calibration matrix. Radial distortion function L : R → R is a function of radius r = x 2 i + y 2 i that models the radial displacement of the ideal image point position from the center of the radial distortion as

xi ŷi = L(r) x i y i . (1)
The function L(r) is only defined for r > 0 and L(0) = 1, L(r) > 0. For the purposes of demonstration of the proposed shape optimization procedure, we will use L(r) defined as follows

L(r) = f (r) g(r) = 1 + k 1 r + k 2 r 2 + k 3 r 3 1 + k 4 r + k 5 k 2 + k 6 r 3 , (2)
where k = (k 1 , k 2 , . . . , k 6 ) is the vector of model parameters. This definition accommodates several models already proposed in the literature [START_REF] Ma | Rational radial distortion models of camera lenses with analytical solution for distortion correction[END_REF]. However, we will see that the shape optimization procedure holds for any rational function.

Extrapolation issues of radial distortion calibration

Let us motivate the need for the radial distortion shape optimization by demonstrating two examples of extrapolation issues arising while using polynomial and rational distortion models. However, the recovered polynomial quickly drifts away elsewhere (red circles depict the distances of the projections of the image corners). In green, a polynomial recovered by the method proposed in this paper is shown. Here, the negativity of the first and the second derivatives of the polynomial on the whole field of view was enforced. This caused the model to fit the original distortion much closer on the whole field of view. Figure 1(b) shows a synthetic checkerboard image (the upper left corner) and the same image distorted by the original barrel distortion (the upper right corner). In the lower left corner, the image is undistorted back using the polynomial recovered by [START_REF] Zhang | A flexible new technique for camera calibration[END_REF]2]. In the lower right corner, the image successfully undistorted by the polynomial recovered using the proposed shape optimization method is shown.

Let us consider a similar calibration scenario to the one from the previous paragraph, this time with a lens causing a mustache type radial distortion, see Figure 1(c). If the radial distortion model is used, L(r) = f (r) g(r) , the classical calibration approach [START_REF] Zhang | A flexible new technique for camera calibration[END_REF]2] is able to correctly recover the original shape. However, the polynomials f (r) and g(r) share a common root (red dash-dot lines), which causes a numerical instability presented as a sharp spike in L(r) around the common root-an issue we will call the zero-crossing problem. When the nonnegativity of g(r) is enforced using the proposed approach, not only is the correct shape recovered, but since there is now no root in the field of view interval (green dash-dot lines), the spike in L(r) is also gone. Figure 1(d) shows a similar arrangement as Figure 1(b), now with only the upper left part of the checkerboard shown. The numerical instability of L(r) is presented as a notable ringing in the upper left corner of the checkerboard. One can argue that the common root is a consequence of the fact that the degrees of f (r), g(r) are higher that needed and that a model with fewer coefficients should be used. This may be true in some cases, however, we observed just as many situations where the lower degree polynomials resolved the zero-crossing problem only at the cost of a considerably higher reprojection error.

Polynomials and PMI Programming

In this section, we present a minimal theoretical background needed for the proposed shape optimization procedure.

Polynomials and polynomial matrices

An univariate polynomial p(x) ∈ R n [x] of degree n ∈ N is a real function defined as

p(x) = p n x n + p n-1 x n-1 + • • • + p 1 x + p 0 = p ⊤ ψ n (x),
where p = (p 0 , p 1 , . . . , p n ) ⊤ ∈ R n+1 is the vector of coefficients with a nonvanishing coefficient p n and

ψ n (x) = (1, x, x 2 , . . . , x n ) ⊤ is the canonical basis. Let q(x) ∈ R 2n [x]. A symmetric matrix Q ∈ R n ′ ×n ′ , Q = (q i,j
), where n ′ = n + 1, is called Gram matrix associated with q(x) and the basis ψ n (x) [START_REF] Choi | Sums of squares of real polynomials[END_REF] if

q(x) = ψ ⊤ n (x) Q ψ n (x). (3)
Generally, there is more than one Gram matrix associated with a polynomial q(x) and we will denote the set of such matrices as G(q(x)).The polynomial q(x) can be expressed in the elements of Q by simply expanding the right hand side of Equation 3 and by comparing the coefficients.

Let x = (x 1 , x 2 , . . . , x d ) ∈ R d be a real vector and α = (α 1 , α 1 , . . . , α d ) ∈ N d an integer vector. A monomial of degree n = α i is defined as x α = n i=1 x α i i . A multivariate polynomial p(x) ∈ R n [x] of degree n ∈ N is a mapping from R d to R defined as a linear combination of monomials up to degree n, p(x) = |α|≤n p α x α = |α|≤n p α x α 1 1 x α 2 2 • • • x α d d = (p α ) ⊤ |α|≤n (x α ) |α|≤d = p ⊤ ψ n (x),
where p ∈ R m is the vector of coefficients and ψ n (x) is the canonical basis of m = d+n d monomials up to degree n. By a polynomial matrix we will understand a symmetric matrix whose elements are polynomials. In the next, S n (R[x]) will denote the set of n × n symmetric polynomial matrices. The degree of P = (p i,j (x)) ∈ S n (R[x]) is the largest degree of all the polynomial elements of P, deg P = max i,j deg p i,j (x).

Besides parameterizing polynomials by the associated Gram matrices, we will also need to "linearize" them, i.e., to substitute every monomial x α by a new variable y α ∈ R. To do this, we define the Riesz functional ℓ y : R n [x] → R[y], a linear functional that for a d-variate polynomial of degree n, p(x) = α p α x α , returns an m-variate polynomial of degree one, ℓ y (p(x)) = α p α y α , m = d+n d . With a slight abuse of notation, we will also use ℓ y as a matrix operator acting on

S n (R[x]): if P ∈ S n (R[x]), then P ′ = ℓ y (P) if and only if p ′ i,j (y) = ℓ y (p i,j (x)).

Polynomials positive on finite intervals

The shape optimization procedure presented in this paper is based on enforcing nonnegativity of certain polynomials. Since most of the real cameras have limited fields of view, we only need to control the behavior of L(r) for values r ∈ [0, r], where r is the maximal distance between the center of the radial distortion and an (undistorted) image point.

For this, we need to characterize the set of univariate polynomials nonnegative on finite intervals. In [START_REF] Nesterov | Squared functional systems and optimization problems[END_REF], based on Markov-Lukacs theorem, Nesterov showed how to characterize such a set using positive semidefinite Gram matrices:

Theorem 1 Let α < β, p(x) ∈ R[x] and deg p(x) = 2n. Then p(x) ≥ 0 for all x ∈ [α, β] if and only if p(x) = s(x) + (x -α)(β -x)t(x), where s(x) = ψ ⊤ n (x) S ψ n (x), t(x) = ψ ⊤ n-1 (x) T ψ n-1 (x)
, such that S, T 0 (i.e., S ∈ G(s(x)), T ∈ G(t(x)) are positive semidefinite Gram matrices of polynomials s(x) and t(x), respectively).

If deg p(x) = 2n + 1, then p(x) ≥ 0 for all x ∈ [α, β] if and only if p(x) = (x -α)s(x) + (β -x)t(x), where s(x) = ψ ⊤ n (x) S ψ n (x), t(x) = ψ ⊤ n (x) T ψ n (x), such that S, T 0.
Even though Theorem 1 is an equivalence, we will only use it as an implication: as long as we will have matrices S, T that are positive semidefinitive, Theorem 1 guarantees that a polynomial p(x) constructed using these matrices will be nonnegative on a given interval.

Polynomial Matrix Inequalities

According to Theorem 1, a polynomial is nonnegative on an interval as long the matrices S, T are positive semidefinite. By combining these constraints with a polynomial cost function, we get a problem of polynomial matrix inequalities (PMI) programming. A PMI program can be formally defined as follows:

Problem 1 (Polynomial matrix inequalities program) minimize p(x) subject to G i (x) 0, i = 1, . . . , m, where p(x) ∈ R[x], G i ∈ S n i (R[x]).
In general, Problem 1 is a hard non-convex problem. Note however, that if the cost function p(x) and the matrices G i (x), i = 1, . . . , m have degree one, then Problem 1 reduces to a linear matrix inequality (LMI) program and as such is a semidefinite program (SDP) solvable by any available SDP solver. In fact, most of the time the shape optimization problems in this paper lead to such a program.

Sometimes still, G i (x) will not be linear. In such cases, we will use the relaxation approach suggested by Henrion and Lasserre [START_REF] Henrion | Convergent relaxations of polynomial matrix inequalities and static output feedback[END_REF]. In [START_REF] Henrion | Convergent relaxations of polynomial matrix inequalities and static output feedback[END_REF], the authors proposed a hierarchy of LMI programs P 1 , P 2 , . . . that produces a monotonically non-decreasing sequence of lower bounds p(x * 1 ) ≤ p(x * 2 ) ≤ . . . on Problem 1 that converges to the global minimum p(x * ). Practically, the series converges to p(x * ) in finitely many steps, i.e., there exists j ∈ N, such that p(x * j ) = p(x * ). The authors also showed how this situation can be detected and how the value of x * can be extracted from the solution of the relaxation by the tools of linear algebra.

Let us show here how to construct P δ , i.e., the LMI relaxation of Problem 1 of order δ; see [START_REF] Henrion | Convergent relaxations of polynomial matrix inequalities and static output feedback[END_REF] for the technical justification of this procedure. Let G ∈ S n (R[x]), n = m i=1 n i denote a block diagonal matrix with matrices G i on it's diagonal. Since (∀i : G i (x) 0) ⇔ G(x) 0, we can replace the PMI constraints G i (x) 0 with one PMI constraint G(x) 0. Next, we construct the so-called moment matrix M δ (y) and localizing matrix M δ (G, y) of G, defined as

M δ (y) = ℓ y (ψ δ (x)ψ ⊤ δ (x)), M δ (G, y) = ℓ y ((ψ δ (x)ψ ⊤ δ (x)) ⊗ G),
where ⊗ denotes the Kronecker product [START_REF] Henrion | Convergent relaxations of polynomial matrix inequalities and static output feedback[END_REF].

Let γ = 1 if deg G ≤ 2, γ = ⌈deg G⌉ 2
otherwise. Now, we can formally write the relaxation P δ as Problem 2 (LMI relaxation P δ of order δ)

minimize ℓ y (p(x)) subject to M δ-γ (G, y) 0, M δ (y) 0.
As the Riesz functional ℓ y was used to "linearize" both the cost function and the constraints, we can easily see that Problem 2 is an LMI program.

Shape optimization for radial distortion calibration

In this section, we show how to combine the results presented in Section 3 into the radial distortion shape optimization procedure. Technically, the procedure consists of minimization of a polynomial cost function in the vector of radial distortion parameters k subject to PMI constraints enforcing nonnegativity of certain polynomials in the radius r. Such a minimization problem is a PMI program that can be dealt with using the approach from Section 3.3.

As mentioned in Section 3.2, we only need to control the shape of L(r) on the interval [0, r]. Note, that r is the maximal distance between the center of the radial distortion and undistorted image points, i.e., the value of r is not known prior to the actual calibration. The value of r is therefore a user supplied parameter. Fortunately, the proposed method is not very sensitive to the value of this parameter and even a gross overestimate yields minima identical to the ground truth value.

Unconstrained radial distortion calibration

There are several ways how to determine the vector of parameters k of the distortion function L(r) [START_REF] Hartley | Multiple view geometry in computer vision[END_REF][START_REF] Szeliski | Computer vision: algorithms and applications[END_REF]. All we need for our shape optimization approach is a polynomial cost function. Here, we will define and use one of such possible cost functions. Let us rewrite Equation 1 using L(r) from Equation 2as

g(r) xi ŷi -f (r) x i y i = g(r) xi -f (r) x i g(r) ŷi -f (r) y i = 0.
By factoring out the vector of parameters k and by denoting

A i = -r x i -r 2 x i -r 3 x i xi r xi r 2 xi r 3 -r y i -r 2 y i -r 3 y i ŷi r ŷi r 2 ŷi r 3 , b i = x i -xi y i -ŷi ,
we get a linear system

A i k = b i . Now, we can stack A = (A ⊤ 1 , A ⊤ 2 , . . . , A ⊤ n ) ⊤ , b = (b ⊤ 1 , b ⊤ 2 , . . . b ⊤ n ) ⊤
and estimate the radial distortion parameters k = (k 1 , k 2 , . . . , k 6 ) as a solution to an overdetermined system Ak = b in the least square sense, i.e., by minimizing Akb 2 . Note that for polynomial model, i.e., g(x) = 1, this corresponds to the minimization of the reprojection error.

Let us now express the minimization of Akb 2 as an LMI program. By expanding

A i k -b i 2 = (A i k -b i ) ⊤ (A i k -b i ) = k ⊤ A ⊤ i A i k -2b ⊤ i A i k + b ⊤ i b i and by denoting M = n i=1 A ⊤ i A i , m = -2 n i=1 A ⊤ i b i , c = n i=1 b ⊤ i b i ,

we can write the polynomial form of the cost function as

Ak -b 2 = k ⊤ M k + m ⊤ k + c. (4)
As expected, Equation 4 is a quadratic polynomial in k and by construction M 0, i.e., M is a positive semidefinite matrix. Even though the cost function is quadratic, it can be converted into a linear function using the Schur complement trick [START_REF] Boyd | Convex Optimization[END_REF]:

F = I Lk k ⊤ L ⊤ -m ⊤ k -c + γ 0 ⇔ k ⊤ L ⊤ Lk + m ⊤ k + c -γ ≤ 0.
By decomposing M as M = L ⊤ L, e.g., using the Cholesky or the spectral decomposition [START_REF] Golub | Matrix computations[END_REF] (recall that M 0), we can rewrite the minimization of Equation 4as the following LMI program:

Problem 3 (Unconstrained radial distortion calibration) minimize γ subject to F = I Lk k ⊤ L ⊤ -m ⊤ k -c + γ 0.

Barrel distortion and the polynomial model

As we can see from the example of barrel radial distortion in Figure 1(a), this type of distortion can be characterized by the negativity of the first and the second derivatives:

∀r ∈ [0, r] : L ′ (r) ≤ 0 & L ′′ (r) ≤ 0, (5)
where [0, r] spans the field of view of the camera. If we consider the polynomial model L(r) = f (r), the constraints above mean that we need to enforce nonnegativity of polynomials

-f ′ (r) = -k 1 -2k 2 r -3k 3 r 2 , -f ′′ (r) = -2k 2 -6k 3 r on the interval [0, r]. According to Theorem 1, -f ′ (r) ≥ 0 for ∀r ∈ [0, r] iff -f ′ (r) = -k 1 -2k 2 r -3k 3 r 2 = ψ 1 (r) ⊤ S 1 ψ 1 (r) + r (r -r) T 1 , (6)
where

S 1 = s 11 s 12 s 12 s 13 0, T 1 = (t 11 ) 0.
By expanding the right hand side of Equation 6and by comparing the polynomial coefficients, we get a parameterization of k in the elements of S 1 and T 1 :

-k 1 = s 11 -2k 2 = 2s 12 + r t 11 -3k 3 = s 13 -t 11    ⇒ k = (-s 11 , -s 12 -1 2 rt 11 , 1 3 (t 11 -s 13 ), 0, 0, 0). (7)
Let's apply Theorem 1 to -f ′′ (r) to get the following constraint:

-f ′′ (r) = -2k 2 -6k 3 r = r S 2 + (r -r)T 2 , S 2 = (s 21 ) 0, T 2 = (t 21 ) 0. ( 8
)
By combining Equations 8 and 7, we can express the entries of S 2 and T 2 in the entries of S 1 , T 1 :

-2k 2 =
Now, we have four PMI constraints on the shape of L(r). If we combine these constraints along with the parameterization of k from Equation 7with Problem 3, we get a radial distortion calibration problem that enforces a barrel type distortion shape of the resulting distortion model:

Problem 4 (Barrel distortion calibration) minimize γ subject to F 0, S 1 0, T 1 = (t 11 ) 0, S 2 = 1 r (2s 12 + 2r s 13 -rt 11 ) 0, T 2 = 2 r (s 12 + 1 2 rt 11 ) 0.
Problem 4 is a PMI program in 5 variables γ, s 11 , s 12 , s 13 , t 11 . Since both the cost function and the PMI constraints have degree one, Problem 4 is in fact an SDP problem. Once it is solved, the unknown distortion parameters k can be easily recovered using Equation 7.

Pincushion distortion and the division model

Let us make an analogous analysis for the pincushion distortion shape and the division model L(r) = 1 g(r) . This type of distortion is characterized by the nonnegativity of the first and the second derivatives of L(r) on the field of view of the camera [0, r]. From the first derivative we get the following constraint on the polynomial denominator g(r):

L ′ (r) = -g ′ (r) g 2 (r) ⇒ L ′ (r) ≥ 0 ⇔ -g ′ (r) ≥ 0.
The second derivative yields a bit more complicated constraint:

L ′′ (r) = g(r)h(r) g 4 (r) = h(r) g 3 (r) ⇒ L ′′ (r) ≥ 0 ⇔ (g(r) ≥ 0 & h(r) ≥ 0) ∨ (g(r) ≤ 0 & h(r) ≤ 0),
where h(r) = 2(g ′ (r)) 2 -g(r)g ′′ (r). However, since we know that L(r) > 0 by definition, we only need to consider the constraints g(r) ≥ 0, h(r) ≥ 0. Let us start with the constraint g(r) ≥ 0. According to Theorem 1, g(r) ≥ 0 for ∀r ∈ [0, r] iff

g(r) = 1 + k 4 r + k 5 r 2 + k 6 r 3 = ψ 1 (r) ⊤ S 1 ψ 1 (r) + (r -r) ψ 1 (r) ⊤ T 1 ψ 1 (r), (10)
where

S
t 11 = 1 r (11)
By applying Theorem 1 to the constraint -g ′ (r) ≥ 0, we get

-g ′ (r) = -k 4 -2k 5 r -3k 6 r 2 = ψ 1 (r) ⊤ S 2 ψ 1 (r) + r (r -r) T 2 , ( 12
)
where

S 2 = s 21 s 22 s 22 s 23 0, T 2 = (t 21 ) 0.
As in the case of the barrel distortion optimization, we can express the entries of S 2 and T 2 in the entries of S 1 , T 1 . This time, however, we have more variables than equations and we have to set one of the entries free-we chose s 22 :

-3k

The final constraint is the most complicated because of the quadratic monomials in k: h(r) > 0 for ∀r ∈ [0, r] iff

h(r) = (6k 6 r 2 + 4k 5 r + 2k 4 )(3k 6 r 2 + 2k 5 + k 4 ) - -(2k 5 + 6k 6 r)(k 6 r 3 -k 5 r 2 + k 4 r + 1) (14) = ψ 2 (r) ⊤ S 3 ψ 2 (r) + (r -r) ψ 1 (r) ⊤ T 3 ψ 1 (r),
where

S 3 =   s 31 s 32 s 33 s 32 s 34 s 35 s 33 s 35 s 36   0, T 3 = t 31 t 32 t 32 t 33 0.
Equation 14 gives us 5 constraints on 9 entries of S 3 and T 3 . We chose to set free variables s 32 , s 34 , s 36 , t 32 ; System 15 shows the form of the remaining 5 variables. Finally, we can combine these 6 PMI constraints, Problem 3 and the parameterization of k from Equation 11 into a radial distortion calibration problem that enforces a pincushion type distortion shape:

Problem 5 (Pincushion distorion calibration) minimize γ subject to F 0, S 1 0, T 1 0, S 2 0, T 2 0, S 3 0, T 3 0. Since S 3 and T 3 are polynomial matrices of degree 2, Problem 5 has to be dealt with using the relaxation scheme from Section 3.3.

           ⇒                                  s 31 = 4t

Zero-crossing problem of the rational model

Also the zero-crossing problem of the rational model L(r) = f (r) g(r) can be dealt with using the proposed shape optimization technique. A sufficient condition for avoiding a common root of the polynomials f (r) and g(r) on the interval [0, r] is to force at least one on them to have no root. Here, we decided on enforcing the constraint ∀r ∈ 0, r : g(r) -p ≥ 0, where p > 0.

(

) 16
Since Theorem 1 guarantees only nonnegativity of a polynomial, we need a strictly positive parameter p to enforce strict positivity of g(r). Even though parameter p must be user supplied, the method is not overly sensitive to its value; in our experiments, we set p = 0.1. By applying Theorem 1 to the above constraint and the interval [0, r], we get

g(r) -p = 1 -p + k 4 r + k 5 r 2 + k 6 r 3 = ψ 1 (r) ⊤ S 1 ψ 1 (r) + (r -r) ψ 1 (r) ⊤ T 1 ψ 1 (r),
where

S

Shape optimization in Camera Calibration Procedure

All of the calibration problems presented in this paper expect the projection coordinates x i , y i , xi , and ŷi to be known, see Equation 1. This assumes a known calibration target X i ∈ R 3 as well as known camera parameters R ∈ SO(3), t ∈ R 3 , and the calibration matrix K ∈ R 3×3 . A straightforward idea how to fold the shape optimized radial distortion calibration into the camera calibration procedure is to first perform "classical" camera calibration [START_REF] Roger | An efficient and accurate camera calibration technique for 3d machine vision[END_REF][START_REF] Zhang | A flexible new technique for camera calibration[END_REF][START_REF] Hartley | Parameter-free radial distortion correction with center of distortion estimation[END_REF][START_REF] Tardif | Calibration of cameras with radially symmetric distortion[END_REF], including radial distortion estimation. Once the projection coordinates are known, the shape optimized radial distortion calibration can be performed to replace the radial distortion parameters estimated by a classical method. One might argue that the quality of such a solution could be compromised, since different error functions may be considered by the camera and the shape optimization calibration methods.

To mitigate this problem, we suggest an alternating approach to "shape-optimize" the results of the classical camera calibration: first, the shape optimization procedure is performed, followed by a bundle adjustment [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF] step where the radial distortion parameters are fixed. This can be repeated in a loop for a fixed number of times, or until desired convergence is reached.

Experiments

To validate the proposed approach, this section presents several experimental results on synthetic as well as real world datasets. We implemented Problems 4, 5, and 6 in Matlab using Yalmip toolbox [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF] with SeDuMi [START_REF] Jos | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF] as the underlying SDP solver. Yalmip toolbox is a modeling language that can be used to solve LMI as well as PMI programs, which it automatically translates into LMI relaxations using the scheme presented in Section 3.3. All of the resulting SDP programs were solved under a second on an Intel i7 3.50GHz based desktop computer running Linux and 64bit Matlab.

Synthetic experiment. In the synthetic experiment, we studied the performance of the proposed method with respect to the image noise. We generated a synthetic 16×16 planar calibration target. A scene consisted of 9 random 640×480 pixel cameras randomly positioned on a hemisphere around the target and rotated to face its center. The focal length was set to approx. 540 px and the distances of the camera centers from the target were set up so that the target (calibration data point set) covered only the middle part of the field of view, approx 50%. For each of the three model-shape problem combinations, we generated 100 scenes and corrupted the projections of the calibration target by an increasing amount of Gaussian image noise in 5 levels, standard deviation σ ∈ [0, 2] px in 1 /2 px steps. We calibrated all scenes with OpenCV [2] made to disregard the radial distortion component. We compare three methods: the first method (BA) is the bundle adjustment method that included the respective radial distortion model performed together with the OpenCV calibration results, the second method (SO) is the respective shape-optimization method performed after the BA step, and the last method (ASO) is the alternating approach from Section 4.5, fixed to 10 iterations.

Barrel distortion. First, we experimented with the barrel distortion and the polynomial model L(r) = f (r). Figure 2(a) shows the mean of the reprojection errors on the calibration data point set for methods BA, SO, and ASO using Matlab function boxplot. The methods show identical performance, however when a validation data set of points covering the whole field of view is used, see Figure 2(a), we see both SO and ASO outperforming the classical BA approach.

Pincushion distortion. Next, Figures 2(c-d) show the analogous measure for the pincushion distortion and the division model L(r) = 1 g(x) . Here, both BA and shape-optimization methods perform significantly better on the validation data point set. Still, we can see superior performance of SO and ASO as the noise increases.

Zero-crossing. Finally, we experimented with the rational model L(r) = f (r) g(r) and the mustache type distortion. Figure 2(d) shows identical performance on the calibration dataset. On the other hand, we can see poor performance on the validation data point set even if no noise is present, Figure 2(e). This is caused by the fact that too few calibration points were on the outer parts of the field of view where the convexity of the distortion Real experiment. In the real experiment, we calibrated a 2 MPix camera from Point Grey's Ladybug 3 system [1] using 12 images of a known 28×20 planar target. Calibration using BA method and the rational model introduced quite noticeable zero-crossing problem. As expected, calibration using ASO method does not suffer from this type of problem. In this experiment, we set r = 4 and p = 0.1. Figure 3(a) shows the upper left corner of a rectified calibration image using k provided by methods BA and ASO, respectively. Figure 3(b) shows the shape of the BA calibration function in red and the ASO calibration in green.

Conclusion

The aim of this work was not to argue for a specific radial distortion model, but to point out extrapolation problems inherent to all polynomial and rational models. We solved these problems by enforcing a predetermined shape of the distortion function. For most shapes and models, the proposed approach leads to small semidefinite programming problems that can be solved fast and globally optimally. We also showed how to deal with shapes and models that lead to PMI problems using a LMI relaxation scheme. We showed experimentally that in terms of the reprojection error on the known data points the proposed approach provides radial distortion models that are equivalent to those provided by the classical bundle adjustment approach, yet with the added value of having the correct shape that mollifies or completely removes all extrapolation issues.

Figure 1 :

 1 Figure 1: Calibration issues. Examples of issues arising while using polynomial and rational function for radial distortion calibration. See text for details.

Figure 2 :

 2 Figure 2: Image noise experiment. Methods BA, SO, and ASO in red, green, and blue, respectively, on calibration and validation data point sets. (a-b) barrel distortion,(c-d) pincushion distortion (e-f) zero-crossing problem. radial distortion component. We compare three methods: the first method (BA) is the bundle adjustment method that included the respective radial distortion model performed together with the OpenCV calibration results, the second method (SO) is the respective shape-optimization method performed after the BA step, and the last method (ASO) is the alternating approach from Section 4.5, fixed to 10 iterations.

Figure 3 :

 3 Figure 3: Real experiment. Correction of the zero-crossing problem of the rational model. function changes. Again, we see better performance of SO and ASO methods.

  1 = s 11 s 12 s 12 s 13 0, T 1 = t 11 t 12 t 12 t 13 0. This leads to the following parameterization of k as well as to a constraint on the variable t 11 : 1 = r t 11 k 4 = s 11 -t 11 + 2rt 12 k 5 = 2s 12 -2t 12 + rt 13 k 6 = s 13 -t 13 11 -t 11 + 2rt 12 , 2s 12 -2t 12 + rt 13 , s 13 -t 13 )

	
	 	 	k = (0, 0, 0, s
		
	⇒
	 	 
	

  4 = s 21 -2k 5 = 2s 22 + rt 21 -3k 6 = s 23 -t 21

		
		
		
	⇒
		
		

s 21 = t 11 -s 11 -2rt 12 s 23 = -1 r (s 12 + 2s 22 -4t 12 + r(3s 13 -t 13 )) t 21 = -1 r (2s 12 + s 22 -2t 12 + rt 13 )

  = 2s 33 +s 34 -t 31 +2rt 32 6k 4 k 5 -6k 6 = 2s 32 +rt 31

	12k 2 6	= s 36 -t 33
	16k 5 k 6 = 2s 35 -2t 32 +rt 33
	6k 2 5 +6k 4 k 6 2k 2 4 -2k 5 = s 31

  12 -s 12 -2rt 31 +2(s 11 -t 11 +2rt 12 ) 2 s 33 = -1 2r (6s 13 +2s 32 -6t 13 +rs 34 -6r(2s 12 -2t 12 +rt 13 ) 2 +2r 2 t 32 --6(s 11 -t 11 +2rt 12 ) (2s 12 -2t 12 +rt 13 +rs 13 -rt 13 )) s 35 = t 32 -2 r s 36 +6r(s 13 -t 13 ) 2 + +8(s 13 -t 13 )(2s 12 -2t 12 +rt 13 ) t 31 = -2 r (3s 13 +s 32 -3t 13 -3(2s 12 -2t 12 +rt 13 )(s 11 -t 11 +2rt 12 )) t 33 = s 36 -12(s 13 -t 13 ) 2 Problem 5 is a PMI program in 11 variables γ, s 11 , s 12 , s 13 , t 12 , t 13 , s 22 , s 32 , s 34 , s 36 , and t 32 .

	(15)

  , k 2 , k 3 , s 11 -t 11 + 2rt 12 , 2s 12 -2t 12 + rt 13 , s 13 -t 13 ) S 1 0, T 1 0. Problem 6 is an LMI program in 9 variables γ, s 11 , s 12 , s 13 , t 12 , t 13 , k 1 , k 2 , k 3 .

	Problem 6 (Zero-crossing distortion calibration)
		minimize γ
		subject to F 0,
	1 =	s 11 s 12 s 12 s 13			0, T 1 =	t 11 t 12 t 12 t 13	0.
	This yields a parameterization of k as well as a constraint on t 11 :
	1 -p = r t 11 k 4 = s 11 -t 11 + 2rt 12 k 5 = 2s 12 -2t 12 + rt 13 k 6 = s 13 -t 13	      	⇒	    	r k = (k 1 t 11 = 1-p	(17)
	Again, by combining the two PMI constraints with Problem 3 and the parameterization
	of k from Equation 17, we get a radial distortion calibration problem that eliminates the
	zero-crossing problem:

Faculty of Electrical Engineering, Czech Technical University in Prague, CZ-16627 Praha 6, Tech-

nická 2, Czech Republic. 2 CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France.[START_REF] Boyd | Convex Optimization[END_REF] Université de Toulouse; F-31400 Toulouse; France.