Introduction

Many characteristics of biological tissues are reflected in their optical properties. Differences in birefringence, i.e. the speed of light through that medium depending on polarisation, may change in healthy and diseased tissues. Mueller polarimetric imaging (MPI) non-invasively measures these optical properties, providing micro-structural features of a sample without contrast media [1,2,3]. Intensity images of the superficial layer are acquired by shining light at different polarisation states. The back-scattered light is then captured by an optical sensor operating in reflection [4] configuration. As per the polarimetric Stokes–Mueller formalism [5], the tissue-specific Mueller coefficients are derived from the acquired intensities by solving a linear system. Polarimetric parameters, including retardance (which characterises the anisotropy of the refractive index of a sample), diattenuation and depolarisation, are determined via different decompositions of the Mueller matrix [6,7,8]. In diagnostic clinical applications, MPI identified disease progression by revealing morphological tissue changes ex vivo [9, 10]. In [11, 12], polarised light first estimated the neuronal fibre bundle orientation in histological sections of formalin-fixed human brain, towards a tractographic reconstruction of the white matter as in diffusion weighted MRI. In [13], a wide-field MPI system showed white matter fibre tracts on fresh and formalin-fixed samples of different specimens, paving the way towards label-free neurooncological visualisations. In the aforementioned studies, real-time performance was not initially sought, since high image quality was prioritised to accurately characterise the samples’ properties. The accurate estimation of the azimuth of the optical axis, indicative of the orientation of fibre bundles within the imaging plane, is key for neurosurgery [14]. Under the assumption that lesions alter the organised arrangement of neuronal fibres, directional cues in axonal pathways may guide the resections of neoplastic tissue in white matter. Neurosurgeons would be informed on tumour boundaries and surrounding healthy tissues, irrespective of tissue ablation and displacement, beyond available navigation systems based on preoperative image planning. High accuracy is also key to discriminate among different tissues, such as neoplastic types and grades showing different degrees of infiltration [1], and to perform tissue classification tasks leveraging artificial intelligence (AI) techniques [15]. With this view, MPI acquisition noise and latency represent two main bottlenecks for the accurate image processing, quantitative analysis, and ultimate neurosurgical feedback. The enhancement of polarimetric image contrast traditionally requires long-time, multi-shots, averaging techniques to reduce the acquisition noise of the optical system and sensor camera [16]. Furthermore, computational noise propagates from unfiltered acquisitions throughout the derivation of polarimetric parameters in the cascade of numeric Mueller decompositions [8, 17]. An optimal MPI enhancement embedded in a performance-optimised image processing pipeline is a key-enabling technology for revolutionising computer-aided neurosurgery, and currently stand as an open translational challenge.

Contribution and outline

Aiming to tackle the aforementioned challenges, in this feasibility study, we introduce an AI-based framework integrated with performance-optimised polarimetric image processing tools, to simultaneously denoise low-quality single-shot polarimetric acquisitions and to boost the performance and the estimation of relevant parameters in an end-to-end pipeline. A denoising diffusion network tailored for polarimetric intensity data is introduced in section ‘Polarimetric denoising diffusion network’, and performance-optimised polarimetric image processing tools are described in section ‘Efficient calculation of polarimetric biomarkers’. Experiments validate the proposed framework for real polarimetric images of human brain tissues in different conditions in section ‘Experiments and results’. Observations on neurosurgical MPI are discussed in section ‘Discussion and conclusions’.

Polarimetric denoising diffusion network

