1. Introduction

In inverse scattering problems, one looks for a quantitatively accurate description of the electrical and geometrical characteristics of an investigated region given a set of incident fields and measures of the corresponding total or scattered fields (both in amplitude and in phase) on a generic surface lying outside the region under test. Due to their wide range of potential applications, the development of accurate and reliable techniques for solving this kind of problems is nowadays a still important challenge [1-11].

By leaving aside peculiar characteristics of each approach, one of the main common drawbacks resides in the need of measuring both amplitude and phase of the scattered fields. As a matter of fact, in several areas of applied science it can be very difficult or very expensive or even not possible at all to perform a measure of the phase of a field. In particular, an accurate knowledge of the phase distribution involves sophisticated measurement equipments, which are increasingly expensive as the working frequency grows, so that phaseless measurements are indeed mandatory at optical frequencies. In addition, existence of minimally invasive (only amplitude) probes strongly suggests adoption of phaseless techniques also at microwave frequencies. In fact, these probes considerably simplify the electromagnetic scenario with respect to classical (amplitude and phase) probes, since they avoid multiple interactions and do not require probe compensation. For these reasons, several approaches for solving inverse scattering problems from only amplitude distribution of the scattered [7-9] and of the total field [10,11] have been proposed in the literature.

Within this framework, an approach based on only amplitude measurements of the total field has been recently proposed by the authors [1,2], first with reference to the case of measures taken on a closed curve surrounding the domain under test [1] and then to that of transmitters and receivers placed over two truncated lines somehow ‘enclosing’ the investigated domain [2]. This approach, which opposite to other contributions [7,8] tackles the inverse problem in its full nonlinearity, splits the problem into two steps. In the first one, the scattered field is estimated from the measurement of the square amplitude of the total field, while the second step is aimed at estimating the unknown dielectric properties from the estimated scattered field. Notably, the separation of the problem into two steps allows a better control of the overall non-linearity with respect to single step procedures [10,11]. In fact, as shown in [1,2] and briefly recalled in the following, a convenient exploitation of theoretical results on the inversion of quadratic operators [12] and on field properties and representations [13] allows to successfully solve the first step, while all the available knowledge about “standard” inverse scattering problems is exploited in the second step. On the other side, this imaging technique [1,2] still requires the knowledge of both amplitude and phase of the incident field on the measurement curve and within the investigation domain, so that, in a sense, is not a strictly phaseless method, as it happens instead for [10,11].

To overcome this limitation, in this paper we show the possibility of estimating the required (complex) incident field from its amplitude distribution, as measured on the measurement curve, through the solution of a traditional Phase Retrieval (PR) problem [12,14-17]. This allows us to introducing a new three-step procedure wherein the huge amount of knowledge about PR problems [12,14-17] can be exploited in the first step, while the advantages of the original approach [1,2] are preserved in the following two steps, consisting in estimating the scattered field and the permittivity profile, respectively. In particular, different from [1], this latter task is pursued by making use of the Contrast Source - Extended Born (CS-EB) inversion method [19,20]. This recently introduced approach relies on a simple rewriting of the scattering equations which has been to some extent inspired by the work by Habashy and co-workers on the so-called Extended Born Approximation (EBA) [21]. However, unlike the EBA, the CS-EB formulation does not introduce any approximation on the scattering equations. Hence, the CS-EB model describes inverse scattering problems in their full non-linearity and it has no theoretical restriction on the class of scatterers which can be dealt with. In addition, analytical studies [19] and experimental comparisons [20] have demonstrated that the CS-EB inversion method results to be less prone to false solutions (i.e., less dependent from the adopted initial estimate) than other inversion schemes dealing with the inverse scattering problem in its full non linearity. Therefore, we make use of the CS-EB inversion method in the final step of the proposed approach in order to further improve its overall effectiveness.

The actual feasibility and performance of the new phaseless non-linear inversion method which is proposed has been tested against experimental data provided by the Institute Fresnel of Marseille, France, recently proposed as a benchmark test for inverse scattering methods [18].

2. Mathematical formulation and rationale

For the sake of simplicity, we refer in the following to the 2D scalar geometry shown in Fig. 1. The region under test D, which without loss of generality is assumed to be a circle, is embedded in a background medium of dielectric permittivity ε_b and encloses one or more homogeneous objects of complex dielectric permittivity equal to ε_r(r)ε_b. The magnetic permeability is everywhere equal to μ _o. As usual, χ(r)=[ε_r(r)-1] defines the contrast function.

 

