Matlab Routines

openairsar.m

Source: Airsar Integrated Processor Documentation (Data Formats)
`Understanding Synthetic Aperture Radar Images' by C. Oliver and S. Quegan, 1998
Direct evaluation in Mathematica - compressed_stokes_simple.nb

openairsar_new.m

Source: Airsar Integrated Processor Documentation (Data Formats)
`Understanding Synthetic Aperture Radar Images' by C. Oliver and S. Quegan, 1998
Direct evaluation in Mathematica - compressed_stokes_simple.nb

posdef_swap.m

Source: see www.mathworld.com for the conditions of positive definite matrices.

coherence_distribution.m

intensity_distribution.m

multiple_polsep.m

power_separation.m

phase_distribution.m

Returns the phase-difference probability distribution function between a pair of channels given the complex correlation coefficient rho. Applicable to both polarimetric and interferometric radar if the underlying statistics are assumed to conform to a Complex Wishart distribution. In general the phase of a single scattered wave has a uniform distribution which provides no useful information. When two channels are present, dependent on the correlation rho, the phase difference between the channels can provide additional classification power (polarimetric) or an interferometry map.

single_polsep.m

cov_channels.m

Generates N samples from channels with normalised correlated complex Gaussian noise of arbitrary correlation, covariance matrix rho. If rho is a scalar then two channels are output with correlation rho, apart from this corrcoef(S) should match rho. The covariance matrix input should be hermitian positive definite.

cov_matrices.m

Where rho is a p x p covariance matrix, this outputs p x p x N sample covariance matrices drawn from the Wishart distribution for L looks. This is the general form of polarimetric and interferometric radar statistics. If rho is a scalar then 2x2 matrices are output as for a correlation co-efficient to exist there must be at least two channels, apart from this mean(S,3) should match rho. The covariance matrix input should be hermitian positive definite. Note the determinant of single look samples should equal zero, but this means they are not strictly positive-definite.
See cwishart_variates.m for an alternative method.

cwishart_variate.sm (not recommended - use cov_matrices.m)

This is an attempt at an order N method for generating N sample covariance matrices from the complex Wishart distribution. For low L looks it appears to be too sensitive to precision and too many outputs are not positive-definite. Note the determinant of single look samples should equal zero, but this means they are not strictly positive-definite. See the source code for more information.
See cov_matrices.m for a more robust order N X L method.

RT_cdf.c

P is the integral from 0 to param1 of distribution cdftype - see dcdflib.fdoc for usage. Q is (1-P) which maintains full double precision when P approaches 1. cdftype = [1..12] as beta(1), binomial(2), chisquare(3), non-central chisquare(4), F(5), non-central F(6), gamma(7), negative binomial(8), normal(9), poisson(10), student T(11), non-central student T(12). Use RT_cdf(cdftype) to display required parameters.

RT_expcdf.m

Returns the cumulative distribution function for the exponential distribution , i.e. C=Integral_0^x[P(x|mu)dx]. This is applicable to the intensity distribution of single look speckle with a mean of mu . Both x and mu can be matrices with compatible sizes.

RT_expinv.m

Returns the inverse cumulative distribution function for the exponential distribution , i.e. Integral_0^x[P(x|mu)dx]=P limited to input 0<P<1. This is applicable to the intensity distribution of single look speckle with a mean of mu . Both x and mu can be matrices with compatible sizes.

RT_exppdf.m

Returns the probability distribution function for the exponential distribution , i.e. P(x|mu)dx. This is applicable to the intensity distribution of single look speckle with a mean of mu . Both x and mu can be matrices with compatible sizes.

RT_exprnd.m

Given scalar parameter mu, this routine returns a matrix of exponential variates of size [rows,cols]. This is applicable to the intensity distribution of single look speckle with a global mean of mu. Alternatively a speckle matrix with a varying underlying mean can be generated by passing this as the matrix mu and omitting rows and cols.

RT_gamcdf.m