The denoising of a short-time, low-quality, single-shot polarimetric acquisition is performed with a deep-learning implementation of a diffusion probabilistic model [18, 20]. Diffusion probabilistic models can generate images as the composition of many small denoising steps. Our polarimetric denoising diffusion network (PDDN) is a tailored adaptation of this framework, leveraged to literally reduce the acquisition noise from polarisation state intensities with high performance. In a probabilistic diffusion process, an input image \(\textbf{x}_{0}\) is gradually corrupted with additional Gaussian noise over a series of T time-points. Such degradation process, i.e. forward diffusion, determines for each \(T-1\) adjacent time-step a pair of images: a degraded instance and less-degraded one. In the reverse diffusion, the ill-posed recovery of clean data from noisy instances is achieved for each time-step by reversing the degradation process with a neural network for conditional inference. Here, the noise level of short-time, single-shot polarimetric scans is assumed comparable to the degradation at specific time-points in a forward diffusion, so that the image restoration is parametrised with a progressively small recovery, as in the reverse diffusion. As in [18], the diffusion is modelled with a Markov chain of T time-points, where the image state at each time-point t only depends on the image state at \(t-1\). For a clean input image \(\textbf{x}_0\) at \(t=0\), the full degrading trajectory to \(\textbf{x}_T\) in Fig. 1 is formulated as the sequential product of the posterior probability \(q(\textbf{x}_t|\textbf{x}_{t-1}) =\mathcal {N}(\textbf{x}_t; \varvec{\mu }_t=\sqrt{1-\beta _t}\textbf{x}_{t-1}, \varvec{\Sigma }_t=\beta _t\textbf{J})\) for each degrading time-step, with \(\varvec{\mu }_t\) and \(\varvec{\Sigma }_t\) the tensorial mean and variance of the Gaussian noise with scalar time-dependent variance \(\beta _t\), and \(\textbf{J}\) the identity matrix. To efficiently sample any degradation state \(\textbf{x}_t\) from \(\textbf{x}_0\), a closed form is obtained by re-parametrising the scalar variance \(\beta _t\) leveraging a canonical Gaussian noise \(\varvec{\epsilon } \sim \mathcal {N}(\varvec{0},\textbf{J})\), as in [18]. In the reverse diffusion, the estimation of the distribution \(q(\textbf{x}_{t-1}|\textbf{x}_t)\) as in Fig. 1 is approximated assuming an underlying Markov chain of \(T\rightarrow \infty \) time-points, and additional Gaussian noise with a small scalar variance \(\beta _t\) at each time-point. The parametric formulation of the reverse trajectory therefore approximates the distribution as \(p_{\theta }(\textbf{x}_{t-1}|\textbf{x}_t) = \mathcal {N}(\textbf{x}_t;\varvec{\mu }_{\theta } (\textbf{x}_t,t),\varvec{\Sigma }_{\theta }(\textbf{x}_t,t))\), where a learning paradigm regresses the tensorial mean \(\varvec{\mu }_{\theta }(\textbf{x}_t,t)\) and variance \(\varvec{\Sigma }_{\theta }(\textbf{x}_t,t)\) from pairs of image states, for each time-point, optimising a loss derived from the evidence lower bound, in the form of a negative log-likelihood. As in [18], a simplified formulation of the loss accounts for the re-parametrisation of the additional noise variance, assumed identical in each tensorial dimension, with further conditioning on the input image. The joint distribution of the reverse diffusion is translated into an encoder–decoder coupling as in a U-Net [21], where denoising kernels are learned in a self-supervised fashion. As in Fig.  1, the PDDN is trained on unpaired high-quality intensities obtained from long-time averaged acquisitions. At convergence, the model denoises short-time, low-quality, single-shot polarimetric images for few terminal steps, at inference.

Fig. 1
figure 1

Denoising diffusion model. a Forward and reverse diffusion as in [18]. Degraded states \(\textbf{x}_t\) of the forward diffusion (black arrow) in the modelled Markov chain of T time-points. Inferential sampling of the reverse diffusion (blue arrow) over the parametrised distribution \(p_{\theta }(\textbf{x}_{t-1}|\textbf{x}_t)\). b Schematic diagram of the AI-based denoising polarimetric framework. PDDN builds on a time-point recursive U-Net for the reverse diffusion, as in [19]

Full size image

Network implementation details The PDDN implements a time-point recursive U-Net [21], as in [19]. Four deep layers of wide ResNet blocks, group normalisation and self-attention blocks are employed, with pooling and upsampling scheme of (1, 2, 4, 8), each with \(3 \times 3\) convolutional kernels, unitary stride and sigmoid linear units (SiLU) activation function, alternated by skip connections between the encoder–decoder branches. A total of \(T=1000\) time-points were considered. The training used an L1-loss with an Adam optimiser, learning rate \(l_{r} = 1e-4\) and 100k epochs, with batches of 32 sampled and augmented patches of data. Full model memory footprint: 250 MB. Polarimetric intensities were arranged in tensor patches of size \(128\times 128\times 16\), with the first two dimensions encompassing spatial extent and the third the fixed polarimetric measurement states. Data augmentation included random rotation, flip and cropping with mirroring padding. Supra-threshold intensity reflections were masked to avoid spurious artefacts and hallucinations. Intensities were linearly re-scaled within \([-1,1]\).

