[SciPy-User] one-sided gauss fit -- or: how to estimate backgound noise ?
Thu Jul 15 17:13:28 CDT 2010
> In image analysis one is often faced with (often unknown) background
> levels (offset) + (Gaussian) background noise.
> The overall intensity histogram of the image is in fact often Gaussian
> (Bell shaped), but depending on how many (foreground) objects are
> present the histogram shows a positive tail of some sort.
> So, I just got the idea if there was a function (i.e. mathematical
> algorithm) that would allow to fit only the left half of a Gaussian
> bell curve to data points !?
> This would have to be done in a way that the center, the variance (or
> sigma) and the peak height are free fitting parameters.
For this task, I usually use some form of robust estimator for the
mean and std, which is designed to ignore noise in the tails.
Below I've pasted code that I use for an "minimum covariance
determinant" estimate, which is translated from some matlab code I
For large images, it's slow and you'll probably want to randomly
sample pixels to feed to the MCD estimator instead of using the entire
image. And there are probably many simpler, faster robust estimators
(like cutting off the tails, etc.) that are out there.
import scipy.stats as stats
"""unimcd(y, h) -> subset_mask
unimcd computes the MCD estimator of a univariate data set. This
given by the subset of h observations with smallest variance. The
estimate is then the mean of those h points, and the MCD scale
their standard deviation.
A boolean mask is returned indicating which elements of the input
in the MCD subset.
The MCD method was introduced in:
Rousseeuw, P.J. (1984), "Least Median of Squares Regression,"
Journal of the American Statistical Association, Vol. 79, pp.
The algorithm to compute the univariate MCD is described in
Rousseeuw, P.J., Leroy, A., (1988), "Robust Regression and Outlier
Detection," John Wiley, New York.
This function based on UNIMCD from LIBRA: the Matlab Library for
Analysis, available at: http://wis.kuleuven.be/stat/robust.html
y = numpy.asarray(y, dtype=float)
ncas = len(y)
length = ncas-h+1
if length <= 1:
return numpy.ones(len(y), dtype=bool)
indices = y.argsort()
y = y[indices]
ind = numpy.arange(length-1)
ay = numpy.empty(length)
ay = y[0:h].sum()
ay[1:] = y[ind+h] - y[ind]
ay = numpy.add.accumulate(ay)
sq = numpy.empty(length)
sq = (y[0:h]**2).sum() - ay2
sq[1:] = y[ind+h]**2 - y[ind]**2 + ay2[ind] - ay2[ind+1]
sq = numpy.add.accumulate(sq)
ii = numpy.where(sq==sqmin)
Hopt = indices[ii:ii+h]
ndup = len(ii)
slutn = ay[ii]
initmean=slutn[numpy.floor((ndup+1)/2 - 1)]/h
# calculating consistency factor
More information about the SciPy-User