Reconstruction of Optical Vector-Fields With Applications in Endoscopic Imaging

We introduce a framework for the reconstruction of the amplitude, phase, and polarization of an optical vector-field using measurements acquired by an imaging device characterized by an integral transform with an unknown spatially variant kernel. By incorporating effective regularization 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 for device 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 and biological tissue samples in an experimental setting, where the measurements of an optical vector-field are acquired by a multicore fiber endoscope, and observe that indeed the recovered Fourier coefficients are useful in distinguishing healthy tissues from tumors in early stages of oesophageal cancer.


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 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 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 1 -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.

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 A : x ∈ C p → y ∈ C n , leading to a transmission matrix A ∈ C n×p which is then characterised through calibration.The individual calibration inputs are arranged into the columns of matrix X cal ∈ C p×m and the corresponding calibration outputs into the columns of matrix Y cal ∈ C n×m .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 X cal is unitary.Then, by assuming that A is also unitary, images can be recovered using phase conjugation.In this approach a (generalised) inverse of the transmission matrix is calculated as X cal Y * cal , where • * denotes the conjugate transpose, and a representation of x with respect to the calibration inputs is recovered as X cal Y * cal y [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, 1 -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.

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 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 vector-field case in Section 2.3.

Imaging model and reconstruction problem
In imaging through fibres or other scattering media, an input optical vector-field F is related to its corresponding output optical vector-field F through an integral transformation with a spatially-varying kernel G, also called Green's function or point-spread function.Specifically, such transformation can be modelled by where F : S → C 2 is a complex-valued vector-field representing the unknown optical field on the input plane S ⊆ R 2 , F : R 2 → C 2 is a complex-valued vector-field on the output plane which can be sampled, and G : R 2 × R 2 → C 2×2 is some unknown bounded matrix-valued function1 [31].In particular, we consider the input field F to be an object with infinite resolution, and thus, we model F as an element of an infinite-dimensional function-space, such as the L 2 -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, F denotes an image, which is observed indirectly at the input imaging plane S 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 S inside the body to the proximal facet outside the body where the imaging sensor directly observes F at the output imaging plane.The question then becomes how to recover the unknown, infinite-dimensional optical field F from the acquired samples of F.
More concretely, given the pointwise measurements of the output vector-field F collected at the imaging sensor where y n ∈ R 2 and N ∈ N is the resolution of the imaging sensor, the goal is to recover the unknown function F 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 G 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 F from finitely many samples of F in scenarios where neither G 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 E m , m = 1, . . ., M , and to measure the corresponding output fields Ẽm , 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 F. It is noted that while the output fields Ẽm are sampled at an output imaging sensor of resolution N , the calibration input fields E m 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 M .Therefore, as for the input F, we model the inputs E m as elements of an infinite-dimensional function-space.Thus, the representation of F 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.

Reconstruction of scalar-fields
We begin approaching the general problem of recovering the complex vector-field F, 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 F, F , E m , Ẽm are scalar valued functions that take values in C rather than C 2 , and accordingly G takes values in C rather than C 2×2 .We highlight this difference by writing these quantities with non-bold symbols.
Our approach starts by considering all fields on the input imaging plane S as elements of the same function-space F, such as the L 2 -space of square-integrable scalar-valued functions supported on S, with inner product defined as E, H := S E(x)H * (x) dx, for E, H ∈ F. We proceed by recovering F ∈ F at resolution K ∈ N in terms of some desired representation system {H k } K k=1 in F, using only the available data in ( 2) and ( 3).Specifically, we aim to estimate the coefficients f = f 1 , . . ., f K ∈ C K of the K-term approximation of F given as Before turning to the computation of f k in (4), it is insightful to work through special cases of S and {H k } K k=1 that are particularly useful in practice.For instance, if the aim is to recover a Fourier representation of F and x 2 ) ∈ S, and (4) specialises to More generally, {H k } K k=1 may contain the first K elements of a Riesz basis in F, such as B-spline wavelets for example, with its corresponding biorthogonal sequence denoted by { Hk } K k=1 , in which case (4) specialises to F K (x) = K k=1 F, Hk H k , i.e., f k = F, Hk .Moreover, since for the reconstruction framework we do not require an explicit form of the coefficients f k , 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 f k in the general K-term approximation F K (4) from the given measurements ( 2)-( 3), we write each H k in terms of the calibration functions {E m } M m=1 as for some coefficients h m,k ∈ C, whose computation we discuss below, and for some error term δ k .Since ( 1) is a linear transformation of F , by substituting in ( 1) and writing F K in terms of ( 4) and ( 6), we have (7) By evaluating equation ( 7) at the measurement points {y n } N n=1 , we obtain the following linear system where g ∈ C N is the vector having its n-th entry equal to F (y n ), E ∈ C N ×M is the matrix having its (n, m)-th entry equal to Ẽm (y n ), H ∈ C M ×K is the matrix with its (m, k)-th entry equal to h m,k and ε ∈ C N is an error term containing the last two terms in the right-hand-side of (7).In addition, the error term ε ∈ C N can be seen also as encapsulating measurement error noise incurred when measuring F (y n ) and Ẽm (y n ) in ( 2) and (3), respectively.We then opt to define the solution of ( 8) as the minimisation problem where • 2 denotes the Euclidean norm on C N , while the regularisation term R and its parameter λ ≥ 0 are described below.Once the coefficients f = f1 , . . ., fK ∈ C K are computed through (9), then in line with (4) we define the reconstruction of F as the approximation given by FK (x) : To obtain the explicit solution defined in (10), it remains to describe the procedure for computing the coefficients of matrix H and to define the regularisation term R.
First, observe that if the same system is used for calibration and reconstruction, then H = I.Otherwise, we can estimate H as follows.Using (6), we write The matrix featuring this linear system is known as the Gram matrix of {E m } M m=1 , which takes the form of the identity matrix when {E m } M m=1 are orthonormal.We note that the accuracy of such estimation of matrix H and its condition number depend on the gap between the function-spaces spanned by {H k } K k=1 and {E m } M m=1 as well as on the conditioning of the Gram matrix.In general, a good estimation requires that {E m } M m=1 forms a good approximation for {H k } K k=1 .We now discuss the choice of the regularisation term R in (9).If R is absent from (9), i.e. if λ = 0, then the solution of this minimisation problem is equivalent to the least-squares solution f := ((EH) * EH) −1 (EH) * g.If the regularisation term is given by R(f ) := f 2 , then ( 9) is known as Tikhonov regularisation and its solution is given by f := ((EH) * EH + λI) −1 (EH) * g.However, if EH is badly conditioned, ε > 0 and, additionally, it is known a priori that only a few elements of {H k } K k=1 are sufficient to represent F well, then R(f ) := f 0 is an appropriate choice of the regularisation term.This is termed as the 0 -regularisation, where the 0 -norm of f , i.e. f 0 , is defined as the number of non-zero coordinates in f .The 0 -regularisation bypasses the ill-conditioning by imposing sparsity in the solution FK with respect to {H k } K k=1 .In practice, solving the minimisation problem with such a non-convex 0 -term is computationally difficult, so typically an 1 -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 R, 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 F − FK = (F − F K ) + (F K − FK ), which depends on several factors.The magnitude of the first term F − F K depends of how well F can be represented by its K-term approximation with respect to {H k } K k=1 , and thus it is expected to decrease with increasing K. On the other hand, the magnitude of the second term F K − FK depends on the conditioning of the matrix EH and the error term ε in (8), and thus, it is expected to increase with increasing K when M and N are fixed.In other words, if the resolution K at which we reconstruct is increased, we also need to increase M and N .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 F −F K and δ k are small or they lie in the span of those eigenfunctions corresponding to a small singular value.Thus, for small ε it is required that {H k } K k=1 and {E m } M m=1 form a good approximation for F 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 E may become worse if the span of {E m } M m=1 includes too many eigenfunctions including those corresponding to a small singular value.Loosely speaking, the calibration functions {E m } M m=1 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 {E m } M m=1 .Thus, the use of regularisation in solving (8) becomes crucial in order to obtain a robust solution.

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 → F := F h F v be the complex-vector-valued functions related as in equation (1), where the superscripts h and v correspond to the horizontal and vertical polarisations of the optical vector-field, respectively.The goal is to recover both polarisations F h and F v , which are scalar-valued functions, by using the measurements (2) and (3).Since each of F h and F v depends on both F v and F h , we cannot consider the reconstructions of F v and F v 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 F. 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 m = 1, . . ., M , b := e βi for a fixed β ∈ (0, 2π) and Bv m are some scalar-valued functions.To make clear the motivation to sample according to (11), observe that if {E h m } M m=1 and {E v m } M m=1 are representation systems for F h and F v respectively, then we have that {A m − B m } M m=1 and {A m − b * B m } M m=1 are representation systems for 0, F v and F h , 0 respectively.In other words, can be used to represent the complex vector-valued function F.
Mimicking the reasoning of the previous subsection, we proceed by recovering approximations of F h and F v with respect to some desired representation systems {H h k } K k=1 and {H v k } K k=1 .Namely we aim to recover where we first write these representation systems in terms of the calibration functions as 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 where .
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 f = f h 1 , . . ., , we can define the reconstructions of F h and F v as Observe that using 1 -regularisation in this case imposes sparsity in the reconstructions F h K and F v K with respect to {H h k } K k=1 and {H v k } K k=1 , respectively.Also similarly as before, we note that matrices H h and H v are identities when the reconstruction functions H h k , H v k and the calibration functions E h m , E v m 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 A m and B m given in (11), it may also be possible to measure where c = b ensures that C m = B m and, as we will see shortly below, c must be carefully chosen so that b + c = 2.With this extra information in hand, rather than using (12), the following functions are used to represent the complex vector fields where a := 1/ (1 − b/2 − c/2) is finite given that c was chosen above in such a way that b + c = 2.We then proceed as above but in place of ( 14) obtain where we now have the additional matrix C ∈ C 2N ×M containing the outputs Cm .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.

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.

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 F (j) (x) = R(x) exp(iP (j) (x)), x ∈ I, j = 1, . . ., 8, with the same amplitude R but with different phases P (j) defined on the compact interval I := [−1/2, 1/2].In particular, for illustration purposes we take R(x) := exp(−x 2 ) and P (j) (x) := τ (j) sin(20x), where 0 < τ (1) < • • • < τ (8) < 2π, 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 I is given by {e 2πikx } k∈Z , then for each F (j) we compute its first 20 Fourier coefficients as where k = −10, . . ., 9, and approximate its Fourier transform by the classical Whittaker-Shannon interpolation formula k f (j) k sinc(w − k), w ∈ R. The amplitude of the approximated Fourier transform of each F (j) is shown in the second panel of Figure 1.Finally, we quantify the decay of the Fourier coefficients by the standard deviation σ (j) of a Gaussian function a (j) exp(−(w − c (j) ) 2 /(2(σ (j) ) 2 )) fitted to the amplitude of the approximated Fourier transform on interval w ∈ [−10, 10).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 τ (j) results in an increased standard deviation σ (j) , 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 σ (j) would increase as well.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.

