Help generating non-Gaussian noise

Question

I'm trying to add non-Gaussian noise to an image. The purpose is to synthesize realistic camera noise. I am doing the following:

1. Obtain a noisy image sample by subtracting a low ISO photo from a high ISO of the same photo.

2. Find a flat 128x128 section of image with no edges by randomly cropping and checking if the standard deviation is small enough.

3. Zero mean the patch.

4. Normalise the patch by the patch standard deviation.

5. Compute PSD (power spectral density) by windowing and then using fast Fourier transform (fftn in Pytorch). Followed by squaring the frequency domain representation.

6. Generate Gaussian noise of same standard deviation as cropped noise patch.

7. Compute frequency representation and PSD of noise as above.

8. Do the following:

   H = gaussian_noise_psd / image_psd 
   output = image_frequency_domain * H
   

9. Lastly I convert back to spatial domain by reverse FFT and reversing the normalization step.

My output doesn't look quite right. The generated noise seems to be a pattern and doesn't look very natural. I'll put a sample below. I've added 0.5 to the output as the noise is zero mean.

Cropped area from noisy sample

Generated Gaussian noise of same std as above

Output of experiment

While the first two images have a std of 2 in this case, the output image has a much higher std of 7.5. You may need to open the images in full screen as they are quite small. How do I replicate the noise properly?

Answer 1

You can create noise with the same frequency-domain characteristic as your example noise image simply by multiplying a white noise frequency spectrum by the magnitude of the frequency spectrum you are trying to emulate. White noise is expected to have equal power at all frequencies (though this will not be exactly the case), so the multiplication will "copy the shape".

Setting the standard deviation of the noise can be done after the fact.

This is example code in MATLAB using DIPimage (disclosure: I'm an author). But I think it's simple enough that you can consider it pseudo-code and replicate it in any other language with any other image-processing library.

a = readim('https://i.sstatic.net/JaeTF.png');
a = a - mean(a);                     % must be zero mean
b = noise(newim(a),'gaussian', 1);   % create white Gaussian noise image
c = real(ift(abs(ft(a)) * ft(b)));   % multiply in frequency domain
c = c * (std(a) /  std(c));          % scale to get the right std. dev.

Here are images a, b and c, stretched for display:

---

PS: Note that there's a difference between "Gaussian noise" (noise created by sampling a Gaussian distribution) and "white noise" (where each value is uncorrelated to other values). You can do the above with other distributions too, try starting with uniform distributed noise instead of Gaussian, for example. b is white noise, c is not. But c is still Gaussian!