How to calculate mean and standard deviation of a set of images

Question

I would like to know I to calculate the mean and the std of a given dataset of RGB images.
For example, with imagenet we have imagenet_stats: ([0.485, 0.456, 0.406], [0.229, 0.224, 0.225].
I tried:

rgb_values = [np.mean(Image.open(img).getdata(), axis=0)/255 for img in imgs_path]
np.mean(rgb_values, axis=0)
np.std(rgb_values, axis=0)

I am not sure that the values I get are correct.
Which could be a better implementation?

Answer 1

Two solutions:

• The first solution iterates over the images. It is MUCH slower than the second solution, and it uses the same amount of memory because it first loads and then stores all the images in a list. So it is strictly worse than the second solution, unless you will change how your images are loaded - load and process them one by one from disc.

• The second solution needs to hold all images in memory at the same time. It is MUCH faster, because it is fully vectorized.

---

First solution (iterating over the images):

For each channel: R, G, B, here is how to calculate the means and stds of all the pixels in all the images:

Requirement:

Each image has the same number of pixels.

If this is not the case - use the second solution (below).

images_rgb = [np.array(Image.open(img).getdata()) / 255. for img in imgs_path]
# Each image_rgb is of shape (n, 3), 
# where n is the number of pixels in each image,
# and 3 are the channels: R, G, B.

means = []
for image_rgb in images_rgb:
    means.append(np.mean(image_rgb, axis=0))
mu_rgb = np.mean(means, axis=0)  # mu_rgb.shape == (3,)

variances = []
for image_rgb in images_rgb:
    var = np.mean((image_rgb - mu_rgb) ** 2, axis=0)
    variances.append(var)
std_rgb = np.sqrt(np.mean(variances, axis=0))  # std_rgb.shape == (3,)

Proof

... that the mean and std will be same if calculated like this, and if calculated using all pixels at once:

Let's say each image has n pixels (with values vals_i), and there are m images.

Then there are (n*m) pixels.

The real_mean of all pixels in all vals_is is:

total_sum = sum(vals_1) + sum(vals_2) + ... + sum(vals_m)
real_mean = total_sum / (n*m)

Adding up the means of each image individually:

sum_of_means = sum(vals_1) / m + sum(vals_2) / m + ... + sum(vals_m) / m
             = (sum(vals_1) + sum(vals_2) + ... + sum(vals_m)) / m

Now, what is the relationship between the real_mean and sum_of_means? - As you can see,

real_mean = sum_of_means / n

---

Analogously, using the formula for standard deviation, the real_std of all pixels in all vals_is is:

sum_of_square_diffs =  sum(vals_1 - real_mean) ** 2
                     + sum(vals_2 - real_mean) ** 2
                     + ... 
                     + sum(vals_m - real_mean) ** 2
real_std = sqrt( total_sum / (n*m) )

If you look at this equation from another angle, you can see that real_std is basically the average of average variances of n values in m images.

---

Verification

Real mean and std:

rng = np.random.default_rng(0)
vals = rng.integers(1, 100, size=100)  # data

mu = np.mean(vals)
print(mu)
print(np.std(vals))

50.93                 # real mean
28.048976808432776    # real standard deviation

Comparing it to the image-by-image approach:

n_images = 10

means = []
for subset in np.split(vals, n_images):
    means.append(np.mean(subset))
new_mu = np.mean(means)

variances = []
for subset in np.split(vals, n_images):
    var = np.mean((subset - mu) ** 2)
    variances.append(var)

print(new_mu)
print(np.sqrt(np.mean(variances)))

50.92999999999999     # calculated mean
28.048976808432784    # calculated standard deviation

---

Second solution (fully vectorized):

Using all the pixels of all images at once.

rgb_values = np.concatenate(
    [Image.open(img).getdata() for img in imgs_path], 
    axis=0
) / 255.

# rgb_values.shape == (n, 3), 
# where n is the total number of pixels in all images, 
# and 3 are the 3 channels: R, G, B.

# Each value is in the interval [0; 1]

mu_rgb = np.mean(rgb_values, axis=0)  # mu_rgb.shape == (3,)
std_rgb = np.std(rgb_values, axis=0)  # std_rgb.shape == (3,)