Optimized Signal Separation for 3D-Polarized Light Imaging

1. Introduction

1.1. Anatomical connectivity mapping

“Why should we bother about connectivity in the age of functional imaging, at a time when magnets of ever increasing strength promise to detect the location of even the faintest thought? Isn’t it enough to locate cortical areas engaged in deception, introspection, empathy? Do we really have to worry about their connections? The answer is “yes”. In the case of the nervous system, the unit of relational architecture that allows the whole to exceed the sum of the parts is known as large-scale network. Its elucidation requires an elaborate understanding of connectivity patterns” [1]. Despite considerable advances in experimental techniques and in our understanding of animal anatomy over the last decades, the real connectivity of the human brain has essentially remained a mystery. It is the human brain’s multiscale topology that poses a particular challenge to any neuroimaging technique and prevented the neuroscientists from unraveling the connectome so far.

However, it is also the brain’s architecture that allows different morphological entities to be defined at different scales depending on the spatial resolution provided by the available neuroimaging techniques and the scientific objectives. Consequently, a comprehensive description of neuronal networks and their intricate fiber connections requires a multimodal approach based on complementary imaging techniques targeting different levels of organization (microscale, mesoscale, and macroscale) [2,3].

MR-based diffusion imaging is the most frequently used method to visualize fiber pathways in both the living and the postmortem human brain (for a comprehensive introduction to the field cf. [4,5]). Diffusion imaging contributes to the understanding of the macroscopic connectivity (i.e., at the millimeter scale) in the living brain, while postmortem studies already explore the upper mesoscale (i.e., the sub-millimeter scale). Hence, not surprisingly, diffusion imaging is of great appeal to neuroscientists as a method for the visualization of connectivity patterns in both clinical and basic research. However, complex fiber networks and small fiber tracts are difficult to be disentangled reliably at present. Furthermore, the termination of fields of pathways emanating from cortical areas no larger than a few millimeters in size cannot be demonstrated with the required precision.

Conversely, microscopic techniques generate data sets of impressing neuroanatomical detail, but they are limited to small sample sizes (i.e., small areas of interest in a small number of subjects) of postmortem tissue. This substantially restricts their predictive power. In the recent years, anatomical connections in the human postmortem brains were studied with dissection techniques [6,7], in myelin stained sections of adult human brains [8], or of immature brains taking advantage of heterochronic myelination of different fiber tracts during pre- and early postnatal development [9], in lesioned brains using various techniques for staining degenerating fibers [10,11], and using tract-tracing methods for discovering local connections [12,13]. These studies have enriched our knowledge about human brain fiber tracts, but all of them suffer from severe restrictions if the 3D courses of fiber pathways are to be mapped in the adult human brain. In contrast to studies in animals, the tight packing of different fiber tracts in the white substance, and the lack of specific tracers for in vitro tracking of long distance fibers made comprehensive fiber tract mapping impossible in the adult human brain [14].

1.2. 3D-polarized light imaging (3D-PLI)

Recently, a neuroimaging technique referred to as 3D-polarized light imaging (3D-PLI) has been introduced that has opened up new avenues to study human brain regions with complex fiber architecture as well as small cortical fibers at the micrometer level [15–17]. This technique is applicable to microtome sections of postmortem brains and utilizes the optically birefringence of nerve fibers, which is induced by the optical anisotropy of the myelin sheaths surrounding axons [18–20]. 3D-PLI provides a three-dimensional description of the anatomical connectivity in form of a vector field indicating the prevailing fiber orientation per voxel. Depending on the used imaging setup and the section thickness the acquirable spatial resolution is 1.6 – 100 µm. Hence, the method bridges the microscopic and the macroscopic description of the anatomical connectivity of the human brain.

The birefringence of brain tissue is measured by passing linearly polarized light through histological brain sections and by detecting local changes in the polarization state of light using a polarimeter setup (Figure 1a-b). The polarimeter is equipped with a pair of crossed polarizers, a tilting specimen stage, a quarter-wave retarder, an LED light source (with a narrow-band green wavelength spectrum), and a charge-coupled device (CCD) camera. By rotating the optical devices simultaneously around the stationary brain section and by imaging the sample at discrete rotation angles ρ, a sinusoidal variation of the measured light intensity (i.e., the light intensity profile) is observed for each image pixel (Figure 1c). The individual course of a light intensity profile essentially depends on the locally prevailing 3D fiber orientation described by an orientation unit vector which is defined by the in-section direction angle φ and the out-of-section inclination angle α. Basic principles of optics (Snell’s law and Huygens-Fresnel principle) and the Jones calculus [21] mathematically link the measured light intensity profile I to the fiber orientation via

$$I=\frac{I_0}{2}\bullet \left[1+sin\left(2\rho -2\phi \right)\bullet sin\delta \right],$$E1

with $\delta \approx 2\pi \bullet \frac{d\bullet \mathrm{\Delta }n}{\lambda }\bullet {cos}^2\alpha$.

