1.

Introduction

During the past few years there has been growing interest in the use of near-infrared light for biological tissue diagnostics (see a review in Ref. 1). One of the numerous applications of light imaging is the visualization of electrical activity in the heart. ^{2, 3, 4, 5, 6, 7, 8, 9} The commonly applied optical technique makes use of voltage-sensitive dyes, which allow rapid changes to be observed in the transmembrane potential of cardiomyocytes. The technique takes advantage of the dyes’ fluorescence or absorption properties when excited by an external source, and could provide invaluable information on the 3-D organization of the heart’s electrical activity during arrhythmias.

Until today, optical mapping of electrical waves in cardiac tissue has yielded 2-D images obtained at the tissue boundaries. Reconstruction of the 3-D activity is limited by several factors. First, the thickness of the myocardial wall $(>1\mathrm{cm})$ prevents the use of confocal microscopy. Second, the electrical activity is not static and excitation waves can have propagation velocities up to $0.5\mathrm{m}∕\mathrm{s}$ , requiring measurements at a time scale no larger than a few milliseconds. This puts severe temporal constraints on the acquisition equipment and complicates the application of conventional diffuse optical tomography (DOT) (see Refs. 10, 11 for reviews on DOT). The most advanced functional DOT systems used for 3-D imaging of brain activity are still too slow for imaging cardiac excitation.^{12, 13, 14}

We propose a method for the localization of an intramural source of electrical activation from pairs of 2-D optical images obtained in reflection and transillumination mode (see Fig. 1 ). We consider an experimental setting described by Baxter, ⁶ where a slab of tissue is placed between two charge-coupled device (CCD) cameras, measuring fluorescence or absorption from both tissue surfaces in rapid succession. We show that signals derived from these measurements allow for the calculation of the exact position of a single point source. The method can be applied to expanding wavefronts. Specifically, we demonstrate that in addition to exact 3-D localization of the excitation center, we can also estimate the distance of the excitation front from the outer tissue boundaries at any given moment of time. Contemporary specialized cooled CCD cameras can readily reach frame rates as high as $2000\mathrm{Hz}$ , which should be sufficient for the experimental implementation of our approach.

Fig. 1

Experimental setting for the 3-D localization of electrical activity in a slab of cardiac tissue. The external light source illuminates a slab surface, inducing voltage-dependent optical signals. The latter is recorded with two CCD cameras: one is aimed at the illuminated surface (reflection mode) and the other at the opposite surface (transillumination mode). (a) shows epicardial illumination, while (b) illustrates endocardial illumination.

One of the potential applications of the proposed method is the localization of so-called ectopic foci—high-frequency sources of electrical excitation responsible for initiation and maintenance of cardiac arrhythmias.¹⁵ The treatment of such arrhythmias is based on radiofrequency ablation of the source.^{16, 17, 18} It requires the knowledge of the exact position of the source and is less concerned about the 3-D shape of the excitation fronts emanating from the source.

Our method is derived from analytical solutions of the light transport equation in the diffusion approximation,^{10, 19} and is validated computationally. Here we assume Robin-type boundary conditions, which represent physical effects at the boundary and take into account Fresnel reflection due to refractive index mismatch. We use the method of images and a configuration involving four source-image pairs to construct a Green’s function in the slab geometry.^{20, 21, 22, 23, 24, 25} The method was tested with both absorptive and fluorescent voltage-sensitive dyes.

The work is organized as follows. In Sec. 2 we derive an analytical formula for the depth estimation of single electrical point sources based on the ratio of total integrated signals in the reflection and transillumination modes of observation. In Sec. 3 we investigate the effects of noise on the depth estimation. In Sec. 4 we consider illumination that alternates between both sides of the slab to improve the accuracy of the proposed method. In Sec. 5 we assess the method’s accuracy for the localization of excitation centers that produce ellipsoidal waves. The technical details of the derivations are presented in Appendices A, B, and C in Secs. 7, 8, 9.

2.

Localization of a Single Point Source

2.2.

Estimation of the Lateral Coordinates

