diff --git a/datasets/noise_free_test_images/estimate_prnu.py b/datasets/noise_free_test_images/estimate_prnu.py
index 5ee7d91..a12089b 100644
--- a/datasets/noise_free_test_images/estimate_prnu.py
+++ b/datasets/noise_free_test_images/estimate_prnu.py
@@ -24,6 +24,7 @@ datasetPath = 'no_noise_images'
 # Note that contrarily to `datasets/fake/`, here we do not have images being Gaussian with `scale` `1` but actual images with pixel values between 0 and 255.
 # In addition to the range difference, note that the distribution in the first set of images was a Gaussian and here is very different and specific.
 PRNU_FACTOR = 0.01
+NOISE_FACTOR = 0.25
 IMAGE_SIZE_SHAPE = (469, 704)
 
 np.random.seed(0)
@@ -48,7 +49,20 @@ for imageName in os.listdir(datasetPath):
 
         imageWithoutPrnuNpArrayTiles = [imageWithoutPrnuNpArray[x : x + m, y : y + n] for x in range(0, imageWithoutPrnuNpArray.shape[0], m) for y in range(0, imageWithoutPrnuNpArray.shape[1], n)]
         for imageWithoutPrnuNpArrayTile in imageWithoutPrnuNpArrayTiles:
-            imageWithPrnuNpArray = imageWithoutPrnuNpArrayTile + prnuNpArray
+            imageNoise = randomGaussianImage(scale = 255 * NOISE_FACTOR, size = imageWithoutPrnuNpArrayTile.shape)
+            imageWithPrnuNpArray = imageWithoutPrnuNpArrayTile + prnuNpArray + imageNoise
+            #showImageWithMatplotlib(imageWithPrnuNpArray)
+            fig, axes = plt.subplots(1, 2)
+            fig.suptitle('Comparison of an image without and with Gaussian noise and PRNU')
+
+            axes[0].set_title('Image without Gaussian noise and PRNU')
+            axes[0].imshow(imageWithoutPrnuNpArray)
+
+            axes[1].set_title('Image with Gaussian noise and PRNU')
+            axes[1].imshow(imageWithPrnuNpArray)
+
+            plt.show()
+            break
             #assert all([isIn256Range(extreme) for extreme in [imageWithPrnuNpArray.max(), imageWithPrnuNpArray.min()]]), 'Adding the PRNU resulted in out of 256 bounds image'
             imageWithPrnuPil = toPilImage(imageWithPrnuNpArray)
             #imagePrnuEstimatePil = contextAdaptiveInterpolator(imageWithPrnuPil.load(), imageWithPrnuPil)
@@ -56,8 +70,9 @@ for imageName in os.listdir(datasetPath):
             imagePrnuEstimateNpArray = imageWithPrnuNpArray - denoise_tv_chambolle(imageWithPrnuNpArray, weight=0.2, channel_axis=-1)
 
             imagesPrnuEstimateNpArray += [imagePrnuEstimateNpArray]
+        break
 
 cameraPrnuEstimateNpArray = np.array(imagesPrnuEstimateNpArray).mean(axis = 0)
 #rms = rmsDiffNumpy(cameraPrnuEstimateNpArray, prnuNpArray, True)
-showImageWithMatplotlib(cameraPrnuEstimateNpArray, f'Camera PRNU estimate\nRMS with actual one: {rmsDiffNumpy(cameraPrnuEstimateNpArray, prnuNpArray):.4f} (normalized RMS: {rmsDiffNumpy(cameraPrnuEstimateNpArray, prnuNpArray, True):.4f})')
+#showImageWithMatplotlib(cameraPrnuEstimateNpArray, f'Camera PRNU estimate\nRMS with actual one: {rmsDiffNumpy(cameraPrnuEstimateNpArray, prnuNpArray):.4f} (normalized RMS: {rmsDiffNumpy(cameraPrnuEstimateNpArray, prnuNpArray, True):.4f})')