Make mean initialization more general in iterativeMean

Verified with:

```py
myIterativeMean = iterativeMean()
print(myIterativeMean.mean)
myIterativeMean.addElement(2)
print(myIterativeMean.mean)
myIterativeMean.addElement(3)
print(myIterativeMean.mean)
```

Note that cannot simplify as follows in the general case of using images as
`element`s for instance:

```diff
diff --git a/datasets/raise/utils.py b/datasets/raise/utils.py
index fba5523..a116edb 100644
--- a/datasets/raise/utils.py
+++ b/datasets/raise/utils.py
@@ -27,12 +27,9 @@ def denoise(imageNpArray, denoiserName):
     return imageDenoisedNpArray

 class iterativeMean():
-    mean = None
+    mean = 0
     numberOfElementsInMean = 0

     def addElement(self, element):
-        if self.mean is None:
-            self.mean = element
-        else:
-            self.mean = ((self.mean * self.numberOfElementsInMean) + element) / (self.numberOfElementsInMean + 1)
+        self.mean = ((self.mean * self.numberOfElementsInMean) + element) / (self.numberOfElementsInMean + 1)
         self.numberOfElementsInMean += 1
\ No newline at end of file
```

datasets/rafael/240424/plot_dates.py

@@ -0,0 +1,59 @@
+from datetime import datetime
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import matplotlib.dates as mdates
+
+from PIL import Image
+from PIL.ExifTags import TAGS
+
+import os
+
+os.chdir('photos')
+
+names = []
+dates = []
+
+for fileName in sorted(os.listdir()):
+    image = Image.open(fileName)
+    imageExif = image.getexif()
+    dateTimeKey = list(TAGS.keys())[list(TAGS.values()).index('DateTime')]
+    dateTime = imageExif[dateTimeKey]
+    names += [fileName.replace('DSC0', '').replace('.JPG', '')]
+    dates += [dateTime[:(-1 if fileName.endswith('.tif') else 0)]]
+
+# Convert date strings to datetime
+dates = [datetime.strptime(d, "%Y:%m:%d %H:%M:%S") for d in dates]
+
+# Choose some nice levels
+NUMBER_OF_LEVELS = 10
+actualLevels = range(-NUMBER_OF_LEVELS * 2 + 1, NUMBER_OF_LEVELS * 2 , 2)
+levels = np.tile(actualLevels,
+                 int(np.ceil(len(dates)/len(actualLevels))))[:len(dates)]
+
+# Create figure and plot a stem plot with the date
+fig, ax = plt.subplots(figsize=(8.8, 4), layout="constrained")
+ax.set(title="Matplotlib release dates")
+
+ax.vlines(dates, 0, levels, color="tab:red")  # The vertical stems.
+ax.plot(dates, np.zeros_like(dates), "-o",
+        color="k", markerfacecolor="w")  # Baseline and markers on it.
+
+# annotate lines
+for d, l, r in zip(dates, levels, names):
+    ax.annotate(r, xy=(d, l),
+                xytext=(-3, np.sign(l)*3), textcoords="offset points",
+                horizontalalignment="right",
+                verticalalignment="bottom" if l > 0 else "top")
+
+# format x-axis with 4-month intervals
+ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y:%m:%d %H:%M:%S"))
+plt.setp(ax.get_xticklabels(), rotation=30, ha="right")
+
+# remove y-axis and spines
+ax.yaxis.set_visible(False)
+ax.spines[["left", "top", "right"]].set_visible(False)
+
+ax.margins(y=0.1)
+plt.show()

datasets/raise/analyze_bayer_filter_mean_denoiser.py

