Python autocorr

SpeckleSizeCode/PythonAutocorr/autocorr.py

@@ -0,0 +1,108 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import tifffile
+import os
+from scipy.ndimage import gaussian_filter
+import cv2
+
+nb = 13
+fname = rf"20190924-200ms_20mW_Ave15_Gray_10X0.4_{nb}.tif"
+
+p = os.path.dirname(os.path.join(os.getcwd(), "..", ".."))
+path = os.path.join(p, "MATLAB", fname)
+
+prefix = "circularWithPhasesSimulations"
+otherPrefix = "50sims"
+path = "4pixelsCircularWithPhasesSimulations.tiff"
+args = [prefix, path]
+#args = [prefix, otherPrefix, path]
+path = os.path.join(*args)
+
+if path.endswith("tif") or path.endswith("tiff"):
+    img = tifffile.imread(path)
+else:
+    img = cv2.imread(path)
+
+plt.imshow(img)
+plt.show()
+
+imgDivide = gaussian_filter(img, 75)
+
+img = img / imgDivide - np.mean(img)
+plt.imshow(img)
+plt.show()
+
+fft = np.fft.fft2(img)
+ifft = np.fft.ifftshift(np.fft.ifft2(np.abs(fft) ** 2)).real
+ifft /= np.size(ifft)
+
+r = (ifft - np.mean(img) ** 2) / np.var(img)
+plt.imshow(r)
+plt.colorbar()
+plt.show()
+
+xSlice, ySlice = r.shape[1] // 2, r.shape[0] // 2
+
+verticalSlice = r[:, xSlice]
+horizontalSlice = r[ySlice, :]
+
+
+def normalize(verticalSlice, horizontalSlice):
+    verticalSlice = verticalSlice - np.min(verticalSlice)
+    verticalSlice = verticalSlice / np.max(verticalSlice)
+
+    horizontalSlice = horizontalSlice - np.min(horizontalSlice)
+    horizontalSlice = horizontalSlice / np.max(horizontalSlice)
+    return verticalSlice, horizontalSlice
+
+
+plt.plot(verticalSlice)
+plt.title("Autocorrelation profile (x centered)")
+plt.xlabel("Vertical position y [px]")
+plt.ylabel("Normalized correlation coefficient [-]")
+plt.show()
+
+plt.plot(horizontalSlice)
+plt.title("Autocorrelation profile (y centered)")
+plt.xlabel("Horizontal position x [px]")
+plt.ylabel("Normalized correlation coefficient [-]")
+plt.show()
+
+relativeErrorHalfMax = 25
+
+# Y FWHM
+print("Y FWHM")
+halfMax = 0.5
+inferiorBound = halfMax - halfMax * relativeErrorHalfMax / 100
+superiorBound = halfMax + halfMax * relativeErrorHalfMax / 100
+
+# Range of values to compute FWHM
+pointsForFW = np.where((verticalSlice >= inferiorBound) & (verticalSlice <= superiorBound))[0]
+print(pointsForFW)
+middlePoint = np.argmax(verticalSlice)  # Find the middle point to separate the values at the left and at the right
+left = np.mean(pointsForFW[pointsForFW < middlePoint])
+right = np.mean(pointsForFW[pointsForFW > middlePoint])
+print(f"left: {left}")
+print(f"right: {right}")
+FWHM = right - left
+print("Diameter : ", FWHM)
+print("Radius : ", FWHM / 2)
+
+# X FWHM
+print("= = = = = = = = = = = = = = = = = = = = = = = =")
+print("X FWHM")
+halfMax = 0.5
+inferiorBound = halfMax - halfMax * relativeErrorHalfMax / 100
+superiorBound = halfMax + halfMax * relativeErrorHalfMax / 100
+
+# Range of values to compute FWHM
+pointsForFW = np.where((horizontalSlice >= inferiorBound) & (horizontalSlice <= superiorBound))[0]
+print(pointsForFW)
+middlePoint = np.argmax(horizontalSlice)  # Find the middle point to separate the values at the left and at the right
+left = np.mean(pointsForFW[pointsForFW < middlePoint])
+right = np.mean(pointsForFW[pointsForFW > middlePoint])
+print(f"left: {left}")
+print(f"right: {right}")
+FWHM = right - left
+print("Diameter : ", FWHM)
+print("Radius : ", FWHM / 2)