Fourier coefficients of simulated tissue images
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 F (x) := R(x) exp(iP (x)), x ∈ S, where the original space-domain S := [−1/2, −1/2] 2 is discretized into a 700 × 700 grid, while R and P are chosen randomly as we now describe.The phase function P := P (τ,ρ) depends on two given parameters τ and ρ, controlling the amplitude and the frequency of phase oscillations, respectively.Specifically, 800 × 800 pixel-values : Six simulated images with the same amplitude but different phase, which are generated from our model with increasing τ (j) /ρ (j) , j = 1, . . ., 6, so that larger τ (j) /ρ (j) characterises larger phase oscillations.In the scatter plot on the right, we report the sum of parameters σ of the Gaussian fitted to the amplitude of the Fourier transform abs(FT), revealing that increased τ (j) /ρ (j) correlates with larger σ are chosen uniformly at random from [−1, 1], which are then filtered by using MATLAB's function 'imgaussfilt' with the smoothing parameter ρ.Following this step, only 700 × 700 pixels are kept by removing 50 pixels from each boundary and such image is then rescaled so that all phase pixel-values are between [−τ, τ ], τ ∈ [0, π].The amplitude function R is selected as the sum of exp(−50 x 2 2 )/1000 and five additional Gaussian functions exp(− x − c 2 2 /d)/2000 with randomly chosen parameters c and d.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 (ρ (j) , τ (j) ), j = 1, . . ., 6 to create six different functions F (j) , where 0 < τ (1) < • • • < τ (6) ≤ π and 0.025 (6) ) −1 ≤ 0.125 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 2 )) fitted to the absolute value of the Fourier transform approximated from the first 20 × 20 Fourier coefficients, where the Fourier coefficients {f k } k∈I 400 of function F are computed using the formula in (5).In particular, in Figure 2, for each F (j) we report the sum of the standard deviations σ 2 of the fitted Gaussian, thereby observing that increased phase oscillations, i.e. increased τ (j) /ρ (j) , results in slower decay of the corresponding Fourier coefficients, i.e. larger σ Next, in Figure 3, for each of the six different values (ρ (j) , τ (j) ), j = 1, . . ., 6, 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 σ 2 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.
j 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 τ (j) and ρ (j) , j = 1, . . ., 6.For each image, we then approximated its Fourier transform and fitted a Gaussian function whose parameter σ 2 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.

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.

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 abs(F h ) phase(F h ) abs( F h ) phase( F h ) 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 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 1200 × 1200 and they are translated to M = 936 different locations across the input imaging plane.Each output is evaluated at N = 34973 different pixels at the output imaging plane.Considering the two polarisation states, the dimension of the system matrix in ( 16) is therefore 1872 × 69946.