The amplitude of the profile quantifies the phase retardation δ induced to the light wave by the myelin. This phase retardation is a function of the light wavelength λ, the section thickness d, the birefringence Δn of the myelin, and the inclination angle α. The transmittance I₀ denotes the intensity of the incident light modified by local extinction effects. While ρ is the azimuth angle of the transmission axis of the first polarizer, $\phi$ represents the projection of the fiber axis onto the section plane with respect to the starting angle of the polarimeter. As a consequence, the fundamental data structure gained by 3D-PLI is a 3D vector field description of fibers and fiber tract orientations – the basis for subsequent tractography and fiber modeling.

For estimating the sinusoidal profile of the 3D-PLI signal Discrete Fourier Analysis (DFA) can be used to deduce I from Equation (1):

$$I=a_0+a_1\bullet \sin \left(2\rho \right)+ b_1\bullet \mathrm{c}\mathrm{o}\mathrm{s}\left(2\rho \right)$$E2

with $a_0=\frac{I_0}{2}\mathrm{ },\mathrm{ }{a}_1=\frac{I_0}{2}\mathrm{ }\bullet \sin \left(\delta \right)\bullet \cos \left(2\phi \right),\mathrm{ }{b}_1=\frac{I_0}{2}\mathrm{ }\bullet \sin \left(\delta \right)\bullet \sin \left(2\phi \right).$

1.3. Relevance and restoration of the light intensity profile

In 3D-PLI, the light intensity profile (cf. Figure 1) represents a crucial data set in terms of fiber orientation determination, since peak position and signal amplitude of the profile are directly related to the in-section direction and out-of-section inclination of the fiber orientation, respectively. A precise and undisturbed recording of the light intensity passing through a microtome section is therefore mandatory for the reliable reconstruction of high-resolution fiber tracts. The signal quality is however influenced by several conditions. Thermal effects and electrical noise in both the light source and the CCD-electronics deteriorate the PLI signal characteristics at each image pixel. Filter inhomogeneities of the polarizers or retarder may also manipulate the intensity profile. Therefore, a standardized image calibration technique using flat fields is usually applied to all raw images prior to analysis to compensate pixel-wise for inhomogeneities across the field-of-view [22]. Depending on the section thickness the pixel-by-pixel intensity is also influenced by the scatter properties of the investigated object. Another possible source of artifacts is dust on the polarizer. Although the polarimeter should be operated in a shielded construction to prevent for external light and dust particles, dust cannot completely be avoided. As a result of the rotation of the polarizer dust particles will deteriorate the measured light intensity only once within a half circle, if and only if they are located on the rotating system (cf. Figure 1). Hence, the measured intensity at the CCD array is a linear mixture of different light sources.

Recently, signal enhancement and restoration techniques for PLI images utilizing independent component analysis (ICA) have been introduced [22,23]. As a result, component maps corresponding to gray and white matter structures as well as noise and artifacts can be identified automatically using statistical analysis tools [23]. Remarkably, even in the presence of dust on the polarimeter ICA can effectively restore the original sinusoidal signal (Figure 2). After ICA filtering the noise and artifacts are removed and the sinusoidal nature of the PLI signal is restored (Figure 2).

Figure 1.

Polarimetry at a glance. (a) Scheme of the rotating polarimeter with tilting stage (N-North, W-West, E-East, S-South). (b) Scheme of the optical fiber model. The refractive index of a negative uniaxially birefringent medium, such as a myelinated axon, is described by an elliptically shaped oblate surface, the indicatrix (gray mesh). A beam of linearly polarized light (blue trace) interacts locally with the myelin sheath of a single axon (black line) and becomes elliptically polarized. (c) A typical 3D-PLI raw image data set consists of 18 images corresponding to equidistant rotation angles between 0° and 170°. Here, a selection of four images of a coronal section is shown, while the sketched arrow indicates one representative pixel. To obtain the fiber orientation, the measured light intensities are studied pixel-wise as a function of discrete rotation angles. The derived physical model provides a precise mathematical description of the measurement (continuous black line) and relates the sine phase to the direction angle φ and the amplitude to the inclination angle α. The highlighted data points correspond to the selected images.

Figure 2.

