Report

Multispectral Imaging and Unmixing Jürgen Glatz Chair for Biological Imaging www.cbi.ei.tum.de Munich, 06/06/12 Intraoperative Fluorescence Imaging Fluorescence Channel Color Channel Outline Multispectral Imaging Unmixing Methods Exercise: Implementation Multispectral Imaging Multispectral Imaging Unmixing Methods Exercise: Implementation Multispectral Imaging Nature Spectral Resolution Sensitivity Range Spatial Resolution Magnification Technology Spectral Resolution has practically not improved since first camera Spatial Resolution and Magnification are significantly improved Color Vision Anyone feeling hungry? Monochrome image of an apple tree Color image of an apple tree • Color vision helps to distinguish and identify objects against their background (here: fruit and foliage) • Color vision provides contrast based on optical properties Color Vision Spectral sensitivity of the human eye short low light mid long wavelength blue green red blue green red perception blue green • Color receptors (cone cells) with different spectral sensitivity enable trichromatic vision • Limited spectral range and poor resolution red Ultraviolet Visible Limited spectral range Evening primrose Cleopatra butterfly • Human eyes can only see a portion of the light spectrum (ca. 400-750nm) • Certain patterns are invisible to the eye Limited spectral resolution Different chemical composition plastic chlorophyll • Color vision is insufficient to distinguish between two green objects Same color appearance • Differences in the spectra reveal different chemical composition blue green red blue green red Optical Spectroscopy • Absorbance • Fluorescence • Transmittance • Emission • Spectroscopy analyzes the interaction between optical radiation and a sample (as a function of λ) • Provides compositional and structural information Directions of optical Methods Imaging Spectroscopy Currently there are two “directions” in optical analysis of an object A Camera B Provides spatial information Reveals morphological features No information about structure or composition / no spectral analysis Spectrometer Provides spectral information Spectrum reveals composition and structure No information about spatial distribution Imaging Spectroscopy Imaging Spectroscopy Spatial dimension x Spatial information Spatial dimension y Spatial dimension y Imaging Spectroscopy Spectral dimension λ Spectral information Spectral Cube Spatial and spectral information Spectral Cube λ6 λ7 λ λ4 5 λ 3 λ1 λ2 λ8 Pseudo-color image representing the distribution of compounds A and B (chlorophyll and plastic) • Acquisition of spatially coregistered images at different wavelengths • The maximum number of components that can be distinguished equals the number of spectral bands • The accuracy of spectral unmixing increases with the number of bands Multispectral Imaging Modalities • Camera + Filter Wheel • Bayer Pattern • Cameras + Prism • Multispectral Optoacoustic Tomography • etc. Let’s find those apples Multispectral imaging alone is only one side of the medal Appropriate data analysis techniques are required to extract information from the measurements Unmixing Methods Multispectral Imaging Unmixing Methods Exercise: Implementation The Unmixing Problem Unmixing Finding the sources that constitute the measurements For multispectral imaging this means separating image components of different, overlapping spectra Unmixing is a general problem in (multivariate) data analysis Multifluorescence Microscopy Disjoint spectra can be separated by bandpass filtering Overlapping emission spectra create crosstalk Autofluorescence I λ Autofluorescence exhibits a broadband spectrum Only mixed observations of the components can be measured Post-processing to unmix them Forward Modeling What constitutes a multispectral measurement at a certain point and wavelength? Principle of superposition: Sum of individual component emission I r , i I 1 r , i I 2 r , i I k r , i A component‘s emission over different wavelengths λ is denoted by its spectrum, its spatial distribution is still to be defined. Setting up a simple forward problem (1) Two fluorochromes on a homogeneous background Note: We define images as row vectors of length n All components are merged in the (n x k) source matrix O n: Number of image pixels k: Number of spectral components Relative Absorption [%] Setting up a simple forward problem (2) Wavelength [nm] Defining the emission spectra for all components at the measurement points Combining them into the (k x m) spectral matrix k: Number of spectral components m: Number of multispectral measurements m≥k Relative Absorption [%] Setting up a simple forward problem (3) Wavelength [nm] • Two fluorochromes on a homogeneous background • Heavily overlapping spectra • 25 equidistant measurements under ideal conditions Mathematical Formulation M OS (+N) Multispectral measurement matrix Original component matrix Spectral mixing matrix Noise, artefacts, etc. (n x m) (n x k) (k x m) (n x m) Multispectral Dataset Mathematical Formulation • = 10000x25 10000x3 3x25 Linear Regression: Spectral Fitting • M OS Reconstructing O • System generally overdetermined: No direct inverse S-1 • Generalized inverse: Moore-Penrose Pseudoinvere S+ • umx MS • Spectral Fitting: Finding the components that best explain the measurements given the spectra • Minimizing the error: • S arg min e S 2 2 e O umx Spectral Fitting Orthogonality principle: optimal estimation (in a least squares sense) is orthogonal to observation space e span M e null M T Spectral Fitting e span M e null M M eM T O MS 0 T T T T M OSS T T M MS M MS M O T M MS SS M MS SS T S S SS S S T T T T SS T 1 T T T Spectral Fitting Spectral Fitting • Given full spectral information (i.e. about all source components) the data can be unmixed • S R S pinv T SS MS T 1 Blood oxygenation in tumors Multifluorescence Imaging RGB image FITC TRITC Nude mice with two different species of autofluorescence and three subcutaneous fluorophore signals: FITC, TRITC and Cy3.5. (Totally 5 signals) Cy3.5 Autofluorescence Food Composite Spectral Fitting Fast, easy and computationally stable Known order and number of unmixed components Quantitative Requires complete spectral information Crucially depends on accuracy of spectra (systematic errors) Suitable for detection and localization of known compositions Still no apples… ? ? S ? Principal Component Analysis • Blind source separation (BSS) technique • Requires no a priori spectral information • Estimates both O and S from M • Assumption: Cov ( o i , o j ) 0 Sources are uncorrelated, while mixed measurements are not Principal Component Analysis • Unmixing by decorrelation: Orthogonal linear transformation • Transforms the data into a space spanned by the orthogonal PCs • Maximum variance along first PC, maximum remaining variance along second PC, etc. Unmixing multispectral data with PCA • 25 multispectral measurements are correlated • Their entire variance can (ideally) be expressed by only 3 PCs Dimension reduction • Those 3 PCs are the unmixed sources • Note that matrix orientations may vary between different implementations Computing PCA Method 1 (preferred for computational reasons) • Subtract mean from multispectral observations • Covariance Matrix: CM Cov m 1 , m 1 Cov m m , m 1 Cov m 1 , m m Cov m m , m m • Diagonalizing CM: Eigenvalue Decomposition • Eigenvectors of CM are the principal components, roots of the eigenvalues are the singular values • Projecting M onto the PCs: R PCA U M T T Computing PCA with the SVD Method 2 (not suitable for implementation) • Subtract mean from multispectral observations • Singular Value Decomposition: M = UΣVT • U is a (m x m) matrix of orthonormal (uncorrelated!) vectors • Projecting M onto those decorrelates the measurements R PCA U M T T • Singular values in Σ denote how much variance is explained by the respective PC PCA does more than just unmix Mixing R PCA U M T S (UT)-1 = U ≈ S Multispectral data space PCA Original data space • UT U is a (non-quantitative) approximation of the PCs spectra • These can be used to verify a components identity • Σ is the singular value matrix Relatively small singular values indicate irrelevant components T PCA Spectra Principal Component Analysis (PCA) Needs no a priori spectral information Also reconstructs spectral properties Significance measurement through singular values Unknown order and number of components Generally not quantitative Crucially depends on uncorrelatedness of the sources Suitable for many compounds and identification of unknown components Advanced Blind Source Separation Independent Component Analysis (ICA): assumes statistically independent source components, which is a stronger condition than PCA’s orthogonality Non-negative Matrix Factorization (NNMF): constraint that all elements must be positive Commonly computed by iterative optimization of cost functions, gradient descent, etc. Independent Component Analysis • Assumes and requires independent sources: P o i o j P o i P o j • Independence is stronger than uncorrelatedness Independent Component Analysis • Central limit theorem: Sum of non-gaussian variables is more gaussian than the individual variables 3E x • Kurtosis measures non-gaussianity: kurt x E x • Maximize kurtosis to find IC • Reconstruction: R ICA U M T 4 2 2 Practical Considerations • Noise • Artifacts (from reconstruction, reflections, measurement,…) • Systematic errors (spectra, laser tuning, illumination,…) • Unknown and unwanted components Exercise: Implementation Multispectral Imaging Unmixing Methods Exercise: Implementation Forward Problem / Mixing • Define at least 3 non-constant images representing the original components • Plot them and store them in the matrix O • Define an emission spectrum for every component at an appropriate number of measurment points • Plot them and store them in the matrix S • Calculate the measurement matrix as M = OS (and save everything) Relative Absorption [%] Forward Problem / Mixing Wavelength [nm] O S Forward Problem / Mixing Useful MatLab functions • Change matrices into vectors: y=reshape(X,…) or y=X(:) • Plot image from a matrix: imagesc(X) or imshow(X) Spectral Fitting Create an m-file and write a function that • Has M and S as input variables • Calculates the pseudoinverse S+ • Returns the unmixing Rpinv • Test it on your data Spectral Fitting S R S pinv T SS MS T 1 Useful MatLab functions • Functions: function [out] = name([input]) • Regular matrix inverse: y = inv(x) Principal Component Analysis Create an m-file and write a function that • Has M as an input variable • Subtracts the mean from the measurements in M • Computes the covariance matrix CM • Performs an eigenvalue decomposition on CM • Sorts the eigenvalues (and corresponding vectors) by size • Projects M onto the eigenvectors • Returns the projected unmixing, the principal components and their loadings Principal Component Analysis CM Cov m 1 , m 1 Cov m m , m 1 Cov m 1 , m m Cov m m , m m Cov x , y n 1 x n 1 i 1 R PCA U M T Useful MatLab functions • Mean: y = mean(x) • Eigenvalue Decomposition: [e_vec e_val] = eig(X) i x y i y Testing your code • Try fitting and PCA on your mixed data • Try adding different types and amounts of noise to M (e.g. using imnoise) • Simulate systematic errors in your spectra (noise, changing values, offset,…) Independent Component Analysis (voluntary) You can download the FastICA MatLab code from http://research.ics.tkk.fi/ica/fastica/ Type doc fastica for function description Use the fastica function to unmix your simulated data Compare the result to PCA. What are advantages and disadvantages of ICA? Recommended Reading • Shlens, J. – A Tutorial on Principal Component Analysis http://www.cfm.brown.edu/people/gk/APMA2821F/PCA-Tutorial- Intuition_jp.pdf • Garini, Y., Young, I.T. and McNamara, G. – Spectral Imaging: Principles and Applications; Cytometry Part A 69A: p.735-747 (2006) http://dx.doi.org/10.1002/cyto.a.20311 • Stone, J.V. – A brief Introduction to ICA; Encyclopedia of Statistics in Behavioral Science, Vol. 2, p. 907-912 http://jimstone.staff.shef.ac.uk/papers/ica_encyc_jvs4everrit2005.pdf