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 vectorfield 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 polarisationsensitive 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 twodimensional complex vectorfield, where different regularisation terms are permitted and the bases used for image representation and device calibration are allowed to be different and/or nonorthogonal. 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 custombuilt 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, finitedimensional 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 phaseonly 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 phaseconjugation such as leastsquares (MoorePenrose 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 infinitedimensional 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 postprocessing step. However, as a twostage 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 nonorthogonal and so phaseconjugation 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 realtime 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 infinitedimensional functionspaces. 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 infinitedimensional imaging model in Section 2.1. We then consider a simplified scalarvalued 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 vectorfield case in Section 2.3.
2.1 Imaging model and reconstruction problem
In imaging through fibres or other scattering media, an input optical vectorfield is related to its corresponding output optical vectorfield through an integral transformation with a spatiallyvarying kernel , also called Green’s function or pointspread function. Specifically, such transformation can be modelled by
(1) 
where is a complexvalued vectorfield representing the unknown optical field on the input plane , is a complexvalued vectorfield on the output plane which can be sampled, and is some unknown bounded matrixvalued function^{1}^{1}1In general, the kernel is also timedependent 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 infinitedimensional functionspace, such as the space of squareintegrable vectorvalued 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, infinitedimensional optical field from the acquired samples of .
More concretely, given the pointwise measurements of the output vectorfield collected at the imaging sensor
(2) 
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 spatiallyvarying 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 vectorvalued functions related through the infinitedimensional 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
(3) 
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 infinitedimensional functionspace. Thus, the representation of as well as the device calibration can be considered with respect to a wide class of infinitedimensional bases or overcomplete systems that may not be orthogonal.
2.2 Reconstruction of scalarfields
We begin approaching the general problem of recovering the complex vectorfield , 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 nonbold symbols.
Our approach starts by considering all fields on the input imaging plane as elements of the same functionspace , such as the space of squareintegrable scalarvalued 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(4) 
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
(5) 
More generally, may contain the first elements of a Riesz basis in , such as Bspline 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 overcomplete representation systems such as those of overcomplete 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
(6) 
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
(7) 
By evaluating equation (7) at the measurement points , we obtain the following linear system
(8) 
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 righthandside 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
(9) 
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
(10) 
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 functionspaces 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 leastsquares 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 nonzero coordinates in . The regularisation bypasses the illconditioning by imposing sparsity in the solution with respect to . In practice, solving the minimisation problem with such a nonconvex 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 crossvalidation 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 illconditioning 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 vectorfields
We now extend the scalarfield reconstruction framework developed in Subection 2.2 to the more general vectorfield problem presented in Subection 2.1. To begin with, let
be the complexvectorvalued functions related as in equation (1), where the superscripts and correspond to the horizontal and vertical polarisations of the optical vectorfield, respectively. The goal is to recover both polarisations and , which are scalarvalued 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 vectorvalued functions such as . To ensure that this is the case, rather than straightforwardly sampling the vectorvalued calibration inputs from (3), i.e.
we instead sample the following two related forms
(11) 
where , for a fixed and are some scalarvalued 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,
(12) 
can be used to represent the complex vectorvalued 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
(13) 
Since
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
(14) 
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 phaseshifts 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
(15) 
where is finite given that was chosen above in such a way that . We then proceed as above but in place of (14) obtain
(16) 
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 onedimensional simulated example
We first consider a simple onedimensional 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
where
, 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.The takeaway message from this simple onedimensional 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
We now generalise our observations to twodimensional complexvalued 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 spacedomain 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, pixelvalues 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 pixelvalues 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 onedimensional 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.
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 spacedomain, the features that we extract are invariant to the shifts of the tissue images in their spacedomain. Second, the quality of the recovered phase in the spacedomain is dependent on a phase unwrapping procedure and is thus highly sensitive to noise, which means that phase may bear more information in the Fourierdomain than in the spacedomain. Third, once the Fourier coefficients are computed from the available measurements, each image can easily be represented in both the Fourier and the original spacedomain, 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 groundtruth 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 groundtruth 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  
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 Gaussianlike 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  
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:

[label=0.]

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

the leastsquares approach ,

the regularisation , and,

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 

leastsquares 

regularisation 

regularisation 
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,

[label=()]

in the Fourier case, we choose both and to be , where , whereas
We can observe in Figure 7 that leastsquares fails to give a useful estimate due to the illconditioning of the system matrix, conveying that it is crucial to use the regularisation term. Although the leastsquares 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 compactlysupported 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 

leastsquares 

regularisation 

regularisation 
4.2 Reconstruction and analysis of biological images
We now apply the framework developed in Sections 2–3 to reconstruct and analyse images of realworld 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 subimages 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.
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 spacedomain (first two columns) and the Fourierdomain (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 Fourierexponentials to obtain images in the spacedomain, as well as, with respect to sincfunctions to obtain images in the Fourierdomain. In the spacedomain, 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 unwrapper^{2}^{2}2Available online at www.ljmu.ac.uk/research/centresandinstitutes/facultyofengineeringandtechnologyresearchinstitute/geri/phaseunwrapping. 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 spacedomain, 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 Fourierdomain 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.
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 boxplot 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 ttest
[42], showing the significant difference in the decay of Fourier coefficients between healthy and lesion 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 spacedomain 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 twofold. 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 realworld 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 realtime imaging through fibre endoscopes in realistic clinical settings. Specifically, future research is needed to lift the timeindependence assumption present in the kernel of the linear model which in everyday clinical use varies across time with bending and temperature. In practice, the timeindependence 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 clinicallyfeasible recovery procedure that accounts for significant fibre changes remains an important open problem.
Acknowledgements:
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 pumppriming 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łodowskaCurie 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.
References
 [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 infinitedimensional 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. RichardsKortum, “Light scattering from collagen fibre networks: microoptical properties of normal and neoplastic stroma,” Biophysical journal, vol. 92, pp. 3260–74, 2007.
 [7] C. Ba, M. Palmiere, J. Ritt, and J. Mertz, “Dualmodality endomicroscopy with coregistered 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 firstorder 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 largescale sparse reconstruction,” http://www.cs.ubc.ca/labs/scl/spgl1, 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 fibrebased 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 nonadapted 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. RichardsKortum “Light scattering from cells: finitedifference timedomain 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. RichardsKortum, “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 multimode 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 noninvasive diagnosis of cancerous tissues and turbid tissuelike 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, “Exvivo 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. KajdacsyBalla, 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 particleprotein 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, anticorrosive properties and biocompatibility of the micro and nanocrystalline 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.
Comments
There are no comments yet.