In (a) multiple artifacts including dust contamination are present. The trajectory of a dust particle is highlighted by the signal power inside the white circles along all rotation angles. After ICA filtering (b) the signal derogation is greatly removed. At the bottom, PLI raw signals (black) are shown together with corresponding signals after artifact rejection (red). The signals were extracted from the two locations as indicated by the red dots. In comparison to the PLI raw data ICA obviously removed the signal derogation and is capable of restoring the original sinusoidal signal in case of dust (signal #1) or noise (signal #2) contaminations.

2. A new concept for optimal signal decomposition in polarized light imaging

Although blind source separation methods are successfully applied for signal separation in all kinds of neuroimaging modalities [24–29] two major problems remain: i) the method applied must be selected carefully and the separation strategy (i.e., the internal cost function) of the applied method should be optimal for the type of data. This however, is one of the reasons why many different ICA algorithms co-exist. ii) Assuming the measured signals are adequately separated into the underlying source signals (i.e., signal of interest and non-interest), identification of the signal of interest should be performed user independently; preferably in an automatic fashion. For the latter, significant effort has to be invested in order to identify components of interest automatically from data recorded utilizing different neuroimaging techniques [29–33]. In contrast to many other ICA applications, the great advantage in PLI signal decomposition is that all profiles of the basis vectors that correspond to the signal of interest must show a sinusoidal waveform [16,22]. Since these waveforms can only vary in amplitude and phase, an automatic extraction of the signal components can be achieved utilizing this property as demonstrated in [23].

Alternatively, it has been demonstrated that by incorporating prior knowledge, i.e., by imposing temporal or spatial constraints for the source separation task, decomposition can be effectively improved [34–37]. Hesse and James [38], for example, showed that using different types of spatial constraints ICA can be trained to i) identify ocular artifacts automatically and ii) to detect and trace ictal activity. Liu and colleagues recently presented a reference based ICA concept to successfully target a specific genetic variation [36] in real and simulated data sets. For the investigation of signal detection in the imaging system of the retina in cats Barriga and colleagues introduced a constrained ICA algorithm which increases the detection of responses to visual stimuli in cats even for low levels of signal-to-noise ratio [37]. In the next section the basic principle of ICA is described shortly. For a more detailed review, we refer to [39,40].

2.1. A short introduction to independent component analysis (ICA)

For 3D-PLI a linear superposition of light at the CCD camera is assumed, where each elementary signal component refers to a distinct region in space. By applying independent component analysis (ICA) to a set of polarized light images (here a stack of 18 images at different rotation angles is used) the decomposition of the data results in spatially independent components (often called feature or basis vectors) yielding maximally (i.e., statistically) independent spatial component maps.

Let $\mathbf{X}=({{\mathrm{X}}_1,{\mathrm{X}}_2,\dots ,{\mathrm{X}}_{\mathrm{N}})}^T$ be the N-dimensional measured PLI signal mixture, where each ${\mathrm{X}}_i$ is one image mixture (flatted to a one-dimensional image vector with M pixels) detected at a corresponding rotation angle ρ. In spatial ICA $\mathbf{X}$ is considered to be a N × M matrix (as opposed to temporal ICA, where the dimension of the matrix is transposed), with M = number of pixels included in the analysis and N = number of rotation angles reflecting N different instances of the signal. The contribution of each source image varies N times over the angles ρ. Similarly, $\mathbf{S}=({{\mathrm{S}}_1,{\mathrm{S}}_2,\dots ,{\mathrm{S}}_N)}^T$ represents the N-dimensional true source signals. Hence, the linear relationship of mixed sources can be expressed as follows:

$$\mathbf{X}=\mathbf{A}\mathbf{S}$$E3

with A being the unknown mixing matrix of dimension N × N. The key problem in ICA is to find an unmixing matrix W (similar to a pseudo-inverse A^-1, with W≈A^-1) while imposing that the sources in S are statistically independent, that is:

$$\mathbf{C}=\mathbf{W}\mathbf{X}.$$E4

Within ICA the N-dimensional data array $\mathbf{X}$ is transformed into an N-dimensional component space C, where each of the spatial component maps carries a minimum amount of mutual information and thus is maximally independent. These independent components (or spatial maps) are stored in the rows of C, while in the columns of ${\mathbf{W}}^{-1}$ the associated spatially independent “temporal” profiles (i.e., the contrast changing function along the rotation angle, the basis vector) are stored.

Signal restoration, i.e. the cleaning process, is performed by zeroing columns in ${\mathbf{W}}^{-1}$ which reflect signal contributions from unwanted (e.g. artifact) sources. This is identical to zeroing rows in C, where the independent component maps of interest are stored in C'. The cleaning of measured data is performed by back transformation of C', which results in a new set of PLI data $\mathbf{X}$'.

$${\mathbf{X}}^{\mathrm{\text{'}}}=\mathrm{ }\widehat{\mathbf{A}}{\cdot \mathbf{C}}^{\mathrm{\text{'}}}$$E5

with ${\mathbf{W}}^{-1}=\widehat{\mathbf{A}}.$

2.2. Constrained ICA for polarized light imaging (cICAP)