The function $\Gamma (Z,\rho )\bullet w$ describing the photon flux through the surface in Eq. 5 is bell shaped. The coordinates of its maximum correspond to the lateral coordinates $X$ and $Y$ of the source. Experimentally, $X$ and $Y$ can be readily estimated by finding the points with maximal signal on either epicardial or endocardial images (see Fig. 2 and its discussion).

Fig. 2

Simulated fluorescence images produced by a single-point source. (a) and (b) show the epicardial and endocardial images obtained from a point source located $3\mathrm{mm}$ from the illuminated epicardial surface. (c) and (d) show the epicardial and endocardial images obtained from a point source located $7\mathrm{mm}$ from the illuminated epicardial surface. (e) and (f) show the same images (source at $7\text{-}\mathrm{mm}$ depth) when illuminating the endocardial surface. The profiles through the middle of each image are also shown. The insets on top illustrate the geometry in a cross section of the slab: the black circle indicates the position of the source, and the arrows indicate the illumination. The lateral coordinates $(X,Y)$ can be inferred from the image, as shown in (a).

3.

Source Localization in the Presence of Noise

Equation 7 allows for the accurate depth estimation of any source location. Similarly, the $x$ and $y$ coordinates are the coordinates of the maximal signal. In the presence of noise, however, the applicability of the proposed approach is determined by the signal-to-noise ratio, which in turn depends on the depth of the source. Let us start with a computational example that shows simulated raw epicardial and endocardial signals obtained for two different depths of the source (see Fig. 2). To simulate the effects of noise, we added a matrix of normally distributed random numbers to the matrix of computed pixel signals. All matrices were sampled from the same distribution with the mean of zero, so that all the noise values were independent and identically distributed. The absolute level of noise, which was defined as the standard deviation of the distribution, was fixed at ${10}^{-7}$ .

Unlike the ideal case described in the previous section, in the presence of noise, the absolute values of signal amplitude become important. Accordingly, we have to take into account the characteristics of the excitation light. In our computational examples, we assume that excitation of fluorescent dyes occurs at $520\mathrm{nm}$ and we used ${\delta }_e=0.8\mathrm{mm}$ , characteristic for this wavelength.⁹ Again, we assume that all recordings, utilizing both fluorescent and absorptive dyes, were made at $650\mathrm{nm}$ , for which the attenuation length $\delta =1.35\mathrm{mm}$ . The extrapolation distance was set at $d=1\mathrm{mm}$ , which corresponds to $D=0.22\mathrm{mm}$ and $\alpha =1.5$ . The quantum yield $\beta$ of both types of dyes was assumed the same.

Figures 2a and 2b show the case when the source is located at $3\mathrm{mm}$ from the epicardial (illuminated) surface. The epicardial signal Fig. 2a is quite strong, and the position $(X,Y)$ of the source in the epicardial plane can be readily detected from the signal profile shown at the bottom of the panel. The same is true for the endocardial signal Fig. 2b obtained in transillumination mode. Although the signal is significantly noisier due to a relatively weaker signal, one can still clearly see the source. The estimated $x$ and $y$ coordinates of the source obtained from Fig. 2b by picking the ten brightest pixels and averaging their coordinates were $X=0.01\pm 0.23\mathrm{mm}$ and $Y=0.01\pm 0.23\mathrm{mm}$ , respectively, very close to the actual coordinates $X=0.0$ , $Y=0.0$ of the source. The average depth, which was estimated by calculating the total signals on both surfaces and subsequently using Eq. 7, is also very close $\left(2.98\mathrm{mm}\right)$ to its actual location $Z=3\mathrm{mm}$ .

The situation changes, however, when the source is located further away from the illuminated surface. Figures 2c and 2d show the case of a source located $7\mathrm{mm}$ from the illuminated epicardial surface. Now, the source can no longer be identified on the epicardial surface and can only be seen in the transillumination mode Fig. 2d. The effect of noise as a function of source depth is illustrated in Table 1 (see columns 3 and 4). The accuracy is higher when the source is located near the epicardial (illuminated) surface and gradually deteriorates as the source moves toward the endocardial surface.

