\subsection{Mathematical framework of the forward problem}
\label{subsec:tMatRecovery}
If we define the electric field of light incident on the distal facet of the fibre bundle at a point $(u,v)$ as $T_i (u,v)$ and the electric field incident on the camera sensor at a point $(x,y)$ as $T_o (x,y)$, then we can model the relationship between the two using a spatially dependent point-spread function (or \emph{Green's function}), $h(u,v,x,y)$ as:

\begin{equation}
\label{eq:greenFunc}
T_o(x,y) = \iint h(u, v, x, y) T_i(u,v) \mathrm{d}u \mathrm{d}v
\end{equation}

In the case of free-space propagation, $h(u,v,x,y)$ can be set as $\frac{e^{ikr}}{i r\lambda}$ to give the Fresnel diffraction integral.  However, in our case $h(u,v,x,y)$ is a more complex function that is not known \emph{a priori}. A simple approach to approximate this is to discretise the inputs and outputs and model $h(u,v,x,y)$ with a matrix. To achieve this, we let $T_i(u,v)$ be discretised and then lexicographically ordered into an $n \times 1$ column vector, $\mathbf{x}_h \in \mathbb{C}^n$. We then let $T_o(x,y)$ be discretised and lexicographically ordered into an $m \times 1$ column vector, $\mathbf{y}_h \in \mathbb{C}^m$. In this latter case, an obvious choice of sampling scheme would be to use the values recorded by each pixel of the image sensor. Discretising $h(u,v,x,y)$ leads to the matrix equation:

\begin{equation}
\mathbf{y}_h = \mathbf{A}_{hh} \mathbf{x}_h
\end{equation}

\noindent where $\mathbf{A}_{hh}$ is an $m \times n$ transfer matrix mapping horizontally polarised light at the input to horizontally polarised light at the output.  Repeating the discretisation for the vertically polarised input and output fields, then combining we find:

\begin{equation}
\mathbf{y} =  \left[ \begin{array}{c}
								\mathbf{y}_h \\
								\mathbf{y}_v \end{array} \right]
\end{equation}

\begin{equation}
\mathbf{x} =  \left[ \begin{array}{c}
								\mathbf{x}_h \\
								\mathbf{x}_v \end{array} \right]
\end{equation}

\begin{equation}
\mathbf{A} =  \left[ \begin{array}{cc}
								\mathbf{A}_{hh} & \mathbf{A}_{hv}\\
								\mathbf{A}_{vh} & \mathbf{A}_{vv} \end{array} \right]
\end{equation}

\begin{equation}
\label{eq:forwardEq}
\mathbf{y} = \mathbf{A} \mathbf{x}
\end{equation}

\noindent where $\mathbf{A}$ is the full $2m \times 2n$ transfer matrix.  In practice we record images of the proximal facet of the MCF in the plane of the camera ($\mathbf{y}$) and use these to reconstruct the input fields in the distal facet plane ($\mathbf{x}$).  Equation \ref{eq:forwardEq} can therefore be considered a linear inverse problem: for a known $\mathbf{y}$ and $\mathbf{A}$, the goal is to find $\mathbf{x}$.  One method of solving this problem is to find a pseudo-inverse of $\mathbf{A}$ such that:

\begin{equation}
\mathbf{x} = \mathbf{A}^{-1} \mathbf{y}
\end{equation}

However, we do not typically know $\mathbf{A}^{-1}$ in advance and so we must first determine this before $\mathbf{y}$ can be recovered.  Consider now combining $p$ vectors representing different input fields at the distal facet to form a matrix, $\mathbf{X}$:

\begin{equation}
\label{eq:xMat}
\mathbf{X} =  \left[ \begin{array}{cccc}
\mathbf{x}_1 & \mathbf{x}_2 & \cdots &\mathbf{x}_p
 \end{array} \right]
\end{equation}

For each of these $p$ vectors, the field exiting the distal facet is recorded at the camera plane, creating another matrix, $\mathbf{Y}$:

\begin{equation}
\label{eq:yMat}
\mathbf{Y} =  \left[ \begin{array}{cccc}
\mathbf{y}_1 & \mathbf{y}_2 & \cdots &\mathbf{y}_p
 \end{array} \right]
\end{equation}

Collectively, $\mathbf{X}$ and $\mathbf{Y}$ are termed the \emph{calibration measurements}.  We can then write the following:

\begin{equation}
\mathbf{X} = \mathbf{A}^{-1} \mathbf{Y}
\end{equation}

\begin{equation}
\label{eq:recoveryEq}
\mathbf{X}^T = \mathbf{Y}^T (\mathbf{A}^{-1})^T
\end{equation}

By solving for each of the $p$ columns of $(\mathbf{A}^{-1})^T$ using the respective column of $\mathbf{X}^T$, we can reduce the problem of finding the inverse transmission matrix to $p$ linear inverse problems. Solving each of these provides a single row of $\mathbf{A}^{-1}$:

\begin{equation}
\label{eq:recoveryEq_sing}
\mathbf{x}_r^T = \mathbf{Y}^T \mathbf{a}_r^T
\end{equation}

\noindent where $\mathbf{x}_r$ is the $r^{th}$ row of $\mathbf{X}$ and $\mathbf{a}_r$ is the $r^{th}$ row of $\mathbf{A}^{-1}$.

\subsection{Experimental measurements}
\label{subsubsec:parallelisation}
The full process of producing appropriate sets of inputs and outputs, $\mathbf{X}$ and $\mathbf{Y}$, for recovering $\mathbf{A}^{-1}$ is detailed in Supplementary Figure \ref{fig:matrixRecoveryAlg}. We control the vectors $\mathbf{x}_i$ of Equation \ref{eq:xMat} for $i=1\cdots p$ by displaying an appropriate hologram on SLM1. In this case, each $\mathbf{x}_i$ represents one of 81 translations (over a $9\times9$ grid) in 9$\mu$m steps (limiting spatial resolution) of an array of 12 equispaced spots (Supplementary Figure \ref{fig:holos}a) in one of 3 elliptical polarisation states, resulting in a total of $p=81 \times 3 = 243$ inputs.  Elliptical states are used so as to maximise power transmitted through the fibre, thereby minimising the effect of noise compared to linear states with only half the power.  A minimum of two states are required to resolve the 2 orthogonal polarisation components but further states can reduce noise.  Using 3 states produced a $\sim 33\%$ reduction in noise, but increasing to 4 states provided little improvement as the noise performance was observed to be limited by approximation noise in the basis pursuit algorithm and non-flatness of the hologram.  3 states are therefore selected as a balance between noise reduction and experimental time reduction.  The translation step size is chosen so that every second fibrelet in each dimension is utilised.  This offers a balance between reducing experimental time ($<$1 hour) while still offering sufficient resolution for examining small lesions (size $>50\mu$m).

The choice of 12 spots ensures individual spots are spaced sufficiently far from each other that there is negligible interference between them at the proximal facet of the fibre bundle (even after been scrambled by the matrix $\mathbf{A}$) and they can be considered to behave independently. In a pixel basis representation, the resulting fibre bundle transmission matrix is both \emph{sparse}, and in some sense \emph{local} (technically, a \emph{band matrix}). Hence, not only are most entries zero, but the non-zero entries are clustered together. Consequently, a single image containing 12 spots can be split up to produce 12 separate characterisation images. Each image simply uses the relevant subsection of the full image and sets all other pixels to zero. This means that there are effectively $12 \times 243 = 2916$ characterisation images from only 243 actual measurements. Each of the 81 spatial translations effectively provides 12 pixels of resolution at the distal facet (one per spot).  Therefore, in total there are $12 \times 81 = 972$ pixels that can be spatially resolved at the distal facet.


When using a phase retrieval algorthm (see Online Methods), the phase recovered for each image has an arbitrary global offset. If measured na\"{\i}vely, the columns of $\mathbf{Y}$ in Equation \ref{eq:yMat} would have an arbitrary phase offset between them. To compare different phase images, a fixed phase reference is created during the characterisation measurements by displaying a superposition of two holograms on SLM1: the first hologram creates the translatable array of spots and the second creates a single spot that is kept in a fixed position. The fixed position is a set of fibrelets at the edge of the MCF, as shown in Supplementary Figure \ref{fig:matrixRecoveryAlg}a. Recovered phase values are normalised to the phase of this spot at the proximal facet of the fibre. This is effectively a form of common-path interferometry. Such a global phase offset is also termed \emph{piston} and can be thought of as a phase surface polynomial of order 0 \cite{Loterie2015a}. Here, we observed experimentally that in addition to this 0$^{th}$ order correction, a 1$^{st}$ order time-dependent phase correction,termed \emph{phase tilt}, was also required. The origin of the time-dependent phase tilt is likely due to temperature, affecting the SLM and optical mounts \cite{Garcia-Marquez2011}. Air currents could also play a role, however the experiment is entirely contained in an optical box so their effect should be negligible. The 1$^{st}$ order phase correction term is tracked by measuring a reference image at 60 second intervals during characterisation. The relative phase tilt between reference images is computed and interpolation with respect to time is then used to estimate the correction to be applied to characterisation measurements, $\mathbf{y_i}$.

\subsection{Solving the inverse problem}
The region of the image sensor used for data collection is 1200$\times$1200 pixels and so the inverse transmission matrix, $\mathbf{A}^{-1}$, has 2$\times$1200$\times$1200 =  2.88 million columns when the two polarisation states are considered. Similarly, the 81 translations of 12 spots in three polarisation states at the proximal facet gives a total of 2916 rows for $\mathbf{A}^{-1}$. To store this 2916 $\times$2880000 matrix at single floating point precision would require $>$30Gb of memory. This can make recovery of the matrix and subsequent operations infeasibly slow.

However, because the transmission matrix is sparse and the few non-zero elements are grouped together (i.e. a banded structure as discussed above), it is possible to reconstruct only the parts of $\mathbf{A}^{-1}$ where the non-zero elements are located. Consider again Equation \ref{eq:recoveryEq_sing} in which we seek to solve for row $r$ of $\mathbf{A}^{-1}$, $\mathbf{a}_r$: 

\begin{equation}
\mathbf{x}_r^T = \mathbf{Y}^T \mathbf{a}_r^T
\end{equation}

Here, $\mathbf{x}_r$ is a single row of $\mathbf{X}$ and thus corresponds to a single spatial position across all input calibration fields at the distal facet.  Similarly, $\mathbf{a}_r$ is a single row of $\mathbf{A}^{-1}$ and therefore represents the coupling to the single spatial position represented by $\mathbf{x}_r$ from all spatial positions at the proximal facet, represented columns of $\mathbf{Y}$.  However, because we know $\mathbf{A}^{-1}$ is sparse and banded, the only contributions to $\mathbf{x}_r$ will come from points in a confined region on the distal facet.  The distal and proximal facet coordinate axes are the same as a result of an affine transform (see Online Methods) and so if the coordinates of the spatial position of $\mathbf{x}_r$ are $(a,b)$, then the relevant points on the distal facet will be neighbouring $(a,b)$.  We therefore select a set of spatial points on the distal facet (i.e. columns $\mathbf{Y}$) that are located within a square centred on $(a,b)$ with side length 68.2$\mu$m.  This length is exactly equal to the spacing between adjacent spots in the array used to characterise the fibre (Supplementary Figure \ref{fig:holos}a), which was in turn set so as to avoid interference of spots at the distal facet (Online Methods). This process is illustrated in Supplementary Figure \ref{fig:matrixRecoveryAlg}.  Equation \ref{eq:recoveryEq_sing} then becomes:

\begin{equation}
\label{eq:subEq}
\mathbf{x}_{r}^T = \mathbf{Y}_{sub}^T \mathbf{a}_{r,sub}^T
\end{equation}

\noindent where $\mathbf{Y}_{sub}$ is the selected subset of rows of $\mathbf{Y}$ and 
 $\mathbf{a}_{r,sub}$ is the corresponding subset of columns of $\mathbf{A}^{-1}$. Typically, $\mathbf{a}_{r,sub}$ has around 60,000 columns over 2 polarisations, a 50-fold reduction from the original 2.88 million.  However, given that $\mathbf{x}_r$ has only $\sim$3000 columns, the problem is underdetermined by a factor of 10.  The problem can still be solved by exploiting the prior knowledge that each $\mathbf{a}_{r,sub}$ ought to be sparse because spatial positions towards the edge would likely couple little or no power directly from the central positions (in agreement with observations). This is achieved using an $L1$ regulariser so that the solution to Equation \ref{eq:subEq} becomes \cite{Zhang2015}:
 
\begin{equation}
\min ||\mathbf{a}_{r,sub}||_1 \mathrm{~~subject~to~~} ||\mathbf{x}_r^T - \mathbf{Y}_{sub}^T \mathbf{a}_{r,sub}^T||_2 < \delta
\label{eq:solver}
\end{equation}

\noindent where $\delta$ is a parameter that can be adjusted to favour a more or less sparse solution.  Equation \ref{eq:solver} is solved numerically using a basis pursuit denoising algorithm implemented in the SPGL1 package \cite{spgl1:2007}.

\pagebreak
\noindent
\Large\textbf{Supplementary Note 2:}\\
\Large\textbf{Evaluation of polarisation parameters}\\
\normalsize
\setcounter{section}{2}
\setcounter{subsection}{0}
\subsection{Polarisation model}
\label{subsec:polModel}
When polarised light interacts with the turbid medium of tissue, it may simultaneously encounter: depolarisation, arising from multiple scattering of sub-cellular structures; birefringence and diattenuation, arising from anisotropic muscle fibres and structural proteins; and circularity caused by chiral molecules and structures. 

The propagation of polarised light and its interaction within a sample can be described by either Jones (field-based) or Stokes-Mueller (intensity-based) calculus \cite{Simon1982}. Stokes-Mueller is more commonly used for tissue polarimetry as it can handle partial polarisation states, represented by the depolarisation parameter  \cite{Ghosh2011,Qi2016,Yao1999,Pierangelo2013a}, though Jones calculus has found use in polarisation-sensitive optical coherence tomography of tissue \cite{Park2009,Li2015}.  However, if the coherence length of the light source is substantially longer than possible scattering paths in the tissue, then the light always has a well-defined polarisation state \cite{Simon1982}.  This corresponds to a light source with a high degree of temporal coherence, i.e. monochromatic light. The laser used here has a coherence length of around 10 m compared to the typical transport mean free path for tissue of $<$ 1 mm \cite{Ntziachristos2010} and so this condition is satisfied. In addition to this temporal coherence, the light field must also exhibit spatial coherence in order to have a well-defined polarisation state.  This can be violated, for example, if the integration area of each camera pixel contains  different polarisation states.  However, in our case the light field from the sample is first passed through the fibre bundle where the single-mode propagation in each fibrelet admits only a single polarisation state -- a filter of sorts.  Each fibrelet is imaged on the camera using $>$10 pixels so it can therefore be assumed that the polarisation state is uniform within each pixel area. Depolarisation is therefore considered to be negligible in this work.

In order to examine the properties of birefringence, diattenuation (also called \emph{linear dichroism}) and circularity, we choose a generalised decomposition of the Jones matrix into an arbitrary elliptical retarder followed by a partial linear polariser. The elliptical retarder accounts for birefringent retardance and circularity. The partial linear polariser accounts for the fact that anisotropic molecules might absorb light preferentially along some axis. This is illustrated in Supplementary Figure \ref{fig:polModel}. It is noted that this model does not account for some polarimetric properties, notably circular dichroism. Circular dichroism has applications in spectroscopic analysis of chiral molecules (e.g. many proteins and biological molecules) but produces weak signals so is typically only effective for molecules in concentrated solution \cite{Fasman1996}. Further, little or no work has been produced to date linking circular dichroism at optical wavelengths to imaging of cancer. Therefore, for simplicity it is neglected in this work.

Having measured a set of $n$ output polarisation vectors from the sample $\mathbf{V}(x,y) = \left[\mathbf{v}_1(x,y) \cdots \mathbf{v}_n(x,y)\right]$ and knowing the corresponding $n$ sample illumination polarisation vectors $\mathbf{U}(x,y) = \left[\mathbf{u}_1(x,y) \cdots \mathbf{u}_n(x,y)\right]$ the 2$\times$2 complex Jones matrices $\mathbf{J}(x,y)$ representing how the sample maps input polarisation state to output polarisation state at each pixel can be determined through:

\begin{equation}
\label{eq:forwardModelJones}
\mathbf{V}(x,y) = \mathbf{J}(x,y) \mathbf{U}(x,y)
\end{equation}

\noindent where $\mathbf{J}(x,y)$ is the Jones matrix of a particular pixel at $(x,y)$. We express our decomposition of $\mathbf{J}(x,y)$ as:

\begin{equation}
\label{eq:jonesMat}
\mathbf{J}(x,y) = \mathbf{A}_{pol} \mathbf{A}_{ret}
\end{equation}

\begin{equation}
\label{eq:polMat}
\mathbf{A}_{pol} = \left( \begin{array}{cc}
\cos \theta_D & \sin \theta_D \\
-\sin \theta_D & \cos \theta_D \end{array} \right)
\left( \begin{array}{cc}
\sqrt{1+D} & 0 \\
0 & \sqrt{1-D} \end{array} \right)
\left( \begin{array}{cc}
\cos \theta_D & -\sin \theta_D \\
\sin \theta_D & \cos \theta_D \end{array} \right)
\end{equation}

\begin{equation}
\label{eq:retMat}
\mathbf{A}_{ret} = \left( \begin{array}{cc}
1 & 0 \\
0 & e^{i\xi}\end{array} \right)
\left( \begin{array}{cc}
\cos \theta_{\eta} & \sin \theta_{\eta} \\
-\sin \theta_{\eta} & \cos \theta_{\eta} \end{array} \right)
\left( \begin{array}{cc}
e^{i\eta/2} & 0 \\
0 & e^{-i\eta/2} \end{array} \right) \cdot \left( \begin{array}{cc}
\cos \theta_{\eta} & -\sin \theta_{\eta} \\ \sin \theta_{\eta} & \cos \theta_{\eta} \end{array} \right)
\left( \begin{array}{cc}
1 & 0 \\
0 & e^{-i\xi}\end{array} \right)
\end{equation}

\noindent allowing us to recover the following 5 polarisation properties: diattenuation, $D$, diattenuation axis orientation, $\theta_D$, retarder circularity, $\xi$, retardance, $\eta$, and retardance axis orientation, $\theta_{\eta}$.  From these equations, it is observed that the inferred parameters are all real numbers limited to the following ranges:

\begin{equation}
\label{eq:paramRanges}
\theta_{\eta} \in (-\pi/2, \pi/2],~~~~~~~\eta \in (-\pi, \pi],~~~~~~~\xi \in (-\pi, \pi], ~~~~~~~D \in [-1, 1],~~~~~~~\theta_D \in (-\pi/2, \pi/2]
\end{equation}

\subsection{Ordering of polarisation elements}
\label{subsec:polOrdering}
The model for inference of polarimetric parameters, introduced in Section \ref{subsec:polModel} and illustrated in Supplementary Figure \ref{fig:polModel}, comprises an elliptical retarded followed by a partial polariser. In the Jones formalism these are each represented by $2\times2$ matrices, which are multiplied together to obtain overall behaviour.  However, matrix multiplication is non-commutative and so it is necessary to consider an alternative model in which the order of these two elements is reversed.  By adapting Equation \ref{eq:jonesMat} we get an alternative composite Jones matrix:

\begin{equation}
\label{eq:jonesMat_rev}
\mathbf{\hat{J}}(x,y) = \mathbf{A}_{ret} \mathbf{A}_{pol} 
\end{equation}

If we then expand both Equation \ref{eq:jonesMat} and \ref{eq:jonesMat_rev} using Equations \ref{eq:polMat} and \ref{eq:retMat}, we get:

\begin{equation}
\label{eq:polOrder1}
\mathbf{J}(x,y) = \left(\begin{array}{c c}
K & L \\
M & N 
\end{array}\right)
\end{equation}

\noindent and

\begin{equation}
\mathbf{\hat{J}}(x,y) = \left(\begin{array}{c c}
\hat{K} & \hat{L} \\
\hat{M} & \hat{N} 
\end{array}\right)
\end{equation}

\noindent where

\begin{equation}
\begin{array}{c}
K = \left({\sin\left(\theta_D\right)}^2 \sqrt{1 - D} + {\cos\left(\theta_D\right)}^2 \sqrt{D + 1}\right) \left(\cos\left(\frac{\eta}{2}\right) + \cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right) i\right) \\
- \frac{1}{2} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \sin\left(2\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\sin\left(\xi\right) + i\cos\left(\xi\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
\hat{K} = \left({\sin\left(\theta_D\right)}^2 \sqrt{1 - D} + {\cos\left(\theta_D\right)}^2 \sqrt{D + 1}\right) \left(\cos\left(\frac{\eta}{2}\right) + \cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right) i\right) \\
- \frac{1}{2} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \sin\left(2\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(- \sin\left(\xi\right) + i\cos\left(\xi\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
L = i \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2\theta_{\eta}\right) \left({\sin\left(\theta_D\right)}^2 \sqrt{1 - D} + {\cos\left(\theta_D\right)}^2 \sqrt{D + 1}\right) \\
- \frac{1}{2} \sin\left(2\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) - i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
\hat{L} = i\mathrm{e}^{i\xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \left({\sin\left(\theta_D\right)}^2 \sqrt{D + 1} + {\cos\left(\theta_D\right)}^2 \sqrt{1 - D}\right) \\
- \frac{1}{2}\sin\left(2 \theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
M = \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \left({\sin\left(\theta_D\right)}^2 \sqrt{D + 1} + {\cos\left(\theta_D\right)}^2 \sqrt{1 - D}\right) \\
- \cos\left(\theta_D\right) \sin\left(\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
\hat{M} = \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \left({\sin\left(\theta_D\right)}^2 \sqrt{1 - D} + {\cos\left(\theta_D\right)}^2 \sqrt{D + 1}\right) \\
-  \cos\left(\theta_D\right)\sin\left(\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
N = \left({\sin\left(\theta_D\right)}^2 \sqrt{D + 1} + {\cos\left(\theta_D\right)}^2 \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)\\
 -\frac{1}{2}\sin(\frac{\eta}{2})\sin\left(2\theta_D\right)  \sin\left(2\theta_{\eta}\right)\left( \sqrt{D + 1} - \sqrt{1 - D}\right) \left(-\sin(\xi)+i\cos(\xi) \right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
\hat{N} = \left({\sin\left(\theta_D\right)}^2 \sqrt{D + 1} + {\cos\left(\theta_D\right)}^2 \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)\\
 - \frac{1}{2}\sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_D\right) \sin\left(2 \theta_{\eta}\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\sin\left(\xi\right) + i\cos\left(\xi\right)\right)
\end{array}
\end{equation}

Taking the difference between these two matrices, we find that:

\begin{equation}
\label{eq:polDiffMat}
\mathbf{J}(x,y) - \mathbf{\hat{J}}(x,y) = \mathbf{J'}(x,y) = \left(\begin{array}{c c}
K' & L' \\
M' & N' 
\end{array}\right)
\end{equation}

\noindent where

\begin{equation}
\begin{array}{c}
K' = K - \hat{K} = - \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \sin\left(2\theta_D\right) \sin\left(\xi\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
L' = L - \hat{L} = \left(\sqrt{D + 1} - \sqrt{1 - D}\right)\left( i \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2\theta_{\eta}\right)\cos(2\theta_D)+ i\sin\left(2\theta_D\right) \left(\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
M' = M - \hat{M} = -\left(\sqrt{D + 1} - \sqrt{1 - D}\right) \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \cos\left(2\theta_D\right) 
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
N' = N - \hat{N} = \left( \sqrt{D + 1} - \sqrt{1 - D}\right) \sin(\frac{\eta}{2})\sin\left(2\theta_D\right)  \sin\left(2\theta_{\eta}\right) \sin(\xi)
\end{array}
\end{equation}

The error terms $K'$, $L'$, $M'$ and $N'$ are all proportional to $\left( \sqrt{D + 1} - \sqrt{1 - D}\right)$, and so as $D \to 0$, so these error terms become negligible.  To investigate the range of $D$ for which the error between these two models is negligible, we performed a Monte-Carlo simulation over the broadest possible range of parameters:

\begin{equation*}
\theta_{\eta} \sim U(-\pi/2,\pi/2)
\end{equation*}

\begin{equation*}
\eta \sim U(-\pi,\pi)
\end{equation*}

\begin{equation*}
\xi \sim U(-\pi,\pi)
\end{equation*}

\begin{equation*}
\theta_D \sim U(-\pi/2,\pi/2)
\end{equation*}

\begin{equation*}
D \sim N(0,\sigma_D)
\end{equation*}

100,000 samples were drawn from these distributions for different values of $\sigma_D$ using a MATLAB script.  For each sample drawn, Jones matrices are computed for the two orderings ($\mathbf{J}$ and $\mathbf{\hat{J}}$ respectively) of the polariser and retarder, enabling computation of the error matrix $\mathbf{J}'$ of Equation \ref{eq:polDiffMat}.  The elements are then normalised and summed giving an error metric, $Q$, defined (with reference to Equations \ref{eq:polOrder1} and \ref{eq:polDiffMat}) as:

\begin{equation*}
Q = \frac{1}{4} \left( \left|\frac{K'}{K}\right| + \left|\frac{L'}{L}\right| + \left|\frac{M'}{M}\right| + \left|\frac{N'}{N}\right| \right)
\end{equation*}

\begin{equation}
Q = \frac{1}{4} \left( \left|\frac{K - \hat{K}}{K}\right| + \left|\frac{L - \hat{L}}{L}\right| + \left|\frac{M - \hat{M}}{M}\right| + \left|\frac{N - \hat{N}}{N}\right| \right)
\end{equation}

This represents the average absolute proportional deviation introduced by reversing polarisation elements.  The measured cumulative distribution function for different values of $\sigma_D$ is shown in Supplementary Figure \ref{fig:polMonteCarlo}.  In the samples of mouse oesophageal tissue presented here, the true value of $\sigma_D$ is estimated to be 0.0081.  As this is less than 0.01, it is observed from the cumulative density function that in $>$99.5\% of possible cases the error introduced by a different ordering of polarisation elements will be $<$5\%. If a 10\% error threshold is considered acceptable, then a value of $\sigma_D = $0.05, 5-times the measured value, would still be sufficiently accurate in $>$96\% of cases.  Of course, if samples are expected to have larger values of $D$ then the two orderings would give different results and it would be necessary to determine which ordering had the stronger biological underpinning and best explained the measured data.  Such a `model comparison' could be achieved through an extension of the Bayesian inference approach of Supplementary Section \ref{subsec:bayesInf}.

\subsection{Bayesian inference}
\label{subsec:bayesInf}
The polarisation model is fitted to the measured data in order to extract the 5 polarisation properties ($\theta_D$, $D$, $\xi$, $\theta_{\eta}$, $\eta$) using a Bayesian inference approach, which has been found to be more robust to noise than clustering or pseudo-inverse approaches \cite{Zallat2008,Kasaragod2014}. In this work maximum likelihood estimates are determined during this inference and are used to plot images.  These could be computed more quickly by deriving analytical expressions via eigenvalue decomposition of the Jones matrix \cite{Makita2010, Lu1996}.  However, such an approach does not easily permit estimation of joint posterior probabilities of neighbouring pixels, which is essential for compensating beam misalignment (see Section \ref{subsubsec:beamMisalignment}). Therefore, we numerically infer joint maximum likelihood values.

To fit the model, we evaluate the \emph{posterior} probability of the parameter set as a function of the data using Bayes' theorem:

\begin{align*}
p\left[ D, \theta_D, \eta, \theta_{\eta}, \xi | \mathbf{U}(x,y), \mathbf{V}(x,y)\right. & \left.\right] \propto \\
p\left[\mathbf{V}(x,y) = \right. & \left. \mathbf{J}(D, \theta_D, \eta, \theta_{\eta}, \xi,x,y) \mathbf{U}(x,y)\right] \cdot p(D, \theta_D, \eta, \theta_{\eta}, \xi)
\end{align*}


The values in $\mathbf{V}(x,y)$ represent measured complex quantities and so are assumed to be independent and drawn from complex Gaussian distributions as:

\begin{equation}
\label{eq:complexNormal}
\mathbf{v}(x,y) \sim \mathcal{CN}\left[\mathbf{J}(D, \theta_D, \eta, \theta_{\eta}, \xi,x,y) \mathbf{u}(x,y), \sigma^2 \mathbf{I}\right]
\end{equation}

\noindent where $\mathbf{v}(x,y)$ is a column of $\mathbf{V}(x,y)$, $\mathbf{u}(x,y)$ is a column of $\mathbf{U}(x,y)$, $\sigma^2 \mathbf{I}$ is the covariance matrix, and $\mathcal{CN} (\mathbf{\mu},\mathbf{\Sigma})$ is a 2-D complex Gaussian distribution of mean $\mathbf{\mu}$ and covariance $\mathbf{\Sigma}$. The noise standard deviation, $\sigma$ is an additional parameter that can be inferred from the data.  However, after initial processing of several samples, the value was fixed at its inferred maximum likelihood value of 0.4 for subsequent calculations in order to increase speed.  We wish to use a broad prior to minimise any bias in the results.  The broadest possible prior would be a 5-dimensional uniform distribution over all the polarisation parameters (ignoring $\sigma$ by assuming it is now fixed).  In this case, the prior distributions for each individual parameter are independent so that we can write:

\begin{equation}
p(D, \theta_D, \eta, \theta_{\eta}, \xi) = p(D)p(\theta_D)p(\eta)p(\theta_{\eta})p(\xi)
\end{equation}

It is assumed that in the case of a slightly more restrictive prior, this independence assumption will still be valid to a good approximation and so it used for the rest of the analysis presented here.  Priors for parameters representing angles, i.e. circular quantities, are modelled with \emph{von Mises} distributions. The von Mises distribution is a continuous probability distribution defined on a circular domain of $[-\pi,\pi]$ that wraps around, giving a probability density function of $\frac{e^{\kappa \cos (x-\mu)}}{2 \pi I_0(\kappa)}$ where $\mu$ is the mean, $1/\kappa$ is analogous to variance, and $I_0(\cdots)$ is a modified Bessel function of order 0. For $\theta_D$, $\theta_{\eta}$ and $\eta$, the $\kappa$ is set to $\infty$ providing the broadest possible prior. For $\xi$, a restricted value of $\kappa = 1$ is used.  This is because when performing inference of known test targets, it was found that this parameter is prone to over-fitting and that a slightly restrictive prior substantially prevents this.  Further, previous studies that have shown that even in samples with high concentrations of chiral molecules such as glucose, the measured circularity is still relatively small \cite{Guo2007}.  It is therefore not expected that this parameter will be very large for the phantoms or tissue samples measured. Diattenuation, $D$, is the only non-circular parameter (i.e. it is not an angle) and so its prior is a truncated Gaussian over its domain $[-1,1]$ with mean 0 and broad variance.  Again, the broad (as opposed to infinite) variance is used to prevent overfitting and is consistent with previous work that has found the values of $D$ to be small (typically $<$0.2) in real tissue samples \cite{Qi2016,Pierangelo2013a}.  It is noted that the ability to find more accurate solutions in the presence of noise or incomplete information using such prior knowledge is a major benefit of the Bayesian inference approach.

The inference of parameters from data was performed numerically using the STAN package \cite{StanMan}, which produces samples from \emph{a posteriori} distributions using a Markov Chain Monte Carlo simulation engine.  It can also provide maximum likelihood estimates through an optimisation engine, which is used here on account of its high speed, and can also perform variational inference.  

\subsection{Compensation for beam misalignment}
\label{subsubsec:beamMisalignment}
The holographic endoscope design exploits two orthogonally polarised arms that are split then recombined before illuminating the sample. This is achieved using a custom housing for the mirrors and polarising beam splitters, which provides a compact way of transferring the beam path between two optical axes at different heights. As a result of this dual arm approach, any slight misalignment of the two polarisation arms leads to a spatially varying illumination polarisation state at the sample, illustrated in Supplementary Figure \ref{fig:waveplateImaging}a. For samples with small retardance this misalignment is not apparent because it is present for both the characterisation of the fibre transmission matrix and the sample imaging steps, and so the relative change in polarisation state is correct.

However, for samples that exhibit large values of retardance, $\eta$, and in which the retardance axes are not aligned with the polarisation axes of the illumination (i.e. ($\theta_{\eta} \neq 0,\pm\pi/2$), light is coupled between the slightly misaligned horizontally and vertically polarised beams. This produces, in the illumination polarisation basis, wavefronts with different tilts in the two polarisation arms.  This manifests as an apparent spatially varying retardance, illustrated in Supplementary Figure \ref{fig:waveplateImaging}a. This apparent retardance does not appear when the polarisation axes are aligned with the retardance axes. Therefore, the effect is removed by re-expressing Equation \ref{eq:forwardModelJones} in a linear polarisation basis at angle $\hat{\theta}_{\eta}$ to the illumination polarisation axes:

\begin{equation}
\label{eq:rotatedPolAxes}
\mathbf{V}(x,y) = \left. \mathbf{J}(x,y) \right|_{\theta_{\eta}=0}  
\left( \begin{array}{cc}
\cos \hat{\theta}_{\eta} & \sin \hat{\theta}_{\eta} \\
-\sin \hat{\theta}_{\eta} & \cos \hat{\theta}_{\eta} \end{array} \right)
\mathbf{U}(x,y)
\end{equation}


We then apply Bayesian inference in this new, rotated basis with $\hat{\theta}_{\eta}$ taking the place of $\theta_{\eta}$ as the inferred orientation angle of the retardance axes. To find the $\hat{\theta}_{\eta}$ that minimises the spatial tilt of the retardance (confirming the new basis is aligned with the retardance axes) we compute the joint probability for a small set of $R$ neighbouring pixels surrounding $(x,y)$.  For computational tractability, we employ a technique from variational Bayesian inference that approximates the posterior probability density function as a product of distributions \cite{MacKay2003}:

\begin{align*}
p[ D(x_1,y_1)\cdots D(x_R, y_R), &\theta_D(x_1,y_1)\cdots, \eta(x_1,y_1)\cdots, \hat{\theta}_{\eta}(x_1,y_1)\cdots, \xi(x_1,y_1)\cdots | \mathbf{U}(x_1,y_1)\cdots, \mathbf{V}(x_1,y_1)\cdots] \\
&= \prod_{r=1}^R p[D(x_r,y_r),\theta_D(x_r,y_r),\eta(x_r,y_r), \hat{\theta}_{\eta}(x_r,y_r), \xi(x_r,y_r) | \mathbf{U}(x_r,y_r), \mathbf{V}(x_r,y_r)]
\end{align*}

\noindent where the term inside the product is computed using equations \ref{eq:complexNormal} and \ref{eq:rotatedPolAxes}. This approximation implies independence of the conditional probabilities of polarisation parameters of neighbouring pixels.  This is a reasonable assumption because the conditional probability distributions (that is, the probability distributions given the experimentally measured Jones matrices) are effectively accounting for noise, which is assumed to be produced by an independent random process at each pixel.

The maximum likelihood probability is maximised when neighbouring pixels have similar conditional distributions for their polarisation parameters. This implicitly biases towards values of $\hat{\theta}_{\eta}$ where the spatial variation of retardance is minimised, i.e. when the means of the retardance distributions across neighbouring pixels are most similar.

It is not known \emph{a priori} whether the phantoms and tissue samples measured will exhibit significant birefringence.  However, the advantage of the above approach is that it will still provide the best estimate of parameters even if no birefringent tilt artefact is present. This is because the probability distribution of $\hat{\theta}_{\eta}$ tends towards being uniform when $\eta \xrightarrow{}$0, which is the same behaviour expected for $\theta_{\eta}$ (i.e. without the correction).  Therefore, the correction is applied to all measurements taken and will have more effect where it is naturally required.

The number of adjacent pixels used for the correction is set as $R$ = 15.  For optimal resolution, this value should be as small as possible because it effectively reduces resolution by low-pass filtering the image.  On the other hand, larger values can correct for smaller tilt artefacts.  $R$ = 15 was found experimentally to be the smallest value that was sensitive enough to correct for the degree of tilt artefact typically observed for the samples measured (which is in turn dependent on the degree of system misalignment). The outcome of this correction is validated by rotating a half-waveplate in the sample plane (Supplementary Figure \ref{fig:waveplateImaging}).

\subsection{Degeneracy of solutions}
We are decomposing a single matrix into the product of 8 matrices (Equations \ref{eq:jonesMat}, \ref{eq:polMat} and \ref{eq:retMat}), each comprising sinusoidal or complex exponential terms, hence there will invariably be degeneracy of the solutions. For example, consider the matrix $\mathbf{A}_{ret}$.  The 5 matrices of equation \ref{eq:retMat} can be multiplied to give:

\begin{equation}
\mathbf{A}_{ret} = \left( \begin{array}{cc}
e^{i\eta/2} \cos^2 \theta_{\eta} + e^{-i\eta/2} \sin^2 \theta_{\eta} & \left(e^{i\eta/2} - e^{-i\eta/2}\right) e^{-i \xi} \cos \theta_{\eta} \sin \theta_{\eta} \\
\left(e^{i\eta/2} - e^{-i\eta/2}\right) e^{i \xi} \cos \theta_{\eta} \sin \theta_{\eta} & e^{i\eta/2} \sin^2 \theta_{\eta} + e^{-i\eta/2} \cos^2 \theta_{\eta}  \end{array} \right)
\end{equation}

By substituting $\theta_{\eta} \rightarrow \theta_{\eta} + \pi/2$ and $\eta \rightarrow -\eta$, the first row, first column element of $\mathbf{A}_{ret}$ becomes:

\begin{equation*}
e^{-i\eta/2} \cos^2 (\theta_{\eta}+\pi/2) + e^{i\eta/2} \sin^2 (\theta_{\eta}+\pi/2)= e^{i\eta/2} \cos^2 \theta_{\eta} + e^{-i\eta/2} \sin^2 \theta_{\eta} 
\end{equation*}

The element of $\mathbf{A}_{ret}$ is therefore unchanged; the other 3 elements are also unchanged under this substitution.  Therefore, the parameter sets  $(\theta_{\eta}, \eta)$ and $(\theta_{\eta} + \pi/2, -\eta)$ are degenerate.  

Considering the ranges of all parameters, outlined in Equation \ref{eq:paramRanges}, there is in fact an 8-fold ambiguity for each point in the 5-D parameter space, comprising the following parameter sets:

\begin{equation}
\begin{array}{ccccccc}
(&\theta_{\eta}, &\eta,&\xi,&D,&\theta_D&)\\
(&\theta_{\eta}+\pi/2,&-\eta,&\xi,&D,&\theta_D&)\\
(&-\theta_{\eta},&\eta,&\xi+\pi,&D,&\theta_D&)\\
(&-\theta_{\eta}+\pi/2,&-\eta, &\xi+\pi,&D, &\theta_D&)\\
(&\theta_{\eta},&\eta, &\xi, &-D, &\theta_D + \pi/2&)\\
(&\theta_{\eta}+\pi/2, &-\eta, &\xi, &-D, &\theta_D + \pi/2&)\\
(&-\theta_{\eta}, &\eta, &\xi+\pi, &-D, &\theta_D + \pi/2&)\\
(&-\theta_{\eta}+\pi/2, &-\eta, &\xi+\pi, &-D, &\theta_D + \pi/2&)
\end{array}
\end{equation}

\noindent where parameters are defined to wrap around their ranges, listed in Equation \ref{eq:paramRanges}, via a modulus operation. If we were to randomly select one of the 8 degenerate solutions for each pixel, it could give the misleading impression that polarisation properties exhibit a high degree of variation.  To avoid this scenario, we select a point in 5-D parameter space (here we use $(\theta_{\eta} = \pi/4, \eta = \pi, \xi = 0, D = 1, \theta_D = \pi/4)$) and choose the degenerate solution closest to that point.  This ensures the maximum degree of smoothness in the resulting polarisation images. As a result, observed variations can be confidently classed as structural tissue features rather than artefacts. 

\pagebreak
\noindent
\Large\textbf{Supplementary Note 3:}\\
\Large\textbf{Preparation and characterisation of tissue mimicking phantoms}\\
\normalsize
\setcounter{section}{3}
\setcounter{subsection}{0}

\subsection{Preparation of optically scattering and birefringent phantoms}
All chemicals were purchased from Sigma Aldrich unless otherwise stated. Birefringent phantoms were made from acrylamide/bis-acrylamide (30\% concentration) solution as the base material. Ammonium persulfate (APS) and tetramethylethylenediamine (TEMED) were used as catalysts for the co-polymerisation of bis-acrylamide and acrylamide gels. Phantoms were fabricated inside a ventilated fumehood cupboard by mixing 30 mL of Acrylamide/Bis-acrylamide with 600 $\mu$L of APS and 40 $\mu$L of TEMED under thorough vortexing.  The mixture was poured immediately into petri dishes and left for approximately 30 minutes until the polymerisation was complete. Birefringent phantoms having thickness $\sim$2 mm were fabricated by pouring 5 mL of the mixture into the petri dishes. The set phantoms were then extracted from the petri dishes and were mounted on to the mechanical stretcher assembly as shown in Supplementary Figure \ref{fig:stretcher}.

Agar-based scattering phantoms were fabricated using 1.5 w/v\% agar solution as the base material as described elsewhere \cite{Joseph2017}. Briefly, agar powder was dissolved in de-ionised water and the solution was heated up to 100$^o$C inside a microwave oven. Pre-warmed intralipid and nigrosin dye solutions were added into the cooling agar solution to introduce scattering and absorption. A vortex mixer was used to homogenise the phantom solution and prevent separation of its components. The mixture was pipetted into petri dishes and the solutions were allowed to cool under room temperature to form phantoms having $\sim$2 mm thickness. 

\subsection{Characterisation of optical absorption and scattering}
A double integrating sphere (DIS) system\cite{Pickering1993} was used to determine the reduced scattering and absorption coefficients of the phantoms. The DIS set-up uses two highly reflective spheres (Labsphere), with Lambertian surfaces made from polytetrafluoroethylene (PTFE), to capture and quantify the amount of transmitted and reflected light. A broadband tungsten halogen lamp (AvaLight-HAL-MINI, Avantes) and two fiber-optic spectrometers (AvaSpec-ULS2048, Avantes) served as the light source and detectors respectively. The inverse-adding doubling algorithm (IAD) was used to compute the optical properties of the phantoms from the reflectance, transmittance and reference measurements \cite{Prahl1993}.

\subsection{Characterisation of stretching deformation and associated birefringence}
Deformation contours across the birefringent phantoms mounted in the stretcher were assessed by drawing a grid pattern onto the polyacrylamide gel using a marker pen. Qualitative observation of the contours after stretching revealed a macroscopic non-linearity in the deformation, however, over the small field of view of the endoscope the stretch was approximately linear. Due to the constant volume of the material, the change in the thickness upon stretching must be accounted for when plotting retardance, by considering the change in area between the two clamps and dividing the measured retardance by the inferred width.


\pagebreak
\noindent
\Large\textbf{Supplementary Note 4:}\\
\Large\textbf{Calculation of entropy and expectation in endoscopic images}\\
\normalsize
\setcounter{section}{4}
\setcounter{subsection}{0}
\subsection{Distribution fitting}
To identify lesions in phase and polarisation images, two image filters/operators are used (see main manuscript): spatial entropy, $H()$, and spatial expectation, $E()$.  These can both determined under a unified framework that estimates the probability distribution of spatial points and then provides the differential entropy of the inferred distribution as $H()$ and the mean as $E()$.  This framework is also easily extended to joint distributions, e.g. for co-occurrence metrics or for multivariate entropy.

The distribution fitting process is detailed in Supplementary Figure \ref{fig:entropyCalcAlg}.  The first step is to defined a spatial region (a set of $M$ pixels, $\mathbf{S}$) to which a distribution will be fitted. It must be decided whether or not the pixel values $\mathbf{S}$ are assumed to be statistically independent.  In the most general case, the $M$ pixels of $\mathbf{S}$, will have a joint distribution density function $p(\mathbf{z}_1, ..., \mathbf{z}_m, ..., \mathbf{z}_M)$, where $\mathbf{z}_m$ is a $P$-dimensional vector containing some subset of available quantities (amplitude, phase, retardance etc.).  Under this model, all $M$ pixels together effectively represent a single sample from an $M\cdot P$-dimensional distribution.  It is not typically feasible to fit such a distribution based on only a single sample and so some assumptions must be made to make the problem tractable.

A widely used assumption is that points separated by some fixed spatial vector $(a,b)^T$ are drawn from a joint distribution, $p(\mathbf{z}(x,y), \mathbf{z}(x+a,y+b))$.  This only has $2\cdot P$ dimensions and there are now $M$ samples for fitting (one for each point in $\mathbf{S}$). This assumption is regularly used in grayscale image amplitude analysis where it is termed the \emph{gray-level co-occurrence matrx (GLCM)} \cite{Albregtsen2008}.  In this work we therefore refer to quantities generated using this assumption as being \emph{co-occurrence} metrics.

A further relaxation from the co-occurrence case is to assume that each pixel, $\mathbf{z}$ in the area $\mathbf{S}$ is drawn independently from a distribution, $p(\mathbf{z})$. This, too, offers $M$ points for fitting.  However, as discussed in the main manuscript, it is found that there is relatively little difference between these two cases for the samples measured and so the latter is typically preferred.

If the quantitites chosen for fitting are all circular (i.e. represent angles), then a \emph{generalisd multivariate von Mises distribution} is used for fitting \cite{Mardia2008}. If only $A$ is chosen a chi distribution is fitted ($A$ is modelled as the absolute value of a Gaussian) and if only $D$ is chosen, a truncated Gaussian on the range $[-1,1]$ is fitted.  If a mixture of circular and non-circular quantities is used it is difficult in general to define a multivariate distribution with appropriate marginal distributions for all quantities. One solution is to scale the domain of the non-circular quantities such that their expected probability mass lies close to 0 on the $(-\pi,\pi]$ range.  In this case, as we let $\kappa \rightarrow 0$, the von Mises distribution starts to approximate a truncated Gaussian distribution which can be fit to the non-circular quantities.
 
Fitting is performed by evaluating predetermined analytic expressions for maximum likelihood parameters of distributions based on measured quantities. Such expressions are known for Gaussian and multivariate von Mises distributions \cite{Mardia2008}.

\subsection{Computation of entropy}
\label{subsec:entropyCalc}
Fitting distributions to points in a spatial region gives estimates for the mean and variance of a quantity in that region.  Given these, the \emph{entropy} of the distribution can then be computed, which gives a measure of spatial heterogeneity that can be indicative of disordered microstructure. There are many different definitions of entropy (e.g. Shannon entropy, thermodynamic entropy), although on some level they are fundamentally analogous. The definition used here is the \emph{differential entropy}, which is a property of the fitted distributions.  This metric can be derived from the Kullback-Leibler divergence, which measures the similarity of probability distributions $P$ and $Q$ with density functions $p(x)$ and $q(x)$ (the 1-D case) respectively:

\begin{equation}
\label{eq:KLdiv}
D_{KL}(P||Q) = \int_{-\infty}^{\infty} p(x) \log \frac{p(x)}{q(x)} dx
\end{equation}

If $Q$ is a uniform distribution then $D_{KL}(P||Q)$ represents how similar to a uniform distribution $P$ is -- that is, how 'spread out' $P$ is.  This is precisely what we are seeking: a measure of heterogeneity.  Setting $Q$ to be a uniform distribution and introducing a scale factor of -1, Equation \ref{eq:KLdiv} becomes the differential entropy, $H$:

\begin{equation}
\label{eq:shannonEnt}
H = - \int_{-\infty}^{\infty} p(x) \log p(x) dx
\end{equation}

This can be extended to multivariate distributions simply by integrating over the additional variables.  However, computation of this integral can be computationally expensive.  Fortunately, for many distributions analytical expressions for the entropy in terms of distribution parameters are known. This is true of the multivariate Gaussian and univariate von Mises distributions \cite{Mardia2008}, shown in Supplementary Figure \ref{fig:entropyCalcAlg}.  Where analytical expressions are not known (e.g. multivariate von Mises) Monte-Carlo sampling methods are used for efficient integral evaluation \cite{MacKay2003}

By computing entropy for a sliding window over the entire raw image, we produce an \emph{entropy filtered} image of the desired quantity. Entropy filtering has been successfully applied in the field of amplitude-only texture analysis for classification of diseased tissues \cite{Yogesan1996,DeArruda2013,Ganeshan2007}.

\subsection{Comparison of entropy and expectation}
If instead of plotting entropy as described in the previous section, we plot the mean of the fitted distributions at each spatial point, this gives images of the spatial expectation, $E$, effectively an averaging filtering on the same length scale as entropy. As shown in Supplementary Figure \ref{fig:CNRall} it was found experimentally that phase entropy is a good indicator of early tumourigenesis, whereas for polarisation properties, expectation offers better performance.

In the case of retardance, $\eta$, this could be because diseased tissue has a higher \emph{average} retardance due to, for example, increasingly dense collagen networks associated with early malignant transformation. The spatial expectation would then represent this change better than entropy.  Further, the anisotropic nature of polarising molecules such as collagen may mean they are preferentially aligned along some axis.  Even if this alignment is weak, we might then expect to see an increase in the average value of the retardance axis orientation, $\theta_{\eta}$ or diattenuation axis orientation, $\theta_D$.  This is consistent with the results of Supplementary Figure \ref{fig:CNRall} that show more contrast in the spatial expectation of $\theta_{\eta}$ and $\theta_{\eta}$ quantities than in the entropy.



\pagebreak
\noindent
\Large\textbf{Supplementary Note 5:}\\
\Large\textbf{Limitations of current design and routes to real-time application}\\
\normalsize
\setcounter{section}{5}
\setcounter{subsection}{0}


\subsection{Experimental imaging acquisition speed}
Currently, the acquisition time for an amplitude, phase and polarisation image is 8.3 s. This is because 11 images are required: 7 for establishing the phase of the horizontally polarised image and 4 for determining the phase difference between the horizontally and vertically polarised images. To infer polarimetric properties of a tissue sample, 3 such images are recorded (see Online Methods and Supplementary Note \ref{subsec:bayesInf}) taking a total of 24.9 s.  For real time operation, this would need to be reduced substantially to $<$0.2 s.  

The following modifications could be made to progress the operation towards real-time imaging. Using a camera with an automated high dynamic range imaging procedure, region of interest readout, and consequently higher frame rate could offer a >15 -fold increase in speed. Next, using an appropriately designed phase-mask in place of the current Fresnel lens, it would be possible to use only 2 images for phase retrieval and 3 for polarimetry, offering a further 2-fold speed up \cite{Horisaki2015,Zhang2007}.  Using only 2 illumination states instead of 3 offers a further reduction. Together these improvements would bring sample imaging time down to below the 0.2 s target. Further speed increases could be achieved through parallelisation of phase/polarisation retrieval by subdividing the camera sensor area, offering a trade-off with resolution \cite{Millerd2004}. At this point the switching speed of the liquid crystals in the SLM (60Hz frame rate) would become the limiting factor but this can be addressed by using high-speed (20kHz frame rate) digital micromirror devices (DMDs) \cite{Dremeau2015, Mitchell2016}.

\subsubsection{Transmission matrix recovery speed}
Although the fibre transmission matrix is typically stable for many months on an optical bench, in reality perturbations such as fibre bending would require periodic re-measurement of the transmission matrix. Our parallel spot characterisation system (described in Online Methods) already offers a speed-up of 12 over conventional spot or plane-wave scanning systems for recording fibre characterisation measurements. As opposed to methods used for MMFs, this approach scales to higher resolutions without requiring additional experimental measurements.

Currently, calibration measurements take 50.8 minutes but by making use of the optimisations discussed in the previous section this could be reduced to around 22 seconds. Re-characterisation may only need to be performed periodically at longer intervals, meaning it should not drastically affect total imaging time.  Further, parallelised compressed imaging techniques could be exploited with an array of high-speed photodiodes to move towards real-time characterisation if required \cite{Macfaden2016}.

\subsection{Image reconstruction speed}
To recover phase and polarisation information from the raw intensity images produced by the camera, a phase retrieval algorithm is used (see Online Methods). This currently takes around 67 s per single image (67 $\times$ 3 = 201 for full polarimetric imaging). The iterative algorithm used here (based on the Gerchberg-Saxton algorithm) requires around 200 iterations to converge. If we were to instead directly solve transport of intensity equations (TIE) this could be reduced to $<$0.1 s per single image \cite{Allen2001}. 

For transmission matrix recovery, it is necessary to solve 293 phase retrieval problems (81$\times$3 characterisation $+$ 50 phase stabilisation), which currently takes 5.5 hours. Again by using a TIE approach this could be reduced to 29.5 s.  Because these problems are highly parallelisable, this could be sped up another $\sim$100-fold by using state-of-the-art processors (e.g. Intel Xeon server with 8 CPUs with 24 cores each).  Further parallelisation could be achieved by using dedicated graphics processing units (GPUs). The actual recovery of the matrix from these images requires solving 81$\times$12$\times$2=1944 inverse problems -- one per row of the desired $\mathbf{A}^{-1}$ (Supplementary Note \ref{subsec:tMatRecovery}).  This currently takes 9.5 hours.  However, by downsampling the camera images (e.g. sampling one pixel for each fibrelet) and increasing the degree of parallelisaiton from the current factor of 8 to $>$100, it should be possible to reduce this to $\sim$30 s. This could be further reduced by adopting solution methods that directly reconstruct input fields as weighted sums of the calibration fields rather than inverting the full transmission matrix \cite{Choi2014}.

Bayesian inference of sample polarimetric properties (described in Supplementary Note \ref{subsec:bayesInf}) is performed using the STAN package maximum likelihood optimiser, which requires 1.6 seconds to perform inference on a single pixel.  A full endoscope frame then takes 24 minutes to reconstruct.  Again, this problem is highly parallelisable and so could benefit from a 100-fold increase in speed on a state-of-the-art CPU. Further, by deriving analytical expressions for the joint maximum likelihood parameters of neighbouring pixels and using these instead of optimisation, inference time could be reduced to $<$0.2 s per image.

Computing the spatial entropy and expectation of an image, which involves fitting distributions as described in Supplementary Note \ref{subsec:entropyCalc}, currently takes 5 s per endoscope frame. Again, by utilising greater parallelisation this could be reduced to $<$0.1 s.
 
Finally, there are solution methods that do not require an explicit phase-retrieval step and can retrieve the desired transfer matrix based purely on intensity-only measurements. For the case where the matrix is known to comprise independent, complex Gaussian distributed entries a variational inference approach can be used \cite{Dremeau2015}. For the more generalised case in which the matrix is sparse, as it is here, algorithms such as \emph{generalised approximate message passing (GAMP)} may be applied directly to amplitude image sets \cite{Rangan2010, Schniter2015}. These could potentially reduce reconstruction time, as well as requiring minimal experimental measurements.

\subsection{Reflection mode operation}
For \emph{in vivo} clinical use, the endoscope must operate entirely in reflection mode rather than its current transmission mode.  In principle, the transmission mode illumination of samples can be simulated in reflection by using oblique back-scattering illumination \cite{Ford2012}. However, the imaging fibre will invariably experience random bending and temperature variations when in use and so the transmission matrix will need to be re-characterised without access to the distal facet.


One approach to reflection mode operation is to attach a calibration beacon to the distal facet, which has been demonstrated to work for creating a single focussed spot in confocal imaging but does not in the same form support wide-field imaging \cite{Farahi2013}. Speckle correlation imaging does not require any calibration of fibres, but relies on the assumption that images are of real, positive targets with finite support, which is not true in phase imaging where objects are generally complex \cite{Porat2016}. Another approach is to attach a known reflective plate to the distal facet and use back reflected light to make dynamic adjustments to a pre-recorded transmission matrix \cite{Gu2015}, however, this has not yet been demonstrated experimentally. Initial work using highly-spaced MCF has shown promise for two-photon imaging \cite{Warren2016}. Experimental investigation of these approaches is underway in our laboratories as a next step towards clinical translation.




\input{supRefs.bbl}
\bibliographystyle{naturemag}
\end{document}