SpeckleSizeCode/PythonAutocorr/autocorrelation.py

@@ -0,0 +1,132 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import tifffile
+from scipy.ndimage import gaussian_filter, median_filter
+from scipy.signal import correlate2d
+import cv2
+
+
+class FileReader:
+    # FIXME: This is only until pydcclab is ok
+
+    @staticmethod
+    def readFile(path: str):
+        if path.endswith(".tif") or path.endswith(".tiff"):
+            pixels = tifffile.imread(path)
+        else:
+            pixels = cv2.imread(path)
+        return pixels.T  # we want shape[0] as the width and shape[1] as the height
+
+
+class Autocorrelation:
+    # FIXME: Use pydcclab when ok
+
+    def __init__(self, imagePath: str):
+        self.__image = FileReader.readFile(imagePath)
+        self.__original = self.image
+        self.__autocorrelation = None
+        self.__slicesObj = None
+
+    @property
+    def image(self):
+        return self.__image.copy()
+
+    @property
+    def autocorrelation(self):
+        if self.__autocorrelation is None:
+            return None
+        return self.__autocorrelation.copy()
+
+    def getSlices(self, indices: tuple = None):
+        if self.__slicesObj is None:
+            raise ValueError("Please compute the autocorrelation to access its slices.")
+        if indices is None:
+            return self.__slicesObj.middleSlices()
+        return self.__slicesObj.slicesAt(indices)
+
+    def computeAutocorrelation(self, gaussianFilterStdDev: float = 75, medianFilterSize: int = 3):
+        self._gaussianNormalization(gaussianFilterStdDev)  # First, gaussian filter normalization
+        self._medianFilter(medianFilterSize)  # Then, median filter to remove noise
+        self._autocorrelationWithFourierTransform()  # Compute the autocorrelation
+
+    def _gaussianNormalization(self, filterStdDev: float = 75):
+        filteredImage = gaussian_filter(self.__image, filterStdDev)
+        self.__image = self.__image / filteredImage - np.mean(self.__image)
+
+    def _autocorrelationWithFourierTransform(self):
+        fft = np.fft.fft2(self.__image)
+        ifft = np.fft.ifftshift(np.fft.ifft2(np.abs(fft) ** 2)).real
+        ifft /= np.size(ifft)
+        self.__autocorrelation = (ifft - np.mean(self.__image) ** 2) / np.var(self.__image)
+        self.__slicesObj = AutocorrelationSlices(self.__autocorrelation)
+
+
+    def _medianFilter(self, filterSize: int = 3):
+        self.__image = median_filter(self.__image, filterSize)
+
+    def showImage(self):
+        plt.imshow(self.__image)
+        plt.show()
+
+    def showAutocorrelation(self, showColorbar: bool = True):
+        if self.__autocorrelation is None:
+            raise ValueError("No autocorrelation computed.")
+        plt.imshow(self.__autocorrelation)
+        if showColorbar:
+            plt.colorbar()
+        plt.show()
+
+    def showAutocorrelationSlices(self, indices: tuple = None, showHorizontal: bool = True, showVertical: bool = True):
+        vSlice, hSlice = self.getSlices(indices)
+        if showHorizontal and showVertical:
+            fig, (ax1, ax2) = plt.subplots(2, sharey="col")
+            fig.suptitle("Autocorrelation slices")
+
+            ax1.plot(hSlice)
+            ax1.set_title(f"Horizontal slice (at index {indices[0]})")
+            ax1.set_xlabel("Horizontal position $x$ [pixel]")
+
+            ax2.plot(vSlice)
+            ax1.set_title(f"Vertical slice (at index {indices[1]})")
+            ax1.set_xlabel("Vertical position $y$ [pixel]")
+
+            ylabel = "Normalized autocorrelation coefficient [-]"
+            fig.text(0.06, 0.5, ylabel, ha='center', va='center', rotation='vertical')
+            plt.subplots_adjust(hspace=0.32)
+            plt.show()
+        elif showVertical:
+            plt.plot(vSlice)
+            plt.title(f"Vertical slice (at index {indices[1]})")
+            plt.xlabel("Vertical position $y$ [pixel]")
+            plt.ylabel("Normalized autocorrelation coefficient [-]")
+            plt.show()
+        elif showHorizontal:
+            plt.plot(hSlice)
+            plt.title(f"Horizontal slice (at index {indices[0]})")
+            plt.xlabel("Horizontal position $x$ [pixel]")
+            plt.ylabel("Normalized autocorrelation coefficient [-]")
+            plt.show()
+
+
+class AutocorrelationSlices:
+
+    def __init__(self, autocorrelation: np.ndarray):
+        if not isinstance(autocorrelation, np.ndarray):
+            raise TypeError("The autocorrelation parameter must be a numpy array.")
+        if not autocorrelation.ndim == 2:
+            raise ValueError("The autocorrelation must be in 2D.")
+        self.__autocorrelation = autocorrelation
+
+    def slicesAt(self, indices: tuple):
+        if len(indices) != 2:
+            raise ValueError("There must be 2 indices of slicing, one horizontal and one vertical.")
+        # Assumes indices[0] is the width and indices[1] is the height
+        xSlice, ySlice = indices[0], indices[1]
+        verticalSlice = self.__autocorrelation[:, ySlice]
+        horizontalSlice = self.__autocorrelation[xSlice, :]
+        return verticalSlice, horizontalSlice
+
+    def middleSlices(self):
+        # Assumes the shape of the autocorrelation is (width, height) or (nb columns, nb lines)
+        middleX, middleY = self.__autocorrelation.shape[0] // 2, self.__autocorrelation.shape[1] // 2
+        return self.slicesAt((middleX, middleY))