The observed trend is the result of the overall reduction of the surface signal with increasing $Z$ . When the source moves away from the illuminated boundary, the random component of the signals recorded from both surfaces becomes dominant. As a result, the average ratio of integrals approaches 1 (the noise is identical for both signals), yielding the estimated depth close to $L∕2=5\mathrm{mm}$ . Indeed, one can see that at $Z=5\mathrm{mm}$ , the estimated depth reaches $5\mathrm{mm}$ and practically does not change as $Z$ increases all the way up to $9\mathrm{mm}$ .

Such underestimation is less pronounced in the case of absorptive dyes (column 4): the largest SD was $0.67\mathrm{mm}$ for a source located at $Z=9.0\mathrm{mm}$ (compare to $1.05\mathrm{mm}$ in the case of fluorescent dyes). This is the result of the longer illumination wavelength in the case of absorptive dyes (650 versus $520\mathrm{nm}$ ). At $650\mathrm{nm}$ , the penetration depth of light is larger (1.35 versus $0.8\mathrm{mm}$ , see earlier) and the signal is less attenuated. Consequently the effect of noise is less pronounced.

4.

Alternating Illumination

The accuracy of the depth detection in the presence of noise can be increased by utilizing alternating illumination when the direction of illumination switches between epicardium [see Fig. 1a] and endocardium [see Fig. 1b], producing two pairs of images. By combining the depths estimated from each pair, one can significantly reduce the depth dependence of the localization accuracy in the presence of noise.

The basic idea of alternating illumination can be readily understood from a simple example. From the previous section, we already know that our method is least accurate when the source is located far from the illuminated surface. One such example is shown in Figs. 2c and 2d. In this case, the epicardial signal [Fig. 2c] is buried in noise and the accuracy of depth detection is not very high. Figures 2e and 2f show the images of the same source after the illumination was switched to the other side of the slab. The situation has dramatically improved [compare Figs. 2c with 2e and 2d with 2f]. Now the signal is clearly identifiable on both surfaces. Note that Figs. 2e and 2f, if flipped, look almost identical to Figs. 2a and 2b, except for the noise components, which are different statistical samples with the same distribution. This resemblance results from the fact that the depth of $7\mathrm{mm}$ for endocardial illumination is equivalent to a depth of $3\mathrm{mm}$ for epicardial illumination (with the slab thickness $L=10\mathrm{mm}$ ).

By using alternating illumination, one can effectively achieve accurate depth detection throughout the whole thickness of the slab. Table 2 shows the estimated depths for image pairs obtained during epicardial $\left(Z_{\mathrm{epi}}\right)$ and endocardial $\left(Z_{\mathrm{endo}}\right)$ illumination. We define both quantities as distances from the epicardium, i.e., $Z_{\mathrm{endo}}=L-Z$ , with $Z$ given by Eq. 7 during endocardial illumination. One can see that $Z_{\mathrm{epi}}$ becomes less accurate as the distance from the epicardial surface increases (column 2), whereas $Z_{\mathrm{endo}}$ has an opposite trend. Qualitatively similar results were obtained for absorptive dyes (compare columns 4 and 5), although the accuracy shows less depth dependence.

Table 2

Effect of endocardial versus epicardial illumination. The absolute noise level was set to 5.00e−6 .

The algorithm of selection between $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ is quite straightforward. One must compare image pairs obtained during epicardial and endocardial illumination and select the dataset with the highest signal-to-noise ratio. If the signal-to-noise ratios of both pairs are comparable, the best accuracy can be achieved by calculating the source coordinates independently from each dataset, i.e., by averaging the results.

5.

Ellipsoidal Waves

Our next step was to determine how well our approach works in more realistic situations when, instead of a single point source, we have an expanding excitation front. To achieve this goal we have extended our model to account for the specific features of cardiac electrical propagation.

The excitation wave in the heart has a very sharp front and a well-defined plateau whose amplitude can be considered constant with a high degree of accuracy. Accordingly, in our simulations, the potential $V$ in Eq. 2 for each voxel was set to 1 (excited state) or 0 (resting state).

Cardiac tissue is anisotropic for electrical propagation: wavefronts have a larger conduction velocity along than across fibers. Moreover, the orientation of the fibers rotates through the heart wall.²⁷ Hence, point stimulation results in ellipsoidal waves and the orientation of their $x$ and $y$ axes changes as they propagate deeper in the tissue.