Distal end
Proximal end abs(E h m ) phase(E h m ) abs( Ãh m ) phase( Ãh m ) Figure 5: Amplitude and phase of the horizontal polarisation of one calibration input (A m 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 H h = H v = I, by inverting the linear system in four different ways: 1. the naive inversion f := E * g, which corresponds to the principle of phase conjugation in that it assumes E * E = I, 2. the least-squares approach f := (E * E) −1 E * g, 3. the 2 -regularisation f := (E * E + λI) −1 E * g, and, 4. the 1 -regularisation f := argmin f ∈C 2M g − Ef 2 + λ f 1 using the iterative solver [10].
In particular, we see in Figure 6 that 1 -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 F h and F v are sparse with respect to m=1 -1 -regularisation successfully removes significant noise while preserving the image details.The amount of noise that is removed by 1 -regularisation depends on the strength of the regularisation parameter λ.
Finally, in Figure 7, we reconstruct the holographic image with respect to different representation systems, namely we solve (16) where both H h and H v correspond to the  representation system of a Fourier or a wavelet basis with cardinality K = 1024.Specifically, (i) in the Fourier case, we choose both ii) in the wavelet case, we choose tensor-products of √ K one-dimensional boundarycorrected Daubechies wavelets with four vanishing moments (DB4) from [17].
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 K M , small K 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 M .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