@@ -0,0 +1,35 @@
+import numpy as np
+from utils import Color
+import os
+from tqdm import tqdm
+import rawpy
+import matplotlib.pyplot as plt
+
+os.chdir('flat-field/NEF')
+
+firstBayerFilterOccurrenceImages = []
+
+for fileName in tqdm(os.listdir()):
+    with rawpy.imread(fileName) as raw:
+        colorDesc = raw.color_desc.decode('ascii')
+        assert colorDesc == 'RGBG'
+        assert np.array_equal(raw.raw_pattern, np.array([[0, 1], [3, 2]], dtype = np.uint8))
+        # RG
+        # GB
+        firstBayerFilterOccurrenceImage = raw.raw_image_visible.copy()[:2, :2]
+        firstBayerFilterOccurrenceImages += [firstBayerFilterOccurrenceImage]
+
+firstBayerFilterOccurrenceImages = np.array(firstBayerFilterOccurrenceImages)
+print(rawImageVisible)
+
+NUMBER_OF_COLORS = 3
+HEX_COLOR = '#' + '%02x' * NUMBER_OF_COLORS
+COLOR_BASE = 256
+
+def getColor(colorIndex):
+    return HEX_COLOR % tuple((255 if colorIndex == colorIndexTmp else 0) for colorIndexTmp in range(NUMBER_OF_COLORS))
+
+for colorIndex, (colorY, colorX) in enumerate([(0, 0), (0, 1), (1, 1)]):
+    X = firstBayerFilterOccurrenceImages[:, colorY, colorX] - np.mean(firstBayerFilterOccurrenceImages, axis = 0)[colorY, colorX]
+    plt.hist(X, bins = len(set(X)), color = getColor(colorIndex), alpha = 0.3)
+plt.show()

datasets/raise/extract_noise.py

@@ -1,6 +1,5 @@
 #!/usr/bin/env python
 
-import skimage.restoration
 from skimage import img_as_float
 import numpy as np
 from PIL import Image
@@ -8,23 +7,19 @@ import os
 from tqdm import tqdm
 import csv
 import rawpy
-from utils import Color
+from utils import Color, denoise, iterativeMean
+import matplotlib.pyplot as plt
+from scipy.ndimage import gaussian_filter
 
 imagesFolderPath = 'rafael/arw'
 imagesFolderPathFileName = imagesFolderPath.replace('/', '_')
 # Among:
-# - `wavelet`
-# - `bilateral`
-# - `tv_chambolle`
-# - `mean`
+# `denoise` possible denoisers and `mean`.
 denoiser = 'mean'
 
 raiseNotFlatFields = False
-# `[Color.RED, Color.GREEN_RIGHT, ...]` or `Color`.
-colors = Color
-
-if denoiser != 'mean':
-    denoise = getattr(skimage.restoration, f'denoise_{denoiser}')
+# `[Color.RED, Color.GREEN_RIGHT, ...]` or `Color` or `[None]` for not raw images.
+colors = [None]
 
 imagesFileNames = os.listdir(imagesFolderPath + ('/png' if raiseNotFlatFields else ''))
 
@@ -79,83 +74,78 @@ def getImageNpArray(imageFileName, computeExtremes, color):
 
             match color:
                 case Color.RED:
-                    colorRawImageVisible = redRawImageVisible
+                    imageNpArray = redRawImageVisible
                 case Color.GREEN_RIGHT:
-                    colorRawImageVisible = greenRightRawImageVisible
+                    imageNpArray = greenRightRawImageVisible
                 case Color.GREEN_BOTTOM:
-                    colorRawImageVisible = greenBottomRawImageVisible
+                    imageNpArray = greenBottomRawImageVisible
                 case Color.BLUE:
-                    colorRawImageVisible = blueRawImageVisible
-
-            if computeExtremes:
-                colorRawImageVisibleMin = colorRawImageVisible.min()
-                colorRawImageVisibleMax = colorRawImageVisible.max()
-                if minColor is None or colorRawImageVisibleMin < minColor:
-                    minColor = colorRawImageVisibleMin
-                if maxColor is None or colorRawImageVisibleMax > maxColor:
-                    maxColor = colorRawImageVisibleMax
-                return
-            imageNpArray = (colorRawImageVisible - minColor) / (maxColor - minColor)
-            # Pay attention to range of values expected by the denoiser.
-            # Indeed if provide the thousands valued raw image, then the denoiser only returns values between 0 and 1 and making the difference between both look pointless.
+                    imageNpArray = blueRawImageVisible
     else:
         imagePil = Image.open(imageFilePath)
         imageNpArray = img_as_float(np.array(imagePil))