In our computational examples, the semi-axis ratio of the ellipsoid was 6:2:15 (in directions $x$ , $y$ , and $z$ , respectively). In each run we simulated six consecutive frames; three of them (frames 2, 4, and 6) are shown in Figs. 3 (fluorescence) and 4 (absorption). The number of voxels in the excited state varied from 30 (frame 1) to 1000 (frame 6). The wave is initiated in a voxel in frame 1 at depth $6.1\mathrm{mm}$ . Ellipsoidal waves produced considerably stronger signals than a single point source. To see the effects of noise in these cases, we increased the absolute noise level by a factor of 50. We repeated the simulations ten times for different random samples from the same noise distribution, and calculated the mean depths and SD.

Fig. 3

Fluorescence images of an expanding hollow ellipsoidal wavefront. The wave was initiated in one voxel at depth $6.1\mathrm{mm}$ from the epicardium. An $x\text{-}z$ cross section of the excited voxels is shown on the top of each frame. A black dotted line in one of the four source images indicates its actual $xy$ position.

Fig. 4

Absorption images of an expanding hollow ellipsoidal wavefront shown in Fig. 3. There is a significantly higher signal-to-noise ratio compared to Fig. 3 due to increased signal magnitude.

Figure 3 depicts the optical images obtained for an expanding hollow ellipsoidal wavefront originating near the endocardial surface in the case of fluorescent dyes. We illustrate different stages of the propagating wave by showing the signals for frames 2, 4, and 6. At an early stage (frame 2), the wave can only be distinguished on the endocardial surface, only during endocardial illumination. As the wave expands and moves closer to the epicardial surface, the signal starts to exceed the noise background in the other images: first in transillumination mode (frame 4) and then, eventually, in all four images (frame 6). The images of the excitation front in frames 4 and 6 acquire a pronounced elliptical shape. One can also see the counterclockwise rotation of the axes of the ellipsoid, simulating the effect of fiber rotation as described before.

Figure 4 shows the optical signals for the same ellipsoidal wave, but in the case of absorptive dyes. The results were qualitatively similar to the ones obtained with fluorescent dyes (Fig. 3). However, absorptive dyes give overall brighter signals, leading to less noticeable effects of noise (compare corresponding frames in Figs. 3 and 4).

Figure 5 shows the depths $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ estimated during epicardial and endocardial illumination in all six frames. As the front expands, the values of $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ steadily start to diverge. However, their average (solid line) remains nearly constant and accurately predicts the position of the excitation center, even at very late stages of the wavefront expansion. The maximal standard deviation (indicated by error bars) was observed at the earliest stages of the excitation, primarily due to the small volume of excited voxels and thus the lower magnitude of the signal. However, the error rapidly decreases as the front becomes larger. Computer simulations show the same effect for absorptive dyes.

Fig. 5

Accuracy of localizing the excitation center of the expanding hollow ellipsoidal wave using the method of total signal ratios. The actual position of the excitation center $\left(6\mathrm{mm}\right)$ is indicated by a dashed line, dotted lines show the calculated depths $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ , and the solid line shows the calculated position of the excitation center. The SD due to noise is indicated for each data point. (a) and (b) are for fluorescence and absorption imaging, respectively.

It is interesting that $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ predict not only the origin of the excitation front by their average, but also the $z$ coordinates of the excitation front above and below the excitation center. Figure 6 shows $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ as well as the $z$ coordinates of the front for fluorescent [Fig. 6a] and absorptive dyes [Fig. 6b] in the absence of noise in the case of hollow and solid ellipsoidal waves. In the fluorescent case, the deviation between the predicted and actual position of the front is less than $0.9\mathrm{mm}$ . In the absorption case, it was higher, up to $1.6\mathrm{mm}$ , yet still rather close.

Fig. 6

Comparison of $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ with the actual positions of the front above and beneath the excitation center. The real positions are indicated by solid lines, dashed lines show the calculated positions for the hollow wave, and dotted lines are for the calculated positions for the solid wave. (a) and (b) are for fluorescence and absorption imaging, respectively. (c) and (d) show the results after the depth estimates $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ were corrected for the form factor of the source distribution (see Appendix C in Sec. 9).