Efficient calculation of polarimetric biomarkers

Divide-and-conquer approaches are employed together with linear algebra vectorisation and parallel computing to boost the performance of Mueller matrix decomposition and the extraction of accurate polarimetric parameters. In the developed performance-optimised image processing tools, polarimetric derivations, as in Fig. 2, are formulated as a system of linear equations [5]. For a MPI system in reflection configuration, the Mueller matrix is derived from a 16-channel tensor of 2D intensities of size \(H \times W\) pixels, for all \((4 \times 4)\) polarisation states as

$$\begin{aligned} \textbf{M}=\textbf{A}^{-1}\,\textbf{I}\,\textbf{G}^{-1}, \end{aligned}$$
(1)

where \(\textbf{M}\) is the unknown full \(4 \times 4\) Mueller matrix of each pixel, \(\textbf{G}\) and \(\textbf{A}\) are the pixel-wise \(4 \times 4\) matrices of the polarisation state generator and analyser, respectively, determined at calibration, and \(\textbf{I}\) is the \(4 \times 4\) tensor of real-valued intensities measured by the camera for the considered pixel, where each component accounts for a different combination of the elicited polarisation states. The linear system in Eq. (1) can be solved in closed form for each Mueller coefficient as sum of scalar products and represents an explicit vectorisation of the solution for tensors of arbitrary dimensions. Such formulation pixel-wise solves for the Mueller matrix coefficients, where parallel computing enables high-dimensional vectorised data processing with arbitrary hardware capacity. Following [7], the Mueller matrix is decomposed as the matricial product of three optical components, i.e. the diattenuator, the retarder and the depolariser, as \(\textbf{M}=\textbf{M}_{\Delta }\,\textbf{M}_{R}\,\textbf{M}_{D}\). Scalar maps of polarimetric parameters are pixel-wise derived from decomposed polarimetric tensors, accounting for total diattenuation (D), total depolarisation (\(\Delta \)), scalar retardance (R), and the azimuth of optical axis (\(\varphi \)) as

$$\begin{aligned} D= & {} \sqrt{M_{D_{12}}^{2} + M_{D_{13}}^{2} + M_{D_{14}}^{2}} ~~ \text {and} ~~ \Delta = 1-\frac{1}{3}\vert tr(\textbf{M}_{\Delta })\vert , \nonumber \\ R= & {} \cos ^{-1} \left( \sqrt{ (M_{R_{22}}+M_{R_{33}})^2+(M_{R_{32}}-M_{R_{23}})^2}-1\right) \nonumber \\{} & {} \quad \text {and }\varphi = \frac{1}{2}\tan ^{{-}1} \left( \frac{M_{R_{24}}}{M_{R_{43}}}\right) . \end{aligned}$$
(2)

All vectorisations are implemented with compiled routines, and all derivations are wrapped in scripting languages for high-level development of AI designs [22].

Fig. 2
figure 2

Wide-field MPI system in reflection configuration. a Schematic: light source, polarisation state generator \(\textbf{G}\), biological sample, polarisation state analyser \(\textbf{A}\), and detecting camera. b Our instrumentation. c Acquired polarisation states intensity image \(\textbf{I}\). d Derived full Mueller matrix \(\textbf{M}\). e Lu-Chipman decomposition: \(\textbf{M}_{\Delta }\), \(\textbf{M}_{R}\), and \(\textbf{M}_{D}\) matrices. f Derived scalar parameters: D, \(\Delta \), R, and \(\varphi \). Enhanced contrast of \(\textbf{M}\) with a sigmoid mapping

Full size image

Experiments and results

