Reconstruction of optical vector-fields with applications in endoscopic imaging

We introduce a framework for the reconstruction of the amplitude, phase and polarisation of an optical vector-field using calibration measurements acquired by an imaging device with an unknown linear transformation. By incorporating effective regularisation terms, this new approach is able to recover an optical vector-field with respect to an arbitrary representation system, which may be different from the one used in calibration. In particular, it enables the recovery of an optical vector-field with respect to a Fourier basis, which is shown to yield indicative features of increased scattering associated with tissue abnormalities. We demonstrate the effectiveness of our approach using synthetic holographic images as well as biological tissue samples in an experimental setting where measurements of an optical vector-field are acquired by a fibre endoscope, and observe that indeed the recovered Fourier coefficients are useful in distinguishing healthy tissues from lesions in early stages of oesophageal cancer.



There are no comments yet.


page 13

page 15

page 16

page 17

page 18

page 19

page 21


Automatic Differentiation for All Photons Imaging to See Inside Volumetric Scattering Media

Imaging through dense scattering media - such as biological tissue, fog,...

Imaging dynamics beneath turbid media via parallelized single-photon detection

Noninvasive optical imaging through dynamic scattering media has numerou...

Tattoo tomography: Freehand 3D photoacoustic image reconstruction with an optical pattern

Purpose: Photoacoustic tomography (PAT) is a novel imaging technique tha...

GANPOP: Generative Adversarial Network Prediction of Optical Properties from Single Snapshot Wide-field Images

We present a deep learning framework for wide-field, content-aware estim...

Interferometry-based modal analysis with finite aperture effects

We analyze the effects of aperture finiteness on interferograms recorded...

Increasing Linear Dynamic Range of Commercial Digital Photocamera Used in Imaging Systems with Optical Coding

Methods of increasing linear optical dynamic range of commercial photoca...

Fourier reconstruction for diffraction tomography of an object rotated into arbitrary orientations

In this paper, we study the mathematical imaging problem of optical diff...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Recently, there has been a significant interest in developing new types of optical fibre endoscopes for medical imaging applications [15, 25, 36, 11, 7]. Typically, these new endoscopes aim to be thinner, and therefore less invasive, and/or use different properties of light than conventional white light endoscopes making them more sensitive for detecting diseases such as cancer [41]. A full optical vector-field reflected from a tissue consists of amplitude, phase and polarisation information. Phase and polarisation have recently shown promise as diagnostic indicators, but are discarded by conventional white light endoscopes which record amplitude information only. Phase is highly sensitive to surface scattering that arises due to microstructural tissue changes in early cancer, creating distorted reflected wavefronts [20, 37, 38, 34]. This effect has been utilised in phase contrast and quantitative phase microscopy to predict recurrence of prostate cancer [33]. Similarly, polarisation information can indicate the formation of dense collagen networks [6], and the concentration of other polarisation-sensitive compounds, such as glucose, linked with early cancer [24, 4]. This has found use in the diagnosis of colon [26, 3] and gastric cancers [39]. Currently, there are no commercial phase and polarisation endoscopes but a number of prototype devices have been demonstrated [32, 15, 27, 36, 40].

To achieve phase and polarisation imaging in fibre endoscopes, it is necessary to characterise the underlying transformation of the optical fibre. In realistic clinical settings, this transformation changes frequently due to bending and temperature fluctuations and it is therefore important that this characterisation is efficient and accurate. For the characterisation, typically a set of known fields that form some kind of a basis are input into one end of the fibre and the resulting outputs are recorded at the other end, a procedure termed calibration. The task then becomes to recover a representation of the optical field reflected from a tissue given the calibration measurements and the samples of the output field measured by an imaging sensor outside of the fibre.

In this paper, we investigate the following questions: (i) is there a particularly useful representation of the full optical field reflected from a tissue that can be used for detecting optical aberrations associated with early cancer, and (ii) how can such a representation be recovered by an efficient and reliable algorithm from raw endoscopic measurements, namely from the calibration measurements and the samples of the output optical field?

To address these questions, we show that a Fourier representation recovered directly from the raw measurements has the statistical power to distinguish healthy tissues from tumours, and we provide a general reconstruction framework that can perform such recovery efficiently and stably.

More concretely, after reviewing related previous works in Section 1.1, in Section 2 we introduce a general reconstruction framework for the recovery of a full two-dimensional complex vector-field, where different regularisation terms are permitted and the bases used for image representation and device calibration are allowed to be different and/or non-orthogonal. In Section 3 we demonstrate that it is possible to extract informative features for detecting cancer from images of simulated tissue samples by projecting them onto a Fourier basis and observing the decay of their respective Fourier coefficients. Finally, in Section 4, we apply our new approach to experimental data acquired using a custom-built fibre endoscope [21] and recover synthetic holographic images as well as images of mouse oesophageal tissue containing small tumours (lesions). In particular, by recovering images of a biological tissue with respect to a Fourier basis using -regularisation, we observe that the corresponding Fourier coefficients are indicative of differences between lesions and healthy tissues and demonstrate their potential for medical diagnostic applications. We conclude with a discussion of our results and directions for future research in Section 5.

1.1 Relation to previous work

Before precisely formulating our new approach in the next section, it is useful to give a brief overview of some existing reconstruction techniques used in the applications of interest and differentiate them from the reconstruction framework developed in this paper.

In imaging through optical fibres or other scattering media, typical recovery procedures use the same, finite-dimensional basis for calibration and image representation in conjunction with standard inversion techniques. They start by discretizing the mathematical operator of the fibre as a mapping between pixels at different ends of the fibre , leading to a transmission matrix which is then characterised through calibration. The individual calibration inputs are arranged into the columns of matrix and the corresponding calibration outputs into the columns of matrix . Most existing systems use a full orthogonal basis of the discretized input space as the calibration inputs, e.g. a set of tilted plane waves (a Fourier basis) [13] or a Hadamard basis generated using a phase-only spatial light modulator [28]. The orthogonality of such bases ensures that matrix is unitary. Then, by assuming that is also unitary, images can be recovered using phase conjugation. In this approach a (generalised) inverse of the transmission matrix is calculated as , where denotes the conjugate transpose, and a representation of with respect to the calibration inputs is recovered as [28, 15, 14]. Although simple and straightforward to compute, the unitary assumptions in this approach are typically violated in practice [12]. In the context of imaging through scattering media, the inversion of the transmission matrix was also performed through alternative approaches to phase-conjugation such as least-squares (Moore-Penrose pseudoinverse) or Tikhonov regularisation [29, 30]. In particular, -regularisation has not been previously explored in these applications.