+
+    if computeExtremes:
+        colorRawImageVisibleMin = imageNpArray.min()
+        colorRawImageVisibleMax = imageNpArray.max()
+        if minColor is None or colorRawImageVisibleMin < minColor:
+            minColor = colorRawImageVisibleMin
+        if maxColor is None or colorRawImageVisibleMax > maxColor:
+            maxColor = colorRawImageVisibleMax
+        return
+
+    if (imageFileName.endswith('.NEF') or imageFileName.endswith('.ARW')) and denoiser != 'mean':
+        imageNpArray = (imageNpArray - minColor) / (maxColor - minColor)
+        # Pay attention to range of values expected by the denoiser.
+        # Indeed if provide the thousands valued raw image, then the denoiser only returns values between 0 and 1 and making the difference between both look pointless.
     return imageNpArray
 
 # `color` is the actual color to estimate PRNU with.
 def treatImage(imageFileName, computeExtremes = False, color = None):
-    global mean, numberOfImagesInMean
+    global estimatedPrnuIterativeMean
     imageNpArray = getImageNpArray(imageFileName, computeExtremes, color)
     if imageNpArray is None:
         return
-    match denoiser:
-        case 'wavelet':
-            imageDenoisedNpArray = denoise(imageNpArray, rescale_sigma=True)
-        case 'bilateral':
-            imageDenoisedNpArray = denoise(imageNpArray, sigma_color=0.05, sigma_spatial=15)
-        case 'tv_chambolle':
-            imageDenoisedNpArray = denoise(imageNpArray, weight=0.2)
-        case 'mean':
-            imageDenoisedNpArray = imageNpArray - means[color]
-    imageNoiseNpArray = imageNpArray - imageDenoisedNpArray
-    if mean is None:
-        mean = imageNoiseNpArray
+    if denoiser != 'mean':
+        imageDenoisedNpArray = denoise(imageNpArray, denoiser)
     else:
-        mean = ((mean * numberOfImagesInMean) + imageNoiseNpArray) / (numberOfImagesInMean + 1)
-    numberOfImagesInMean += 1
+        imageDenoisedNpArray = means[color]
+    imageNoiseNpArray = imageNpArray - imageDenoisedNpArray
+    estimatedPrnuIterativeMean.add(imageNoiseNpArray)
 
-# Assuming same intensity scale across color channels.
-for imageFileName in tqdm(imagesFileNames, 'Computing extremes of images'):
-    for color in colors:
-        treatImage(imageFileName, computeExtremes = True, color = color)
+if (minColor is None or maxColor is None) and denoiser != 'mean':
+    # Assuming same intensity scale across color channels.
+    for imageFileName in tqdm(imagesFileNames, 'Computing extremes of images'):
+        for color in colors:
+            treatImage(imageFileName, computeExtremes = True, color = color)
 
-# To skip this step next time.
-# Maybe thanks to `rawpy.RawPy` fields, possibly stating device maximal value, can avoid doing so to some extent.
-print(f'{minColor=}')
-print(f'{maxColor=}')
+    # To skip this step next time.
+    # Maybe thanks to `rawpy.RawPy` fields, possibly stating device maximal value, can avoid doing so to some extent.
+    print(f'{minColor=}')
+    print(f'{maxColor=}')
+
+def saveNpArray(fileName, npArray):
+    with open(f'{fileName}.npy', 'wb') as f:
+        np.save(f, npArray)
 
 if denoiser == 'mean':
     means = {}
-    for color in Color:
-        colorMean = None
-        numberOfImagesInColorMean = 0
-        for imageFileName in tqdm(imagesFileNames, f'Computing mean of {color} images'):
+    for color in colors:
+        colorIterativeMean = iterativeMean()
+        for imageFileName in tqdm(imagesFileNames, f'Computing mean of {color} colored images'):
             imageNpArray = getImageNpArray(imageFileName, False, color)
-            if colorMean is None:
-                colorMean = imageNpArray
-            else:
-                colorMean = ((colorMean * numberOfImagesInColorMean) + imageNpArray) / (numberOfImagesInColorMean + 1)
-            numberOfImagesInColorMean += 1
-        means[color] = colorMean
+            imageNpArray = gaussian_filter(imageNpArray, sigma = 5)
+            colorIterativeMean.add(imageNpArray)
+        means[color] = colorIterativeMean.mean
+        fileName = f'mean_{imagesFolderPathFileName}_{color}'
+        # Then use `merge_single_color_channel_images_according_to_bayer_filter.py` to consider all color channels, instead of saving this single color channel as an image.
+        saveNpArray(fileName, colorIterativeMean.mean)
 
 for color in colors:
-    mean = None
-    numberOfImagesInMean = 0
+    estimatedPrnuIterativeMean = iterativeMean()
 
     for imageFileName in tqdm(imagesFileNames, f'Denoising images for color {color}'):
         treatImage(imageFileName, color = color)
 
-    npArrayFilePath = f'mean_{imagesFolderPathFileName}_{denoiser}_{color}.npy'
-    with open(npArrayFilePath, 'wb') as f:
-        np.save(f, mean)
+    npArrayFilePath = f'mean_{imagesFolderPathFileName}_{denoiser}_{color}'
+    saveNpArray(npArrayFilePath, estimatedPrnuIterativeMean.mean)

datasets/raise/flat-field_pointer.py

@@ -0,0 +1,92 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.backend_bases import MouseButton
+import matplotlib.image as mpimg
+import os
+from tqdm import tqdm
+
+# `Zoom to rectangle` shortcut is `o`.
+
+os.chdir('flat-field/TIF')
+
+fileNames = sorted(os.listdir())[41:]
+fileNameIndex = 0
+xys = []
+progressBar = tqdm(total = len(fileNames))
+
+def displayImage():
+    global fileNameIndex
+    if len(fileNames) > fileNameIndex:
+        fileName = fileNames[fileNameIndex]
+        image = mpimg.imread(fileName)
+
+        fig, ax = plt.subplots()
+
+        ax.imshow(image)
+        fileNameIndex += 1
+        plt.connect('button_press_event', onClick)
+
+        mng = plt.get_current_fig_manager()
+        mng.full_screen_toggle()
+
+        fig.canvas.toolbar.zoom()
+        fig.show()
+
+def onClick(event):
+    global xys
+    button = event.button
+    if button is MouseButton.RIGHT:
+        xy = [event.xdata, event.ydata]
+        xys += [xy]
+    if button is MouseButton.MIDDLE or button is MouseButton.RIGHT:
+        if button is MouseButton.MIDDLE:
+            print(f'Skipped {fileName}')
+        plt.close()
+        displayImage()
+    progressBar.update(1)
+
+displayImage()
+
+##
+
+xys = [[3061.83568617998, 842.2814890347822], [3048.9553647053262, 891.7951806771933], [3109.6923734795832, 878.8472318972372], [3095.1992301129835, 859.646383417094], [3044.950680354029, 830.531447782855], [3034.039239345914, 775.7097977790371], [3034.0562409262707, 753.1780473232427], [2991.7499740162116, 742.1417475753151], [3021.5541233968024, 737.7119154247914], [3022.412107608028, 703.0028101114289], [3094.4565209481857, 680.7339885543926], [3101.840188177543, 692.6886985125186], [3150.09819940581, 686.3913997808281], [3148.4955568040136, 694.3747788064451], [3165.156239230552, 703.9346483213995], [3158.7235775979357, 710.9522942666542], [3156.5302112965014, 697.5993396871294], [3106.8119050930513, 698.4561327048264], [3151.4728367372877, 694.3380683877142], [3137.283014559081, 662.9543422436452], [3008.1987322850605, 692.0180565561739], [3036.647953172549, 705.1188029786949], [2988.9166660643627, 684.0851408200217], [3001.5070402052334, 684.1944057536488], [2967.282534077184, 680.6072888791263], [2983.5140619626286, 693.469904886339], [2979.972210442175, 702.9033995969894], [2974.7535915953795, 709.6086279669086], [2972.55952140686, 700.4689248964266], [2967.1769510144627, 703.4861098128929], [2968.4714535085923, 708.5924315970815], [2978.725084963369, 725.0064286729727], [2992.7472737250696, 718.9682686788141], [2999.106406229902, 713.0463041988097], [2981.5746364633296, 692.5739107725338], [3014.22885383495, 675.1638989177214], [3004.351989366563, 635.6738446427469], [2964.9109157636794, 580.2416938434527], [2976.672441933737, 599.3922555819147], [2972.9623072548534, 578.6218471929612], [2694.5865522760287, 972.5526873692775], [2924.1211630176217, 1115.9816839025543], [2910.2861978522596, 1095.5172128924614], [2917.4530137143342, 1067.9680470058072], [2919.3469031337613, 1068.44539073963], [2959.4757277956496, 1105.6227950624348], [2936.370033263817, 1157.9142039933472], [2965.3568614486203, 1145.3186319170795], [2951.1991473505136, 1070.8989245366433], [2946.108094501549, 1045.891549058405], [2947.19213475732, 1074.4303801863052], [2943.690837961984, 1071.4536958497956], [3021.083385372544, 1116.881810271317], [3051.7205580454297, 1103.0992679894152], [3107.0873956690143, 972.8919609441568], [2861.932349690137, 1131.4847600231471], [2789.583044951935, 964.6512841164608]]
+
+x, y = [[xy[dimensionIndex] for xy in xys] for dimensionIndex in range(2)]
+
+plt.title('Center wall marker locations')
+
+plt.scatter(x, y, c = range(len(x)))
+for xyIndex, xy in enumerate(xys):
+    plt.text(xy[0], xy[1] + 10, str(xyIndex), horizontalalignment = 'center')
+
+plt.show()
+
+##
+
+greatestDistances = []
+
+#for xy in xys:
+for xy, otherXy in zip(xys, xys[1:]):
+    greatestDistance = None
+    #for otherXy in xys:
+    if xy != otherXy:
+        distance = sum([(xy[dimensionIndex] - otherXy[dimensionIndex]) ** 2 for dimensionIndex in range(2)]) ** 0.5
+        if greatestDistance is None or distance > greatestDistance:
+            greatestDistance = distance
+    greatestDistances += [greatestDistance]
+
+fig, ax = plt.subplots()
+plt.title('Distance between a wall marker location and the next one')
+ax.boxplot(greatestDistances)
+
+ys = []
+
+for percentile in [25, 50, 75]:
+    y = np.percentile(greatestDistances, percentile)
+    ys += [y]
+    ax.axhline(y = y)
+
+ax.set_yticks(list(ax.get_yticks())[1:-1] + ys)
+
+ax.set_xticks([], [])
+fig.show()

