1.

Introduction

Photoacoustic (PA) tomography is a multiple wave imaging method that seeks primarily to obtain the optical absorption coefficient of heterogeneously absorbing media, by the resolution of an inverse problem. Biological tissues are both heterogeneously optically absorbing and diffusing media. The accurate quantification of the optical absorption coefficient can be obtained provided the optical diffusion process is accounted for,¹ as well as experimental factors (e.g., transducers response, illumination profile). Multiple source illumination^2,3 is a method to discriminate optical absorption and diffusion where the medium is actively probed, locally, under various incidences, to assess its response.

Another advantage, not exploited so far, of performing local illumination is that it generates a PA wave near the interface. Usually, this signal is time-filtered or not measured (transducers oriented in specific directions). However, this PA wave, originating from the external boundary of the object [which we refer to as the boundary PA (BPA) wave], generally follows a longer ray path within the object. Therefore, it is more likely to be scattered by any of the acoustical heterogeneities than PA waves. In addition, knowing the exact position of the BPA source makes it possible to perform scattering measurements that lead to the separation of the scattered field from the incident field. While the incident field is representative of the BPA source (power, directivity) and gives information about the mean speed of sound along the ray path,⁴ the scattered field carries information on the acoustical heterogeneities (localization, geometry, and impedance contrast). In addition, the pressure amplitude of the BPA wave can eventually be increased by considering an optically absorbing coating on the external boundary of the object probed, or, as proposed in Refs. 4 and 5, by placing a localized absorber in the matching medium at a distance from the sample. This last situation is even more advantageous for tissue safety and for discriminating (in time) between the different sorts of signals. Thus, with local source illumination, the interesting point is that two different echo sequences coexist that can be separated with respect to the physical processes from which they originate. The two waves (BPA and PA) can be exploited separately. Indeed, with BPA, both the time of generation and the localization of its point source can be used as prior knowledge. As a consequence, the BPA wave brings spatial information that can be introduced afterward in the PA signal processing to retrieve the optical properties.

Usually, these secondary echoes are filtered by considering reflectivity (half-time) measurements^6,7 or weighted in the reconstruction process in order to penalize long-time arriving signals.⁸ Otherwise, they produce artifacts in the reconstructions of the initial pressure distribution maps. In this work, we propose a method that utilizes these echoes of purely acoustic nature to retrieve prior information in the locations of the abnormalities in the tissues. Even if the quantification focuses on the reconstruction of the optical parameters, having access to acoustic information is of primary importance to access to higher accuracy in the quantification. A number of multimodality systems^9–12 have been developed to access this information. They show the improvement that measuring the acoustic properties can bring.

In this paper, we want to take advantage of the pure acoustic information contained in the PA measured signal, by additionally taking into consideration the fact that a biological abnormality, such as a tumor, usually presents both optical and acoustic contrasts. First, the phenomenon is described through examples of measurements performed on tissue-mimicking phantoms, showing the feasibility of such measurements without the use of an external or superficial absorber. The information content of the acquired signals is presented and discussed. The signal processing protocol is detailed and the motivation for using the acoustic echoes measurement as prior knowledge is illustrated through reconstructions performed on synthetic phantoms.

2.

Photoacoustic Signal in the Presence of Acoustic and Optical Heterogeneities

The principle of PA probing consists of heating the medium with a time varying optical source and collecting the acoustic waves generated by thermal expansion with conventional acoustic detectors (transducers). Sound and light velocities being so different, the optical fluence can be deposited instantaneously within the whole medium according to the spatial distribution of the absorption and the diffusion coefficients. The deposited energy is dissipated into heat, and microdilations are produced that give birth to pressure waves that propagate and probe the medium.