When compared with these conventional techniques, we emphasise that our new framework is able to recover a representation of the unknown optical field with respect to any particular infinite-dimensional basis which is allowed to be different from the one used for calibration, directly from the raw measurements. If an image representation with respect to a particular basis (such as Fourier) is desired, alternatively to our new approach one could in principle use the conventional techniques to recover an approximation to such a representation as we now describe. One could calibrate the fibre with respect to a Fourier basis and use standard recovery techniques to reconstruct images with respect to the same basis. However, in high resolution imaging, calibration with respect to a Fourier basis may become prohibitively slow in practice and it may be preferable to use different, more efficient systems for calibration, as we do in this paper. Another possible approach would be to first recover the image with respect to the calibration basis and then approximate its Fourier coefficients in a post-processing step. However, as a two-stage procedure, such approach is inherently less efficient and suffers from greater error than the approach proposed in this paper, which can recover Fourier coefficients directly from the raw measurements.

In the earlier work [21], a fibre endoscope was developed to produce images of phase and polarisation for early cancer detection. In this case, a set of calibration inputs was chosen so as to greatly speed up experimental characterisation measurement time by enabling parallelized calibration that exploits the localised confinement of light of the underlying fibre structure. However, this input basis is non-orthogonal and so phase-conjugation cannot be naively applied. Using the framework presented in this paper, we are able to reconstruct phase and polarisation images with respect to a diagnostically relevant representation system, while simultaneously preserving the benefits of efficient experimental calibration achieved with a system tailored to the fibre structure. Moreover, this new approach significantly decreases the reconstruction time compared to the previously implemented reconstruction technique [21], providing an important advance towards real-time image reconstruction.

Finally, we mention that changing representation systems between image recovery and sampling has previously been applied to inverse problems arising in various image and signal processing applications (see [1, 2] and references therein). Typically, it is assumed that the imaging device of a known linear transformation provides image samples with respect to a specified sampling system, while the aim is to recover a representation of the image with respect to a different system chosen so that a good approximation of the image is obtained or the number of required samples is decreased. As in this paper, the systems considered are modelled by Riesz bases or frames of infinite-dimensional function-spaces. By contrast, in this paper the imaging device produces pixel samples of a transformed image where the underlying transformation is unknown and is characterised through a calibration procedure.

2 Reconstruction framework

In this section we introduce our reconstruction framework. We start by presenting an infinite-dimensional imaging model in Section 2.1. We then consider a simplified scalar-valued setting in Section 2.2, where we derive a linear system of equations and its regularised solution while providing flexibility in choosing different systems for calibration and image representation. We then extend our framework to the vector-field case in Section 2.3.

2.1 Imaging model and reconstruction problem

In imaging through fibres or other scattering media, an input optical vector-field is related to its corresponding output optical vector-field through an integral transformation with a spatially-varying kernel , also called Green’s function or point-spread function. Specifically, such transformation can be modelled by


where is a complex-valued vector-field representing the unknown optical field on the input plane , is a complex-valued vector-field on the output plane which can be sampled, and is some unknown bounded matrix-valued function111In general, the kernel is also time-dependent as it depends on many factors such as bending of the fibre and temperature. In this paper, we account for significant measurement noise but only consider imaging at a single time point; c.f. Section 5. [31]. In particular, we consider the input field to be an object with infinite resolution, and thus, we model as an element of an infinite-dimensional function-space, such as the -space of square-integrable vector-valued functions.

In endoscopic imaging, we are interested in capturing a full optical field (i.e., amplitude, phase and polarisation) reflected from a human tissue inside the body, which is also called a wavefront, and which in this paper, we refer to as an image. In the terminology above, denotes an image, which is observed indirectly at the input imaging plane located at the end of the fibre placed inside the body, termed as the distal facet of the fibre. The fibre then transports light from the distal facet inside the body to the proximal facet outside the body where the imaging sensor directly observes at the output imaging plane. The question then becomes how to recover the unknown, infinite-dimensional optical field from the acquired samples of .

More concretely, given the pointwise measurements of the output vector-field collected at the imaging sensor


where and is the resolution of the imaging sensor, the goal is to recover the unknown function via equation (1). It is important to note that these measurements will also contain noise introduced by the measurement procedure.

This linear inverse problem is especially challenging because both the spatially-varying kernel

as well as the eigenfunctions associated with the underlying integral transform (

1) are unknown. Such eigenfunctions are termed modes of the fibre and their analytic form is available only for some limited fibres such as parabolic graded index multimode fibres [35].

To be able to recover from finitely many samples of in scenarios where neither nor the eigenfunctions are known, one strategy that may be considered is to employ a calibration procedure. Concretely, it is possible to design calibration input fields , , and to measure the corresponding output fields , which in line with the notation above are vector-valued functions related through the infinite-dimensional model given in (1). The advantage of calibration is that we now have access not only to the data given in (2) but also to the calibration data


which forms additional information with which to recover .

It is noted that while the output fields are sampled at an output imaging sensor of resolution , the calibration input fields can be evaluated on a discretised grid whose resolution does not depend on any physical limitation imposed by the fibre or by the sensor collecting the transmitted image; it only depends on the resolution of the sensors used for calibration, which may be much larger than . Therefore, as for the input , we model the inputs as elements of an infinite-dimensional function-space. Thus, the representation of as well as the device calibration can be considered with respect to a wide class of infinite-dimensional bases or over-complete systems that may not be orthogonal.

2.2 Reconstruction of scalar-fields

We begin approaching the general problem of recovering the complex vector-field , described in Section 2.1, by first solving a simpler but related problem, which once solved will provide us with the methodology and insights necessary to tackle the problem in its full generality in the next Section 2.3. Specifically, we assume in this subsection that are scalar valued functions that take values in rather than , and accordingly takes values in rather than . We highlight this difference by writing these quantities with non-bold symbols.