Fig. 1. Reference geometry.

Download Full Size | PDF

Let us consider as incident fields TM-polarized unitary plane or cylindrical waves coming from a curve Γ_s encircling the targets and assume that measurement probes are located on circle Γ_s of radius R _o. According to the CS-EB model, the scattering problem at hand can be exactly formulated through the system of coupled integral equations [19]:

wherein the variables θ_i and θ_o denote the generic incidence angle and the generic observation angle, respectively and g(∙) is the free-space Green’s function. For each θ_i, E_inc,i is the incident field in D, J is the contrast source [5] induced by this latter in D, E_inc,e and E_tot are the incident and total field on Γ_o, respectively and E_s is the scattered field from the unknown object. A_imod and A_e denote the integral operators relating the contrast source to the scattered fields outside and inside the domain under test, respectively, while the auxiliary function p is defined as [19]:

When considering the inverse scattering problem in the framework of the CS-EB model given by equations (1.a) and (1.b), the unknowns quantities are then given by the contrast source J and the auxiliary function p, in which the “usual” contrast unknown χ is embedded. It is worth to remark that, despite the CS-EB model is just a simple rewriting of the traditional Contrast Source model [5], it has proved to be a more effective tool to formulate and solve both forward and inverse scattering problems [19]. Notably, while its derivation was inspired by some mathematical and physical considerations related to presence of losses in the host medium and/or in the targets [19,21], processing of experimental data has shown that accurate and reliable results can be achieved also for lossless inhomogeneous targets in free space [20].

Traditionally, in standard inverse scattering problems, one assumes the knowledge of the total fields in both amplitude and phase. As incident fields are supposed to be known, this is equivalent to the knowledge of the scattered fields for each illumination. The problem we want herein to solve consists instead in retrieving the dielectric characteristics of a region under test from measurements of square amplitude distribution of the total and incident fields on the measurements domain Γ_o. Accordingly, no phase measurements are required.

Formally, our aim is to determine the contrast function χ(r), or equivalently the auxiliary function p(r), starting from the knowledge of the square amplitude distribution of the total and incident field on Γ_o. In order to do that, under actual conditions, we have first to come to a convenient finite dimensional approximation of the problem. This approximation should be such to lead to a finite number of experiments and measurements, while possibly not excluding any useful information. In this respect, it proves fruitful to recall properties and possible representations of both scattered and incident fields, which will also allow to develop and discuss the three steps procedure outlined in the Introduction.

With reference to the geometry of Fig. 1, it can be proved that the scattered field corresponding to a given incident field can be accurately represented with a number of Fourier harmonics given by 2βa, a being the radius of the minimum circle enclosing the target [17,22]. As a Fourier series can be turned into a Dirichlet sampling series [17], 2βa samples uniformly spaced in angle accurately represent each scattered field as well. By using reciprocity, one can also prove that the number of non-superdirective independent incident fields impinging on the domain under test is 2βa as well [23]. Hence, by excluding superdirective sources, 2βa plane waves uniformly spaced in angle form a complete family of independent incident fields. Therefore, as a function of θ_i and θ_o, the scattered field can be accurately represented by a number of samples given by (2βa)×(2βa)=(2βa)². Note that, by virtue of reciprocity, only one half of these samples is actually independent [23].

As far as the incident fields as measured on Γ_o are concerned, a different result holds true. In fact, by paralleling the above reasoning on the representation of the incident field in D, one can prove that each incident field on Γ_o can be accurately represented by 2βb Dirichlet samples (b being the radius of the circle Γ_o, see Fig. 1), and that 2βb (non-superdirective) independent incident fields (constituted by plane waves uniformly spaced in angle) exist therein. Therefore, as discussed for the scattered field, the incident field on Γ_o as a function of both θ_i and θ_o can be accurately represented by a number of samples given by (2βb)×(2βb)=(2βb)². Note that also in this case only one half of these samples is actually independent by virtue of reciprocity [23].

Let us explicitly observe that, when considering the square amplitude patterns of the above fields, the number of samples required for a faithful representation becomes four times larger (with respect to amplitude and phase measurements) as the sampling step has to be halved along each of the two coordinates.

