1 files changed, 123 insertions, 0 deletions
diff --git a/utils/raspberrypi/ctt/ctt_noise.py b/utils/raspberrypi/ctt/ctt_noise.py
new file mode 100644
index 00000000..b84cf0ca
--- /dev/null
+++ b/utils/raspberrypi/ctt/ctt_noise.py
@@ -0,0 +1,123 @@
+# SPDX-License-Identifier: BSD-2-Clause
+#
+# Copyright (C) 2019, Raspberry Pi (Trading) Limited
+#
+# ctt_noise.py - camera tuning tool noise calibration
+
+from ctt_image_load import *
+import matplotlib.pyplot as plt
+
+"""
+Find noise standard deviation and fit to model:
+
+    noise std = a + b*sqrt(pixel mean)
+"""
+def noise(Cam,Img,plot):
+    Cam.log += '\nProcessing image: {}'.format(Img.name)
+    stds = []
+    means = []
+    """
+    iterate through macbeth square patches
+    """
+    for ch_patches in Img.patches:
+        for patch in ch_patches:
+            """
+            renormalise patch
+            """
+            patch = np.array(patch)
+            patch = (patch-Img.blacklevel_16)/Img.againQ8_norm
+            std = np.std(patch)
+            mean = np.mean(patch)
+            stds.append(std)
+            means.append(mean)
+
+    """
+    clean data and ensure all means are above 0
+    """
+    stds = np.array(stds)
+    means = np.array(means)
+    means = np.clip(np.array(means),0,None)
+    sq_means = np.sqrt(means)
+
+
+    """
+    least squares fit model
+    """
+    fit = np.polyfit(sq_means,stds,1)
+    Cam.log += '\nBlack level = {}'.format(Img.blacklevel_16)
+    Cam.log += '\nNoise profile: offset = {}'.format(int(fit[1]))
+    Cam.log += ' slope = {:.3f}'.format(fit[0])
+    """
+    remove any values further than std from the fit
+
+    anomalies most likely caused by:
+    > ucharacteristically noisy white patch
+    > saturation in the white patch
+    """
+    fit_score = np.abs(stds - fit[0]*sq_means - fit[1])
+    fit_std = np.std(stds)
+    fit_score_norm = fit_score - fit_std
+    anom_ind = np.where(fit_score_norm > 1)
+    fit_score_norm.sort()
+    sq_means_clean = np.delete(sq_means,anom_ind)
+    stds_clean = np.delete(stds,anom_ind)
+    removed = len(stds) - len(stds_clean)
+    if removed != 0:
+        Cam.log += '\nIdentified and removed {} anomalies.'.format(removed)
+        Cam.log += '\nRecalculating fit'
+        """
+        recalculate fit with outliers removed
+        """
+        fit = np.polyfit(sq_means_clean,stds_clean,1)
+        Cam.log += '\nNoise profile: offset = {}'.format(int(fit[1]))
+        Cam.log += ' slope = {:.3f}'.format(fit[0])
+
+    """
+    if fit const is < 0 then force through 0 by
+    dividing by sq_means and fitting poly order 0
+    """
+    corrected = 0
+    if fit[1] < 0:
+        corrected = 1
+        ones = np.ones(len(means))
+        y_data = stds/sq_means
+        fit2 = np.polyfit(ones,y_data,0)
+        Cam.log += '\nOffset below zero. Fit recalculated with zero offset'
+        Cam.log += '\nNoise profile: offset = 0'
+        Cam.log += ' slope = {:.3f}'.format(fit2[0])
+        # print('new fit')
+        # print(fit2)
+
+    """
+    plot fit for debug
+    """
+    if plot:
+        x = np.arange(sq_means.max()//0.88)
+        fit_plot = x*fit[0] + fit[1]
+        plt.scatter(sq_means,stds,label='data',color='blue')
+        plt.scatter(sq_means[anom_ind],stds[anom_ind],color='orange',label='anomalies')
+        plt.plot(x,fit_plot,label='fit',color='red',ls=':')
+        if fit[1] < 0:
+            fit_plot_2 = x*fit2[0]
+            plt.plot(x,fit_plot_2,label='fit 0 intercept',color='green',ls='--')
+        plt.plot(0,0)
+        plt.title('Noise Plot\nImg: {}'.format(Img.str))
+        plt.legend(loc = 'upper left')
+        plt.xlabel('Sqrt Pixel Value')
+        plt.ylabel('Noise Standard Deviation')
+        plt.grid()
+        plt.show()
+    """
+    End of plotting code
+    """
+
+    """
+    format output to include forced 0 constant
+    """
+    Cam.log += '\n'
+    if corrected:
+        fit = [fit2[0],0]
+        return fit
+
+    else:
+        return fit