Our approach starts by considering all fields on the input imaging plane as elements of the same function-space , such as the -space of square-integrable scalar-valued functions supported on , with inner product defined as , for . We proceed by recovering at resolution in terms of some desired representation system in , using only the available data in (2) and (3

). Specifically, we aim to estimate the coefficients

of the -term approximation of given as


Before turning to the computation of in (4), it is insightful to work through special cases of and that are particularly useful in practice. For instance, if the aim is to recover a Fourier representation of and is the unit square, then is the -dimensional Fourier basis where , , , and (4) specialises to


More generally, may contain the first elements of a Riesz basis in , such as B-spline wavelets for example, with its corresponding biorthogonal sequence denoted by , in which case (4) specialises to , i.e., . Moreover, since for the reconstruction framework we do not require an explicit form of the coefficients , the notion of basis can be further relaxed to over-complete representation systems such as those of over-complete frames [16].

Returning to the key issue of approximating the coefficients in the general -term approximation (4) from the given measurements (2)–(3), we write each in terms of the calibration functions as


for some coefficients , whose computation we discuss below, and for some error term . Since (1) is a linear transformation of , by substituting with in (1) and writing in terms of (4) and (6), we have


By evaluating equation (7) at the measurement points , we obtain the following linear system


where is the vector having its -th entry equal to , is the matrix having its -th entry equal to , is the matrix with its -th entry equal to and is an error term containing the last two terms in the right-hand-side of (7). In addition, the error term can be seen also as encapsulating measurement error noise incurred when measuring and in (2) and (3), respectively. We then opt to define the solution of (8) as the minimisation problem


where denotes the Euclidean norm on , while the regularisation term and its parameter are described below. Once the coefficients are computed through (9), then in line with (4) we define the reconstruction of as the approximation given by


To obtain the explicit solution defined in (10), it remains to describe the procedure for computing the coefficients of matrix and to define the regularisation term .

First, observe that if the same system is used for calibration and reconstruction, then . Otherwise, we can estimate as follows. Using (6), we write , , and therefore, provided , we have

The matrix featuring this linear system is known as the Gram matrix of

, which takes the form of the identity matrix when

are orthonormal. We note that the accuracy of such estimation of matrix and its condition number depend on the gap between the function-spaces spanned by and as well as on the conditioning of the Gram matrix. In general, a good estimation requires that forms a good approximation for .

We now discuss the choice of the regularisation term in (9). If is absent from (9), i.e. if , then the solution of this minimisation problem is equivalent to the least-squares solution . If the regularisation term is given by , then (9) is known as Tikhonov regularisation and its solution is given by . However, if is badly conditioned, and, additionally, it is known a priori that only a few elements of are sufficient to represent well, then is an appropriate choice of the regularisation term. This is termed as the -regularisation, where the -norm of , i.e. , is defined as the number of non-zero coordinates in . The -regularisation bypasses the ill-conditioning by imposing sparsity in the solution with respect to . In practice, solving the minimisation problem with such a non-convex -term is computationally difficult, so typically an -relaxation is considered instead. The corresponding relaxed minimisation problem can then be solved by fast iterative algorithms [9, 8]. In addition to the choice of the form of the regularisation term , the strength of the regularisation is controlled by the parameter , which can be chosen by cross-validation techniques [18].

We conclude this subsection by a discussion on the accuracy and robustness of the solution defined in (10).

The reconstruction error can be quantified by the magnitude of , which depends on several factors. The magnitude of the first term depends of how well can be represented by its -term approximation with respect to , and thus it is expected to decrease with increasing . On the other hand, the magnitude of the second term depends on the conditioning of the matrix and the error term in (8), and thus, it is expected to increase with increasing when and are fixed. In other words, if the resolution at which we reconstruct is increased, we also need to increase and . However, it may be possible to attain higher resolutions if some form of regularisation is used when solving (8).

As previously noted, the error term contains the measurement error as well as the last two terms of the right hand side in (7), which can be disregarded provided and

are small or they lie in the span of those eigenfunctions corresponding to a small singular value. Thus, for small

it is required that and form a good approximation for or for the eigenfunctions with large singular values. However, if the singular values of the underlying integral operator accumulate at zero, the conditioning of the matrix may become worse if the span of includes too many eigenfunctions including those corresponding to a small singular value. Loosely speaking, the calibration functions should form a good representation for the span of the eigenfunctions, excluding those eigenfunctions corresponding to small singular values if they exist. However, as we do not have access to the true eigenfunctions, we do not have control over the ill-conditioning that we are introducing by using a particular choice of . Thus, the use of regularisation in solving (8) becomes crucial in order to obtain a robust solution.

2.3 Reconstruction of vector-fields

We now extend the scalar-field reconstruction framework developed in Subection 2.2 to the more general vector-field problem presented in Subection 2.1. To begin with, let

be the complex-vector-valued functions related as in equation (1), where the superscripts and correspond to the horizontal and vertical polarisations of the optical vector-field, respectively. The goal is to recover both polarisations and , which are scalar-valued functions, by using the measurements (2) and (3). Since each of and depends on both and , we cannot consider the reconstructions of and as two independent problems. However, as it will be demonstrated below, we can still use the reconstruction framework introduced earlier in Section 2.2, as long as the calibration inputs can form a representation system for vector-valued functions such as . To ensure that this is the case, rather than straightforwardly sampling the vector-valued calibration inputs from (3), i.e.

we instead sample the following two related forms


where , for a fixed and are some scalar-valued functions. To make clear the motivation to sample according to (11), observe that if and are representation systems for and respectively, then we have that and are representation systems for and respectively. In other words,


can be used to represent the complex vector-valued function .

Mimicking the reasoning of the previous subsection, we proceed by recovering approximations of and with respect to some desired representation systems and . Namely we aim to recover

where we first write these representation systems in terms of the calibration functions as and , for some coefficients and some error terms . It thus follows that



by applying (1) to (13), we obtain