3. A three steps “phaseless” inverse scattering approach

3.1 First step: estimating the incident filed

The first step of the approach aims at estimating the incident field given its amplitude distribution on the measurement domain (and possibly in some additional points). Note that this information plays a key role in both the following steps, as the knowledge of the incident field on the measurement domain is required to estimate the scattered field, while the knowledge of the incident field within the domain under test is needed in the inverse scattering procedure.

First, let us assume that the primary source adopted to generate the different incident field is always the same for each fixed θ_i, so that a single estimation problem has to be dealt with. Note this assumption entails no loss of generality, as 2βb different estimation problems of the same kind would have instead to be solved in case of different primary sources. Then, by exploiting an expansion into cylindrical harmonics, the field radiated by this primary source (i.e., the incident field) can be accurately represented as:

wherein ŕ and $\acute{\theta }$ are the polar coordinates of the generic observation points in a reference system centred on the primary source, and

wherein c is the radius of the minimum circle enclosing the source, and ΔN controls the approximation error of the representation (3), which decays exponentially fast with increasing ΔN [24].

Since, provided ŕ is larger than half a wavelength [24], expansion (3) is valid for any pair (ŕ,$\acute{\theta }$ ), from the knowledge of the {a_n} coefficients we can obtain the required information, i.e., the incident field in both amplitude and phase on the measurement domain Γ_o as well as within the domain under test D. Hence, the problem at hand is now that of estimating the 2P+1 coefficients {a_n, n=-P, …,-1,0,1, …P} from only amplitude measurements of the incident field itself on Γ_o (or part of it).

The above problem can be regarded as a PR problem, as the knowledge of the incident field phase and the knowledge of the coefficients {a_n} are each other equivalent.

PR problems have received great deal of attention in different fields, as testified by the many contributions present in the literature [14-17]. In particular, it has been shown that PR of radiated (or scattered) fields under actual measurement conditions is indeed possible [25-27]. These results rely on the circumstance that (radiated or scattered) fields are essentially band-limited functions with respect to observation variables and, as a consequence, their squared amplitude distributions are band-limited as well. Therefore, provided data are properly sampled, there is no loss of information (but for very negligible approximations comparable to noise) associated to the fact that data are sampled. Then, the presence of noise on data can be conveniently dealt with by defining a suitable ‘generalized solution’ as the global minimum of a cost functional [25]. Accordingly, the coefficients of the representation (3) can be determined by iteratively minimizing the squared distance between the computed square amplitude distribution of the incident field on Γ_o and the measured one M ²(r, θ), i.e.,

However, from the huge amount of results available on 1D PR problems, it is known that the problem of estimating the coefficients of (3) has not generally a unique solution when using a single measurement surface [12]. Therefore the problem arises to achieve a “convexification” of (5), meant to lead the optimization process to the actual solution, i.e., the actual phase distribution of the incident field.

In order to overcome this not trivial difficulty, two different strategies can be devised. The first one consists in providing a reasonable starting guess for the optimization process and/or to exploit a priori information about the primary source used in the experimental setup. Of course, this “software” strategy strictly depends on the availability of a ‘off line’ characterization of the source as much reliable and accurate as possible. Alternatively, a more robust “hardware” strategy can be devised in which two sets of independent measurements of the amplitude distribution of the incident field are exploited [12]. Hence, the PR problem is solved without the need of a priori information, [but for the size and shape of the source, which allow to fixing the number of terms to be retained in (3)]. On the other hand, in order to achieve independent information, the two scanning surfaces should be sufficiently spaced in terms of the wavelength and their distance must be known with wavelength accuracy [12].

This latter requirement, which may seem problematic, is indeed not surprising if one takes into account that a sub-wavelength reconstruction accuracy is required at the end of the overall reconstruction process.

As a final comment, let us explicitly note that, while the present paper is devoted to 2D inverse scattering problems, thus dealing with 1D phase retrieval problems in this first step, actual 3D inverse scattering problems, which are the final goal, will require to deal with 2D phase retrieval problems, wherein convexification of functional (5) is much simpler [12].

3.2 Second step: estimating the amplitude and phase of the scattered field

