diff --git a/algorithms/image_utils/image_utils.py b/algorithms/image_utils/image_utils.py
new file mode 100644
index 0000000..4316f95
--- /dev/null
+++ b/algorithms/image_utils/image_utils.py
@@ -0,0 +1,10 @@
+from PIL import Image
+import numpy as np
+from matplotlib import pyplot as plt
+
+def randomGaussianImage(scale, size):
+    return np.random.normal(loc = 0, scale = scale, size = size)
+
+def showImageWithMatplotlib(npArray):
+    plt.imshow(npArray)
+    plt.show()
\ No newline at end of file
diff --git a/datasets/fake/generate_dataset.py b/datasets/fake/generate_dataset.py
index eac391d..1e94df3 100644
--- a/datasets/fake/generate_dataset.py
+++ b/datasets/fake/generate_dataset.py
@@ -10,6 +10,10 @@ from rms_diff import rmsDiffPil, rmsDiffNumpy
 sys.path.insert(0, '../../algorithms/context_adaptive_interpolator/')
 
 from context_adaptive_interpolator import contextAdaptiveInterpolator
+
+sys.path.insert(0, '../../algorithms/image_utils/')
+
+from image_utils import showImageWithMatplotlib, randomGaussianImage
 from tqdm import tqdm
 
 IMAGE_SIZE = 64
@@ -21,13 +25,9 @@ PRNU_FACTOR = 0.1
 IMAGE_SIZE_SHAPE = (IMAGE_SIZE, IMAGE_SIZE)
 
 # Generate PRNUs and images of phones.
-# Is such `np.maximum` probabilistically correct with our theoretical method? See #19.
-def randomGaussianImage(scale):
-    return np.random.normal(loc = 0, scale = scale, size = IMAGE_SIZE_SHAPE)
+imagesWithoutPrnu = [[randomGaussianImage(scale = 1, size = IMAGE_SIZE_SHAPE) for _ in range(NUMBER_OF_IMAGES_PER_PHONE)] for phoneIndex in range(NUMBER_OF_PHONES)]
 
-imagesWithoutPrnu = [[randomGaussianImage(scale = 1) for _ in range(NUMBER_OF_IMAGES_PER_PHONE)] for phoneIndex in range(NUMBER_OF_PHONES)]
-
-prnus = [randomGaussianImage(scale = PRNU_FACTOR) for _ in range(NUMBER_OF_PHONES)]
+prnus = [randomGaussianImage(scale = PRNU_FACTOR, size = IMAGE_SIZE_SHAPE) for _ in range(NUMBER_OF_PHONES)]
 
 imagesWithPrnu = [[imageWithoutPrnu + prnus[phoneIndex] for imageWithoutPrnu in imagesWithoutPrnu[phoneIndex]] for phoneIndex in range(NUMBER_OF_PHONES)]
 
@@ -41,10 +41,6 @@ def showImageWithPil(npArray):
     npArray = (npArray / npArray.max()) * 255
     Image.fromarray(npArray).show()
 
-def showImageWithMatplotlib(npArray):
-    plt.imshow(npArray)
-    plt.show()
-
 plt.title('RMS between actual PRNU and the mean of the first $N$ images with PRNU (i.e. estimated PRNU)')
 plt.xlabel('$N$ first images with PRNU')
 plt.ylabel('RMS')