Signal separation in constrained ICA is based on incorporating prior knowledge about the underlying signals. In the study of Barriga and colleagues [37], a priori information is incorporated into one column of the estimated mixing matrix of the spatial decomposition via an off-on constraint, meaning that the visual response signal used in the study has to follow the stimulus signal. For this, the authors filled one column of the mixing matrix with ones for the time of the stimulation, while zero values reflect off-stimulation. Similar to the concept introduced by Barriga and colleagues in 2011, it is possible to incorporate the a priori information from 3D-PLI, but with the major difference that for birefringent signals an analytical expression exists to model the expected signal of interest, which only differs in amplitude and phase. Therefore, a new algorithm motivated by the results described in [23] has been developed, where signal separation and the identification of the underlying 3D-PLI source signals are combined into one method. The new algorithm (referred to as constrained ICA for polarized light imaging, cICAP) is based on a modified version of the Infomax principle [41] by means of incorporating prior information to the cost function of Infomax; this is similar to the work reported in [34,37]. In cICAP a theoretically expected function $f_i\left(u_{\rho }\right)$ is implemented into the Infomax cost function to control for the results of the decomposition. Recall that if the i-th column in $\widehat{\mathbf{A}}$ has a sinusoidal profile, the corresponding component map in C represents a signal of interest (cf. Equation (4-5)). The expected function $f_i\left(u_{\rho }\right)$ can be derived from Equation (2), where the parameters a₀, a₁ and b₁ are fitted involving all rotation angles ρ [16].

Starting from Equation (4) the natural-gradient version in Infomax is used to update the weight matrix W during the ICA learning procedure [42,43]:

$$∆\mathbf{W}=ϵ[\mathbf{I}+\left(1-2\mathrm{y}\right){{\mathbf{C}}^{\mathit{*}}}^{\mathit{T}}]\mathbf{W}$$E6

with $\mathrm{y}=1/(1+e^{-{\mathbf{C}}^{\mathit{*}}})$.

In Equation (6) I and $ϵ$ denote the identity matrix and the learning rate, respectively, as it is used in the standard Infomax algorithm [41]. C^* denotes the intermediate component matrix which is identical to C when the learning procedure is finished. In common with the approach of [37], the mixing matrix $\widehat{\mathbf{A}}$ is estimated by $\widehat{\mathbf{A}}≔(\mathbf{P}\bullet {\mathbf{W})}^{-1}$, with $\mathbf{P}$ being the sphering matrix of the decorrelation process within Infomax. By incorporating the a priori information at each iteration step, $\widehat{\mathbf{A}}$ is updated by:

$${\widehat{\mathbf{a}}}_{\mathit{i}}^{\mathbf{*}}=\left(1-\mathit{c}\right)\bullet \mathit{ }{\widehat{\mathbf{a}}}_{\mathit{i}}+\mathit{ }\mathit{c}\mathit{ }\bullet \mathit{ }{\mathit{f}}_{\mathit{i}}\left({\mathit{u}}_{\mathit{\rho }}\right)$$E7

where ${\widehat{\mathrm{a}}}_i$ refers to the i-th column of the estimated mixing matrix $\widehat{\mathbf{A}}$ and c is a confidence parameter ranging from 0 to 1. The confidence parameter c in Equation (7) indicates the percentage of weighting of the prior knowledge by means of the expected function $f_i\left(u_{\rho }\right)$. In the limit of c→0 the original Infomax weight update is used (cf. Equation (6)). It is important to note that for each iteration step only one column in $\widehat{\mathbf{A}}$ is updated as expressed in Equation (7) since we do not want to influence the decomposition of the remaining components during the unmixing process. The column in $\widehat{\mathbf{A}}$ which entries are most similar to the corresponding expectation function is selected and modified according to Equation (7). As a measure for similarity the kurtosis of the deviation function $d_i= {\widehat{a}}_i- f_i\left(u_{\rho }\right)$ is used:

$$\mathit{k}\mathit{u}\mathit{r}\mathit{t}\left({\mathit{d}}_{\mathit{i}}\right)=\mathit{ }\frac{1}{\mathit{N}}\mathit{ }\bullet \mathit{ }\sum _{\mathit{\rho }}^{\mathit{N}}{\left(\frac{{\mathit{d}}_{\mathit{i},\mathit{\rho }}-{\stackrel{-}{\mathit{d}}}_{\mathit{i}}}{{\mathit{\sigma }}_{\mathit{i}}}\right)}^4-3,$$E8

where ${\stackrel{-}{d}}_i$ is the expected value of the deviation function of the i-th component. In the ideal case d_i will be zero (or close to zero); for any deviation at one or more corresponding angles $\rho$ the measure $kurt\left(d_i\right)$ will largely be increased. For a subsequent selection of columns in $\widehat{\mathbf{A}}$ the column profile with smallest kurtosis value as expressed by Equation (8) is used first. Moreover, if the mean square error $\left(MSE\right)$ between the updated column ${\widehat{a}}_i^{*}$ and the corresponding theoretically expectation function $f_i\left(u_{\rho }\right)$ is less than a predefined tolerance value t the entries of ${\widehat{a}}_i^{*}$ are fixed. Thus, the $MSE$ is expressed by:

$$\mathit{M}\mathit{S}\mathit{E}\left({\widehat{\mathit{a}}}_{\mathit{i}}^{\mathit{*}},\mathit{ }{\mathit{f}}_{\mathit{i}}\left({\mathit{u}}_{\mathit{\rho }}\right)\right)=\mathit{ }\frac{1}{\mathit{N}}\mathit{ }\bullet \mathit{ }\sum _{\mathit{\rho }}^{\mathit{N}}{\left({\widehat{\mathit{a}}}_{\mathit{i}}^{\mathit{*}}\left(\mathit{\rho }\right)-\mathit{ }{\mathit{f}}_{\mathit{i}}\left({\mathit{u}}_{\mathit{\rho }}\right)\right)}^2.$$E9

The basic idea here is that a very small MSE (i.e., MSE ≤ t) in Equation (9) indicates that ${\widehat{a}}_i^{*}$ and $f_i\left(u_{\rho }\right)$ are very similar, which means that the $i$-th column in $\widehat{\mathbf{A}}$ then represents a source of interest. This approach is repeated for the remaining columns in $\widehat{\mathbf{A}}$ until no further components of interest are found. As a stopping criterion for the iteration, a predefined threshold, ε, for the MSE (with $\epsilon \gg t$) is used to prevent the algorithm of running into an endless loop. This differentiation of updating ${\widehat{a}}_i^{*}$ is expressed in Equation (10):

$${\widehat{\mathit{a}}}_{\mathit{i}}^{\mathit{*}}=\left\{\begin{array}{cc}\left(1-\mathit{c}\right)\bullet {\widehat{\mathit{a}}}_{\mathit{i}}+\mathit{c}\bullet {\mathit{f}}_{\mathit{i}}\left({\mathit{u}}_{\mathit{\rho }}\right)& ,\mathit{i}\mathit{f}\mathit{ }\mathit{k}\mathit{u}\mathit{r}\mathit{t}\left({\mathit{d}}_{\mathit{i}}\right)<\epsilon \wedge MSE\left({\widehat{\mathit{a}}}_{\mathit{i}},\mathit{ }{\mathit{f}}_{\mathit{i}}\left({\mathit{u}}_{\mathit{\rho }}\right)\right)>t\\ {\widehat{\mathit{a}}}_{\mathit{i}}& ,\mathit{e}\mathit{l}\mathit{s}\mathit{e}\end{array}\right.$$E10

For the weight update described in Equation (10) the choice of the confidence value c and the tolerance values t and ε should be performed in a representative but independent data set, which is described later.

2.2.1. Evaluation of the signal enhancement

For measuring the noise reduction and signal enhancement in a set of 3D-PLI data after ICA application we use the weighted reduced chi-squared statistic as introduced in [23]. The reduced chi-squared statistic χ² involves the variance σ² of the observation, where the statistical output is weighted based on the measurement error

$${\mathit{\chi }}^2=\frac{1}{\mathit{\nu }}\mathit{ }\bullet \mathit{ }\sum _{\mathit{\rho }}^{\mathit{N}}\frac{{({\mathit{I}}_{\mathit{\rho }}-\mathit{f}({\mathit{u}}_{\mathit{\rho }}\left)\right)}^2}{{\mathit{\sigma }}_{\mathit{\rho }}^2}.$$E11

Here, $I_{\rho }$ denote the intensity measured at angle ρ. The variance σ² used in Equation (11) was obtained from 100 flat field images, which were independently generated for the image calibration process [22]. The weighted sum over N angles is further normalized by the degrees of freedom υ. The normalization ensures a good fit when the reduced chi-square value equals one, which is achieved when the squared difference between the measurement and the expected function resemble the variance of the measurement. In order to measure the signal improvement (and the reduction of noise) through the ICA process the test statistic is determined before $\left({\chi }_{raw}^2\right)$ and after $\left({\chi }_{ICA}^2\right)$ ICA. Using the ratio of these values and introducing a weight factor $\omega$ that penalizes the goodness of fit, whenever a signal component is missing [23]. The weighted test statistic reads as follows:

$$\mathit{w}\mathit{r}\mathit{G}\mathit{O}\mathit{F}\left(\mathit{x},\mathit{y}\right)=\mathit{ }\frac{{\mathit{\chi }}_{\mathit{r}\mathit{a}\mathit{w}}^2}{\mathit{\omega }\mathit{ }\bullet \mathit{ }{\mathit{\chi }}_{\mathit{I}\mathit{C}\mathit{A}}^2}\mathit{ },$$E12

with $\mathit{\omega }=\mathit{ }\frac{1}{\mathit{\nu }}\mathit{ }\bullet \mathit{ }\sum _{\mathit{\rho }}^{\mathit{N}}\frac{{({\mathit{f}}_{\mathit{r}\mathit{a}\mathit{w}}\left({\mathit{u}}_{\mathit{\rho }}\right)-\mathit{ }{\mathit{f}}_{\mathit{I}\mathit{C}\mathit{A}}({\mathit{u}}_{\mathit{\rho }}\left)\right)}^2}{{\mathit{\sigma }}_{\mathit{\rho }}^2}\ge 1.$

The weight factor ω increases when the squared difference between the two expectation functions ${\mathit{f}}_{\mathit{r}\mathit{a}\mathit{w}}\left({\mathit{u}}_{\mathit{\rho }}\right)$ and ${\mathit{f}}_{\mathit{I}\mathit{C}\mathit{A}}\left({\mathit{u}}_{\mathit{\rho }}\right)$ (derived from the raw and the ICA filtered PLI data sets, respectively) is large. In case of missed components of interest, the signal strength of the ICA filtered 3D-PLI signal would be largely reduced at corresponding pixel locations. Consequently the expectation function ${\mathit{f}}_{\mathit{I}\mathit{C}\mathit{A}}\left({\mathit{u}}_{\mathit{\rho }}\right)$ would differ in terms of its amplitude, while the shape of the waveform may be still sinusoidal. The assumption here is that the signal power (across the rotation angles) of the signal of interest is larger compared to noise. The restriction of ω ≥ 1 is needed in order to not artificially improve the goodness-of-fit in cases where the two expectation functions ${\mathit{f}}_{\mathit{r}\mathit{a}\mathit{w}}\left({\mathit{u}}_{\mathit{\rho }}\right)$ and ${\mathit{f}}_{\mathit{I}\mathit{C}\mathit{A}}\left({\mathit{u}}_{\mathit{\rho }}\right)$ are very similar (e.g., at pixel locations with low noise).

2.2.2. Finding the optimal parameters

The optimal parameters c, t and ε for the weight update as described above (cf. Equation (10)), should be determined using a representative but independent data set. In the test data from the same post-mortem brain were used, but from a different microtome section which is not included in the subsequent analysis. The determination of the optimal parameters for cICAP requires the maximization of the weighted goodness-of-fit (wrGOF) values across all pixel locations (≈$3.41\bullet {10}^5\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}\mathrm{s}$). As a control measure for this requirement the minimal wrGOF value was checked for all trajectories (i.e., the intensity profiles at pixel location x,y) within the brain section. In total we performed 1225 (=35²) cICAP calculations while varying both the confidence value c and the tolerance value t accordingly to the range from 0.0 to 0.7 and ${10}^{-3}$ to ${10}^{-7}$, respectively. Figure 3a shows the results of this benchmark run. The largest improvement (cf. blue dot in Figure 3a) was found for the confidence and tolerance values of c=0.16 and $t=1.07\bullet {10}^{-6}$, respectively. Once both parameters c and t are fixed the stopping criterion ε, which controls the number of iterations, can be determined. For this task, we investigated the smallest and largest error found in the base functions of $\widehat{\mathbf{A}}$ by means of the MSE value as expressed by Equation (9). As a reference for the definition of signal components, the user independent component selection tool was used, which is based on the statistically analysis of the base functions as reported in [23]. Throughout all iterations in cICAP the maximal MSE across the columns in matrix $\widehat{\mathbf{A}}$ that were identified as signal components was found to be well below 0.01. Figure 3b shows that after about 35 iterations the MSE value converges to almost a constant for both the signal (red) and noise components (blue), where both types of components can clearly be identified by means of the corresponding MSE value.