The second step of the procedure is aimed at retrieving the amplitude and phase of E_s from the measurements of ∣E_tot(θ_o,θ_i)∣². As ∣E_tot∣² = ∣E_inc,e∣² + ∣E_s∣² + 2ℜe[E_s(E_inc,e)*], where ℜe[°] denotes the real part of its argument, and the incident field E_inc,e is now available on the observation domain Γ_o (in amplitude and phase), the problem can be pursued by simply imposing the equality amongst the measured values of the square amplitude of the total field and the one corresponding to the current value of the scattered field. While enforcing such an equality, it is convenient that both terms are deprived from ∣E_inc,e∣² [1], as this latter choice allows to avoid the oversampling which would result in the discretization of the initial equation because of the properties of ∣E_inc,e∣² Accordingly (see also [1]), by defining the quadratic operator B as:

the problem of estimating E_s is formulated as the inversion of the nonlinear equation:

wherein ∣E_tot(θ_o,θ_i)∣² and ∣E_inc,e(θ_o,θ_i)∣² are the measured amplitudes of the total and incident fields, respectively. Then, one can estimate the complex scattered field by minimizing a functional given by the norm of the discrepancy between the measured amplitude distribution pattern, given by the right hand side of (7) and the one calculated using (6). In particular, the unknown scattered field is represented by means of a finite dimensional expansion (Fourier harmonics in the following) and the coefficients of the representation become the unknowns of the problem in the second step.

By addressing the reader to [1,2] for a more detailed discussion, let us here recall the differences of the above estimation problem with the one faced in the first step. First, note that in this case both the unknown, i.e., E_s(θ_o,θ_i), and the data, i.e., ∣E_tot(θ_o,θ_i)∣², of the problem are functions of two variables (i.e., θ_i and θ_o). Therefore, this latter procedure is eventually related to 2D PR problems which, unlike 1D ones, have got a unique solution, but for a zero measure set of cases and trivial ambiguities [12]. In this respect, as discussed with more details in [1], it is important to remark that a key point in the correct setting of the scattered field retrieval procedure is indeed that of correctly sampling the (possibly noisy) data. When these sampling rules are fulfilled, and a generalized solution has been defined in order to take into account noise on data, one can expect that, within the required accuracy, the continuous problem and its discrete counterpart are equivalent (as discussed at the end of previous section), so that only the possible occurrence of false solutions has to be avoided [1]¹.

As a second obvious distinction, note problem (7) differs from the usual PR one [12] as the operator B to be inverted has both a non-linear term (due to the square amplitude of the scattered field, i.e. the unknown of the problem) and a linear one (given by the interference between the incident field and the unknown scattered field). As a consequence, whenever the intensity of the scattered field does not overcome that of the incident one, the linear term is dominant in the operator to be inverted, so that, by using this information, solution of (7) is relatively simple. Unfortunately, this hypothesis is not generally satisfied, so that the problem at hand is non-linear. Therefore, besides requiring some form of regularization to tackle ill-posedness, the inversion of operator B is subject to possible false solutions corresponding to local minima of the cost functional whose global minimum defines the solution. Luckily, a second difference with respect to standard PR problems comes into play. As a matter of fact, the interference between the incident and scattered fields, by virtue of previously discussed properties, makes the essential dimension of the space of the data larger as compared to the case wherein only the square amplitude of the unknown function is given [1,2]. This circumstance comes out to be of the outmost importance, as it has been proved that in solving quadratic inverse problems of this kind false solutions can be avoided provided that a sufficiently large ratio between the number of (linearly) independent data and the number of unknowns is available [12]. Accordingly, the presence of a linear term (the interference) is of help in successfully reconstructing the unknown scattered field, as by properly choosing the incident field one can greatly enlarge the number of independent data, thus avoiding at all any false solution problem.

For the geometry at hand (see Fig. 1), implications arising from this latter circumstance have been extensively discussed in [1], leading to the following important results:

Therefore, provided the above rules are fulfilled, solution of the inverse problem (7) can be safely performed, i.e., avoiding occurrence of false solutions. Note that, by properly taking into account the different geometry, similar kind of rules can be derived also in the case of aspect limited data [2].

In summary, the main idea of the method proposed in [1] and adopted in the following is indeed that of taking advantage of the above mentioned sampling criteria together with results on quadratic inverse problems to determine which the proper number of data which have to be collected. As shown in [1,2] and in the results of the present paper, it seems that the proposed method (and the corresponding rules) allow robust and reliable reconstructions of the scattered field amplitude and phase patterns.