The spatial distribution of the light fluence depends on the optical properties of the medium. In the wavelength range of interest (red to near infrared), biological tissues are highly scattering media, with a mean free path less than $100\mu \mathrm{m}$, but also absorbing (absorption coefficient of skin varying in this wavelength ranges from 0.01 to $1{\mathrm{cm}}^{-1}$).¹³ When considering large organs, such as the breast, the fluence will be mostly deposited near the illuminated surface. If the medium is homogeneous in terms of both optical and acoustic properties, a single acoustic pressure wave signal is going to be measured due to the deposition of energy. Hence, its time-of-flight corresponds mostly to the distance between the illuminated surface and the detector position. This initial pressure wave, i.e., the BPA wave, propagates and potentially probes the medium as a conventional acoustic wave (Fig. 1, signal in black). If acoustic heterogeneities are present in the medium, for example, because of a different density as in tumors (mammary tumors are visible with ultrasound echography), then this BPA wave is going to be perturbed, and echoes related to the presence of acoustic heterogenities, which we refer to as E-BPA, will be sensed at longer times-of-flight (Fig. 1, signal in blue). If optically absorbing heterogeneities are present within the medium, more energy will be absorbed locally giving birth to the conventional PA echoes (Fig. 1, signal in green) arriving at a shorter time-of-flight compared to the above-mentioned BPA echo generated at the interface. Schematically, in Fig. 1, on the one hand, the PA signal (in green) localizes an abnormality $T$ belonging to a circle $\mathcal{C}$ of radius $\overline{TD}=t_{\mathrm{id}}/\nu$, $\nu$ being the speed of sound, $t_{\mathrm{id}}$ the measured time-delay of the wave propagating from the inclusion to the detector, and additional measurements are needed to localize the object without ambiguity. The E-BPA signal (in blue), on the other hand, localizes this abnormality belonging to an ellipse $\mathcal{E}$ of focii the position of the detector $D$ and the projection of the laser beam at the surface of the probed organ such that $\overline{ST}+\overline{TD}=(t_{\mathrm{si}}+t_{\mathrm{id}})/\nu$, $t_{\mathrm{si}}$ being the time-delay representative of the wave propagating from the source to the inclusion. Hence, the object $T$ can be localized without ambiguity at the intersection of $\mathcal{C}$ and $\mathcal{E}$ that occurs inside the probed sample area. For multiple object localization, one has to multiply the measurements and plot all possible occurrences.

Fig. 1

Schema illustrating the principle: (a) a sample containing an optical and acoustic abnormality centered at position $T$ is locally illuminated at the boundary in $S$ by a pulsed laser (pulse width approximately few nanoseconds) that generates two simultaneous PA waves; (b) three echoes can be probed with a transducer located at position $D$: the PA and BPA signals, and E-BPA resulting from the perturbation of the BPA by the presence of the acoutic abnormality; (c) the time-delay of PA signal localizes the abnormality on a circle while the one of the E-BPA localizes it on an ellipse.

2.1.

Illustration on Synthetic Data

This principle can be illustrated first through simulations preformed by solving the forward softly coupled model describing the propagation of the time evolution of PA wave fields $p(\mathbf{r},t)$ (equations of linear acoustics for isotropic media, irrotational pressure flow, shear waves neglected)¹⁴

Eq. (1)