provided the two last terms in (13) are small or they become small after applying (1). By using the pointwise measurements from (2) and (11), this leads to the following linear system


where , is such that its -th entry is , is such that its -th entry is , is the error term, and , , are defined as

We propose to solve the linear system (14) in a similar manner to that used in (9). Finally, once (14) is solved for the coefficients , we can define the reconstructions of and as

Observe that using -regularisation in this case imposes sparsity in the reconstructions and with respect to and , respectively. Also similarly as before, we note that matrices and are identities when the reconstruction functions and the calibration functions are the same, otherwise they can be computed approximately from the equations

Finally, we note that in order to reduce the impact of noise it may be possible to include measurements of additional phase-shifts of the calibration functions. In particular, in addition to the calibration inputs and given in (11), it may also be possible to measure

where ensures that and, as we will see shortly below, must be carefully chosen so that . With this extra information in hand, rather than using (12), the following functions are used to represent the complex vector fields


where is finite given that was chosen above in such a way that . We then proceed as above but in place of (14) obtain


where we now have the additional matrix containing the outputs . As we will see below in Section 4, augmenting the calibration data in such a way is indeed an effective manner to decrease the influence of the measurement noise on the reconstruction.

3 Fourier coefficients as informative features

Since inhomogeneities on a cellular scale caused by cancer result in increased scattering of an optical field reflected from a tumourous tissue [19], it is expected that they also result in higher spatial frequencies of the reflected optical field. Hence, we propose that by representing such optical fields in a Fourier basis and by inspecting the corresponding Fourier coefficients it is possible to detect the increased scattering of an optical field associated with the tumourous tissue and thereby gain insight into the disease status of the tissue.

In this section, we focus exclusively on the merits of the Fourier coefficients as indicative features of increased phase scattering. By using simulated data, we show how increased changes in phase result in a slower decay of the corresponding Fourier coefficients and how this effect can be quantified. In the next section, we confirm the findings of this section on real biological data, where we use the reconstruction framework developed in Section 2 to recover tissue images directly in a Fourier basis and demonstrate that the recovered Fourier coefficients are indeed useful for detecting cancer.

3.1 Fourier coefficients of a one-dimensional simulated example

We first consider a simple one-dimensional example to illustrate the effect of increased phase oscillations on the decay of the corresponding Fourier coefficients. In this example, we measure the decay of the Fourier coefficients associated with different functions , , , with the same amplitude but with different phases defined on the compact interval . In particular, for illustration purposes we take and , where , so that different phase functions exhibit different degrees of oscillations. These phase functions are shown in the first panel of Figure 1. Since the Fourier basis in the domain is given by , then for each we compute its first 20 Fourier coefficients as


, and approximate its Fourier transform by the classical Whittaker–Shannon interpolation formula

, . The amplitude of the approximated Fourier transform of each is shown in the second panel of Figure 1

. Finally, we quantify the decay of the Fourier coefficients by the standard deviation

of a Gaussian function fitted to the amplitude of the approximated Fourier transform on interval . The fitted Gaussian functions are shown in the third panel of Figure 1. From the forth panel of Figure 1, we observe that an increased magnitude of the phase oscillations results in an increased standard deviation , which can be used to quantify the decay of the Fourier coefficients. It is important to note that although in this example the zeros of the different phase functions coincide, the same effect is observed even if this is not the case. Also, if the frequency of the phase oscillation is increased while their magnitude is kept constant, then would increase as well.

[width=1trim=.10pt .000pt .090pt .00pt,clip]intro-eps-converted-to.pdf

Figure 1: Higher phase oscillations result in a slower decay of the Fourier coefficients.

The takeaway message from this simple one-dimensional example is that representing a signal with respect to a Fourier basis is especially useful to identify variations in oscillating phase, and that the decay of the corresponding Fourier coefficients is sensitive to variations in phase scattering in a manner that can be easily identified. As we will see in the remainder of the paper, these observations remain true also in higher dimensional practical examples.

3.2 Fourier coefficients of simulated tissue images

Phase                         abs(FT)                   Fitted Gaussian 1 [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]1ph-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]1ft-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]1gs-eps-converted-to.pdf 2 [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]2ph-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]2ft-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]2gs-eps-converted-to.pdf 3 [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]3ph-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]3ft-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]3gs-eps-converted-to.pdf 4 [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]4ph-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]4ft-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]4gs-eps-converted-to.pdf 5 [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]5ph-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]5ft-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]5gs-eps-converted-to.pdf 6 [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]6ph-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]6ft-eps-converted-to.pdf [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]6gs-eps-converted-to.pdf Amplitude [scale=0.5,trim=.050pt .050pt .050pt .000pt,clip]amp-eps-converted-to.pdf [scale=0.49,trim=.000pt .020pt .020pt .020pt,clip]plot-eps-converted-to.pdf
Figure 2: Six simulated images with the same amplitude but different phase, which are generated from our model with increasing , , so that larger characterises larger phase oscillations. In the scatter plot on the right, we report the sum of parameters and of the Gaussian fitted to the amplitude of the Fourier transform abs(FT), revealing that increased correlates with larger .

We now generalise our observations to two-dimensional complex-valued functions. In particular, we perform a simulation study to demonstrate that higher and more frequent changes in phase result in a slower decay of the Fourier coefficients of the corresponding function. For this purpose, we create a model mimicking tissue samples with a different level of phase oscillations, which we then use to generate images and compute their Fourier coefficients.

In our model, we use randomness to achieve certain variability across different samples and two different parameters to control the degree of phase oscillations. In detail, the model that we use in our simulation study corresponds to a complex function , , where the original space-domain is discretized into a grid, while and are chosen randomly as we now describe. The phase function depends on two given parameters and , controlling the amplitude and the frequency of phase oscillations, respectively. Specifically, pixel-values are chosen uniformly at random from , which are then filtered by using MATLAB’s function ‘imgaussfilt’ with the smoothing parameter . Following this step, only pixels are kept by removing pixels from each boundary and such image is then rescaled so that all phase pixel-values are between , . The amplitude function is selected as the sum of and five additional Gaussian functions with randomly chosen parameters and .