Returns the cumulative distribution function for the gamma distribution as P. The complementary cumulative distribution is returned as Q, but calculated separately. Q is (1-P) which maintains full double precision when P approaches 1. shape can be non-integer, but integer values correspond to multi-look speckle.

RT_gaminv.m

Returns the inverse cumulative distribution function for the gamma distribution. shape can be non-integer, but integer values correspond to multi-look speckle. 0 < p < 1. Input the complementary value Q if p > (1-eps) such that P + Q = 1.

RT_gampdf.m

Returns the probability distribution function for the gamma distribution. shape can be non-integer, but integer values correspond to multi-look speckle.

RT_gamrnd.m

Returns a matrix of gamma variates of size [rows,cols]. shape can be non-integer, but integer values correspond to multi-look speckle.

RT_invcdf.c

X is returned so that the integral from 0 to X of distribution cdftype gives P, where Q = (1-P). See dcdflib.fdoc for usage. cdftype = [1..12] as beta(1), binomial(2), chisquare(3), non-central chisquare(4), F(5), non-central F(6), gamma(7), negative binomial(8), normal(9), poisson(10), student T(11), non-central student T(12). Use RT_invcdf(cdftype) to display required parameters.

RT_Kcdf.m

Returns the cumulative probability distribution function for the K distribution of mean mu. The distribution approaches the exponential distribution as the shape parameter nu goes to infinity. Returns UERR when the besselk function has problems, pass noflag!=0 to suppress this.

RT_Kcdf_noise.m

Returns the cumulative probability distribution function for the K distribution of mean mu in the presence of thermal noise of variance rho. The overall mean is (mu + rho). The distribution approaches the exponential distribution as the shape parameter nu goes to infinity or if rho>>mu. Attempts to achieve an absolute accuracy of tol (default sqrt(eps)) , and will produce a warning if the relative accuracy is less than 0.1%. Some effort was required to get this numerically stable and it may fail for extreme nu values - manually use the limiting distribution when this is the case. For reasonable nu, the limiting distribution is returned based on the Clutter(mu) to Noise(rho) Ratio. CNR>+30dB implies pure K-distribution, and CNR<-30dB implies exponential distribution.

RT_Kpdf.m

Returns the probability distribution function for the K distribution of mean mu. The distribution approaches the exponential distribution as the shape parameter nu goes to infinity. Returns UERR when the besselk function has problems, pass noflag!=0 to suppress this.

RT_Kpdf_noise.m

Returns the probability distribution function for the K distribution of mean mu in the presence of thermal noise of variance rho. The overall mean is (mu + rho). The distribution approaches the exponential distribution as the shape parameter nu goes to infinity or if rho>>mu. Attempts to achieve an absolute accuracy of tol (default sqrt(eps)) , and will produce a warning if the relative accuracy is less than 0.1%. Some effort was required to get this numerically stable and it may fail for extreme nu values - manually use the limiting distribution when this is the case. For reasonable nu, the limiting distribution is returned based on the Clutter(mu) to Noise(rho) Ratio. CNR>+30dB implies pure K-distribution, and CNR<-30dB implies exponential distribution.

RT_Krnd.m

Returns a matrix of K variates of size [rows,cols]. S has an expected mean of mu and a shape parameter of nu. The compound K formulation is used but Independent Identically Distributed variates are returned.

RT_Krnd_noise.m

Returns a matrix of K variates in the presence of thermal noise of size [rows,cols]. S has an expected mean of (mu + rho) and a (noise contaminated) shape parameter of nu. The variance of the zero mean Gaussian noise is rho. The compound K formulation is used but Independent Identically Distributed variates are returned.

RT_weibcdf.m

Returns the cumulative distribution function for the Weibull distribution as P. The form used is such that the mean value is mu and larger a is closer to a Gaussian distribution. a = 1 is equivalent to an exponential distribution. The complementary cumulative distribution is returned as Q, but calculated separately. Q is (1-P) which maintains full double precision when P approaches 1.