Data Polarimetric data \(\textbf{I}\) were acquired with a wide-field imaging Mueller polarimeter as in [16] Fig.  2, at 550 nm wavelength, with a CCD camera (Stingray F080B, Allied Vision, Germany) \(512\times 384\) pixels (\(20\times 24\) mm FoV, resulting in \(\approx 50\ \upmu \)m resolution). Matrices \(\textbf{G}\) and \(\textbf{A}\) in Eq. (1) were determined at calibration.

Training PDDN was trained on 200 high-quality (HQ) images from multi-shots averaged (\(n=8\)) acquisitions of fresh human brain tissues from neurosurgical resections and post-mortem examinations. Portions of grey and white matter were resected from cortical regions involving eloquent areas of the brain. These included tumour-free and neoplastic samples (\(\approx 50\%\) ratio) of different types (gliomas, meningiomas, metastases), at varying degrees of severity and infiltration.

Validating 50 HQ images of mixed fresh and formalin-fixed brain tissues from healthy human and animal specimens exhibited similar contrast and minor biological heterogeneity (cortical and deep-brain structures) for the model optimisation.

Testing A different set of 200 rigidly co-registered paired images of mixed fresh and formalin-fixed human brain samples including tumour-free and neoplastic tissues (\(\approx 50\%\) ratio) was acquired at low quality (LQ), i.e. short time, noisy, single shot (\(n=1\)), and at HQ, respectively. Representative annotated data were acquired also at super high-quality (SHQ) multi-shot (\(n=16\)) averaged acquisitions, for a case study.

Evaluation The accuracy and the performance of the denoising framework were compared to traditional approaches, alternative methods and the state-of-the-art. The evaluated polarimetric instances comprised denoised \(\textbf{I}\) and derived \(\textbf{M}\), as well as the scalar maps D, \(\Delta \), R, and \(\varphi \). Image quality scores including the root-mean-squared error (RMSE), the normalised peak signal-to-noise ratio (nPSNR) and the structural similarity index (SSIM) were pixel-wise computed for the paired test data, as in [23]. Values and deviations of angular data (R and \(\varphi \in [0, \pi ]\)), were computed with circular statistics and reported in degrees. Scores were evaluated within a region of interest (ROI) separating tissues from background. Reflection artefacts were excluded from the analysis. Significant differences were assessed with a pairwise Wilcoxon rank sum test [24].

Denoising polarimetric intensities of human brain tissues

Low-quality intensity images of polarimetric states of the testing set were denoised with a single-pass (\(t=1\)) PDDN filtering step and with a set of traditional and AI-based methods. Traditional averaging schemes yielded HQ images as reference ground-truth. State-of-the-art polarimetric convolutional denoising networks (PCNN) [23, 25] were trained on paired intensity instances of the training set. Also, deterministic denoising algorithms were considered as baseline: the median filter (MEDF: 3-kernel) [26], the Gaussian blur (GBLR: 5-kernel) [27], and the gradient anisotropic diffusion (GRAD: 5 steps, 1 conductance) [28]. Qualitative results are shown in Figs. 3, 4 for representative tumour-free and neoplastic cases. In general, minor and subtle changes are observed in the processed intensities of polarisation states and derived Mueller matrices. Conversely, the effect of denoising was predominant on derived polarimetric parameters of clinical relevance. Successful denoising is obtained for PDDN, PCNN, and GBLR with improved rejection of acquisition noise compared to LQ acquisitions. Limited denoising is found for MEDF and GRAD, where noisy patterns remained visible, and the physical characterisation of the underlying sample remained unclear or partially altered. Overall, polarimetric parameters showed high sensitivity to acquisition noise, to the computational error propagation, and to the denoising method of choice, where D, R, and \(\varphi \) exhibited major deviations between noisy and processed instances. Image quality scores in Eq. (1) supported the qualitative analysis, where the proposed PDDN reported best values in all cases for all indices, followed by PCNN and GBLR. This suggests the early and optimal rejection of subtle acquisition noise with PDDN is effective to reduce further error propagation in the computational cascade of polarimetric parameters. Major deviations were found against LQ instances, where oriented patterns of white matter fibres in \(\varphi \) can only be clearly observed after denoising. The significant improvements obtained with PDDN in Eq. (1) suggest the learned filtering kernels in the proposed polarimetric denoising diffusion network are suitable for enhancing the image quality with minimal deviations compared to reference HQ data using a fast, single-pass, filtering step (Table 1).