In Figure 2 we demonstrate how changing phase parameters and while keeping amplitude fixed changes the decay of Fourier coefficients. Specifically, we use six different values , to create six different functions , where and are increasing logarithmically. Similarly to the one-dimensional example of Figure 1, the decay of corresponding Fourier coefficients is measured by standard deviation of a Gaussian function fitted to the absolute value of the Fourier transform approximated from the first Fourier coefficients, where the Fourier coefficients of function are computed using the formula in (5). In particular, in Figure 2, for each we report the sum of the standard deviations of the fitted Gaussian, thereby observing that increased phase oscillations, i.e. increased , results in slower decay of the corresponding Fourier coefficients, i.e. larger .

Next, in Figure 3, for each of the six different values , , chosen as in Figure 2, we generate a hundred different images using our model (with each image having a different phase and a different amplitude) and we report the value of the fitted Gaussian. We observe the same trend in the decay of the Fourier coefficients in Figure 3 as in Figure 2, but now across 600 different images.

[width=0.8trim=.000pt .020pt .020pt .020pt,clip]boxplot-eps-converted-to.pdf
Figure 3: For each of the six categories, we generated a hundred images with different phase and amplitude from the tissue model with fixed parameters and , . For each image, we then approximated its Fourier transform and fitted a Gaussian function whose parameter is reported.

We conclude this section by noting that the features extracted from Fourier coefficients as we described in this section have three additional useful properties. First, since the amplitude of the Fourier transform is invariant to the shifts of the corresponding complex function in its space-domain, the features that we extract are invariant to the shifts of the tissue images in their space-domain. Second, the quality of the recovered phase in the space-domain is dependent on a phase unwrapping procedure and is thus highly sensitive to noise, which means that phase may bear more information in the Fourier-domain than in the space-domain. Third, once the Fourier coefficients are computed from the available measurements, each image can easily be represented in both the Fourier and the original space-domain, allowing for additional flexibility.

4 Experimental results

Having established the utility of Fourier coefficients in quantifying phase scattering using simulated data in Section 3, we now apply the reconstruction framework developed in Section 2 to measurements obtained experimentally using the prototype fibre endoscope developed in [21], which can measure optical phase, polarisation and amplitude. In Section 4.1, we first demonstrate the recovery of a synthetic holographic image with a known ground-truth that can be used for validation. Next, in Section 4.2, we apply our reconstruction framework to biological images of tissue samples taken from mice and demonstrate that reconstruction with respect to a Fourier basis can be used as a diagnostic indicator of early tumorigenesis.

4.1 Reconstruction of a synthetic holographic image

To demonstrate our reconstruction algorithm on experimental data we first reconstruct a synthetic holographic image. Since in this case we have access to the ground-truth image at the distal end of the fibre, we can visually assess the quality of our proposed imaging methodology. Specifically, in this subsection, using the raw output of the synthetic holographic image shown in Figure 4, we test our general reconstruction framework in combination with different representation systems as well as different regularisation terms.

Holographic image at the distal end Raw output at the proximal end
Figure 4: Amplitude and phase of the horizontal and vertical polarisation of the ground-truth synthetic holographic image at the distal end and the corresponding output at the proximal end.

The fibre is calibrated using input and output pairs such as those shown in Figure 5. For the purposes of reconstruction, each calibration input and output is separated from the others by evaluating each of them only over a circular region around the centre of the corresponding Gaussian-like spot. In particular, the calibration inputs in Figure 5 are evaluated on a grid with resolution and they are translated to different locations across the input imaging plane. Each output is evaluated at different pixels at the output imaging plane. Considering the two polarisation states, the dimension of the system matrix in (16) is therefore .

Distal end Proximal end
Figure 5: Amplitude and phase of the horizontal polarisation of one calibration input ( of Eq. (15)) at the distal end and at the proximal end.

In Figure 6, we recover the amplitude and the phase of the horizontal and vertical polarisations of the holographic image from raw endoscopic measurements using different inversion techniques while reconstructing with respect to the calibration coefficients. In particular, we solve (16) where , by inverting the linear system in four different ways:

  1. [label=0.]

  2. the naive inversion , which corresponds to the principle of phase conjugation in that it assumes ,

  3. the least-squares approach ,

  4. the -regularisation , and,

  5. the -regularisation using the iterative solver [10].

In particular, we see in Figure 6 that -regularisation performs well when compared with the other approaches. In fact, since our holographic image is sparse with respect to the calibration inputs – namely, since and are sparse with respect to and -regularisation successfully removes significant noise while preserving the image details. The amount of noise that is removed by -regularisation depends on the strength of the regularisation parameter .

Horizontal polarisation Vertical polarisation
Amplitude Phase Amplitude Phase

naive approach




Figure 6: Reconstructed amplitude and phase of the horizontal and vertical polarisations of the holographic image from Figure 4 with respect to the calibration functions such as those in Figure 5, using naive, least-squares, and approaches. The regularisation parameter in the and -regularisation is and , respectively.

Finally, in Figure 7, we reconstruct the holographic image with respect to different representation systems, namely we solve (16) where both and correspond to the representation system of a Fourier or a wavelet basis with cardinality . Specifically,

  1. [label=()]

  2. in the Fourier case, we choose both and to be , where , whereas

  3. in the wavelet case, we choose tensor-products of

    one-dimensional boundary-corrected Daubechies wavelets with four vanishing moments (DB4) from


We can observe in Figure 7 that least-squares fails to give a useful estimate due to the ill-conditioning of the system matrix, conveying that it is crucial to use the regularisation term. Although the least-squares could still be used in the case where , small does not necessarily lead to a good approximation of the image, and so to achieve the desired resolution one would need to increase the number of calibration measurements . This is undesirable as it would incur additional experimental time. On the other hand, given that our holographic image is sparse with respect to compactly-supported wavelets, we see from Figure 7 that -regularisation performs quite well in combination with DB4 wavelets even though .

Fourier exponentials DB4 wavelets
Amplitude Phase Amplitude Phase

naive approach




Figure 7: Reconstructed images (horizontal polarisation) with respect to the different bases using different inversion approaches. We used Fourier exponentials and DB4 wavelets. Regularisation parameter in the and -regularisation is and , respectively.