datasets/raise/generate_histogram.py

@@ -0,0 +1,41 @@
+import os
+from PIL import Image
+import matplotlib.pyplot as plt
+import numpy as np
+from tqdm import tqdm
+
+os.chdir('flat-field/TIF')
+
+NUMBER_OF_COLORS = 3
+HEX_COLOR = '#' + '%02x' * NUMBER_OF_COLORS
+COLOR_BASE = 256
+
+def getColor(colorIntensity, colorIndex):
+    return HEX_COLOR % tuple((colorIntensity if colorIndex == colorIndexTmp else 0) for colorIndexTmp in range(NUMBER_OF_COLORS))
+
+def getHistogram(fileName):
+    image = Image.open(fileName)
+    #image = image.crop((1, 1, 2, 2))
+    #print(image.size)
+
+    histogram = image.histogram()
+
+    colors = [histogram[COLOR_BASE * colorIndex:COLOR_BASE * (colorIndex + 1)] for colorIndex in range(NUMBER_OF_COLORS)]
+    return colors
+
+def plotHistogram(colors):
+    for colorIndex, color in enumerate(colors):
+        for colorIntensity in range(COLOR_BASE):
+            plt.bar(colorIntensity, color[colorIntensity], color = getColor(colorIntensity, colorIndex), alpha = 0.3)
+
+fileNameColors = []
+
+for fileName in tqdm(os.listdir()):
+    colors = getHistogram(fileName)
+    fileNameColors += [colors]
+
+meanFileNameColors = np.mean(fileNameColors, axis = 0)
+
+plotHistogram(meanFileNameColors)
+
+plt.show()

datasets/raise/merge_single_color_channel_images_according_to_bayer_filter.py

@@ -41,4 +41,4 @@ for color in tqdm(Color, 'Color'):
 
             multipleColorsImage[newX, newY] = pixel
 
