How to calculate median and standard deviation from histogram?

Question

I have a 256 bin histogram of an 8 bit image.

Is it possible to calculate the Median, Variance and Standard deviation of the image from the given histogram?

If yes, then what is the procedure?

Thank you.

Answer 1

You can't calculate any of them exactly because all you have is the interval of values that they belong to and not their exact values. It is the mode and not the median that is in the tallest bin. You can determine which bin the median is in and thus know the two end points of its bin are values that it falls betweem. To find where the median is you just total the number of data points in each bin starting from the left unit the get to the integer equal to (n-1)/2 when n is odd and (n+1)/2 if (n-1)/2 and (n+1)/2 are in the same bin. If (n+1)/2 is in a higher bin then (n-1)/2 you can't be sure which bin the median is in but you know it is near the boundary separating the two adjacent bins.

You can calculate grouped mean and grouped variances which may be rough approximations to the actual sample means and variances but not exact.

Answer 2

Since you have exactly all the pixel values and their occurrences, you could reorder them in an array in O(n) – with n the number of pixels – and simply pick the value in the middle of the array.

int* histogram(Image& img)
{
    int hist[256] = {0};
    int i,j;
    for (i=0; i<img.height; i++)
    {
        for (j=0; j<img.width; j++)
        {
            hist[img.at(i,j)]++;
        }
    }
}

double median(int* hist, int n)
{
    double med;
    int *values = (int*) malloc(n, sizeof(int));
    double position = double(n+1)/2.0;
    int _position_ = (int) floor(position);
    int i,j,k;
    // Order pixel values
    for (i=0; i<256; i++)
    {
        for (j=0; j<hist[i]; j++)
        {
            values[k++] = i;
        }
    }
    // Get median value
    if (position == (double)_position_)
    {
        // Odd number of pixels
        med = (double) values[_position_];
    }
    else
    {
        // Even number of pixels
        med = double(values[_position_] + values[_position_+1])/2.0;
    }
    free(values);
    return med;
}

int main(int argc, char* argv[])
{
    Image img;
    int *hist;
    double med;

    // Load image
    img = imread(argv[argc-1]);

    // Compute histogram
    hist = histogram(img);

    // Compute median
    med = median(hist, img.height*img.width);

    printf("The median value is: %f", med);
    return 0;
}

Answer 3

Adding to Michaels answer: You can compute grouped mean and std.dev, but they are approximations. While there are no systematic error in the approximation to the mean, the st.dev will be underestimated in this way (since you are calculating as if all observations falling in the same bin are in fact equal, so you miss the contribution to variance from variation within bins). Long time ago corrections was developed for this, based on the fact that the variance of a uniform distribution on $[a,b]$ is $\frac{ (b-a)^2}{12}$. They are Sheppards corrections: http://mathworld.wolfram.com/SheppardsCorrection.html

Answer 4

I think that if you normalize your histogram in the sense that the total area under that is $1$, then you can take the mode, the mean, and the median. So if you have a not-normalize histogram $\hat H$, where $\Delta x_k$ is the width of each interval and the height of each interval is $\hat f(x_k)$, the area $A(\hat H)$ can be define as

$$x_k^*\in\Delta x_k $$

$$A(\hat H):=\sum_{k=1}^n \Delta x_k\cdot \hat f(x_k^*)$$

for $n$ intervals in $\hat H$, and the normalize histogram $H$ is $\hat H/A(\hat H)$. Then, for giving a continuos meaning to $H$, you use the Heaviside's step function $u(x)$ define as a map $u:\Bbb R\to\{0,1\}$ and

$$u(x):= \begin{cases} 0 & x<0\\[2ex] 1 & x\ge0\end{cases}$$

Then $H$ can be defined as

$$\sup_x\Delta x_k=x_k^s \qquad \inf_x\Delta x_k=x_k^i \qquad x_{k+1}^i=x_k^s $$

$$H(x):=\sum_{k=1}^{n} f(x_k^*)(u(x-x_k^i)-u(x-x_k^s))$$

Where $f(x_k^*)$ is the normalize height of each interval. Now $A(H)$ can be calculated as

$$A(H)=\int_{-\infty}^\infty H(x)\ dx =\int_{-\infty}^\infty \sum_{k=1}^{n} f(x_k^*)(u(x-x_k^i)-u(x-x_k^s)) dx$$