4.2 Reconstruction and analysis of biological images

We now apply the framework developed in Sections 23 to reconstruct and analyse images of real-world biological samples. We imaged ex vivo samples of mouse oesophagus from healthy controls and from carcinogen treated animals with induced oesophageal tumours (lesions) using the model presented in [5]. More concretely, 3 control mice (6 healthy areas analysed) and 6 mice with induced tumours (6 distinct lesions analysed) were used. Each sample was segmented into areas of healthy and lesion tissue by an expert using DAPI fluorescence images.

For clarity in the examples below, we index different areas by , where the first six are healthy and the rest are lesions. Due to the limited field of view of the endoscope () relative to the sample size (mm), each of the 9 mice produce 6–20 individual images corresponding to different areas of the same sample that may overlay by up to 15%. We therefore also introduce index to denote individual sub-images within a larger area on a given sample. Each individual sample in the data set thus has associated index , where and , for in the range 6–20.

Space domain                                         Fourier domain
       Amplitude          Unwrapped phase            Amplitude             Fitted Gaussian
LESION [scale=0.6]L1_ABS-eps-converted-to.pdf[scale=0.6]L1_PH-eps-converted-to.pdf[scale=0.6]L1_FT-eps-converted-to.pdf[scale=0.6]L1_GS-eps-converted-to.pdf

Figure 8: Healthy and lesion tissue reconstructed with respect to Fourier coefficients by solving the linear system (14) with -regularisation with . Only the horizontal polarisation in shown due to space limitations.

In Figure 8, we show the reconstruction of the horizontal polarisation of one healthy image indexed as (first row) and one lesion image indexed as (second row), in both the space-domain (first two columns) and the Fourier-domain (last two columns). Specifically, we reconstruct Fourier coefficients per polarisation by solving the linear system (14) with -regularisation. We then expand these coefficients with respect to Fourier-exponentials to obtain images in the space-domain, as well as, with respect to sinc-functions to obtain images in the Fourier-domain. In the space-domain, we show the amplitude (first column) and unwrapped phase (second column) of the reconstructed image, where for the unwrapping we used an efficient 2D phase unwrapper222Available on-line at In the Fourier domain, we show the amplitude (third column) of the reconstructed Fourier transform as well as a Gaussian fit to the amplitude of the Fourier transform (fourth column), where we used the procedure explained in Section 3. While the difference between the healthy and the lesion sample is not so apparent from the amplitude and phase in the space-domain, the difference becomes more pronounced in the Fourier domain; specifically, it can be observed that the Fourier coefficients decay slower in the lesion than in the healthy tissue, where the decay is quantified by the standard deviation of the fitted Gaussian.

To see if the standard deviation of the fitted Gaussian can be used to discriminate between healthy and lesion samples in general, we present in Figure 9 the variation of the Fourier-domain information across different tissue samples. In particular, for each reconstructed image , , we show the amplitude of the Fourier transform and the associated Gaussian function. We can see in Figure 9 that the discriminative behaviour of the standard deviation of the fitted Gaussian between one particular pair of healthy and lesion samples seen in Figure 8, holds more generally throughout the dataset.

HEALTHY [scale=0.61]FT1-eps-converted-to.pdf[scale=0.61]FT2-eps-converted-to.pdf[scale=0.61]FT3-eps-converted-to.pdf[scale=0.61]FT4-eps-converted-to.pdf[scale=0.61]FT5-eps-converted-to.pdf[scale=0.61]FT6-eps-converted-to.pdf
LESION [scale=0.61]FT7-eps-converted-to.pdf[scale=0.61]FT8-eps-converted-to.pdf[scale=0.61]FT9-eps-converted-to.pdf[scale=0.61]FT10-eps-converted-to.pdf[scale=0.61]FT11-eps-converted-to.pdf[scale=0.61]FT12-eps-converted-to.pdf
HEALTHY [scale=0.61]GS1-eps-converted-to.pdf[scale=0.61]GS2-eps-converted-to.pdf[scale=0.61]GS3-eps-converted-to.pdf[scale=0.61]GS4-eps-converted-to.pdf[scale=0.61]GS5-eps-converted-to.pdf[scale=0.61]GS6-eps-converted-to.pdf
LESION [scale=0.61]GS7-eps-converted-to.pdf[scale=0.61]GS8-eps-converted-to.pdf[scale=0.61]GS9-eps-converted-to.pdf[scale=0.61]GS10-eps-converted-to.pdf[scale=0.61]GS11-eps-converted-to.pdf[scale=0.61]GS12-eps-converted-to.pdf

Figure 9: (a): Amplitude of the Fourier transform of different healthy and lesion tissues, computed from the Fourier coefficients of the horizontal polarisation. (b): Gaussian function fitted to the amplitude of the Fourier transforms shown in (a).

Finally, in Figure 10 we perform a statistical test which confirms that the standard deviation of the fitted Gaussians is an informative feature to distinguish between healthy tissues and lesion tissues. In particular, for each individual sample in the data set, we compute of the fitted Gaussian with parameters and . For each tissue sample , we then compute the average value and show a box-plot of these twelve values grouped according to their class label ‘healthy’ or ‘lesion’, for each polarisation as well as for both polarisations combined. We also compute the

-value of Welch’s t-test

[42], showing the significant difference in the decay of Fourier coefficients between healthy and lesion samples.

[scale=0.75,trim=.00pt .00pt .050pt .000pt,clip]box_hor-eps-converted-to.pdf [scale=0.75,trim=.00pt .00pt .050pt .000pt,clip]box_ver-eps-converted-to.pdf [scale=0.75,trim=.00pt .00pt .050pt .000pt,clip]box_horver-eps-converted-to.pdf

Figure 10: Standard deviation , , of a Gaussian function fitted to the amplitude of the Fourier transform of 6 healthy and 6 lesion areas across our sample set of mouse oesophagus tissue. In the box-plot on the right, a two-dimensional feature corresponding to two different polarisations is used to characterise each sample, and thus, prior to the t-test, each sample is rescaled by the mean and standard deviation of the total of 12 samples.

We may conclude that the flexibility of the approach developed in Section 2