RT_weibest.m

Estimates the parameters a and mu of a Weibull distribution that best fits the data in X. The form used is such that the mean value is mu and larger a is closer to a Gaussian distribution. a = 1 is equivalent to an exponential distribution. The method uses an initial least squares fit to Weibull paper. It is not the Maximum Likelihood Solution, but a is reasonably accurate (mu is quite poor). To return this estimate, set the fast flag. With fast unset or 0, a numerical MLE solution is generated (slow) but this can fail when a is too small (a << 0.1).

RT_weibinv.m

Returns the inverse cumulative distribution function for the Weibull distribution. The form used is such that the mean value is mu and larger a is closer to a Gaussian distribution. a = 1 is equivalent to an exponential distribution. 0 < P < 1. Input the complementary value Q if P> (1-eps) such that P + Q = 1.

RT_weibpdf.m

Returns the probability distribution function for the Weibull distribution. The form used is such that the mean value is mu and larger a is closer to a Gaussian distribution. a = 1 is equivalent to an exponential distribution.

RT_weibrnd.m

Returns a matrix of Weibull variates of size [rows,cols]. The form used is such that the mean value is mu and larger a is closer to a Gaussian distribution. a = 1 is equivalent to an exponential distribution.

cgamgen.m

Returns the mapping of a correlation Rg under an MNLT for a measured gamma correlation. When a correlated Gaussian is directly mapped to a gamma distribution using a Memoryless Non-Linear Transform (MNLT) the correlation structure is not preserved. It has been shown by Tough and Ward (see below) how to correctly obtain the desired correlation structure while maintaining the distribution. The demo shows how to invert this to obtain a general gamma shape distribution with a desired correlation structure. Many additional functions are associated with this command.

covfullint.c

Fully polarimetric, multifrequency classification algorithm based on expectation maximisation of the covariance matrices. The polarisation inputs (HHHH etc.) are vectors of image segments size [Ndata, Nfreq] this requires use of the utility reclassrho.c to subsequently map to the original image. Note that the weights are output based on the number of segments, i.e. they are not weighted wrt segment population. The initialisation is based on a dB K-means search over the first frequency. A maximum of seven classes is recommended although 9 may be selected. If Ndata <= 10*Nclasses the initialisation algorithm usually cannot proceed.

nleU.m

Returns the normalised log estimate U of the data, a reasonable parameter estimate for long-tailed distributions. For an image pass data(:) otherwise it will process column-wise.

nleU2nu.m

Converts the U estimate to the K distribution parameter nu. Inverts U=log(nu) - psi(nu) - psi(1) by an iterative scheme. Does not use a numerical search so it's very fast.

nu2nleU.m

Converts the K distribution parameter nu to the equivalent normalised measure U. U=log(nu) - psi(nu) - psi(1)

RT_cgaussian_iir.m

Returns Nsamples of correlated zero-mean, unit-variance gaussian variates with an approximate autocorrelation function of acf. Note that Nsamples does not have to be the same size as acf. It uses ac2poly from the Signal Processing toolbox. This function is not necessarily stable as it uses an IIR filter, however this can make it more accurate. err can be returned to give the overall 'prediction error'. The normalised, one sided form for acf should be used.

RT_cgaussian_fft.m

Returns up to two innovations of correlated zero-mean, unit-variance gaussian variates with an approximate autocorrelation function of acf. It shapes the power spectrum via a standard FFT method and so the number of samples is limited by the length of acf. The normalised, one sided form for acf should be used.

segfullX.c

