What jmch didn't say is, if sqrt( C'+A'-B'-D'/K - (mean*mean) )
is not how you compute the standard deviation from the integral image, then how you do it?
First, let me switch to Python / numpy code, so we get a modicum of notation consistency and expressions are easier to check. Given a sample array X, say:
X = array([random() * 10.0 for i in range(0, 9)])
The uncorrected sample standard deviation of X
can be defined as:
std = (sum((X - mean(X)) ** 2) / len(X)) ** 0.5 # 1
Applying the binomial theorem to (X - mean(X)) ** 2
we get:
std = (sum(X ** 2 - X * 2 * mean(X) + mean(X) ** 2) / len(X)) ** 0.5 # 2
Given the identities of the summation operation, we can make:
std = ((sum(X ** 2) - 2 * mean(X) * sum(X) + len(X) * mean(X) ** 2) / len(X)) ** 0.5 # 3
If we make S = sum(X)
, S2 = sum(X ** 2)
, M = mean(X)
and N = len(X)
we get:
std = ((S2 - 2 * M * S + N * M ** 2) / N) ** 0.5 # 4
Now for an image I
and two integral images P
and P2
calculated from I
(where P2
is the integral image for squared pixel values), we know that, given the four edge coordinates A = (i0, j0)
, B = (i0, j1)
, C = (i1, j0)
and D = (i1, j1)
, the values of S
, S2
, M
and N
can be calculated for the range I[A:D]
as:
S = P[A] + P[D] - P[B] - P[C]
S2 = P2[A] + P2[D] - P2[B] - P2[C]
N = (i1 - i0) * (j1 - j0)
M = S / N
Which can then be applied to equation (4) above yielding the standard deviation of the range I[A:D]
.
Edit: It's not entirely necessary, but given that M = S / N
we can apply the following substitutions and simplifications to equation (4):
std = ((S2 - 2 * M * S + N * M ** 2) / N) ** 0.5
std = ((S2 - 2 * (S / N) * S + N * (S / N) ** 2) / N) ** 0.5
std = ((S2 - 2 * ((S ** 2) / N) + (S ** 2 / N)) / N) ** 0.5
std = ((S2 - ((S ** 2) / N)) / N) ** 0.5
std = (S2 / N - (S / N) ** 2) ** 0.5 # 5
Which is quite close to the equation remi gave, actually.