3.3 Third step: estimating the permittivity profile

Given the scattered field E_s and incident field E_inc,i estimated in the previous steps, the last step of the procedure consists in the reconstruction of the unknown dielectric profile.

As mentioned in the introduction, in this step we make use of the CS-EB inversion method [19,20]. This method, which has proved to be an effective and robust tool to solve inverse scattering problems in their full-nonlinearity [19,20], consists of two parts. First, one aims at determining the auxiliary function p from the (reconstructed) scattered field data, then, from the knowledge of this latter one achieves the dielectric profile, i.e., the contrast function χ, by inverting point-wise equation (2).

As far as the first part of the inversion method is concerned, this is formulated as the problem of finding, for a fixed frequency, the global minimum of the cost functional:

wherein the first two terms directly descend from the CS-EB formulation eqs. (1.a) and (1.b), while the third one is the Tichonov regularization [28], whose weighting parameter τ is iteratively updated during the minimization process, being in each iteration equal to the value of the cost functional normalized to the mesh dimension. In this way, the effect of the regularization term is more relevant at the beginning of the minimization process, while its contribution decreases for increasing iterations. Note that, due to the large number of unknowns involved in (8), the minimization is performed using a CG-FFT based scheme [19,20].

In (8), N_v is the total number of illuminations, i.e., the number of scattering experiments. Of course, this latter should be such to collect all the available information about the targets, while being non redundant. From the field properties recalled in Section 2, it immediately descends that 2βa incident fields (plane waves or cylindrical waves) have to be considered at least.

4. Testing on experimental data

The performances of the proposed three-step phaseless approach have been tested using the experimental data provided by the Institute Fresnel of Marseille [18]. In these experiments, measurements are collected under an aspect limited configuration in which, for each position of the primary source θs∈[0°,360°], measurements are gathered over an open arc θo∈[θs + 60°,θs + 300°]. In particular, we have considered the FoamDielIntTM.exp data-set [18] which is concerned with an inhomogeneous target embedded in free space. The target is made of two purely dielectric cylinders, the outer one of relative permittivity 1.45±0.1 and radius 0.04m and the inner one (which is completely embedded in the first one) of relative permittivity 3±0.3 and radius 0.015m. At the working frequency of 4GHz, a square investigated region of side 0.2062m, subdivided in 50×50 pixels has been assumed. The overall number of scattering experiments is N_v=72 and M=61 measurements are collected for each one of them.

As the data-base provides amplitude and phase of the total and incident fields as measured over the observation domain, we can conveniently exploit it to check the accuracy of all the three steps of the proposed approach.

In the first step, by using the measured square amplitude distribution of the incident field, we have solved the problem described in Section 3.1. Since the incident field is the same for all the scattering experiments [18], it is sufficient to solve the estimation problem for the case θs=0° and θo∈[60°, 300°]. The incident field has been modelled using the expansion (3), with P=7 and the pertaining coefficients have been determined by minimizing the cost functional (5). The minimization has been initialized by assuming a field having the measured amplitude and a constant (zero) phase on the measurement surface, i.e., from the blue-dotted pattern of Fig. 2. In Fig. 2 we compare the real (Fig. 2.a) and the imaginary parts (Fig. 2.b) of the measured incident field (red dotted-line) provided by the data-set, the estimated one (green solid-line) and the starting guess (blue dotted-line). As it can be observed, notwithstanding the starting guess is quite different from the actual pattern, a good reconstruction is achieved. As expected, a slightly worse estimation is obtained at the end of the observation arc due the truncation of the measurement domain [17]. Also note that similar results can be obtained by using the field radiated by a line source as starting guess.

On the other side, it has to be noted that other starting guesses (such as for example completely random patterns) do not allow a faithful reconstruction of the incident field. However, a large number of numerical simulations (measured data were not available) has shown that by using two sufficiently different measurement surfaces (which can also be achieved by moving the source) the incident field is accurately retrieved even starting from a completely random initial point (see also [12]).

 

Fig. 2. Real (a) and imaginary (b) part of the measured incident field (red-dotted line), the estimated one (green-solid line) and the starting guess (blue-dotted line).

Download Full Size | PDF