SpeckleSizeCode/PythonAutocorr/gaussianNormalization.py

@@ -0,0 +1,61 @@
+import tifffile
+from scipy.ndimage import gaussian_filter
+import numpy as np
+import cv2
+import os
+
+
+class GaussianNormalization:
+
+    def __init__(self, image: np.ndarray):
+        self.image = image
+        self.normalized = None
+
+    def normalize(self, stdDev: float):
+        gaussianFiltered = gaussian_filter(self.image, stdDev)
+        self.normalized = self.image / gaussianFiltered - np.mean(self.image)
+
+    def saveToOriginalFormat(self, name: str, stdDevNormalization: float = 75):
+        raise NotImplementedError("You must implement this method in format-specific classes.")
+
+
+class TiffFileGaussianNormalization(GaussianNormalization):
+
+    def __init__(self, path: str):
+        image = tifffile.imread(path)
+        super(TiffFileGaussianNormalization, self).__init__(image)
+
+    def saveToOriginalFormat(self, name: str, stdDevNormalization: float = 75):
+        if self.normalized is None:
+            self.normalize(stdDevNormalization)
+        if not name.endswith(".tif") and not name.endswith(".tiff"):
+            name += ".tif"
+        tifffile.imwrite(name, self.image)
+
+
+class PNGFileGaussianNormalization(GaussianNormalization):
+
+    def __init__(self, path: str):
+        image = cv2.imread(path)
+        super(PNGFileGaussianNormalization, self).__init__(image)
+
+    def saveToOriginalFormat(self, name: str, stdDevNormalization: float = 75):
+        if self.normalized is None:
+            self.normalize(stdDevNormalization)
+        if not name.endswith(".png"):
+            name += ".png"
+        cv2.imwrite(name, self.image)
+
+
+if __name__ == '__main__':
+    # This is where you can normalize and save files.
+    # This is an example that normalizes avery tif file in MATLAB folder.
+    for nb in range(1, 32):
+        fname = r"20190924-200ms_20mW_Ave15_Gray_10X0.4_{}.tif".format(nb)
+        p = os.path.dirname(os.path.join(os.getcwd(), "..", ".."))
+        path = os.path.join(p, "MATLAB", fname)
+        newName = path[:-4] + "_GaussianNorm_75StdDev.tif" # We don't want the previous .tif in the new name
+
+        tif = TiffFileGaussianNormalization(path)
+        tif.saveToOriginalFormat(newName)
+        print(f"{nb} / 31 processed")