Fig. 3
figure 3

Denoising polarimetric intensities of human brain tissues: gallery of polarimetric instances. (top) Tumour-free sample of the testing set. \(\text {I}_{11}\) component shown with the evaluation ROI contour. (bottom) Details of polarimetric parameters in a centre-cropped area. Images in each row have the same range of values as in the colour-bar. R and \(\varphi \) reported in degrees

Full size image
Table 1 Image quality scores: human brain tissues
Full size table
Fig. 4
figure 4

Denoising polarimetric intensities of human brain tissues: gallery of polarimetric instances. (top) Neoplastic lesion of the testing set. \(\text {I}_{11}\) component shown with the evaluation ROI contour. (bottom) Details of polarimetric parameters in a centre-cropped area. Images in each row have the same range of values as in the colour-bar. R and \(\varphi \) reported in degrees

Full size image
Fig. 5
figure 5

Azimuth \(\varphi \) variations and distributions: intensity and circular standard deviation csd(\(\varphi \)) in degrees. Low values of csd(\(\varphi \)) for homogeneous directional patterns, whereas high values for high disruption of fibres orientation or change of directional patterns. (top) Tumour-free sample: variability of fibres orientations, csd distributions in annotated white (WM) and grey (GM) matter. (bottom) Neoplastic lesion: variability of fibres orientations in diseased white matter, csd distributions in annotated tumour centre (TC) and infiltration (TI). High rejection of background noise, and visually comparable directional patterns of the fibres after denoising, similarly to high-quality image standards. Boxplots: consistent csd drop in PDDN as in HQ and SHQ. Better PDDN separation in tumour areas

Full size image

Fibre orientation and Azimuth variation

Assuming lesions alter the directional arrangement of fibres in white matter, we focus on PDDN and evaluate the effect of denoising on the derived \(\varphi \) as a sanity check in tumour-free brain tissue and in a glioma lesion. The directional variability is evaluated with the azimuth circular standard deviation, i.e. csd(\(\varphi \)) in a 5x5 image pixel neighbourhood. Low csd(\(\varphi \)) indicates homogeneous directional patterns, whereas high csd(\(\varphi \)) corresponds to degrees of directional disruption. A histological section was annotated by a neuropathologist, delineating the tumour centre and the infiltration area, as well as, grey and white matter in the tumour-free sample. In Fig. 5, the csd(\(\varphi \)) in a fast and single-pass denoising with PDDN is compared against the LQ, HQ and SHQ instances as reference. The denoising reduces the csd(\(\varphi \)) in both tumour-free and diseased samples similarly to HQ and SHQ data. LQ instances show a high level of angular variability due to the intrinsic polarimetric acquisition noise and computational error propagation. More homogeneous directional patterns are found in tumour-free white matter, with lower csd(\(\varphi \)) compared to grey matter, suggesting more organised fibre tracts, while crossing fibres increase the angular deviation. The denoised lesion shows a clear difference between tumour centre and infiltration area, similarly to HQ and SHQ data. A substantially higher csd(\(\varphi \)) is found for the tumour centre compared to the infiltration area, where higher variability suggests a higher degree of disruption of axonal fibres. Conversely, the LQ instance showed limited separation of the tumour core from the background and the neighbouring structures. The qualitative observations are reflected in the box-plot distributions, in line with the underlying tissue classes, even for higher polarimetric image quality standards.

End-to-end computational performance

The end-to-end polarimetric processing pipeline accounts for single-shot denoising and parameters derivation. Traditional multi-shot averaging techniques and derivation algorithms in [16] are considered as reference for HQ MPI. The performance is reported in Eq. (2) for a local patch and a full-scale image. The comparison showed a substantial reduction in total processing time. Our end-to-end pipeline achieved real-time performance (\(<40\) ms) for a tensor patch (size: \(128\times 128\times 16\)) of polarimetric intensities. This suggests that translation to in vivo, real-time MPI can already be achieved by focusing on a smaller field of view, with further optimisation necessary for full-frame images. All computations were performed on a Linux Ubuntu 20.04 laptop, \(16\times \) CPU at 2.6 GHz, 64 GB RAM, NVIDIA RTX A5000 GPU (Table 2).