, which permits the use of different systems for calibration and image representation, is particularly useful when reconstructing images of biological tissues. The possibility to use a Fourier basis for reconstruction provides the opportunity to investigate images in both the space-domain and the Fourier domain, and thus investigate the degree to which the reconstructed Fourier coefficients decay. We showed that this decay quantified by the standard deviation of the fitted Gaussian is a feature with a discriminative power between healthy tissues and lesions, which in the future, in conjunction with a larger data set, could be used to build an automated classifier for distinguishing healthy and lesion samples.

5 Discussion and future research

Towards the development and practical use of fibre endoscopes, the main contributions in this paper are two-fold. Firstly, we demonstrated that a Fourier basis yields a diagnostically relevant representation of the optical field reflected from a tissue, using both simulated and experimental real-world data. Secondly, we provided a general reconstruction algorithm that through regularisation can stably recover such representation directly from the calibration measurements and the measurements of the output optical field transmitted through a fibre, where the system used for calibration is allowed to be different and thus more efficient than a Fourier basis.

Nevertheless several open problems remain. One possible direction for future research, which would require a significantly larger number of biological samples to be tested, relates to learning an ‘optimal’ dictionary (alternative to Fourier) as a means to minimise the classification error between healthy and lesion tissues. More importantly, further work is required to achieve real-time imaging through fibre endoscopes in realistic clinical settings. Specifically, future research is needed to lift the time-independence assumption present in the kernel of the linear model which in everyday clinical use varies across time with bending and temperature. In practice, the time-independence assumption means that the calibration measurements need to be taken often and under similar bending and temperature conditions as when sampling the output optical field, which is difficult to achieve in realistic clinical deployments. Although recently there have been some initial steps in addressing this issue, as for example [22], the development of a clinically-feasible recovery procedure that accounts for significant fibre changes remains an important open problem.