The difference between the actual front coordinates and the respective values of $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ represents a systematic error resulting from the fact that the pairs of images used to calculate $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ represent selectively the portions of the excitation wave on the epicardial and endocardial side, respectively. This systematic error can be expressed in terms of the so-called form factor, a geometric characteristic of the source configuration (see Appendix C in Sec. 9). The less accurate results obtained for “solid” waves in Figs. 6a and 6b are a manifestation of differences in $\Lambda$ for hollow and solid source configurations. Indeed, in the case of solid waves, the inner core makes a contribution to the signal additional to that from the front surface, resulting in the shift of the estimated values to the center of the core.

As shown in Appendix C in Sec. 9, the form factor is determined by the vertical size of the ellipsoid and does not depend on the semi-axes ratios. This property of the systematic error allows for the design of an iterative procedure to correct the values of the depths $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ (Appendix C in Sec. 9). Figures 6c and 6d show the results of such correction. We see that, after correction, the dotted and solid lines match perfectly, i.e., $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ accurately predict the distance of the wavefront to the epicardial and endocardial surfaces, respectively. The corrected depths obtained in the worst case, i.e., in frame 6 with the largest ellipsoid, differ from the actual depths by less than 0.6%.

A qualitative understanding of the physical meaning of $Z_{\mathrm{epi}}$ and $Z_{\mathrm{endo}}$ for source configurations extended in the $z$ direction can be achieved by analyzing the case of two point sources located at two different depths, $Z_A$ and $Z_B$ , as illustrated in Fig. 7 . Indeed, during epicardial illumination, both the epicardial (reflect-epi) and the endocardial (trans-endo) images show only source A located near the epicardial surface. Source B, located near the endocardium, becomes visible only after the illumination is switched to the endocardium (reflect-endo and trans-epi). Note that the depth selectivity is more pronounced in the case of fluorescent dyes (not shown). An important factor that determines the higher depth selectivity in this case is the difference in the attenuation lengths of the excitation and emission light.⁹

Fig. 7

Simulated fluorescence images produced by two point sources. Reflection and transillumination signals were obtained by illuminating either the epicardium (reflect-epi and trans-epi) or endocardium (reflect-endo and trans-endo). The depths $Z_A$ and $Z_B$ of the sources were respectively 2.0 and $8.0\mathrm{mm}$ from the epicardium (see inset on top).

6.

Concluding Remarks

We are proposing a method for the localization of sources of electrical activity inside the myocardial wall. It is based on ratios of the total photon flux signals in reflection and transillumination. We demonstrate that even with a considerable level of noise, the method provides 3-D information on the source position and size. In general, we found that the depth estimate was more accurate when using fluorescent rather than absorptive dyes. However, absorptive dyes provided brighter signals and made the depth estimation method less sensitive to noise and more attractive for experiments.

We show that by using alternating illumination, one can significantly increase the accuracy of the proposed method. Importantly, it can be used not only to localize the intramural sources of excitation, but also to study qualitatively the dynamics of the wavefront expansion. It should be noted, however, that the experimental implementation of this approach might require rather high acquisition speeds. Specifically, to achieve an error due to motion-related smearing of less than $1\mathrm{mm}$ , the whole measurement cycle, including epicardial and endocardial illumination, should be completed in $2\mathrm{ms}$ (in view of the propagation speed in the heart of about $0.5\mathrm{m}∕\mathrm{s}$ ). However, this is well within the reach of contemporary CCD cameras that can readily achieve frame rates of $2000\text{frames}∕\mathrm{s}$ .

Our theoretical approach is based on previous studies of the forward problem of light propagation in cardiac tissue.^{9, 26} A major improvement added to the present study consists in using the more accurate Robin boundary condition. The latter takes Fresnel reflection at the tissue boundaries into account, as opposed to the previously used zero-boundary conditions.