Table 2 End-to-end Computational Performance: pipeline processing time on same hardware
Full size table
Fig. 6
figure 6

Denoised polarimetric tractography: tumour-free brain sample. LQ and denoised azimuth with recursive PDDN\(^{+}\). MPI PTF: parameters mapped into an ellipsoidal model for tractography. Colours code for trace of ellipsoids eigenvalues. Fibres in seeding regions as in neurosurgical probing

Full size image

Discussion and conclusions

In this feasibility study, we introduced a novel polarimetric denoising framework, with the goal of enabling high quality, high-performance MPI for neurosurgery. Developments combined our PDDN, for accurately enhancing images from short-time low-quality acquisitions, with a performance-optimised toolkit to efficiently derive parameters of clinical relevance. The validation reported significantly improved image quality and achieved real-time performance for a local field of view. The denoising accuracy was tested on multiple and diverse instances of human brain samples for different image restoration methods. Our self-supervised PDDN yielded best rejection of the acquisition noise and limited the error propagation in the computational cascade, with comparable values to reference HQ data. While multi-shots averaging [16, 29] produces reference MPI, it is incompatible with in vivo neurosurgery, where real-time feedback is needed. Bypassing time-consuming MPI with enhanced image processing was first proposed in [23, 25], where U-Net-like architectures (PCNN) denoised Mueller matrices derived from noisy, short-time acquisitions. In [23], Mueller coefficients were denoised after training on large, paired, histological data, with inferential performance not yet compatible with real-time applications. In our experiments, PDDN over-performed PCNN and traditional denoising methods. This is likely due to a combination of factors: the different nature of input polarimetric data, the type of noise, and the underlying probabilistic model. While our model is specific to the considered human brain samples, and for the specific polarimetric acquisition conditions, differently from [23, 25], the PDDN denoises source intensities corrupted by acquisition noise, with the Mueller coefficients being subsequently derived. The rationale behind adopting PDDN builds on empirical similarities between measured MPI acquisition noise, i.e. pseudo-Gaussian: symmetric, zero-mean, bell-shaped, with cumulative slightly deviating from the Normal reference, and the additive noise in the probabilistic formulation. As denoising diffusion networks can generalise for complex distributions [18], we aimed to reduce MPI acquisition noise by generalising for contrast variability in biological structures with few filtering steps. The initial calibration mitigated systematic errors in the polarisation states; however, different wavelengths and varying exposure time may introduce a nonlinear intensity bias together with other specific human brain tissue structures (e.g. cortical regions in eloquent areas vs. deep-brain structures of corpus callosum, or other structures of the cerebellum), which may potentially alter the image contrast, structural patterns, values and noise distributions propagated in the Mueller derivations. In this case, our specific model was able to generalise for the considered polarisation states, for the intensity bias, for fresh and formalin-fixed samples, and for tumour-free and neoplastic tissues, by preserving the underlying micro-structure after denoising. Clear cortical white matter fibres orientations and comparable azimuth deviations were observed after denoising with respect to HQ and SHQ data for a representative tumour-free sample and a glioma lesion in Fig. 5. Angular deviations were visible after denoising, in keeping with underlying tissue: the consistent drop in azimuth variability showed higher compression and reduced overlap among pathological regions, better than HQ and SHQ data. Interestingly, PDDN was only trained on HQ images, yet azimuth deviations were similar to SHQ data, suggesting high MPI quality is achievable with AI beyond conventional acquisition paradigms. Prospectively, advanced configurations (PDDN\(^{+}\)) may enable neurosurgical fibre tracking with polarimetric tensor fields (PTF) in Fig. 6. Future analyses will test multi-spectral denoising in vivo, accounting for motion and bleeding artefacts. MPI instrumentation optimisations and image processing developments will be tailored on edge-computing solutions, for real-time wide-field MPI video streams.

Supplementary materials Representative polarimetric brain data employed in the validating dataset are available at: https://osf.io/9ynmf/.