Figure 3.

Parameter benchmark run. In (a) the optimal confidence value c and the tolerance value t was determined through a test run of 1225 signal decomposition steps varying both parameters c and t 35 times each through the range from 0.0 to 0.7 and ${10}^{-3}$ to ${10}^{-7}$, respectively. As a control measure to find the best parameters, where the goodness-of-fit statistics (wrGOF) is increased at all pixel locations, the largest minimal wrGOF value (blue circle) is determined. The best set of parameters was found for c=0.16 and $t=1.07\bullet {10}^{-6}$ with a minimal wrGOF value of about 21. (b) The development of the maximal mean squared error (MSE) is shown for 147 iterations to determine the stopping criterion ε, which controls the number of iterations. The maximal and minimal MSE value is plotted for components reflecting signal of interest (red) and noise components (blue), respectively. After a few iterations the largest MSE value for signal components was found to be well below 1%. Therefore, setting ε to 0.01 both types of components can clearly be identified.

In Figure 4 the progress of changes in the sinusoidal profile of one exemplary basis vector is shown for all 18 iterations (blue) within cICAP. The algorithm converges with a very small fit error (MSE) of $1.07\bullet {10}^{-6}$. The resulting sinusoidal trajectory (black) almost perfectly fits the theoretically expectation function (red).

The process of updating the weight matrix and how a priori knowledge is included in cICAP is shown in Figure 5. Once a basis vector resembles the expectation function under the criteria of MSE(${\widehat{a}}_i^{*}$) < t, this basis vector is not changed in further iterations. Note, throughout the optimization the selection of the underlying signal components is automatically included in cICAP and no further user interaction is required.