MG and SEB were supported by an EPSRC grant EP/N014588/1 for the centre for Mathematical and Statistical Analysis of Multimodal Clinical Imaging. GSDG and SEB received funding from CRUK (C47594/A16267, C14303/A17197, C47594/A21102) and a pump-priming award from the Cancer Research UK Cambridge Centre Early Detection Programme (A20976). The research of FR was performed in part while he was at University of Cambridge, and it was funded in part by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 655282 and in part by the FCT grant SFRH/BPD/118714/2016. The research of AGCPR was conducted during January–June 2017, while visiting University of Cambridge.


  • [1] B. Adcock, A. C. Hansen, B. Roman, and G. Teschke, “Chapter Four - Generalized Sampling: Stable Reconstructions, Inverse Problems and Compressed Sensing over the Continuum,” Advances in Imaging and Electron Physics, vol. 182, pp. 187–279, 2014.
  • [2] B. Adcock, and A. C. Hansen, “Generalized sampling and infinite-dimensional compressed sensing,” Foundations of Computational Mathematics, vol. 16, no. 5, 2016.
  • [3] I. Ahmad, M. Ahmad, K. Khan, S. Ashraf, S. Ahmad, and M. Ikram, “Ex vivo characterization of normal and adenocarcinoma colon samples by Mueller matrix polarimetry,” Journal of Biomedical Optics, vol. 20, no. 5, p. 056012, 2015.
  • [4] S. Alali, and I. A. Vitkin, “Polarized light imaging in biomedicine: emerging Mueller matrix methodologies for bulk tissue assessment,” Journal of Biomedical Optics, vol. 20, no. 6, p. 20209, 2015.
  • [5] M. P. Alcolea, P. Greulich, A. Wabik, J. Frede, B. D. Simons, and P. H. Jones, “Differentiation imbalance in single oesophageal progenitor cells causes clonal immortalization and field change,” Nature Cell Biology, 16(6), 615–22, 2014.
  • [6] D. Arifler, I. Pavlova, A. Gillenwater, and R. Richards-Kortum, “Light scattering from collagen fibre networks: micro-optical properties of normal and neoplastic stroma,” Biophysical journal, vol. 92, pp. 3260–74, 2007.
  • [7] C. Ba, M. Palmiere, J. Ritt, and J. Mertz, “Dual-modality endomicroscopy with co-registered fluorescence and phase contrast,” Biomedical optics express, vol. 7, no. 9, pp. 3403–3411, 2016.
  • [8] S. Becker, J. Bobin, and E. J. Candès, “NESTA: a fast and accurate first-order method for sparse recovery,” SIAM Journal on Imaging Sciences, vol. 4, no. 1, pp. 1–39, 2011.
  • [9] E. van den Berg and M. P. Friedlander, “Probing the Pareto frontier for basis pursuit solutions,” SIAM Journal on Scientific Computing, vol. 31, no. 2, pp. 890–912, 2008.
  • [10] E. van den Berg and M. P. Friedlander, “SPGL1: A solver for large-scale sparse reconstruction,”, 2007.
  • [11] J. A. Carpenter, B. J. Eggleton, and J. Schroeder, “Maximally efficient imaging through multimode fibre,” in CLEO: 2014, (Washington, D.C.), p. STh1H.3, OSA, 2014.
  • [12] J. Carpenter, B. Eggleton, and J. Schröder, “110x110 optical mode transfer matrix inversion,” Optics Express, vol. 22, no. 1, pp. 96–101, 2014.
  • [13] T. Čižmár and K. Dholakia, “Tunable Bessel light modes: engineering the axial propagation,” Optics Express, vol. 17, no. 18, p. 15558, 2009.
  • [14] T. Čižmár and K. Dholakia, “Shaping the light transmission through a multimode optical fibre: complex transformation analysis and applications in biophotonics,” Optics express, vol. 19, pp. 18871–84, 2011.
  • [15] T. Cižmár and K. Dholakia, “Exploiting multimode waveguides for pure fibre-based imaging,” Nature communications, vol. 3, p. 1027, 2012.
  • [16] O. Christensen, “An Introduction to Frames and Riesz Bases,” Birkhäuser, 2002.
  • [17] A. Cohen, I. Daubechies, and P. Vial, “Wavelets on the interval and fast wavelet transforms,” Applied and Computational Harmonic Analysis, vol. 1, no. 1, pp. 54–81, 1993.
  • [18] A. Doostan and H. Owhadi, “A non-adapted sparse approximation of PDEs with stochastic inputs,” Journal of Computational Physics, vol. 230, no. 8, pp. 3015–3034, 2011.
  • [19] R. Drezek, A. Dunn, and R. Richards-Kortum “Light scattering from cells: finite-difference time-domain simulations and goniometric measurements,” Applied Optics, vol. 38, no. 16, 1999.
  • [20] R. Drezek, M. Guillaud, T. Collier, I. Boiko, A. Malpica, C. Macaulay, M. Follen, and R. Richards-Kortum, “Light scattering from cervical cells throughout neoplastic progression: influence of nuclear morphology, DNA content, and chromatin texture,” Journal of biomedical optics, vol. 8, pp. 7–16, 2003.
  • [21] G. S. D. Gordon, J. Joseph, M. P. Alcolea, T. Sawyer, A. J. Macfaden, C. Williams, C. R. M. Fitzpatrick, P. H. Jones, M. di Pietro, R. C. Fitzgerald, T. D. Wilkinson, and S. E. Bohndiek, “Quantitative phase and polarisation endoscopy for detection of early oesophageal tumourigenesis,” Under review
  • [22] R. Y. Gu, R. N. Mahalati, and J. M. Kahn, “Design of flexible multi-mode fibre endoscope,” Optics Express, vol. 23, no. 21, pp. 26905–26918, 2015.
  • [23] M. Kim, W. W. Choi, Y. Choi, C. Yoon, and W. W. Choi, “Transmission matrix of a scattering medium and its applications in biophotonics,” Optics Express, vol. 23, p. 12648, 2015.
  • [24] B. Kunnen, C. Macdonald, A. Doronin, S. Jacques, M. Eccles, and I. Meglinski, “Application of circularly polarized light for non-invasive diagnosis of cancerous tissues and turbid tissue-like scattering media,” Journal of Biophotonics, vol. 8, pp. 317–323, 2015.
  • [25] D. Loterie, S. Farahi, I. Papadopoulos, A. Goy, D. Psaltis, and C. Moser, “Digital confocal microscopy through a multimode fibre,” Optics Express, vol. 23, p. 23845, 2015.
  • [26] A. Pierangelo, A. Benali, M.-R. Antonelli, T. Novikova, P. Validire, B. Gayet, and A. De Martino, “Ex-vivo characterization of human colon cancer by Mueller polarimetric imaging,” Optics Express, vol. 19, p. 1582, 2011.
  • [27] M. Plöschner, T. Tyc, and T. Čižmár, “Seeing through chaos in multimode fibres,” Nature Photonics, vol. 9, pp. 529–535, 2015.
  • [28] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, a. C. Boccara, and S. Gigan, “Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media,” Physical Review Letters, vol. 104, p. 100601, mar 2010.
  • [29] S. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Image transmission through an opaque material,” Nature Communications, vol. 1, pp. 1–5, 2010.
  • [30] S. M. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Controlling light through optical disordered media: Transmission matrix approach,” New Journal of Physics, vol. 13, p. 123021, 2011.
  • [31] S. Rotter and S. Gigan, “Light fields in complex media: Mesoscopic scattering meets wave control,” Reviews of Modern Physics, vol. 89, no. 1, 2017.
  • [32] J. Qi and D. S. Elson, “A high definition Mueller polarimetric endoscope for tissue characterisation,” Scientific Reports, vol. 6, no. April, p. 25953, 2016.
  • [33] S. Sridharan, V. Macias, K. Tangella, A. Kajdacsy-Balla, and G. Popescu, “Prediction of Prostate Cancer Recurrence Using Quantitative Phase Imaging,” Scientific Reports, vol. 5, p. 9976, 2015.
  • [34] J. W. Su, Y. H. Lin, C. P. Chiang, J. M. Lee, C. M. Hsieh, M. S. Hsieh, … K. B. Sung, “Precancerous esophageal epithelia are associated with significantly increased scattering coefficients,” Biomedical Optics Express, 6(10), 3795, 2015.
  • [35] A. W. Snyder and J. D. Love, Optical Waveguide Theory. New York: Chapman and Hall, 1983.
  • [36] A. J. Thompson, C. Paterson, M. a. a. Neil, C. Dunsby, and P. M. W. French, “Adaptive phase compensation for ultracompact laser scanning endomicroscopy,” Optics letters, vol. 36, pp. 1707–9, 2011.
  • [37] H.-D. Wang, C. H. Niu, Q. Yang, and I. Badea, “Study on protein conformation and adsorption behaviors in nanodiamond particle-protein complexes,” Nanotechnology, vol. 22, p. 145703, 2011.
  • [38] Z. Wang, K. Tangella, A. Balla, and G. Popescu, “Tissue refractive index as marker of disease,” Journal of biomedical optics, vol. 16, no. 11, p. 116017, 2011.
  • [39] J. Wang, J. Zhou, H. Long, Y. Xie, X. Zhang, H. Luo, Z. Deng, Q. Wei, Z. Yu, J. Zhang, and Z. Tang, “Tribological, anti-corrosive properties and biocompatibility of the micro- and nano-crystalline diamond coated Ti6Al4V,” Surface and Coatings Technology, vol. 258, pp. 1032–1038, 2014.
  • [40] S. C. Warren, Y. Kim, J. M. Stone, C. Mitchell, J. C. Knight, M. A. A. Neil, C. Paterson, P. M. W. French, and C. Dunsby, “Adaptive multiphoton endomicroscopy through a dynamically deformed multicore optical fibre using proximal detection,” Opt. Express, vol. 24, no. 19, pp. 21474–21484, 2016.
  • [41] D. J. Waterhouse, C. R. M. Fitzpatrick, M. di Pietro, and S. E. Bohndiek, “Emerging optical methods for endoscopic surveillance of Barrett’s oesophagus,” The Lancet Gastroenterology & Hepatology, vol. 3, no. 5, pp. 249–362, 2018.
  • [42]

    B. L. Welch, “The generalization of “Student’s problem when several different population variances are involved,”

    Biometrika, vol. 34, pp. 28–35, 1947.