-plt.imsave(PREFIX + 'multiple_colors.png', multipleColorsImage)
+plt.imsave(PREFIX + 'multiple_colors.png', multipleColorsImage)

datasets/raise/show_mean_noise.py

@@ -3,7 +3,8 @@ import matplotlib.pyplot as plt
 
 fileName = 'mean_flat-field_nef_wavelet_blue'
 npArray = np.load(f'{fileName}.npy')
-print(f'{npArrayMin=}')
+# For other than raw images:
+#npArray = (npArray - npArray.min()) / (npArray.max() - npArray.min())
 plt.imsave(f'{fileName}.png', npArray)
 #plt.imshow(npArray)
 #plt.show()

datasets/raise/split_and_compare_prnus_of_subgroups.py

@@ -2,19 +2,16 @@ from PIL import Image
 import numpy as np
 import matplotlib.pyplot as plt
 from tqdm import tqdm
-from skimage.restoration import denoise_tv_chambolle
+from utils import denoise
 from skimage import img_as_float
 import sys
 
-sys.path.insert(0, '../../algorithms/image_utils/')
-
-from image_utils import showImageWithMatplotlib, toPilImage
-
 sys.path.insert(0, '../../algorithms/distance/')
 
 from rms_diff import rmsDiffNumpy
 
 NUMBER_OF_SUBGROUPS = 1
+DENOISER = 'wavelet'
 
 IMAGES_FOLDER = 'flat-field/TIF'
 imagesFileNames = os.listdir(IMAGES_FOLDER)
@@ -31,7 +28,7 @@ for subgroupIndex in range(NUMBER_OF_SUBGROUPS):
         imagePath = f'{IMAGES_FOLDER}/{imageFileName}'
         imagePil = Image.open(imagePath)
         imageNpArray = img_as_float(np.array(imagePil))
-        imagePrnuEstimateNpArray = imageNpArray - denoise_tv_chambolle(imageNpArray, weight=0.2, channel_axis=-1)
+        imagePrnuEstimateNpArray = imageNpArray - denoise(imageNpArray, DENOISER)
         imagesPrnuEstimateNpArray += [imagePrnuEstimateNpArray]
 
     subgroupPrnuEstimateNpArray = []
@@ -46,7 +43,7 @@ for numberOfImagesIndex, numberOfImages in enumerate(numberOfImagesThresholds):
     rms = rmsDiffNumpy(subgroupsPrnuEstimatesNpArray[0][numberOfImagesIndex], subgroupsPrnuEstimatesNpArray[1][numberOfImagesIndex])
     rmss += [rms]
 
-plt.title('RMS between both subgroups estimated PRNUs for a given number of images among them')
+plt.title(f'RMS between both subgroups estimated PRNUs with {DENOISER} denoiser for a given number of images among them')
 plt.xlabel('Number of images of each subgroup')
 plt.ylabel('RMS between both subgroups estimated PRNUs')
 plt.plot(rmss)

datasets/raise/utils.py

@@ -1,4 +1,5 @@
 from enum import Enum, auto
+import skimage.restoration
 
 class Color(Enum):
     RED = auto()
@@ -8,3 +9,27 @@ class Color(Enum):
 
     def __str__(self):
         return self.name.lower()
+
+# Among:
+# - `wavelet`
+# - `bilateral`
+# - `tv_chambolle`
+def denoise(imageNpArray, denoiserName):
+    skImageRestorationDenoise = getattr(skimage.restoration, f'denoise_{denoiserName}')
+
+    match denoiserName:
+        case 'wavelet':
+            imageDenoisedNpArray = skImageRestorationDenoise(imageNpArray, rescale_sigma=True)
+        case 'bilateral':
+            imageDenoisedNpArray = skImageRestorationDenoise(imageNpArray, sigma_color=0.05, sigma_spatial=15)
+        case 'tv_chambolle':
+            imageDenoisedNpArray = skImageRestorationDenoise(imageNpArray, weight=0.2)
+    return imageDenoisedNpArray
+
+class iterativeMean:
+    mean = None
+    numberOfElementsInMean = 0
+
+    def add(self, element):
+        self.mean = ((self.mean * self.numberOfElementsInMean) + element) / (self.numberOfElementsInMean + 1)
+        self.numberOfElementsInMean += 1