Course Outline¶
Quick tour of the IPython Notebook
Mutlispectral visual/infrared images
Probability distributions
Principal components
Clustering
SAR and polarimetric SAR images
Data acquisition and multilooking
Probability distributions
Speckle filtering
Multispecral change detection and radiometric normalization
Polarimetric SAR change detection
Software and Notebook Installation on Windows¶
Install the Docker Engine from https://docs.docker.com/
Ensure that hardware virtualization is enabled.
Download DockerToobox from https://www.docker.com/docker-toolboxand
and install with all defaults. Be sure to accept all Oracle drivers.
In the Docker Quickstart terminal, run the commands
docker pull mort/zfldocker
docker pull mort/datadocker
docker run -d -v /imagery --name=data mort/datadocker
docker run -d -p 463:8888 --name=zfl --volumes-from=/data mort/zfldockerRun the command
docker-machine ip default
to get the IP address of your virtual machine.Point your browser to
http://ip-address:463Click on the course notebook /home/imagery/ZFL-Course.ipynb,
or open a new notebook with New/Python 2.
!ls -l /home
The /home/imagery directory contains the image data and this notebook:
!ls -l /home/imagery
The /home/imagery/figures directory contains some figures:
!ls -l /home/imagery/figures
which can be displayed within a markdown cell:
SAR imaging geometry¶
We can talk directly to the Python kernel:
import numpy as np
np.random.rand(5,5)
To enable graphics output, we run an IPython "magic"
%matplotlib inline
import matplotlib.pyplot as plt
x = np.linspace(0,10,100)
plt.plot(x,np.sin(x))
plt.show()
Python programs (scripts) are called with the run command:
run /home/dispms.py -h
2. Multispectral Visual/Infrared Images¶
2.1 Scalar (one-band) images¶
The statistical distributions of optical and SAR measurements are fundamentally different: SAR intensity observations are exponentially distributed, optical/infrared observations are normally distributed.
The normal or Gaussian density function of a random variable $Z$ representing an image pixel is
$$ p(z) = {1\over \sqrt{2\pi}\sigma}\exp\left(-{1\over 2\sigma^2}(z-\mu)^2\right),\quad $$where $-\infty<\mu<\infty$ and $\sigma^2>0$. The mean value and variance are given by
$$ \langle Z\rangle = \int_{-\infty}^\infty z\cdot p(z) dz = \mu,\quad $$$$ {\rm var}(Z) = \langle(Z-\langle Z\rangle)^2\rangle = \int_{-\infty}^\infty (z-\langle Z\rangle)^2 p(z) dz = \sigma^2.\quad $$This is commonly abbreviated by writing
$$ Z \sim \mathcal{N}(\mu,\sigma^2).\quad $$Here is a plot of the normal distribution for $\mu=0$ and $\sigma=1$ (the standard normal distribution)
import scipy.stats as st
z = np.linspace(-5,5,2000)
plt.plot(z,st.norm.pdf(z))
plt.show()
Parameter estimation¶
Consider a single-band image and a specific land cover category within it. We might choose $m$ pixels $Z_i$, $i=1\dots m$, belonging to that category and use them to estimate the moments of the underlying distribution. That distribution will be determined not only by measurement noise, atmospheric disturbances, etc., but also by the spread in reflectances characterizing the land cover category itself.
We can estimate the mean and variance of the distribution with two sample functions: The sample mean
$$ \bar{Z} = {1\over m}\sum_{i=1}^m Z_i\quad $$with realizations for actual measurements $z_i$
$$ \bar{z} = {1\over m}\sum_{i=1}^m z_i\quad $$and the sample variance
$$ S = {1\over m-1}\sum_{i=1}^m(Z_i - \bar Z )^2.\quad $$with realizations
$$ s = {1\over m-1}\sum_{i=1}^m(z_i - \bar z )^2.\quad $$These two sample functions are also called unbiased estimators because their average values are equal to the corresponding moments of the distribution, that is,
$$ \langle\bar Z\rangle = \langle Z\rangle\quad ( = \mu ) $$and
$$ \langle S \rangle = {\rm var}(Z)\quad ( = \sigma^2). $$Distributions of the sample mean and variance¶
The sample mean $\bar Z$ is normally distributed with mean $\mu$ and variance $\sigma^2/m$.
$$ p(z) = {1\over \sqrt{2\pi\over m}\sigma}\exp\left(-{1\over 2(\sigma^2/m)}(z-\mu)^2\right),\quad $$For the sample variance $S$ then
$$ (m-1)S/\sigma^2 \quad $$is independent of $\bar Z$ and has the chi-square distribution with $m-1$ degrees of freedom.
The gamma distribution and its relatives¶
A random variable $Z$ is said to have a gamma distribution if its probability density function is given by
$$ p_\gamma(z) = {1\over \beta^\alpha\Gamma(\alpha)}z^{\alpha-1}e^{-z/\beta}, \quad z>0, $$where $\alpha>0$ and $\beta>0$ and where the gamma function $\Gamma(\alpha)$ is given by
$$ \Gamma(\alpha) = \int_0^\infty z^{\alpha-1} e^{-z}dz,\quad \alpha>0. \quad $$The gamma function has the recursive property
$$ \Gamma(\alpha) = (\alpha-1)\Gamma(\alpha-1),\quad \alpha>1, \quad $$and generalizes the notion of a factorial:
$$ \Gamma(n) = (n-1)!\ . $$The gamma distribution has mean and variance
$$ \langle Z\rangle = \alpha\beta, \quad {\rm var}(Z) = \alpha\beta^2. \quad $$A special case of the gamma distribution arises for $\alpha = 1$. Since $\Gamma(1) = \int_0^\infty e^{-z}dz = 1$, we obtain the exponential distribution with density function
$$ p_e(z) = \cases{{1\over \beta}e^{-z/\beta} & for $z>0$\cr \quad 0 & elsewhere,} $$where $\beta >0$. The exponential distribution has mean $\beta$ and variance $\beta^2$. In addition If random variables $Z_1,Z_2\dots Z_m$ are independent and exponentially distributed, then
$$ Z = \sum_{i=1}^m Z_i \quad $$is gamma distributed with $\alpha = m$.
The chi-square distribution with $m$ degrees of freedom is another special case of the gamma distribution. We get its density function with $\beta = 2$ and $\alpha = m/2$, i.e.,
$$ p_{\chi^2;m}(z)=\cases{{1 \over 2^{m/2}\Gamma(m/2)} z^{(m-2)/2} e^{-z/2} & for $z>0$\cr 0 & otherwise.} \quad $$It follows that the chi-square distribution has mean $\mu = m$ and variance $\sigma^2 = 2m$.
Here are plots of the chi-square distribution for different degrees of freedom:
z = np.linspace(1,40,200)
for m in [1,2,3,4,5,10,20]:
plt.plot(z,st.chi2.pdf(z,m))
plt.show()
The gamma and exponential distributions will be essential for the characterization of SAR speckle noise in SAR imagery. The chi-square distribution plays a central role in the iterative change detection algorithm IR-MAD that we will be looking at later:
If the random variables $Z_i$, $i=1\dots m$, are independent and standard normally distributed (i.e., with mean $0$ and variance $1$), then the random variable $Z=\sum_{i=1}^m Z_i^2$ is chi-square distributed with $m$ degrees of freedom.
2.2 Vector (multispectral or polarimetric) images¶
# vector and matrix multiplication
a = np.mat([[1],[2]])
print 'column vector\n %s' %str(a)
print 'row vector\n %s' %str(a.T)
print 'inner product\n %s' %str(a.T*a)
print 'outer product\n %s' %str(a*a.T)
b = np.mat([[1,2],[3,4]])
print 'matrix\n %s' %str(b)
print 'matrix times vector\n %s' %str(b*a)
print 'matrix times matrix\n %s' %str(b*b)
print 'matrix inverse\n %s' %str(b.I)
print 'identity matrix\n %s' %str(b*b.I)
print 'singular matrix\n %s' %str((a*a.T).I)
A measured pixel or observation in an N-band multispectral image corresponds to a random vector
$$ Z = \pmatrix{Z_1\cr \vdots \cr Z_N} \quad $$i.e., a vector the components of which are random variables.
The covariance of two random variables $Z_1$ and $Z_2$ is a measure of how their realizations are dependent upon each other and is defined as
$$ {\rm cov}(Z_1,Z_2) = \left\langle\ (Z_1-\langle Z_1\rangle)(Z_2-\langle Z_2\rangle)\ \right\rangle.\quad $$Useful relations
$$ {\rm cov}(Z_1,Z_1) = {\rm var}(Z_1)\\ {\rm cov}(Z_1,Z_2) = \langle Z_1 Z_2\rangle - \langle Z_1\rangle\langle Z_2\rangle $$The covariance matrix is a convenient way of representing the variances and covariances of the components of a random vector. Let $a = (a_1,a_2)^\top$ be any constant vector. Then the variance of the random variable $ a^\top Z = a_1Z_1+a_2Z_2$ is
$$ {\rm var}(a^\top Z) = {\rm cov}(a_1Z_1+a_2Z_2,a_1Z_1+a_2Z_2) \\ =a_1^2{\rm var}(Z_1)+a_1a_2{\rm cov}(Z_1,Z_2)+a_1a_2{\rm cov}(Z_2,Z_1)+a_2^2{\rm var}(Z_2)\\ =(a_1,a_2) \pmatrix{{\rm var}(Z_1)&{\rm cov}(Z_1,Z_2)\cr{\rm cov}(Z_2,Z_1)&{\rm var}(Z_2)} \pmatrix{a_1\cr a_2} \\ = a^\top\Sigma a. $$Thus:
$$ {\rm var}(a^\top Z) = a^\top\Sigma a. $$The random vector $Z$ has a multivariate normal density function when
$$ p(z) = {1\over (2\pi)^{N/2}\sqrt{|\Sigma|}} \exp\left(-{1\over 2}(z-\mu)^\top \Sigma^{-1} (z-\mu)\right), \quad $$The mean vector is
$$ \langle Z\rangle = \mu \quad $$and $\Sigma$ is the covariance matrix. This is indicated by writing
$$ Z \sim \mathcal{N}(\mu,\Sigma).\quad $$Parameter estimation¶
The vector sample mean is given by
$$ \bar{Z} = {1\over m}\sum_{\nu=1}^m Z(\nu) $$and the sample covariance matrix by
$$ S = {1\over m-1}\sum_{\nu=1}^m( Z(\nu) - \bar Z ) (Z(\nu) - \bar Z)^\top. $$If we can arrange that the mean values are zero, then
$$ S = {1\over m-1}\sum_{\nu=1}^m Z(\nu) Z(\nu)^\top. $$In the limit for large $m$,
$$ \Sigma = \langle Z Z^\top\rangle. $$Principal components¶
The principal components transformation, also called principal components analysis (PCA), generates linear combinations of multi-spectral pixel intensities which are mutually uncorrelated and which have maximum variance.
pca.ipynb
Consider a multispectral image represented by the random vector $Z$ ($Z=$ vector of gray-scale values) and assume that $\langle Z\rangle=0$, so that the covariance matrix is given by
$$ \Sigma=\langle Z Z^\top\rangle. $$Let us seek a linear combination
$$ Y = w^\top Z $$whose variance
$$ {\rm var}(Y) = w^\top \Sigma w $$is maximum. This quantity can trivially be made as large as we like by choosing $w$ sufficiently large, so that the maximization only makes sense if we restrict $w$ in some way. A convenient constraint is
$$ w^\top w=1. $$We can solve this problem by maximizing the Lagrange function
$$ L = w^\top \Sigma w - \lambda(w^\top w - 1). $$This leads directly to the eigenvalue problem
$$ \Sigma w = \lambda w. $$This procedure is coded in the python script pca.py.
ls -l /home/imagery
run /home/pca /home/imagery/20010525
run /home/dispms -f /home/imagery/20010525_pca -p [1,2,3] -e 3
Denote the orthogonal and normalized eigenvectors of $\Sigma$ obtained by solving the above problem by $w_1\dots w_N$, sorted according to decreasing eigenvalue $\lambda_1\ge \dots\ge\lambda_N$. These eigenvectors are the principal axes and the corresponding linear combinations $w_i^\top Z$ are projections along the principal axes, called the principal components of $Z$. The individual principal components
$$ Y_1 = w_1^\top Z,\ Y_2= w_2^\top Z,\ \dots,\ Y_N= w_N^\top Z $$can be expressed more compactly as a random vector $Y$ by writing
$$ Y = W^\top Z, $$where $W$ is the matrix whose columns comprise the eigenvectors, that is,
$$ W = (w_1\dots w_N). $$Since the eigenvectors are orthogonal and normalized, $W$ is an orthonormal matrix:
$$ W^\top W= I. $$If the covariance matrix of the principal components vector $ Y$ is called $\Sigma'$, then we have
$$ \Sigma' = \langle Y Y^\top \rangle= \langle W^\top Z Z^\top W\rangle = W^\top\Sigma W $$or, since $\Sigma w_i = \lambda_i w_i$,
$$ \Sigma'= W^\top \Sigma W = \pmatrix{\lambda_1 & 0 & \cdots &0\cr 0 &\lambda_2& \cdots&0\cr \vdots&\vdots&\ddots&\vdots\cr 0 & 0& \cdots&\lambda_N} $$The eigenvalues are thus seen to be the variances of the principal components, and all of the covariances are zero. The first principal component $Y_1$ has maximum variance ${\rm var}(Y_1)=\lambda_1$, the second principal component $Y_2$ has maximum variance ${\rm var}(Y_2)=\lambda_2$ subject to the condition that it is uncorrelated with $Y_1$, and so on. The fraction of the total variance in the original multispectral image which is accounted for by the first $i$ principal components is $$ {\lambda_1+\dots +\lambda_i\over\lambda_1+\dots +\lambda_i+\dots +\lambda_N}. $$
Let's look at the first and last principal components of the image above:
run /home/dispms -f /home/imagery/20010525_pca -p [1,1,1] -e 4 -F /home/imagery/20010525_pca -P [6,6,6] -E 4
Distribution of sample mean vector and sample covariance matrix¶
In the multivariate case we have a similar situation to the univariate case. The sample mean vector
$$ \bar Z = {1\over m}\sum_{\nu=1}^m Z(\nu) $$is multivariate normally distributed with mean $\mu$ and covariance matrix $\Sigma/m$ for samples
$$ Z(\nu)\sim\mathcal{N}(\mu,\Sigma),\quad \nu = 1\dots m. $$Suppose
$$ Z(\nu)\sim\mathcal{N}(0,\Sigma),\quad\nu=1\dots m. $$Then the sample covariance matrix
$$ S = {1\over m-1}\sum_{\nu=1}^m Z(\nu) Z(\nu)^\top $$is described by a Wishart distribution:
Define
$$ X = (m-1)S = \sum_{\nu=1}^m Z(\nu) Z(\nu)^\top. $$The probability density function of $X$ is the Wishart density with $m$ degrees of freedom:
$$ p_W(x) ={|x|^{(m-N-1)/2}\exp(-{\rm tr}(\Sigma^{-1} x)/2) \over 2^{Nm/2}\pi^{N(N-1)/4}|\Sigma|^{m/2}\prod_{i=1}^N\Gamma[(m+1-i)/2]} $$for $x$ positive definite, and $0$ otherwise.
We write
$$ X \sim \mathcal{W}(\Sigma,N,m). $$Since the gamma function $\Gamma(\alpha)$, is not defined for $\alpha\le 0$, the Wishart distribution is undefined
for $m We don't need the Wishart distribution here, but later when we talk about SAR change detection, we will need its analog for complex multinormal random variables representing polarimetric data: The probability density function of $X$ is
the complex Wishart density with $m$ degrees of freedom: This is denoted
Clustering or unsupervised classification¶
We represent the observed pixels in a multispectral image by the set
$$ \{z(\nu) \mid \nu = 1\dots m\}. $$A given partitioning or clustering may be written in the form
$$ C=[C_1,\dots C_k,\dots C_K], $$where
$$ C_k = \{\ \nu\mid z(\nu)\hbox{ is in class } k \}. $$Let $u_{k\nu}$ be the probability that pixel $\nu$ is class $C_k$ given its observed value $z(\nu)$
$$ u_{k\nu} = {\rm Pr}(k \mid z(\nu)) $$This is called the posterior probability and can be written, from Bayes' Rule, as
$$ u_{k\nu} = {p(z(\nu)\mid k){\rm Pr}(k)\over p(z(\nu))} $$${\rm Pr}(k)$ is called the prior probability for class $C_k$.
Now assume that the pixels in each class are (approximately) multivariate normally distributed with mean $\mu_k$ and covariance matrix $\Sigma_k$. Then
$$ u_{k\nu} = {1\over \sqrt{|\Sigma_k|}}\exp\left[-{1\over 2}(z(\nu)-\mu_k)^\top\Sigma_k^{-1} (z(\nu)-\mu_k)\right]p_k\quad (1) $$apart from a normalization constant. In the above expression the prior probability ${\rm Pr}(k)$ has been replaced by
$$ p_k = {m_k\over m} = {1\over m}\sum_{\nu=1}^m u_{k\nu}. \quad (2) $$To complete the picture we still need expressions for $\mu_k$ and $\Sigma_k$.
$$ \mu_k = {1\over m_k}\sum_{\nu\in C_k}z(\nu) = {\sum_{\nu=1}^m u_{k\nu}z(\nu) \over \sum_{\nu=1}^m u_{k\nu}},\quad k=1\dots K, \quad (3) $$$$ \Sigma_k= {\sum_{\nu=1}^m u_{k\nu}(z(\nu)-\mu_k) (z(\nu)-\mu_k)^\top \over \sum_{\nu=1}^m u_{k\nu}},\quad k=1\dots K.\quad (4) $$Algorithm (Expectation maximization clustering)¶
Choose random values for the initial class memberships $u_{k\nu}$ with $\sum_{k=1}^K u_{k\nu}=1$
Determine the cluster means $\mu_k$ with Equation (3) , then the covariance matrices $\Sigma_k$ for each cluster with Equation (4) and the priors $p_k$ with Equation (2)
Calculate the new class membership probabilities $u_{k\nu}$ with Equation (1), and normalize to ensure that $\sum_{k=1}^K u_{k\nu}=1$
If the $u_{k\nu}$ values have stopped changing significantly, stop, else go to step 2.
run /home/em -h
run /home/em -p [1,2,3,4] -K 8 /home/imagery/20010525_pca
run /home/dispms -c [1,2,3,4,5,6,7,8] -f /home/imagery/20010525_pca_em -e 2