$$\{\begin{array}{l}\frac{1}{{\nu }^2}\frac{{\partial }^{2}p}{\partial t^2}-\rho \nabla ·(\frac{1}{\rho }\nabla p)=0;\\ \frac{\partial p(\mathbf{r},t=0)}{\partial t}=0;p(\mathbf{r},t=0)=p_0(\mathbf{r})=\mathrm{\Gamma }(\mathbf{r}){\mu }_a(\mathbf{r})\phi (\mathbf{r}),\end{array}$$

$\nu$ ($\mathrm{m}{\mathrm{s}}^{-1}$) being the speed of sound and $\rho$ ($\mathrm{kg}{\mathrm{s}}^{-3}$) the mass density, $\mathrm{\Gamma }$ ($-$) being the Grueneisen coefficient, ${\mu }_a$ (${\mathrm{m}}^{-1}$) the optical absorption coefficient, $\phi$ ($\mathrm{J}{\mathrm{m}}^{-2}$) the deposited fluence, and ${\mu }_a\phi$ ($\mathrm{J}{\mathrm{m}}^{-3}$) being the absorbed energy per unit volume transformed into heat at position $\mathbf{r}$. As the medium fulfills, the diffusion approximation assumptions¹⁵ (large medium, scattering predominating absorption), $\phi (\mathbf{r})$ can be determined by solving the diffusion equation numerically, with the finite-element method (FEM) in the present case, in the continuous wave regime

Eq. (2)

$${\mu }_a(\mathbf{r}){\phi }_s(\mathbf{r})-\nabla ·[D(\mathbf{r})\nabla {\phi }_s(\mathbf{r})]=S_s(\mathbf{r}),$$

with optical source $S_s$ located at points belonging to the periphery of the organ, and $D$ (m) the diffusion coefficient. Robin (or mixed) boundary conditions apply at any point ${\mathbf{r}}_s$ of the diffusing/nondiffusing boundary between the organ and the outer region

Eq. (3)

$${\phi }_s({\mathbf{r}}_s)+2AD({\mathbf{r}}_s)\nabla {\phi }_s({\mathbf{r}}_s)·\mathbf{n}=0,$$

where $A=(R-1)/(R+1)$ is an effective parameter that accounts for index of refraction mismatch at the boundary¹⁶ and $\mathbf{n}$ an outward unit vector normal to the boundary. An impedance matching fluid is often used to match the acoustic impedance between the sample and the acoustic detector. Water is the most widely used with an index of refraction ($\sim 1.33$) close to the one of most biological tissues ($\sim 1.4$) leading to a coefficient $A\sim 1$. The optical propagation equation is solved with the FEM (freefem++)¹⁷ that allows calculation of the initial pressure distribution map. The acoustic time evolution propagation equation is then solved by using the MATLAB^® toolbox $k$-wave.¹⁸

Examples of simulations are shown in Fig. 2. Two discs with 6-cm diameter and background physical properties similar to those of soft tissues are considered. Illumination is a point-like source in both situations. Inclusions with 6-mm diameter and 10 times more absorbing inclusions have been inserted, right at the center for simplicity of interpretation of the sinograms. The first one has the same acoustic properties as the background, while in the second the speed of sound is increased to $\nu =1700\mathrm{m}{\mathrm{s}}^{-1}$, corresponding to a higher limit of values in cancer tissues.¹⁹ The corresponding simulated sinograms measured by considering 64 point detectors evenly distributed at the periphery of the tissue are represented in Fig. 2, after subtraction of the homogenous background signal to enhance the visibility of the signals (relative amplitude is $\sim 4$ between BPA and PA, and 2 between PA and E-BPA).

Fig. 2

Illustration on synthetic data. Simulations are performed on discs of diameter 6 cm (${\mu }_a=0.05{\mathrm{cm}}^{-1}$, ${\mu }_s^{\prime }=7{\mathrm{cm}}^{-1}$, $\rho =1000\mathrm{kg}{\mathrm{s}}^{-3}$, $\nu =1500\mathrm{m}{\mathrm{s}}^{-1}$) in the center of which are placed disc-inclusions (10 times more absorbing) of diameter 6 mm. (a) Sinograms obtained by considering 1 point source illumination and 64 detectors evenly distributed; (b) crossplots corresponding to transmission measurement at detector # 32. Top: no contrast in the acoustic properties; bottom: speed of sound of the inclusion is $\nu =1700\mathrm{m}{\mathrm{s}}^{-1}$.

The PA echo arriving at early times is visible in both cases while the E-BPA is visible only in the second situation. These two signals will be more identifiable as the object is more contrasted optically or acoustically. The time-evolution of the two signals is clearly different. In comparison to the PA signal, the E-BPA signal exhibits a time-evolution profile that depends on the source-transducer configuration, the time extension of the E-BPA being shorter in transmission configuration than in reflection. In reflection, the time course of the E-BPA is representative of twice the diameter of the inclusion, while in transmission, the time course is representative of the diameter of the inclusion modified by the speed of sound in the inclusion.

2.2.

Illustration on Experimental Data

Preliminary experimental measurements have been performed on phantoms to verify the measurability of the E-BPA signals. In our setup [Figs. 3(a) and 3(b)], a pulsed laser (Quantel, France, Nd:Yag, 1064 nm, 330 mJ, 5 ns, 10 Hz, frequency doubled by a KTP crystal at 532 nm) and immersion acoustic transducers (diameter 10 mm, focal length 1 in., central frequency 3.5 MHz, Olympus, France) are used. The illumination beam was enlarged in the transverse direction (beam diameter on the surface of the object 4 mm). This is a compromise between having a localized source and stable measurements (less sensitive to the surface roughness, the presence of dusts and air bubbles). The collected signal is then preamplified (preamplifier Olympus, France, gain 40 dB, cutting frequency 10 MHz) and driven into an oscilloscope (Tektronix, 4 channels, bandpass width 200 MHz) for signal analyses.

Fig. 3

(a) Schema of the experimental setup; (b) photograph of the tranducers and phantom holder; (c) phantom containing hair inclusion: photograph, sinogram of the 32 measurements from 180 deg (transmission) to 330 deg, and two extracted crossplots; (d) phantom containing a purely acoustically contrasted inclusion, sinogram of the 64 measurements from 30 deg (reflection) to 330 deg (reflection), and three extracted crossplots.

A series of demonstrative experiments was conducted under the experimental configuration described in Fig. 3(a): cylindrical phantoms were made of an aqueous gel of a mixture of agar-agar (wt. 2%, Jeulin, France) and titanium dioxide (${\mathrm{TiO}}_2$) powder (wt. 0.05%, Sigma-Aldrich, France) to mimic the level of optical scattering of the tissues (estimated reduced scattering coefficient $6{\mathrm{cm}}^{-1}$, corresponding diffusion coefficient $5.5\times {10}^{-2}\mathrm{cm}$, and absorption coefficient $0.05{\mathrm{cm}}^{-1}$ at 532 nm), in the center of which were introduced inclusions of different nature: (i) a thin black hair, typical example of dual heterogeneity (high absorption, high density) inserted in a cylindrical phantom $Ø=2.5\mathrm{cm}$ Fig. 3(c); (ii) a cylindrical inclusion $Ø=0.6\mathrm{cm}$ inserted in a cylindrical phantom $Ø=5\mathrm{cm}$, containing the same mixture of $\text{water}/\text{agar}\text{-}\text{agar}/{\mathrm{TiO}}_2$ to which was added a wt. 30% of orgasol (Orgasol 2001 EXD NAT1, Arkema, France) that modifies essentially the speed of sound (estimated speed of sound $\sim 1700\mathrm{m}{\mathrm{s}}^{-1}$) and not much density ($1.03\mathrm{g}{\mathrm{cm}}^{-3}$ for Orgasol wt. 100%)²⁰ [Fig. 3(d)]. The sinograms of the measurements collected by transducer D positioned around the phantoms, placed at a distance such that the focal volume is close to the centers of the phantoms, are reported in Figs. 3(c) and 3(d), as well as crossplots measurements; the corresponding detector positions are indicated by arrows. The first sinogram exhibits clearly the existence of three kinds of signals: (i) the sinusoidal BPA, low-frequency PA signal generated at the boundary; (ii) the quasi vertical PA signal, with constant intensity and time course ($\sim 10\mu \mathrm{s}$), typical of an optically contrasted object positioned close to the center of the circle described by the different positions of the detector; (iii) the E-BPA signal with variable intensity and time course, corresponding to the presence of an acoustically contrasted object, arriving at a delayed time-of-flight compared to the PA signal corresponding to the path source-inclusion: $t_{\mathrm{si}}\sim 8\mu \mathrm{s}$ corresponding to $\sim 1.2\mathrm{cm}$. In the second situation, almost no PA signal is measured, but one can still clearly measure the E-BPA at a delayed time-of-flight compared to PA $t_{\mathrm{si}}\sim 15\mu \mathrm{s}$ corresponding to $\sim 2.3\mathrm{cm}$.

Hence, exploiting these secondary acoustic echoes due to the presence of buried structures that are inhomogeneous both acoustically and optically, as in the case of tumors, offers the possibility to retrieve spatial information on the location and/or on the acoustic properties of the medium (speed of sound). This information can afterward be included in the reconstruction of the optical parameters as prior knowledge, in order to, alternatively, speed up the reconstructions or improve their qualities. The main difficulty lies in temporally sorting signals of different nature. One has to clearly discriminate the two species of waves, the PA ones from the scattered E-BPA ones. Fortunately, in transmission detection half-space, the corresponding species of echoes belong to two different time spaces delimited by the position of the interface: (i) PA localizes an object $T$ belonging to a circle $\mathcal{C}$ of radius $\overline{TD}=t_{\mathrm{id}}/\nu$, $\nu$ being the sound velocity; (ii) E-BPA localizes an object $T$ belonging to an ellipse $\mathcal{E}$ of focii the position of the detector $D$ and the projection of the laser beam at the surface of the phantom such that $\overline{ST}+\overline{TD}=(t_{\mathrm{si}}+t_{\mathrm{id}})/\nu$. Hence, the object $T$ belongs to the intersection of $\mathcal{C}$ and $\mathcal{E}$ that occurs inside the phantom area. For multiple objects localization, one has to multiply the measurements to build a map of probability of occurrences of the intersections. By redundancy, the map will show most likely regions of intersections.

3.

Reconstructions on Simulated Data Using Intrinsic Prior Knowledge

Hereafter we explain, through a 2-D simulation example, how the measurements of these secondary signals can be exploited and introduced in a reconstruction algorithm to reconstruct the optical parameters. The simulations were done on a homogeneous disc phantom (diameter 5 cm) containing four inhomogeneous regions. The maps of the values of the optical parameters are reported in Fig. 4.

Fig. 4

Left: description of the geometry of the numerical phantom: the phantom is a disc of diameter 5 cm with four inclusions, three discs of diameter 1 cm (a) and one square $0.5\times 0.5{\mathrm{cm}}^2$ (b). It is illuminated by a point source positioned at eight different positions around the phantom (indicated by the green points). Center: absorption coefficient ${\mu }_a$ (${\mathrm{m}}^{-1}$); right: diffusion coefficient $D$ (m).

3.1.

Abnormalities Localization Algorithm

Synthetic pressure measurements are simulated by first computing the initial pressure distribution by considering the speed of sound and density abnormalities of the medium. In the present example, 360 transducers were evenly distributed at the periphery of the object and a set of eight different point sources experiments, also evenly distributed, was run. The final sinograms were composed of $8\times 360$ time-resolved measurements. An example of a sinogram obtained for one source is plotted in Fig. 5, for given speed of sound and density maps, showing the complexity of the signals.

Fig. 5

(a) Sinogram for source #1. (b) Mass density $\rho$ ($\mathrm{kg}{\mathrm{s}}^{-3}$). (c) Speed of sound $\nu$ ($\mathrm{m}{\mathrm{s}}^{-1}$) map.

From the sinograms, an algorithm of peak detections is applied,²¹ and the spatial histogram of probabilities of intersections between $\mathcal{C}$ and $\mathcal{E}$ is calculated (Fig. 6). After filtering, the gradient of the resulting image (Fig. 6) is taken to obtain enhanced contours and a final thresholding allows the identification of a limited number of regions of interest, representing the most favorable regions of locations of tumors.

Fig. 6

Steps of localization: maps of speed of sound and density, followed by the map of probabilities of localizations, gradient of the same map after filtering and labeling of regions of interest.

These localization maps depend on the acoustic properties of the medium. They will be different as the contrast of the acoustic properties differ. However, according to tests performed over a range of acoustic parameters values that can be found in tissues, the regions will be slightly modified. Three examples are presented in Fig. 7, for the following range of acoustic properties: $1000<\rho <1300\mathrm{kg}{\mathrm{m}}^{-3}$ and $1485<\nu <1600\mathrm{m}{\mathrm{s}}^{-1}$.

Fig. 7

Examples of localization maps: (a) speed of sound map; (b) density map; (c) localization map. First row: phantom with homogeneous speed of sound; second row: phantom with homogeneous density; third row: inhomogeneous phantom.

For a reasonable range of acoustic properties values, the regions can be reasonably found. Reconstructions of the unknown physical parameters can then be spatially focused to these regions of interest. However, in vivo, the distribution of the heterogeneities may be more complex and the discrimination between the signals may be more difficult. Strategies consisting of enhancing the BPA signals by using a localized absorber, placed at the surface of the sample or at a distance within the matching medium,⁵ can be adopted in these situations. In all cases, the fluence deposited at the surface and into the tissue should be smaller than the ANSI safety limit.

4.

Results and Discussion

The reconstructions obtained for ${\mu }_a$ only with and without use of the localization map are shown in Fig. 8. The reconstruction results are slightly improved by the use of the localization map, but obtained with a much shorter computation time ($\sim 3$ times faster) as the algorithm converges faster as it can be seen in the map of the norm of the gradient.

Fig. 8

Reconstructions ${\mu }_a$ only (a) crossplot along the horizontal (b) and vertical axis (c). Red line: target values; blue line: reconstructed values. Top row: without prior knowledge; bottom: with prior knowledge.

The norm of the gradient has been plotted as a function of the number of iterations [Fig. 9(a)]. It can be seen that using prior knowledge allows a much faster convergence (approximately a factor 3). The reconstructions were obtained with a much faster computation time, as the number of unknown is reduced with localized reconstructions.

Fig. 9

Norm of the gradient as a function of the number of iterations until convergence: without (blue) and with (red) prior knowledge. (a) reconstruction of ${\mu }_a$ only (first step); (b) simultaneous reconstruction of ${\mu }_a$ and $D$ (second step).

In the second step, the joint minimization of the gradients was run. The results of the reconstruction of ${\mu }_a$ and $D$ are represented in Fig. 10. After convergence, as no noise was considered, the results of the reconstructions without and with prior knowledge on the position of the abnormalities become quite similar but are obtained much faster when prior knowledge is considered (less iterations). These results may be improved by using an adaptive mesh in the forward and adjoint problems resolution by FEM, refined on the regions of interest and coarser outside. Here, the mesh was kept the same for all calculations.

Fig. 10

Reconstructions of ${\mu }_a$ (first row) and $D$ (second row): (a) initial values; (b) reconstructed values after convergence; (c) crossplots along the horizontal; and (d) vertical axis. Red dashed line: target values; blue line: initial values; green line: reconstructed values.

Note that in the initial pressure measurements, the Grueneisen coefficient $\mathrm{\Gamma }$ was kept constant for simplicity of description of the reconstruction process. In the general case, when $\mathrm{\Gamma }$ is spatially varying, one can consider measurements normalized $P_s^m(\mathbf{r})=p_{s_{\ne }0}^m(\mathbf{r})/p_0^m(\mathbf{r})$ to one of them $p_0^m(\mathbf{r})$ that can serve as the reference. A new functional can be defined

The expressions of the gradients can then be derived. After reconstruction of the optical parameters, the reconstruction of $\mathrm{\Gamma }$ can be immediately obtained.

Although the algorithm has not been tested on experimental data, the robustness of the method to noisy data has been examined by adding a random noise to the initial pressure distribution maps (mean value: $p_0^m(\mathbf{r})$; standard deviation: 1%). The results of the dual step reconstruction process are presented in Fig. 11. In the first step, the reconstructions of the absorption coefficient map are degraded, but the convergence of the algorithm is not very affected by the noise. In the second step, it is to be noted that the algorithm has reached convergence only when the localization maps are considered. Nonetheless, the crosstalk already observed in the noise free data reconstructions is enhanced. While it seemed not to impact the reconstruction of ${\mu }_a$, it produces here large artifacts in the reconstruction of $D$, especially in the regions close to the sources. A Thikonov regularization (L2-norm, with regularization factor 0.25) has been introduced to lower the effect of these artifacts.

Fig. 11

Reconstructions of ${\mu }_a$ (first row) and $D$ (second row) with noisy initial pressure distribution map (1% noise) and prior knowledge: (a) initial values (first step reconstruction for ${\mu }_a$); (b) reconstructed values after convergence; (c) crossplots along the horizontal; and (d) vertical axis. Red dashed line: target values; blue line: initial values; green line: reconstructed values.

5.

Conclusion

In conclusion, the method proposed for registration and processing of the PA measurements allows more accurate reconstructions with efficient computation times. The method relies in sorting the different sources of echoes, which is made feasible as primary PA and secondary BPA echoes are clearly identified in time. It has been shown experimentally, through demonstrative experiments performed on phantoms with physical properties close to those of biological tissues that these signals can be measured experimentally through tomographic measurements. The method can be improved by quantitatively exploiting the E-BPA echoes by using proper acoustic model-based reconstruction algorithms.²⁷ The measurement of these echoes may also be facilitated by applying a more absorbing layer at the surface of the sample or by using a time-gated amplification, in order to amplify the signal in a given time-band. Measurements on biological tissues are programmed in the near future. The localization process relies on the redundancy of the measurements. It might result in detection of false regions-of-interest, which will slow the reconstruction process but it does not modify the final reconstruction results. In all cases, the threshold can be adapted. The reconstruction method can also be improved by using an adaptative mesh in the resolution of the forward and the adjoint problems to benefit even more of the prior knowledge. Before applying the reconstruction algorithm to experimental data, the behavior of the reconstruction algorithm is to be studied with noisy data. Preliminary tests show that it is even more efficient in this situation.