Fully polarimetric, multifrequency segmentation algorithm that merges a scene from single pixels. The method is reasonably fast for what it does but it was developed for a computer with 2GB memory so speed was optimised rather than space. Each polarisation input (HH_HH etc.) should have the size [X,Y{,Nfreqs}] for image size X by Y. The pixel connectivity can be 4 or 8, but I recommend only using 8 otherwise you can't segment diagonal lines. The number of looks is input so that the mean variance ratio of the segmented and original co-polarised images can be calculated. By setting no_regions=0 , the algorithm will automatically stop when the variance ratio matches the expected value. Otherwise, the desired no_regions can be input - but be warned that oversegmentation can take a very long time. A positive penalty value can be used to constrain the surface tension of each individual segment, typical values are 0.01. For multilook data such as AirSAR there is no need for any penalty (use 0.0), but for a two look radar the images can be significantly improved. See the theory in Stewart et al. for a good explanation, but check the source code for the exact algorithm used. The routine uses a slightly modified version of Red Black Tree code by kind permission of James S. Plank at http://www.cs.utk.edu/~plank/

Source: Theory explained in `Polarimetric Classification using Expectation Methods', G Davidson et al.. This is a very complicated routine. An undocumented Bayesian input is also possible - see the source code.
`Optimal approach to SAR image segmentation and classification', D. Stewart, D. Blacknell, A. Blake, R. Cook and C. Oliver, IEE Proc. Radar, Sonar and Navigation, 147(3), pp. 134-148, June 2000.

RT_cfar.m

A wrapper function for the cfarmex.c file which implements fast CFAR processors for 1D vector data. Returns the test statistic Z (as defined in [1]) for the following CFAR models: CA - Cell Average, CAGO - Cell Average Greatest Of, CASO - Cell Average Smallest Of, OS - Order Statistic (k1th largest), TM - Trimmed Mean (k1 smallest & k2 largest trimmed). The total window size Nwindow and the number of guard cells Nguard must be even numbers. These are then split equally, centred over the test cell.

RT_fcdf.m

Returns the F cumulative probability density with degrees of freedom of df1(numerator) and df2(denominator) at the values in vector X. The F distribution arises from the ratio of two homogenous multilook speckle fields. Strictly speaking, df1 and df2 should be integers. Q is returned as the complementary value, and is calculated separately from P via dcdflib. This allows the full range of realmin < x < realmax to be processed for inversion via RT_finv.

RT_finv.m

Returns the inverse of the F cumulative probability density with degrees of freedom of df1(numerator) and df2(denominator) for the value P(0 < P< 1). The F distribution arises from the ratio of two homogenous multilook speckle fields. Input the complementary value Q if P> (1-eps) such that P + Q = 1.

RT_fpdf.m

Returns the F probability density with degrees of freedom of df1(numerator) and df2(denominator) at the values in vector X. The F distribution arises from the ratio of two homogenous multilook speckle fields. Strictly speaking, df1 and df2 should be integers.

RT_Pcfar.m

Given a linear threshold T, returns the CFAR probability for models: CA - Cell Average, CAGO - Cell Average Greatest Of, CASO - Cell Average Smallest Of, OS - Order Statistic (k1th largest), TM - Trimmed Mean (k1 smallest & k2 largest trimmed). N is the total background window size which has to be an even number. Using S = 0 returns the Probability of False Alarm. Passing a linear target strength S > 0 returns the Probability of Detection for a Swerling 1 target. T can be an array if size is compatible with S. Use RT_Tcfar.m to invert for threshold T from probability P.

RT_Tcfar.m

Inverts RT_Pcfar.m to obtain the required threshold given a probability goal P. See RT_Pcfar.m for parameters, pass non-zero trace to follow the integration steps.

swerling_pd.m

Gives the Probability of Detection with respect to number of pulses for a class of Swerling target models. Signal to Noise Ratio SNR should be in linear format (not dB). The fluctuating target models are numbered 1 to 4. A non-fluctuating target is commonly referred to as model 0 and 5, both are accepted. Some cases are necessarily an approximation, don't rely on the exact value but in practice it should be accurate enough. Additionally, the approximation error is far outweighed by the variability in real radar data. The Swerling theory assumes an optimum processor, see cfar_theory_demo.m for the losses associated with practical processors. Uses logs for numerical stability.