Figure 4.

Stepwise changes of one exemplary basis function (changing from light blue to black). After 18 iterations using cICAP the algorithm converges with a very small fit error (MSE) of $1.07\bullet {10}^{-6}$. The resulting sinusoidal trajectory (black) almost perfectly fits the theoretically expectation function (red).

Figure 5.

Scheme of the weight matrix update and optimization procedure in cICAP.

3. Signal enhancement in 3D-PLI data sets

3.1. Signal acquisition and preprocessing

The performance of cICAP was tested on 102 histological sections of one adult human brain acquired from the body donor program of the University of Rostock, Germany, in accordance with legal requirements. The brain was fixed in 4% formalin for 6 months, cryo-protected with glycerin and cut after freezing in a cryostat microtome (Polycut CM 3500, Leica, Germany). The brain was sectioned coronally with slice a thickness of 70 µm. The sections were mounted on glass slides with Aquatex© (Merck, Germany) and cover-slipped. Note that all preparation steps are conforming to preserve the integrity of the birefringent myelin sheaths.

The sections were digitized using the polarimeter setup described in section 1.2. Each brain section was imaged at 18 equidistant rotation angles of the polarimeter covering an angle range between 0° and 170° (Figure 1b). The acquired RGB-colored images have a size of 2776 × 2080 pixels with a pixel size of 64 µm × 64 µm. The intensities were sampled with a dynamic range of 14 bits per color channel. Since the light source of the polarimeter is composed of light-emitting diodes (LED) emitting a narrow-band green wavelength spectrum (central wavelength of 525 nm), only the green channel from the RGB color triplet was used for further analysis.

Before the images were decomposed by ICA, a standardized calibration technique using flat fields was applied, in order to compensate pixel-wise for inhomogeneities across the field of view [22]. In addition, the separation of brain tissue from the non-relevant image background was done by means of the interactive learning and segmentation toolkit (ilastik) [44]. Thus, the following calculations and statistical analyses were solely done on basis of brain tissue measurements.

3.2. Performance on signal decomposition using cICAP

To test the performance of the signal decomposition and signal restoration cICAP was applied to the data set described in section 3.1. The mean number of pixels used for ICA across the 102 brain sections was found to be ($3.32\pm 0.95) \bullet {10}^5$. The results of signal decomposition and enhancement gained by cICAP were compared to the automatic component selection routine, which were recently published in [23], whereas the threshold for this data set was adapted to 2.2.

After the determination of the parameters c, t and ε all values were fixed within cICAP and both ICA algorithms were applied to the same set of PLI data for comparison. As shown in Table 1, the identified number of signal components (on average) does not vary dramatically across both algorithms. However, the smallest variability with the fewest number of iterations needed was found in the results of cICAP. The number of iterations needed (on average) in cICAP was 192, while for the Infomax algorithm using the automatic selection routine it was found to be 406. More importantly, the goodness-of-fit criterion wrGOF was found to be largest using cICAP. It is important to note that this criterion expresses the signal enhancement (or degradation) after ICA filtering and was calculated for all trajectories (i.e., at each pixel location). If the wrGOF value is larger than 1, then noise and artifacts were removed (or suppressed) successfully. In cases, where the signal decomposition fails (or is not optimal) the wrGOF value becomes even smaller than 1 (Figure 6). To check the decomposition results for signal degradation after the application of ICA, the minimal wrGOF value was calculated. As illustrated in Table 1 a wrGOF value smaller than 1 was found in 3.57% and 0.35% of all pixels in the results from Infomax and cICAP, respectively.

In general both ICA algorithms produced remarkable decomposition results, where signal enhancement (wrGOF ≥ 1) is evident in almost all pixels. However, we found a much better signal decomposition performance in the results of cICAP. This is expressed by the fact that the wrGOF value was found to be equal or larger than 10 in about 80% of all trajectories, while for the Infomax algorithm it was 39%. For cICAP we still found a goodness-of-fit value of larger or equal than 100 in 56% of all pixels. In comparison to the standard Infomax this number drops down to about 15%.

	I) Infomax	II) cICAP
identified ICs (mean)	$7.7\mathrm{ }\pm 4.75$	3.8$\mathrm{ }\pm 1.19$
number of iterations (mean)	405.5$\pm \mathrm{ }32.3$	191.8$\pm \mathrm{ }23.13$
wrGOF (median)	18.8	1910.5
wrGOF < 1 (%)	3.57 %	0.35 %
wrGOF ≥ 10 (%)	38.8 %	80.0 %
wrGOF ≥ 100 (%)	15.1 %	55.6 %

Table 1.

Statistical analysis of ICA decomposition results. The signal components for the ICA routine I were selected by means of a user independent algorithm (see text).