The present work represents the first step in exploring the feasibility of a 3-D visualization of electric activity inside the cardiac muscle. We focus on the transmural depth of sources of electrical activity and do not attempt to reconstruct the amplitude or 3-D distribution of the transmembrane voltage $V$ . It should be noted, however, that just knowing the location of sources of electrical activity in cardiac tissue might have important applications in cardiac electrophysiology as a tool to understand the mechanisms underlying life-threatening arrhythmias, and toward the development of new antiarrhythmic procedures.

7.

Appendix A: Derivation of a Point-Spread Function in a Thick-Slab Approximation

7.3.

Point-Spread Function in Thick-Slab Approximation

For the slab geometry with the observation surface at $z=0$ and another boundary at $z=L$ , we assume that the slab thickness $L$ is much greater than the attenuation length $\delta$ . Strictly speaking, in the slab geometry, the method of images requires an infinite sequence of mirror sources. However, when $\delta$ is much smaller than $L$ , the contributions from distant image sources decrease exponentially. Thus, to satisfy the boundary condition at $z=L$ , it is found to be sufficient to add two more images to a solution satisfying a boundary condition at $z=0$ (see Fig. 8 ). The function $[{\Phi }_0\left(r_1\right)-{\Phi }_0\left(r_2\right)]-[{\Phi }_0\left(r_3\right)-{\Phi }_0\left(r_4\right)]$ with $r_1={[{(Z-z)}^2+{\rho }^2]}^{1∕2}$ , $r_2={[{(Z+z)}^2+{\rho }^2]}^{1∕2}$ , $r_3={[{(2L-Z-z)}^2+{\rho }^2]}^{1∕2}$ , and $r_4={[{(2L-Z+z)}^2+{\rho }^2]}^{1∕2}$ for the four sources satisfies the condition $\Phi =0$ at $z=0$ . This function also approximately satisfies the condition $\Phi =0$ at the $z=L$ boundary, since the contribution to the photon density from the real source 1 cancels the contribution from source 3, and the contributions from sources 2 and 4 at this border are exponentially small [less than $\exp (-L∕\delta )$ ]. The photon flux (point-spread function) at the observation surface $z=0$ is then given by:

Eq. 13

$$\Gamma (Z,\rho )={\Gamma }_0(Z,\rho )-{\Gamma }_0(2L-Z,\rho ).$$

For the total flux at this surface, we find $J$ $\left(Z\right)=\exp (-Z∕\delta )-\exp [-(2L-Z)∕\delta ]$ . Replacing $Z$ with $L-Z$ gives the signal on the surface $z=L$ .

Fig. 8

The method of images for the solution in slab geometry. The real source located at depth $Z$ inside the tissue (hatched region) is represented by a filled circle. The mirror images are shown as open circles.

We obtain solutions satisfying the Robin boundary conditions by simply substituting $L+2d$ for $L$ (for identical conditions at the two boundaries) and $Z+d$ for $Z$ :

Eq. 14

$$\Gamma (Z,\rho )={\Gamma }_0(Z+d,\rho )-{\Gamma }_0(2L+3d-Z,\rho ),$$

Eq. 15

$$J\left(Z\right)=\exp [-(Z+d)∕\delta ]-\exp [-(2L+3d-Z)∕\delta ]=2\exp [-(L+2d)∕\delta ]\sinh \left(\frac{L-Z+d}{\delta }\right).$$

Here, the function ${\Gamma }_0(\zeta ,\rho )$ is given by Eq. 11. Again, the signal on the opposite surface can be found by replacing $Z$ with $L-Z$ .

8.

Appendix B: Analytical Approximation of the Depth Estimate

By assuming that a point source is far from the boundaries ( $L-Z+d⪢\delta$ and $Z+d⪢\delta$ ), one can obtain from Eq. 7 in the main text an explicit expression for the source position:

9.

Appendix C: Form Factor and Its Calculation

In Sec. 5, we showed computationally that the ratio of total photon fluxes measured in reflection and transillumination modes can be used to estimate the minimal distance of the excitation front from the illuminated surface (see Fig. 6). The systematic error of such estimates can be expressed as a function of the so-called form factor $\Lambda$ , a geometric characteristic of distributed source configurations. Next, we introduce a general definition of $\Lambda$ and present an analytical derivation for ellipsoidal waves.