The antiderivative of $u(x)$ is know as the ramp function defined as the map $R:\Bbb R\to[0,\infty)$ and:

$$R(x):= \begin{cases} 0 & x<0\\[2ex] x & x\ge0\end{cases}$$

$$\int u(x-a)\ dx=R(x-a) $$

So now $A(H)$ is

$$\int_{-\infty}^\infty \sum_{k=1}^{n} f(x_k^*)(u(x-x_k^i)-u(x-x_k^s))\ dx$$

$$=\left.\sum_{k=1}^{n} f(x_k^*)(R(x-x_k^i)-R(x-x_k^s))\right\rvert_{x=a}^b =1$$

where $[a,b]$ is the range of the random variable $X:\Omega\to [a,b]$. Given that $H(x)$ is a linear combination of step function, the mean or expected value of the random variable $\Bbb E[x]$

$$\Bbb E[X]=\int_{-\infty}^\infty xH(x)dx $$

is also a linear combination, but of the integral of the form

$$\int_{-\infty}^x t\cdot u(t-a)\ dt$$

Using integration by parts result in...

$$\int_{-\infty}^x t\cdot u(t-a)\ dt=xR(x-a)-\int_{-\infty}^x R(t-a)\ dt =xR(x-a)-{\left(R(x-a)\right)^2\over 2}$$

Also, note that the antiderivative just calculated is a parabola described by the function $p(x)$, which has the property that when $x=a$, $p(a)=0$, but it's not shifted because the term $xR(x-a)$ is not invariant under translation, so $p(x)$ is not the form $(x-a)^2/ 2$, rather...

$$p(x)={x^2\over 2}+c $$

So for $p(a)=0$, we have

$$p(a)=0={a^2\over 2}+c \Rightarrow c=-{a^2\over 2}$$

and the antiderivative of $x\cdot u(x-a)$ is

$$xR(x-a)-{\left(R(x-a)\right)^2\over 2}={x^2-a^2\over2} $$

Now we can calculate the expected value fairly easy...

$$\Bbb E[X]=\int_{-\infty}^\infty xH(x)dx=\int_{-\infty}^\infty x\left(\sum_{k=1}^{n} f(x_k^*)(u(x-x_k^i)-u(x-x_k^s))\right)\ dx$$

$$=\left.\sum_{k=1}^{n}f(x_k^*)\left({x^2-(x_k^i)^2\over 2}-{x^2-(x_k^s)^2\over 2}\right)\right\rvert^{x=b}$$

$$\sum_{k=1}^{n}f(x_k^*)\left((x_k^s)^2-(x_k^i)^2\over 2\right)={f(x_n^*)(x_n^s)^2-f(x_1^*)(x_1^i)^2\over 2}+\sum_{k=2}^{n}{(f(x_{k-1}^*)-f(x_k^*))(x_k^i)^2\over 2} $$

Defining the norm of the width as $|\Delta x_k|:=x^s_k-x^i_k$ and if $|\Delta x_k|=|\Delta x|\quad \forall k\in{1,2,\dots,n}$ and $|\Delta x|\in\Bbb R$ simplify the result to

$$\sum_{k=1}^{n}f(x_k^*)\left((x_k^s)^2-(x_k^i)^2\over 2\right)=\sum_{k=1}^{n}f(x_k^*)\left(2|\Delta x|(x^i_k)+|\Delta x|^2\over 2\right) $$

Note that taking the limit as $|\Delta x|\to 0$ gives the definition of the Riemann integral for the left of the expected value of a continuos random variable.

For the last values of interest, the mode should be $\sup H(x)$, the median is the solution of $\sum_{k=1}^{n} f(x_k^*)(R(x-x_k^i)-R(x-x_k^s))=0.5$.

For the variance, using the same approach of $\int xu(x-a)dx$ to

$$\int x^2u(x-a) = {x^3-a^3\over 3}$$

and making the same assumtions as before gives...

$$\Bbb E[X^2]=\int_{-\infty}^\infty x^2H(x)\ dx =\sum_{k=1}^{n}f(x_k^*)\left((x_k^s)^3-(x_k^i)^3\over 3\right)=\sum_{k=1}^{n}f(x_k^*)\left(3|\Delta x|^2(x^i_k)+3|\Delta x|(x_k^i)^2+|\Delta x|^3\over 3\right)$$

again, taking the limit as before gives the Riemann integral for the left of $\int x^2f_X(x)dx$