#acl AdminGroup:read,write EditorGroup:read,write All:read
#format wiki
#language en
#pragma section-numbers off


= Random variables =

The measurements that form an image show statistical fluctuations. The image is characterised by a probability density function (PDF). The PDF ''f'' describes the probability of the occurrence of a discrete grey level ''q'' in the range of grey levels ''Q''.

{{{#!latex
$\sum_{q=1}^Qf_q=1$
}}}

== Probability density functions and histograms ==

The PDF of an image can be calculated and displayed using the {{{pylab.hist}}} function or it can be
calculated using the {{{scipy.histogram}}} function.

=== ASAR image of sea ice ===

{{attachment:seaice_intensity.png}} {{attachment:seaice_intensity_db.png}}

The received power (intensity) of a radar system is proportional to the (normalized) radar backscatter coefficient 
{{{#!latex 
$\sigma^0(f,\theta)$
}}} as a function of frequency and incidence angle. The backscatter coefficient describes how much of the transmitted energy is backscattered from the surface media. 

A (calibrated) ASAR image ''img(y,x)=I'' is made of the measured intensities ''I'' which can be directly related to the backscatter coefficient. A value of zero means that no energy is reflected from the surface whereas a value of one means the total reflection.

A useful representation of the intensity is the logarithmic transformation to the [[http://en.wikipedia.org/wiki/Decibel| decibel unit]]
{{{#!python
img_dB=10*log10(img)
}}}

The PDF of the image above can be estimated from the number of occurrence of grey levels in the ''B'' intervals between ''q'' and ''q+dq'' and displayed using
{{{#!python
hist(img,bins=50)
}}}
with the resolution ''B'' for 50 ''bins'' (intervals).

{{attachment:seaice_histogram.png}}

[[/ExampleCode|Source code and data]]

== Estimation of the PDF parameters ==

With assumptions about the statistical processes it is possible to model the probability distribution. In the following example it is assumed that the image is comprised of two surface types. The gamma distribution describes the statistical process that leads to the noise in radar images. A linear superposition of two gamma functions is used as a model for the PDF. A least square optimization is used with a cost/error function to estimate the model shape parameters (Line 9,11) and the linear superpostion coefficients (Line 14). The fitted model parameters {{{p[6]}}} and {{{p[7]}}} yield the fractions of the two surface types which are 47% and 53%, respectively

{{attachment:model_stat.png}}


{{{#!python
from scipy import *
from pylab import *
import scipy.optimize as opti
import scipy.stats as stats


def my_pdf(p,x):
    """A linear mixture of two gamma PDF"""
    pdf1=stats.gamma.pdf(x,p[0],loc=p[1],scale=p[2])
    pdf2=stats.gamma.pdf(x,p[3],loc=p[4],scale=p[5])
    return pdf1*p[6]+pdf2*p[7]

def my_cdf(p,x):
    """A linear mixture of two gamma CDF"""
    cdf1=stats.gamma.cdf(x,p[0],loc=p[1],scale=p[2])
    cdf2=stats.gamma.cdf(x,p[3],loc=p[4],scale=p[5])
    return cdf1*p[6]+cdf2*p[7]

def my_cost(p,x,y):
    """The cost (error) function"""
    f=my_pdf(p,x)
    return y-f

#Read data
img=reshape(fromfile('ASAR_seaice_mixed_20080421_f32_1000x1000.dat',dtype=float32),(1000,1000))

# Calculate relative frequency of occurence
MIN,MAX,N=0.0,0.3,400
h=histogram(img,bins=N,range=[MIN,MAX],normed=True)
y,x=h[0],h[1]
# Cumulative frequency (normalized)
cdf=y.cumsum()*(MAX-MIN)/N

# Fit PDF model parameters to data by least squares optimization
p0=[7.0,0.003,0.005,7.0,0.1,0.005,0.5,0.5 ] # initial guess
p,success = opti.leastsq(my_cost, p0[:], args = (x, y))


figure(1)
subplot(2,1,1)
plot(x,y)
plot(x,my_pdf(p,x))
legend(('Observed PDF','Model PDF'),'upper right')
xlabel('Backscatter coefficient')
ylabel('Probability')
subplot(2,1,2)
plot(x,cdf)
plot(x,my_cdf(p,x))
legend(('Observed CDF','Model CDF'),'lower right')
xlabel('Backscatter coefficient')
ylabel('Cumulative probability')
savefig('model_stat.png',dpi=100)
show()
}}}