In the second step, by using the estimated incident field as computed on Γ_o through (3) and the measured square amplitude distribution of the total field, we have solved the problem (7). The complex scattered field has been modelled as a superposition of Fourier harmonics [1], whose number is related to the electrical dimension of the investigated domain [17]. The coefficients of this representation [i.e., the actual unknowns of problem (7)] are evaluated by minimizing the squared distance between the actual total field and the experimental one, as measured on the 240° arc. A random distribution for the unknowns has been used as starting guess in the minimization procedure.

From the comparison between the retrieved scattered field and the actual one (which is available in the considered data set) for all the illuminations, it can be observed that a good reconstruction for both the amplitude (see Fig. 3) and phase (see Fig. 4) distribution is achieved, although, similarly as before, a slightly worse reconstruction is achieved at the end of the observation arc, due to the truncation of the measurement domain [2,17].

It is interesting to note that since the radius of circle wherein the receivers are located is r_m=1.765m and b_cr=0.2113m, the condition r_m > b_cr [1] holds. In agreement with the analysis in [1], this condition indeed prevents local minima occurrence. In the minimization process the scattered field is represented through a 2D Fourier expansion where 9x9 coefficients represent the actual unknowns of the problem.

 

Fig. 3. Amplitude distribution of the measured scattered field (a) and of the estimated one (b)

Download Full Size | PDF

 

Fig. 4. Phase distribution of the measured scattered field (a) and of the estimated one (b)

Download Full Size | PDF

The third and final step of the proposed phaseless imaging technique deals with a standard inverse scattering problem. Of course, the estimated scattered and incident fields are necessary as input data in this step. The “background” (better to say, an initial contrast function slightly different from zero) has been used as starting guess in the minimization of (8). The very accurate reconstruction of the real part of the contrast function which is achieved is shown in Fig. 5.a. Note that the corresponding imaginary part has not been displayed as, in agreement with the lossless nature of the targets, this latter is indeed negligible with respect to the real part. The total number of iterations required in the minimization process is 176 and the final value of the cost functional is 4.54e-3. The numerical procedure is stopped when the difference between the previous value of the cost functional (8) and the actual one is less than 1e-4. The maximum value of the estimated contrast function is 1.84, which is within the measurement accuracy [18]. The overall computational time is about of 10 minutes with an AMD-Athlon 64 800MHz Processor.

As final remark to the above results, note that the considered experimental data-set is known in the literature as a benchmark for non-linear inverse scattering methods, as linearized approaches are understood to fail when applied to it [18].

It is worth to note that a comparable reconstruction is actually obtained when using the measured (amplitude and phase) scattered fields, thus confirming the possibility of performing a faithful phaseless quantitative imaging without any loss of accuracy.

It has also to be remarked that, while the Marseille data are usually elaborated by using multi-frequency (i.e., broad-band) data [18], the accurate results achieved herein only rely on monochromatic data. Of course, use of multi-frequency information can further improve the final result also in this case. To show this, we have considered data at 4GHz and 8GHz. In particular, after applying the three-step procedure to the 4GHz data, we have repeated the first two steps (i.e., estimating the incident and scattered fields, respectively) with the 8GHz data and then we have tackled the final inverse scattering step exploiting a frequency hopping scheme [20] in which the reconstructed contrast at 4GHz has been used as a starting guess for the inversion at 8GHz. In this latter case, the total number of iterations is 92, the final value of the cost functional is 1.02e-2 and the final maximum value of the contrast function is 1.89. As it can be observed from Fig. 5.b a further improvement of the reconstruction is achieved as the contrast now appears to be sharper, both in terms of maximum value and shape characterization. Also note that the position of the inner cylinder is much more accurate.

 

Fig. 5. Real part of the estimated contrast function at f=4GHz (a), and f=8GHz (b)

Download Full Size | PDF

5. Conclusions

A new three-step imaging technique based on phaseless measurements has been presented and tested against experimental data. The good results achieved with experimental data commonly adopted to benchmark inverse scattering procedures have fully confirmed the possibility of performing a faithful phaseless nonlinear imaging by properly combining different “know-how”, ranging from results on phase retrieval problems to advanced electromagnetic inverse scattering techniques. It is interesting to remark that the processing required to extract the information needed by the inverse scattering procedures (i.e., the incident and the scattered field complex distributions) represents a negligible part of the overall computational time.