In Figure 6a pixel locations of a representative section is shown, where the test statistic (wrGOF) is larger or equal 300. In this example signal enhancement was evident at all pixel locations (i.e., wrGOF > 1; not shown). By focusing on larger values of the test statistic (e.g., wrGOF ≥ 300) the figure demonstrates that the enhancement of the optical signals is ensured in particular in the region of the neo-cortex and at the white to gray matter boundary (Figure 6a). This in general is challenging, since in such regions the 3D-PLI signal has low amplitudes and therefore the measured signal has a poor signal-to-noise ratio.

With the introduction of the weighting factor (or penalty factor) as expressed by Equation (12), the test statistic wrGOF is sensitive to both, changes in the signal-to-noise ratio as well as in reductions of signal strength after ICA. This in particular will be the case when one or more components of interest are not selected for signal reconstruction (e.g., when the selection is performed manually). As a result, at corresponding pixel locations the wrGOF value will be less than 1. The idea behind this is, that the weighted test statistic wrGOF accounts for missing gray and white matter components and improvements in SNR, independently to the metric used for identification of components of interest. In Figure 6b such a testing scenario is demonstrated, where one component of interest was manually excluded before the 3D-PLI signal was reconstructed. For this section, 24% of all pixels have a wrGOF value below 1.

Figure 6.

${\chi }^2$ ratios and pixel locations illustrating goodness-of-fit values. In a) ${\chi }^2$ ratios are shown using all components identified by cICAP, where the goodness-of-fit statistics was larger than 1 at all pixel locations. In the coronal brain section pixels highlighted in red show locations with goodness-of-fit values larger than 300. b)To demonstrate the effect of being sensitive to missing signal of interest, one identified component was manually deselected. As a result, large areas of white matter structure do show a wrGOF value of less than 1 (highlighted in blue).

4. Conclusion

In recent years 3D-polarized light imaging (3D-PLI) has been shown to provide new insights into the organization of the human brain including mapping of the fiber anatomy [15,45,46]. Through advances in the experimental setup of the employed polarimeter, as well as in signal processing, this modality provides unique data sets to explore the 3D fiber architecture in the human brain at a submillimeter resolution [16,47,48]. Though the technique as described here is applicable solely to postmortem brain tissue, the comprehensive description of complex fiber orientations in distinct brain regions (e.g., prevalent fiber crossings) can be used to guide and evaluate fiber tractography algorithms based on diffusion MRI. By this means the fiber orientation maps provided by 3D-PLI might help to optimize the reliability of in-vivo diffusion MRI results. Precise information about the local individual fiber architecture of a patient is of particular interest in case of planning and performing a neurosurgical intervention, for instance.

However, in 3D-PLI the reconstruction of nerve fiber pathways strongly depends on the quality of the measured intensities represented by the so-called light intensity profiles or 3D-PLI signals, respectively. Hence, advanced signal processing tools are required to enable the precise determination of locally prevailing fiber orientations in form of unit vectors defined by the in-section direction angle φ and the out-of-section inclination angle α.

Independent component analysis (ICA) turned out to improve 3D-PLI signals significantly [22,23]. It was shown that ICA is capable of restoring the original birefringent signal by effectively removing noise and artifact components in the measured data. In addition, measures for the qualitative and quantitative evaluation of the 3D-PLI signals before and after the ICA filtering were introduced. In particular, the signal enhancement after ICA based denoising is large at the white to gray matter boundary, where the 3D-PLI signal is weak due to decreasing fiber density when approaching gray matter domains [22].

In cICAP, the “constrained ICA for polarized light imaging”, the signal decomposition is optimized using a weight update, which incorporates the prior knowledge that the expected signal in 3D-PLI can be modeled utilizing the Jones calculus (cf. Equation (1)). This prior knowledge gives rise to a training of the unmixing procedure in cICAP, where the weight matrix update is optimal with respect to the sources of interest (cf. Equation (10)). Utilizing this approach the source separation and identification in cICAP is carried out automatically in one routine, where no further component selection tool is needed for subsequent extraction of the signal of interest.

With the introduction of cICAP an ICA-based source separation method is available which is optimal for the extraction of the sinusoidal trajectory along different rotation angles in a set of spatially mixed 3D-PLI images. The better quality of the source separation in cICAP is reflected by the increased goodness-of-fit statistic $\left(wrGOF\right)$, which not only takes the improvement in the signal-to-noise ratio (SNR) into account, but also accounts for signal restoration [23].

The dedicated signal processing tool cICAP largely contributes to the reliability of 3D-PLI signals and, hence, the accuracy of subsequent fiber tract reconstruction. As a consequence, cICAP plays a key role in the 3D-PLI processing chain and helps to establish the so far missing link between the microscopic and the macroscopic characterization of the anatomical connectivity of the human brain.

Acknowledgments

We thank Prof. Dr. A. Wree, Institute for Anatomy, University of Rostock, for providing us with the postmortem human brain.