Reconstruction and analysis of biological images
We now apply the framework developed in Sections 2-3 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 In Figure 8, we show the reconstruction of the horizontal polarisation of one healthy image indexed as (1, 1) (first row) and one lesion image indexed as (7, 1) (second row), in both the space-domain (first two columns) and the Fourier-domain (last two columns).Specifically, we reconstruct K = 400 Fourier coefficients per polarisation by solving the linear system ( 14) with 1 -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 unwrapper2 .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 (n, i), n = 1, . . ., 12, 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 (n, i) in the data set, we compute σ .For each tissue sample n = 1, . . ., 12, we then compute the average value σ ) and show a box-plot of these twelve values σ 2 grouped according to their class label 'healthy' or 'lesion', for each polarisation as well as for both polarisations combined.We also compute the p-value of Welch's t-test [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.

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 diag-  nostically 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 trans-

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

Figure 2
Figure2: Six simulated images with the same amplitude but different phase, which are generated from our model with increasing τ (j) /ρ (j) , j = 1, . . ., 6, so that larger τ (j) /ρ (j) characterises larger phase oscillations.In the scatter plot on the right, we report the sum of parameters σ

Figure 6 :
Figure 6: Reconstructed amplitude and phase of the horizontal and vertical polarisations of the holographic image from Figure4with respect to the calibration functions such as those in Figure5, using naive, least-squares, 2 and 1 approaches.The regularisation parameter in the 2 and 1 -regularisation is λ = 0.3 and λ = 0.267, respectively.

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

Figure 9 :
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).