diff --git a/.gitignore b/.gitignore index b6e4761..69c085c 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,9 @@ __pycache__/ *.py[cod] *$py.class +# Mac stuff +.DS_STORE + # C extensions *.so diff --git a/Code/.DS_Store b/Code/.DS_Store deleted file mode 100644 index 5008ddf..0000000 Binary files a/Code/.DS_Store and /dev/null differ diff --git a/Code/exceptions.py b/Code/exceptions.py deleted file mode 100644 index a9d5f7a..0000000 --- a/Code/exceptions.py +++ /dev/null @@ -1,36 +0,0 @@ -import numpy as np - -def isANumpyArray(parameter): - if type(parameter) is not np.ndarray: - raise Exception(f"{parameter} does not result in a numpy.ndarray.") - -def isInt(parameter): - if type(parameter) is not int: - raise Exception(f"{parameter} must be an int.") - -def isIntOrFloat(parameter): - if type(parameter) is not int and type(parameter) is not float: - raise Exception(f"{parameter} must be an int or a float.") - -def isString(parameter): - if type(parameter) is not str: - raise Exception(f"{parameter} must be a string.") - -def isTifOrTiff(path): - pathExtention = path.split(".")[-1] - if pathExtention != "tif" and pathExtention != "tiff": - raise Exception("Image must be a .tif of .tiff") - -def areValuesEqual(value1, value2): - if value1 != value2: - raise Exception("Values are different.") - -def isSameShape(image1, image2): - if image1.shape[0] != image2.shape[0]: - raise Exception("Images don't have the same size.") - elif image1.shape[1] != image2.shape[1]: - raise Exception("Images don't have the same size.") - -def isPixelATuple(pixel): - if type(pixel) is not tuple and len(pixel) != 2: - raise Exception("Pixel must be a tuple of two integers (starting from 0)!") diff --git a/Code/functions.py b/Code/functions.py deleted file mode 100644 index b689992..0000000 --- a/Code/functions.py +++ /dev/null @@ -1,555 +0,0 @@ -import tifffile as tiff -import scipy.ndimage as simg -import numpy as np -import math as math -import exceptions as exc -import cmath as cmath -import time -from rich import print -import matplotlib.pyplot as plt - -def createImage(imagepath): - """ - This function is used to open an image in Python by using the path of the image defined by imagepath. - 1. imread is used to open the image in a 2D numpy array. - 2. Convert the image to 16-bit unsigned integer format. - Returns the 2D numpy array with 16-bit unsigned integers. - """ - exc.isTifOrTiff(imagepath) - exc.isString(imagepath) - - image = tiff.imread(imagepath) - if np.min(image) < 0: - imageRescale = (image + abs(np.min(image))).astype("uint16") - else: - imageRescale = image.astype("uint16") - return imageRescale - - -def createDifferenceImage(speckle, uniform): - """ - This function does image_speckle-image_unform to have an image that has the speckle pattern on the object without - the object. The resulting image is called the difference image. - """ - exc.isANumpyArray(speckle) - exc.isANumpyArray(uniform) - exc.isSameShape(image1=speckle, image2=uniform) - - # Convert speckle and uniform images into np.float64 data type, so that the subtraction can be negative. - newSpeckle = speckle.astype(np.float64) - newUniform = uniform.astype(np.float64) - - # Subtraction - difference = newSpeckle - newUniform - - # Find the minimal value - minPixel = np.min(difference) - - # Rescale the data to only have positive values. - difference = difference + abs(minPixel) - - # Convert data in uint16 - differenceInInt = difference.astype(np.uint16) - - return differenceInInt - - -def cameraOTF(image): - """ - Optical transfer function is the fft of the PSF. The modulation transfer function is the magnitude of the - complex OTF. - """ - sizex, sizey = image.shape[0], image.shape[1] - pixels = np.zeros(shape=(sizex, sizey), dtype=np.float64) - - for x in range(sizex): - for y in range(sizey): - x1 = (2 * x - sizex) * np.pi / sizex - y1 = (2 * y - sizey) * np.pi / sizey - if x1 != 0 and y1 != 0: - pixels[x, y] = (math.sin(x1) * math.sin(y1)) / (x1 * y1) - elif x1 == 0 and y1 != 0: - pixels[x, y] = math.sin(y1) / y1 - elif x1 != 0 and y1 == 0: - pixels[x, y] = math.sin(x1) / x1 - elif x1 == 0 and y1 == 0: - pixels[x, y] = 1 - if pixels[x, y] < 0: - pixels[x, y] = 0 - return pixels - - -def imagingOTF(numAperture, wavelength, magnification, pixelSize, image): - sizex, sizey = image.shape[0], image.shape[1] - bandwidth = 2 * numAperture / (wavelength * 1e-9) - scaleUnits = magnification / (pixelSize * 1e-6) / bandwidth - pixel = np.zeros(shape=(sizex, sizey), dtype=np.float64) - - for x in range(sizex): - for y in range(sizey): - pixel[x, y] = scaleUnits * math.sqrt( - ((x + 0.5 - sizex / 2)**2 / (sizex / 2 - 1)**2) + ((y + 0.5 - sizey / 2)**2 / (sizey / 2 - 1)**2) - ) - if pixel[x, y] > 1: - pixel[x, y] = 1 - pixel[x, y] = 0.6366197723675814 * ( - math.acos(pixel[x, y]) - pixel[x, y] * math.sqrt(1 - pixel[x, y] * pixel[x, y]) - ) - if pixel[x, y] < 0: - pixel[x, y] = 0 - return pixel - - -def obtainFFTFitler(image, filteredImage): - """" - returns: - 2D numpy array of complex128 numbers. - """ - exc.isANumpyArray(image) - exc.isANumpyArray(filteredImage) - exc.isSameShape(image1=image, image2=filteredImage) - - fftImage = np.fft.fft2(image) - fftFilteredImage = np.fft.fft2(filteredImage) - fftFilter = np.divide(fftFilteredImage, fftImage) - normfftFilter = np.divide(fftFilter, np.amax(fftFilter))/fftFilter.shape[0]/fftFilter.shape[1] - - return normfftFilter - - -def applyGaussianFilter(sigma, image, truncate=4.0): - """ - Apply a Gaussian filter on the imput image. The higher the value of sigma, the bigger the filter window, the more - intense the filter will be on the image. Another way of saying this is that a grater value of sigma will make the - image blurrier. - In the HiLo algorithm, we use a gaussian filter to increase the defocus dependence of the PSF, which leads to an - improved or "stronger" optical sectioning. Sigma is the parameter that drives the optical sectioning thickness in - the HiLo algorithm. - """ - exc.isANumpyArray(image) - exc.isIntOrFloat(sigma) - imgGaussian = simg.gaussian_filter(image, sigma=sigma, truncate=truncate) - return imgGaussian - - -def valueOnePixel(image, pixelPosition): - """ - Simply extracts the value at a specific position in the image. The only reason why this is in a function is because - we also want to make sure that thing pixel value is not lower than 0. - """ - value = image[pixelPosition[0]][pixelPosition[1]] - if value < 0: - value = 0 - return value - - -def valueAllPixelsInSW(image, pixel, size, absoluteOn): - """ - Generates a list of the values of all elements in image at the center pixel position px in the area of the sampling - window samplingWindow. - """ - n = size // 2 - xpixel, ypixel = pixel[0], pixel[1] - values = [] - - for xvalue in range(xpixel - n, xpixel + n + 1): - for yvalue in range(ypixel - n, ypixel + n + 1): - if not((xvalue < 0) or (yvalue < 0) or (xvalue > image.shape[0] - 1) or (yvalue > image.shape[1] - 1)): - if absoluteOn == True: - values.append(abs(valueOnePixel(image, (xvalue, yvalue)))) - else: - values.append(valueOnePixel(image, (xvalue, yvalue))) - return values - - -def absSum(array, samplingWindow): - """ - returns: - numpy.nparray of the absolute sum of the values included in the sampling window samplingW of all pixels/elements in - array. - """ - sizex, sizey = array.shape[0], array.shape[1] - absSumImage = np.zeros(shape = (sizex, sizey)) - for x in range(sizex): - for y in range(sizey): - absSumImage[x, y] = sum(valueAllPixelsInSW(array, (x, y), samplingWindow, True)) - return absSumImage - - -def squared(array): - """ - returns: - The squared value of the list, numpy.ndarray or number. The type of array is kept throughout the execution of the - function. - """ - if type(array) is np.ndarray: - newArray = array**2 - elif type(array) is list: - newArray = [i**2 for i in array] - else: - newArray = array**2 - - return newArray - - -def stDevOnePixel(image, samplingWindow, pixel): - """ - """ - exc.isANumpyArray(image) - exc.isInt(samplingWindow) - exc.isPixelATuple(pixel) - - valuesOfPixelsInSamplingWindow = valueAllPixelsInSW(image, pixel, samplingWindow, False) - stdDev = np.std(valuesOfPixelsInSamplingWindow) - - return stdDev - - -def meanOnePixel(image, samplingWindow, pixel): - """ - """ - exc.isANumpyArray(image) - exc.isInt(samplingWindow) - exc.isPixelATuple(pixel) - - valuesOfPixelsInSamplingWindow = valueAllPixelsInSW(image, pixel, samplingWindow, False) - meanOfList = np.mean(valuesOfPixelsInSamplingWindow) - - return meanOfList - - -def stdevWholeImage(image, samplingWindow): - """ - """ - exc.isANumpyArray(image) - exc.isInt(samplingWindow) - - sizex, sizey = image.shape[0], image.shape[1] - stDevImage = np.zeros(shape=(sizex, sizey)) - - for x in range(sizex): - for y in range(sizey): - stDevImage[x, y] = stDevOnePixel(image, samplingWindow, (x, y)) - return stDevImage - - -def meanWholeImage(image, samplingWindow): - """ - """ - exc.isANumpyArray(image) - exc.isInt(samplingWindow) - - sizex, sizey = image.shape[0], image.shape[1] - meanImage = np.zeros(shape=(sizex, sizey)) - - for x in range(sizex): - for y in range(sizey): - meanImage[x, y] = meanOnePixel(image, samplingWindow, (x, y)) - return meanImage - - -def noiseInducedBias( - cameraGain, - readoutNoiseVariance, - imageSpeckle, - imageUniform, - difference, - samplingWindowSize, - sigma - ): - """ - """ - exc.isIntOrFloat(cameraGain) - exc.isIntOrFloat(readoutNoiseVariance) - exc.isANumpyArray(imageSpeckle) - exc.isANumpyArray(imageUniform) - exc.isInt(samplingWindowSize) - - meanSpeckle = meanWholeImage(image=imageSpeckle, samplingWindow=samplingWindowSize) - meanUniform = meanWholeImage(image=imageUniform, samplingWindow=samplingWindowSize) - - fftFilter = obtainFFTFitler(image=imageUniform, filteredImage=difference) - sumAbsfftFilter = absSum(array=fftFilter, samplingWindow=samplingWindowSize) - squaredSumAbsfftFilter = sumAbsfftFilter**2 - - # Calculate the noise-induced biased function. - sizex, sizey = imageSpeckle.shape[0], imageSpeckle.shape[1] - bias = np.zeros(shape=(sizex, sizey)) - for x in range(sizex): - for y in range(sizey): - bias[x, y] = ( - (np.sqrt(cameraGain * meanUniform[x, y]) + (cameraGain * meanSpeckle[x, y]) + readoutNoiseVariance) - * squaredSumAbsfftFilter[x, y] - ) - return bias - - -def contrastCalculation(uniform, speckle, samplingWindow, sigma): - """ - This function generates the contrast weighting function that peaks at 1 for in-focus regions and goes to 0 for - out-of-focus regions of the difference image. - TODO : noise function calculation doesn't work - Returns the contrast weighting function that is the 2D size of the raw input images. - """ - exc.isANumpyArray(speckle) - exc.isANumpyArray(uniform) - exc.isSameShape(image1=uniform, image2=speckle) - exc.isInt(samplingWindow) - - # create difference image and apply a gaussian filter on it - differenceImage = createDifferenceImage(speckle=speckle, uniform=uniform) - gaussianImage = applyGaussianFilter(sigma=sigma, image=differenceImage) - - # calculate the noise-induced bias of the image - - #TODO : - noiseFunction = np.zeros(shape=(uniform.shape[0], uniform.shape[1])) - # noiseFunction = noiseInducedBias( - # cameraGain=1, - # readoutNoiseVariance=0.3508935, - # imageSpeckle=speckle, - # imageUniform=uniform, - # difference=gaussianImage, - # samplingWindowSize=samplingWindow, - # sigma=sigma - # ) - - # calculate the stdev and the mean of the speckle and the uniform images - stdevDifference = stdevWholeImage(image=gaussianImage, samplingWindow=samplingWindow) - meanUniform = meanWholeImage(image=uniform, samplingWindow=samplingWindow) - - # subtract the noise-induced bias from the stdev of the filtered difference image - reducedStdevDifference = np.subtract(stdevDifference, noiseFunction) - - # calculate the contrast function - contrastFunction = reducedStdevDifference / meanUniform - - return contrastFunction - - -def lowpassFilter(image, sigmaFilter): - """ - Creates a 2D numpy array the size used as a low-pass Fourier filter. - Sigma is used again in this function to evaluate the size of the center circle, aka the frequencies to get rid of. - """ - exc.isANumpyArray(image) - exc.isIntOrFloat(sigmaFilter) - - x, y = np.meshgrid(np.linspace(-1, 1, image.shape[0]), np.linspace(-1, 1, image.shape[1])) - d = np.sqrt(x**2 + y**2) - sigma = (sigmaFilter * 0.18) / (2 * math.sqrt(2 * math.log(2))) - mu = 0.0 - gauss = (1 / (sigma * 2 * np.pi)) * np.exp(-((d - mu)**2 / (2.0 * sigma**2))) - maxPixel = np.max(gauss) - gaussNorm = gauss / maxPixel - - return gaussNorm - - -def highpassFilter(low): - """ - The high-pass Fourier filter is produced from the low-pass Fourier filter. It is litteraly its inverse. - """ - exc.isANumpyArray(low) - return 1 - low - - -def estimateEta(speckle, uniform, sigma, ffthi, fftlo, sig): - """ - TODO : It is not clear how we should evaluate the eta. - The LO image is always less intense than the HI image after treatments, so we need a scaling function eta to - readjust the intensities of the LO image. This scaling function eta compensates for the fact that the contrast - weighting function never reaches 1 even for the perfectly in-focus regions. It also prevents the discontinuities at - the cutoff frequency. - Returns the function eta in numpy.ndarray. Calculations taken from the java code for the HiLo Fiji plugin and from - personal communication with Olivier Dupont-Therrien at Bliq Photonics - """ - differenceImage = createDifferenceImage(speckle=speckle, uniform=uniform) - gaussianImage = applyGaussianFilter(sigma=sigma, image=differenceImage) - bandpassFilter = obtainFFTFitler(image=uniform, filteredImage=gaussianImage) - - illuminationOTF = imagingOTF(numAperture=1, wavelength=488e-9, magnification=20, pixelSize=6.5, image=uniform) - detectionOTF = imagingOTF(numAperture=1, wavelength=520e-9, magnification=20, pixelSize=6.5, image=uniform) - camOTF = cameraOTF(image=uniform) - # print(f"Ill OTF : {illuminationOTF}{illuminationOTF.dtype}") - # print(f"DET OTF {detectionOTF}{detectionOTF.dtype}") - - # Method 1 : Generate a function eta for each pixels. - # eta = np.zeros(shape=(uniform.shape[0], uniform.shape[1]), dtype=np.complex128) - # x = 0 - # y = 0 - # while x int: + """DOCS + """ + return int(imageFilePath.split("TimeLapse_")[1].split("-")[0]) + + +def sortListOfTimeAcqPerPlaneFiles(imageFilePath: str) -> int: + """DOCS + """ + return int(imageFilePath.split("planeNumber")[1].split("UncorrectedImages")[0]) + + +def sortListOfImageFilesSparq(imageFilePath: str) -> float: + """DOCS + """ + typeOfImage = imageFilePath.split("Sparq_")[1].split(".")[0] + timeOfImage = int(imageFilePath.split("Frame_")[1].split("-")[0]) + if typeOfImage == "Processed": + return 1 + timeOfImage + elif typeOfImage == "Uniform": + return 2 + timeOfImage + elif typeOfImage == "Speckle": + return 3 + timeOfImage + + +def readFilesInDirectory(directoryToRead: str, funcToSort) -> list: + """DOCS + """ + allFilesInDirectory = [files for files in glob.glob(directoryToRead + "*.tif")] + return sorted(allFilesInDirectory, key = funcToSort) + + +def obtainAllImagePath(unifDirPath: str, speckDirPath: str, timeAcquisitionPath: str) -> list: + """DOCS + """ + return [readFilesInDirectory(directory, sortListOfImageFiles) + for directory in (unifDirPath, speckDirPath, timeAcquisitionPath)] + + +def obtainAllImageDataFromPath(allImagePathForHiLoProcessing: list) -> list: + """DOCS + """ + allImageDataForHiLoProcessing = [ + [tf.imread(path) for path in pathLists] for pathLists in allImagePathForHiLoProcessing + ] + return allImageDataForHiLoProcessing + + +def extractDataFromSparqDirectory(sparqDirectoryPath: str) -> list: + """DOCS + """ + pathsInDir = readFilesInDirectory(sparqDirectoryPath, sortListOfImageFilesSparq) + allProcessedPaths = pathsInDir[0::3] + allUnifPaths = pathsInDir[1::3] + allSpeckPaths = pathsInDir[2::3] + return allProcessedPaths, allUnifPaths, allSpeckPaths + + +def extractDataFromDirectoriesForHiLoProcessing(unifDirPath: str, speckDirPath: str, timeAcquisitionPath: str) -> tuple: + """DOCS + """ + allImagePath = obtainAllImagePath(unifDirPath, speckDirPath, timeAcquisitionPath) + allImageData = obtainAllImageDataFromPath(allImagePath) + return (*allImageData,) + + +def removeWasteImagesFromTimeAcquisition( + timeAcquisitionData: list, + imagesPerZstack: int + ) -> list: + """DOCS + """ + indicesToRemove = [idx for idx in range(len(timeAcquisitionData)) if idx % imagesPerZstack == 0] + return np.delete(timeAcquisitionData, indicesToRemove, axis = 0) + + +def mergePlanesOfDifferentZstacks(noWasteImageTimeAcquisitionData: list, imagesPerZstack: int) -> list: + """DOCS + """ + numberOfZstack = len(noWasteImageTimeAcquisitionData) / (imagesPerZstack - 1) + splittedZstack = np.array_split(noWasteImageTimeAcquisitionData, numberOfZstack) + allPlanesMatched = [np.array(i) for i in zip(*splittedZstack)] + return allPlanesMatched + + +def createTifFilesOfPlanesFromTimeAcquisition( + directoryToInputResults: str, + timeAcquisitionData: list, + imagesPerZstack: int + ) -> tuple: + """DOCS + """ + directoryToInputData = directoryToInputResults + "timeAcquisitionInOrderOfPlanes/" + noWasteImageData = removeWasteImagesFromTimeAcquisition(timeAcquisitionData, imagesPerZstack) + mergedPlanesFromZstacks = mergePlanesOfDifferentZstacks(noWasteImageData, imagesPerZstack) + createDirectoryIfInexistant(directoryToInputData) + for index, plane in enumerate(mergedPlanesFromZstacks): + tf.imwrite(directoryToInputData + f"planeNumber{index + 1}UncorrectedImages.tif", plane) + return mergedPlanesFromZstacks, directoryToInputData + + +def createDirectoryIfInexistant(directoryName: str) -> None: + """ + This function creates a directory if it is not yet been created. + + Parameters + ---------- + directoryName : string + The name of the directory to create + + """ + if os.path.exists(directoryName) == False: + os.mkdir(directoryName) + return + + +def setParameterForHiLoComputation() -> dict: + """DOCS + """ + parameters = { + "cameraGain" : 1, + "readoutNoiseVariance" : 0.3508935, + "sigma" : 2, + "waveletGaussiansRatio" : 2, + "windowValue" : 7, + "illuminationAperture" : 1, + "detectionAperture" : 1, + "illuminationWavelength" : 488e-9, + "detectionWavelength" : 520e-9, + "pixelSize" : 4.5, + "magnification" : 20, + "superGauss": 1, + "eta": 0 + } + parameters["sigmaLP"] = 5.2 * parameters["sigma"] * (1 + parameters["waveletGaussiansRatio"] / 2) + return parameters + + +def sortMotionCorrectedImageFiles(imageFilePath: str) -> int: + """DOCS + """ + return int(imageFilePath.split("planeNumber")[1].split("U")[0]) + + +def readFilesInAllCorrectedDirectories(directoryWhereCorrectedAre: str) -> list: + """DOCS + """ + allFilesNotMergedInSubdirectory = [files for files in glob.glob(directoryWhereCorrectedAre + "*/corrected_0.tif")] + allFilesMergedInSubdirectory = [files for files in glob.glob(directoryWhereCorrectedAre + "*/corrected_and_*.tiff")] + sortedListOfNotMergedFiles = sorted(allFilesNotMergedInSubdirectory, key = sortMotionCorrectedImageFiles) + sortedListOfMergedFiles = sorted(allFilesMergedInSubdirectory, key = sortMotionCorrectedImageFiles) + + if len(allFilesMergedInSubdirectory) == 0: + return getDataOfCorrectedFilesInCorrectedDirectories(sortedListOfNotMergedFiles) + + elif len(allFilesNotMergedInSubdirectory) == len(allFilesMergedInSubdirectory): + return getDataOfCorrectedFilesInCorrectedDirectories(sortedListOfMergedFiles) + + else: + sortedListOfMergedAndNotMergedFiles = keepOnlyCorrectedOrMergedFiles( + sortedListOfNotMergedFiles, + sortedListOfMergedFiles + ) + return getDataOfCorrectedFilesInCorrectedDirectories(sortedListOfMergedAndNotMergedFiles) + + +def keepOnlyCorrectedOrMergedFiles(listOfNotMergedFiles: list, listOfMergedFiles: list) -> list: + """DOCS + """ + listOfPlaneNumbersNotMerged = [ + int(file.split("planeNumber")[1].split("U")[0]) for file in listOfNotMergedFiles + ] + listOfPlaneNumbersMerged = [ + (int(file.split("planeNumber")[1].split("U")[0]), file) + for file in listOfMergedFiles + ] + + for planeNumbers, file in listOfPlaneNumbersMerged: + listOfNotMergedFiles[planeNumbers - 1] = file + + return listOfNotMergedFiles + + +def getDataOfCorrectedFilesInCorrectedDirectories(listOfCorrectedFiles: list) -> list: + """DOCS + """ + return [tf.imread(file) for file in listOfCorrectedFiles] + + +def removeDirectory(directoryPath: str) -> None: + """DOCS + """ + return shutil.rmtree(directoryPath) + + +def temporalMeanOfTimeAcquisitionAfterMc( + timeAcqDirectoryToCorrect: str, + directoryToInputResults: str, + nameOfAveragedAndMotionCorrectedFile: str + ) -> list: + """DOCS + """ + nameOfFile = directoryToInputResults + nameOfAveragedAndMotionCorrectedFile + eachCorrectedPlanes = readFilesInAllCorrectedDirectories(timeAcqDirectoryToCorrect) + averagedAndMotionCorrectedPlanes = [np.mean(plane, axis = 0) for plane in eachCorrectedPlanes] + + # removeDirectory(timeAcqDirectoryToCorrect) + + tf.imwrite(nameOfFile, averagedAndMotionCorrectedPlanes) + return averagedAndMotionCorrectedPlanes + + +def temporalMeanOfTimeAcquisitionNoMc( + timeAcqDirectoryToCorrect: str, + directoryToInputResults: str, + nameOfAveragedFile: str + ) -> list: + """DOCS + """ + nameOfFile = directoryToInputResults + nameOfAveragedFile + eachPlanes = readFilesInDirectory(timeAcqDirectoryToCorrect, sortListOfTimeAcqPerPlaneFiles) + averagedPlanes = [np.mean(tf.imread(plane), axis = 0) for plane in eachPlanes] + tf.imwrite(nameOfFile, averagedPlanes) + return averagedPlanes + + +if __name__ == "__main__": + pass diff --git a/HiLoProcessingUtilities/motionCorrection.py b/HiLoProcessingUtilities/motionCorrection.py new file mode 100644 index 0000000..b723562 --- /dev/null +++ b/HiLoProcessingUtilities/motionCorrection.py @@ -0,0 +1,483 @@ +""" +------------------------------------------------------------------------------------------------------------------------ +Script that sends a whole directory to be motion corrected. There are a few word distinctions to make before reading +the code: + +Directory : Main directory that contains single or multiple `.tiff` files to correct + +File : Main `.tiff` file contained in the directory that need motion correction + +Subfile : To correct the main file faster, we split it into multiple smaller files called subfiles that will all + undergo motion correction + +Subdirectory : Sub directories containing the subfiles that will need motion correction +------------------------------------------------------------------------------------------------------------------------ + +Example on the main directory for only one file to correct after the correction : + + Directory + | + Subdirectory 1 ___|___ File 1 + | + Subfile 1 ___|... ___ Subfile N ___ CorrectedSubfile 1 ___... ___ CorrectedSubfile N ___ Correctedfile 1 + +------------------------------------------------------------------------------------------------------------------------ +""" +import calimba.processing.caiman.CaimanPDK as cal +import caiman as cm +import tifffile as tf +import os +import ants +import numpy as np +import globalUtilities as gu + + +def wholeMotionCorrectionOnDirectory(directoryToCorrect: str) -> None: + """ + This function determines the files to correct and sends individual motion correction onto each one. + + Parameters + ---------- + directoryToCorrect : string + Path of the main directory where all the files needing motion correction are located + """ + uncorrectedFiles = cal.find_uncorrected_files(directoryToCorrect) # Find files to correct in directory + print("Files to correct in directory: ", uncorrectedFiles) + for fileToCorrect in uncorrectedFiles: + motionCorrectionOnFile(directoryToCorrect, fileToCorrect) + return + + +def motionCorrectionOnFile(directoryToCorrect: str, fileToCorrect: str) -> None: + """ + This function creates the subdirectories and the subfiles, apply motion correction on each one of those subfile and + finally merges them back together using image registration to align them correctly. + + Parameters + ---------- + directoryToCorrect : string + Path of the main directory where all the files needing motion correction are located + + fileToCorrect : string + Name of the file to apply motion correction onto + """ + subfilesDirectory = createSubdirectoryAppendSubfiles(directoryToCorrect, fileToCorrect) + correctedSubfiles = correctTheSubfilesInSubdirectory(subfilesDirectory) + if len(correctedSubfiles) > 1: + applyRegistrationToFramesOfSubfileAndMerge(correctedSubfiles) + return + +def createSubdirectoryAppendSubfiles( + directoryToCorrect: str, + fileToCorrect: str + ) -> str: + """ + Function that creates the subdirectory and append the subfiles to it. These subfiles all come from the original file + to correct and will each need independant correction. + + Parameters + ---------- + directoryToCorrect : string + Path of the main directory where all the files needing motion correction are located + + fileToCorrect : string + Name of the file to apply motion correction onto + + Returns + ------- + nameOfSubdirectory : string + Name of the subdirectory where all the subfiles, corrected subfiles and corrected file will be + """ + fileAsArray = tf.imread(directoryToCorrect + fileToCorrect) + numberOfSubfilesToCreate = computeNumberOfSubfilesForFileToCorrect(fileAsArray) + + print(f"The file {fileToCorrect} will be splitted into", numberOfSubfilesToCreate, "subfiles to correct.") + + nameOfSubdirectory = splitFileCreateSubdirectoryAddSubfiles( + fileAsArray, + numberOfSubfilesToCreate, + directoryToCorrect, + fileToCorrect + ) + return nameOfSubdirectory + + +def correctTheSubfilesInSubdirectory(nameOfSubdirectory: str) -> list: + """ + This function sends motion correction on each `.tiff` subfile and their corrected corresponding subfile will be + outputed as `.tiff` in the same directory. + + Parameters + ---------- + nameOfSubdirectory : string + Name of the subdirectory where all the subfiles, corrected subfiles and corrected file will be + + Returns + ------- + nameOfCorrectedSubfiles : list of strings + List that contains all the names of the corrected subfiles + """ + applyMotionCorrectionWithCaiman(nameOfSubdirectory) + os.chdir(nameOfSubdirectory) + nameOfCorrectedSubfiles = sorted(getNamesOfCorrectedSubfiles(nameOfSubdirectory)) + return nameOfCorrectedSubfiles + + +def applyRegistrationToFramesOfSubfileAndMerge(nameOfCorrectedSubfiles: str) -> None: + """ + This function applies image registration on the corrected subfiles to align them all together after motion + correction. They will then be merged together to create the final corrected file. The registration is done choosing + a fixed referential subfile on which all the other subfiles will be maped. Here, the first subfile will be chosen as + the fixed referential, and the second subfile will be the moving referential. Then the transformation computed with + the registration will be applied individually onto all the frames of all the subfiles (except the first, chosen as + the fixed referential). It is also important to know that we will use the mean frame of both these fixed and moving + subfile. + + Parameters + ---------- + nameOfCorrectedSubfiles : string + List that contains all the names of the corrected subfiles + """ + correctedSubfilesAsArray = [ + np.array(tf.imread(correctedSubfile), dtype = np.float32) + for correctedSubfile in nameOfCorrectedSubfiles + ] + fixedCorrectedSubfile = correctedSubfilesAsArray[0] + + # Already include the reference corrected subfile which doesnt need registration + correctedAndRegisteredSubfilesIntoSingleArray = [frames.astype(np.float16) for frames in fixedCorrectedSubfile] + + registrationResult, meanFixedFrame = applyRegistrationOnMeanFixedAndMovingSubfiles( + fixedCorrectedSubfile, + correctedSubfilesAsArray + ) + + applyTransformationToCorrectedSubfileFrames( + correctedSubfilesAsArray, + registrationResult, + meanFixedFrame, + correctedAndRegisteredSubfilesIntoSingleArray + ) + + mergeSubfilesIntoOneCorrectedFile("corrected_and_merged.tiff", correctedAndRegisteredSubfilesIntoSingleArray) + return + + +def computeNumberOfSubfilesForFileToCorrect(fileAsArray: np.ndarray) -> int: + """ + This function computes the number of subfiles the principal `.tiff` file will be splitted into to perform motion + correction. It is to be noted that each subfile as a maximum size of 2GB. + + Parameters + ---------- + fileAsArray : numpy array of numpy arrays + Principal `.tiff` file as a single numpy array containing all image frames + + Returns + ------- + numberOfSubfilesToCreate : integer + The number of subfiles the principal file will be splitted into + """ + sizeForWholeFile = fileAsArray.size + numberOfImagesInFile = fileAsArray.shape[0] + GbInBytes = 1e-9 # Conversion factor + maximumSizeForOneSubfile = 2 # Size in GB + + sizePerImageInFile = (sizeForWholeFile / numberOfImagesInFile) * GbInBytes + numberOfImagesInOneSubfile = maximumSizeForOneSubfile / sizePerImageInFile + ratioImagesInFileToImagesInOneSubfile = numberOfImagesInFile / numberOfImagesInOneSubfile + numberOfSubfilesToCreate = ratioImagesInFileToImagesInOneSubfile * maximumSizeForOneSubfile + + if numberOfSubfilesToCreate >= maximumSizeForOneSubfile: + return int(numberOfSubfilesToCreate) + else: + return 1 + + +def splitFileCreateSubdirectoryAddSubfiles( + fileAsArray: np.ndarray, + numberOfSubfilesToCreate: int, + directoryToCorrect: str, + fileToCorrect: str + ) -> str: + """ + This function split the main file as multiple subfiles given the `numberOfSubfilesToCreate` parameters, creates the + subdirectory and appends the subfiles as `.tiff` to it. + + Parameters + ---------- + fileAsArray : numpy array + Principal `.tiff` file as a single numpy array containing all image frames + + numberOfSubfilesToCreate : integer + Number of subfiles that will be create based on the size of the main file + + fileToCorrect : string + Principal `.tiff` file as a single numpy array containing all image frames + + Returns + ------- + nameOfSubdirectory : string + The name of the subdirectory created that will contain the subfiles, the corrected subfiles and the corrected + main file + """ + nameOfSubdirectory = directoryToCorrect + fileToCorrect.split(".")[0] + "/" + splittedFileAsArray = np.array_split(fileAsArray, numberOfSubfilesToCreate) + gu.createDirectoryIfInexistant(nameOfSubdirectory) + buildSubfilesInSubdirectory(splittedFileAsArray, nameOfSubdirectory) + return nameOfSubdirectory + + +def applyMotionCorrectionWithCaiman(nameOfDirectoryToCorrect: str) -> None: + """ + This function applies motion correction to all the `.tiff` files in the specified directory. It uses the CaImAn + software tools via the calimba package. The parameters for the motion correction are all chosen for the particular + correction we want on zebrafishes and their effects are described in the CaImAn documentation. + + Parameters + ---------- + nameOfDirectoryToCorrect : string + Path of the directory the function will motion correct + + See also + ---------- + CaImAn documentation : https://caiman.readthedocs.io/en/master/ + + Calimba private github : https://github.com/PDKlab/Calimba + """ + parametersForMotionCorrection = cal.get_custom_mc_params() + parametersForMotionCorrection['strides'] = (144, 144) + parametersForMotionCorrection['overlaps'] = (72, 72) + + cal.correct_motion_directory(nameOfDirectoryToCorrect, parameters = parametersForMotionCorrection) + return + + +def getNamesOfCorrectedSubfiles(nameOfSubdirectory: str) -> list: + """ + This function searches through the the subdirectory and finds all the files that have been corrected. + + Parameters + ---------- + nameOfSubdirectory : string + The name of the subdirectory created that will contain the subfiles, the corrected subfiles and the corrected + main file + + Returns + ------- + nameOfCorrectedSubfiles : string + List that contains all the names of the corrected subfiles + """ + nameOfCorrectedSubfiles = [ + nameOfSubdirectory + subfile + for subfile in os.listdir(nameOfSubdirectory) + if subfile.split("_")[0] == "corrected" + ] + return nameOfCorrectedSubfiles + + +def applyRegistrationOnMeanFixedAndMovingSubfiles( + fixedCorrectedSubfile: np.ndarray, + correctedSubfilesAsArray: list + ) -> tuple: + """ + This function applies the registration on the chosen mean fixed frame and the chosen mean moving frames. + + Parameters + ---------- + fixedCorrectedSubfile : numpy array + Array containing all the frames of the chosen fixed corrected subfile, in this case, the first corrected + subfile + + correctedSubfilesAsArray : list of numpy arrays + List containing all the numpy arrays corresponding to each corrected subfiles + + Returns + ------- + registrationResult : dictionnary + Dictionnary that contains the result of the registration (see the ANTsPy documentation here: + https://antspy.readthedocs.io/en/latest/index.html for the contents of the dictionnary) + + meanFixedFrame : numpy array + Array containing the chosen fixed corrected subfile mean frame + """ + meanFixedFrame, meanMovingFrame = computeMeanFixedAndMovingImages(fixedCorrectedSubfile, correctedSubfilesAsArray) + + registrationResult = ants.registration( + fixed = ants.from_numpy(meanFixedFrame), + moving = ants.from_numpy(meanMovingFrame), + type_of_transform = "Affine" + ) + return registrationResult, meanFixedFrame + + +def applyTransformationToCorrectedSubfileFrames( + correctedSubfilesAsArray: list, + registrationResult: dict, + meanFixedFrame: np.ndarray, + correctedAndRegisteredSubfilesIntoSingleArray: list + ) -> None: + """ + This function applies the transformation obtained with the registration to all frames of each corrected subfiles and + appends these corrected and transformed frames to a final list used for merging all the subfiles together. + + Parameters + ---------- + correctedSubfilesAsArray : list of numpy arrays + List containing all the numpy arrays corresponding to each corrected subfiles + + registrationResult : dictionnary + Dictionnary that contains the result of the registration (see the ANTsPy documentation here + https://antspy.readthedocs.io/en/latest/index.html for the contents of the dictionnary) + + meanFixedFrame : numpy array + Array containing the chosen fixed corrected subfile mean frame + + correctedAndRegisteredSubfilesIntoSingleArray : list of numpy arrays + List containing all the corrected and registered frames of each corrected subfiles used to create the final + corrected file + + """ + for index, correctedSubfile in enumerate(correctedSubfilesAsArray): + if index != 0: # Skip the first subfile because it is the fixed file + temporaryTransformation = appendTransformsToTemporaryList( + registrationResult, + meanFixedFrame, + correctedSubfile + ) + + appendFramesToCorrectedAndRegisteredSingleArray( + correctedAndRegisteredSubfilesIntoSingleArray, + temporaryTransformation + ) + return + + +def mergeSubfilesIntoOneCorrectedFile( + nameOfCorrectedFile: str, + correctedSubfilesIntoSingleArray: list + ) -> None: + """ + This function merges all the corrected and registered frames into a single final corrected `.tiff` file. + + Parameters + ---------- + nameOfCorrectedFile : string + Name that will be given to the final corrected file + + correctedSubfilesIntoSingleArray : list of numpy arrays + List containing all the corrected and registered frames of each corrected subfiles used to create the final + corrected file + """ + tf.imwrite(nameOfCorrectedFile, correctedSubfilesIntoSingleArray, bigtiff = True, dtype = np.float16) + return + + +def buildSubfilesInSubdirectory(splittedFileAsArray: np.ndarray, nameOfSubdirectory: str) -> None: + """ + This function creates all the `.tiff` subfiles into the subdirectory. + + Parameters + ---------- + splittedFileAsArray : numpy array + Numpy array of numpy arrays where each of them contain the information of one `.tiff` subfile + + """ + os.chdir(nameOfSubdirectory) + for index, subfileAsArray in enumerate(splittedFileAsArray): + tf.imwrite(f"{index}.tif", subfileAsArray) + os.chdir("..") + return + + +def computeMeanFixedAndMovingImages( + fixedCorrectedSubfile: np.ndarray, + correctedSubfilesAsArray: list + ) -> tuple: + """ + This function computes the fixed corrected subfile mean frame and the moving corrected subfile mean frame. + + Parameters + ---------- + fixedCorrectedSubfile : numpy array + Array containing all the frames of the chosen fixed corrected subfile, in this case, the first corrected + subfile + + correctedSubfilesAsArray : list of numpy arrays + List containing all the numpy arrays corresponding to each corrected subfiles + + Returns + ------- + meanFixedFrame : numpy array + Array containing the chosen fixed corrected subfile mean frame + + meanMovingFrame : numpy array + Array containing the chosen moving corrected subfile mean frame + """ + movingCorrectedSubfile = correctedSubfilesAsArray[1] # Chosen to be the second corrected image but could be other + meanFixedFrame = fixedCorrectedSubfile.mean(axis = 0) + meanMovingFrame = movingCorrectedSubfile.mean(axis = 0) + return meanFixedFrame, meanMovingFrame + + +def appendTransformsToTemporaryList( + registrationResult: dict, + meanFixedFrame: np.ndarray, + correctedSubfile: list + ) -> list: + """ + This function applies the registration transformation to all the frames a corrected subfile and appends them to a + temporary list. + + Parameters + ---------- + registrationResult : dictionnary + Dictionnary that contains the result of the registration (see the ANTsPy documentation here + https://antspy.readthedocs.io/en/latest/index.html for the contents of the dictionnary) + + meanFixedFrame : numpy array + Array containing the chosen fixed corrected subfile mean frame + + correctedSubfile : numpy array + Array corresponding to the frames of a corrected subfile + + Returns + ------- + temporaryTransformation : list of ANTs images + Temporary list containing the corrected and transformed frames of a subfile as ANTs images (see the ANTsPy + documentation here https://antspy.readthedocs.io/en/latest/index.html) + """ + temporaryTransformation = [ + ants.apply_transforms( + fixed = ants.from_numpy(meanFixedFrame), + moving = ants.from_numpy(frame), + transformlist = registrationResult["fwdtransforms"] + ) + for frame in correctedSubfile + ] + return temporaryTransformation + + +def appendFramesToCorrectedAndRegisteredSingleArray(toAppendTo: list, temporaryTransformation: list) -> None: + """ + This function appends the registered and transformed frames of a corrected subfile as a numpy array to the + final list that contains the corrected and transformed frames of all corrected subfiles. + + Parameters + ---------- + toAppendTo : list of numpy arrays + final list that contains the corrected and transformed frames of all corrected subfiles + + temporaryTransformation : list of ANTs images + Temporary list containing the corrected and transformed frames of a subfile as ANTs images (see the ANTsPy + documentation here https://antspy.readthedocs.io/en/latest/index.html) + """ + for correctedAndRegisteredSubfileFrame in temporaryTransformation: + toAppendTo.append( + correctedAndRegisteredSubfileFrame.numpy().astype(np.float16) + ) + return + + +if __name__ == "__main__": + pass diff --git a/HiLoProcessingUtilities/processingHiLoImages.py b/HiLoProcessingUtilities/processingHiLoImages.py new file mode 100644 index 0000000..ab95190 --- /dev/null +++ b/HiLoProcessingUtilities/processingHiLoImages.py @@ -0,0 +1,549 @@ +""" +Module containing all the functions to do HiLo image treatement given a speckle and uniform stack of single image +""" +import numpy as np +import math +import matplotlib.pyplot as plt +from scipy.fft import fft2, ifft2, fftshift, fftfreq +from scipy.integrate import simpson +import multiprocess as mp +import glob +import os +import tifffile as tf +import globalUtilities as gu +from rich import print + + +class FiltersForHiLo: + """DOCS + """ + + def __init__(self, sizeOfFilter: tuple, filterType: str, sigma: float) -> None: + """DOCS + """ + self.sizeAlongKSpaceX = sizeOfFilter[0] + self.sizeAlongKSpaceY = sizeOfFilter[1] + self.kSpaceX, self.kSpaceY = np.meshgrid( + range(self.sizeAlongKSpaceX), + range(self.sizeAlongKSpaceY) + ) + self.kVectorDistance = ( + (self.kSpaceX + 0.5 - (self.sizeAlongKSpaceX / 2))**2 + + (self.kSpaceY + 0.5 - (self.sizeAlongKSpaceY / 2))**2 + ) + self.sigma = sigma + self.data = self.findFilterType(filterType) + + + def simpleGaussianFilter(self) -> np.ndarray: + """DOCS + """ + simpleGaussFilter = np.exp(-(self.kVectorDistance / (2.0 * self.sigma**2))**globalParams["superGauss"]) + normSimpleGaussFilter = simpleGaussFilter / np.max(simpleGaussFilter) + return normSimpleGaussFilter + + + def doubleGaussianFilter(self) -> np.ndarray: + """DOCS + """ + firstGaussianConstant = ( + -1.0 + * ((2 * np.pi) / self.sizeAlongKSpaceX) + * ((2 * np.pi) / self.sizeAlongKSpaceY) + * self.sigma**2 + ) + secondGaussianConstant = ( + -1.0 + * ((2 * np.pi) / self.sizeAlongKSpaceX) + * ((2 * np.pi) / self.sizeAlongKSpaceY) + * self.sigma**2 + * globalParams["waveletGaussiansRatio"]**2 + ) + doubleGaussFilter = ( + np.exp(self.kVectorDistance * firstGaussianConstant) + - np.exp(self.kVectorDistance * secondGaussianConstant) + ) + normDoubleGaussFilter = doubleGaussFilter / np.max(doubleGaussFilter) + return normDoubleGaussFilter + + + def findFilterType(self, filterType: str): + """DOCS + """ + filterContainer = { + "lowpass": self.simpleGaussianFilter(), + "highpass": 1 - self.simpleGaussianFilter(), + "doublegauss": self.doubleGaussianFilter() + } + + if filterType in filterContainer.keys(): + return filterContainer[filterType] + else: + raise Exception("The filter you want does not exist, maybe you can create it?") + + + +class ImageForHiLo: + """DOCS + """ + + def __init__(self, samplingWindow: int, **kwargs) -> None: + """DOCS + """ + if "imagePath" in kwargs: + self.data = (createImageFromPath(kwargs["imagePath"])).astype(np.float64) + else: + self.data = kwargs["imageArray"] + + self.sizeAlongXAxis = self.data.shape[0] + self.sizeAlongYAxis = self.data.shape[1] + self.samplingWindow = samplingWindow + self.standardDev = np.ndarray(shape = (self.sizeAlongXAxis, self.sizeAlongYAxis)) + self.mean = np.ndarray(shape = (self.sizeAlongXAxis, self.sizeAlongYAxis)) + + + def fftOnImage(self) -> np.ndarray: + """DOCS + """ + fourierTransform = fft2(self.data) + return fftshift(fourierTransform) + + + def applyFilter(self, theFilter: FiltersForHiLo) -> np.ndarray: + """DOCS + """ + imageInKSpace = self.fftOnImage() + imageFilteredInKSpace = imageInKSpace * theFilter.data + return ifft2(imageFilteredInKSpace) + + + def viewAllPixelsInSamplingWindow(self, pixelCoords: tuple) -> list: + """DOCS + """ + n = self.samplingWindow // 2 + xPixel, yPixel = pixelCoords[0], pixelCoords[1] + pixelValuesInSamplingWindow = [] + + for x in range(xPixel - n, xPixel + n + 1): + for y in range(yPixel - n, yPixel + n + 1): + if not ((x < 0) or (y < 0) or (x > self.sizeAlongXAxis - 1) or (y > self.sizeAlongYAxis - 1)): + pixelValuesInSamplingWindow.append(self.data[x, y]) + return pixelValuesInSamplingWindow + + + def stdDevOfWholeImage(self) -> None: + """DOCS + """ + for x in range(self.sizeAlongXAxis): + for y in range(self.sizeAlongYAxis): + self.standardDev[x, y] = np.std(self.viewAllPixelsInSamplingWindow((x, y))) + return + + + def meanOfWholeImage(self) -> None: + """DOCS + """ + for x in range(self.sizeAlongXAxis): + for y in range(self.sizeAlongYAxis): + self.mean[x, y] = np.mean(self.viewAllPixelsInSamplingWindow((x, y))) + return + + + def showImageInRealSpace(self) -> None: + """DOCS + """ + xAxis, yAxis = range(self.sizeAlongXAxis), range(self.sizeAlongYAxis) + xAxisGrid, yAxisGrid = np.meshgrid(xAxis, yAxis) + fig, ax = plt.subplots() + c = ax.pcolormesh( + xAxisGrid, + yAxisGrid, + self.data, + cmap = "gray", + vmin = np.min(self.data), + vmax = np.max(self.data) + ) + ax.set_xlabel("x") + ax.set_ylabel("y") + fig.colorbar(c, ax = ax, label = "Intensity") + plt.show() + return + + +def contrastCalculation( + uniformImage: ImageForHiLo, + speckleImage: ImageForHiLo, + bandpassFilter: FiltersForHiLo + ) -> np.ndarray: + """DOCS + """ + differenceImage = createDifferenceImage(speckleImage, uniformImage) + differenceImageWithDefocusIncrease = ImageForHiLo( + globalParams["windowValue"], + imageArray = differenceImage.applyFilter(bandpassFilter) + ) + + differenceImageWithDefocusIncrease.stdDevOfWholeImage() + speckleImage.meanOfWholeImage() + + noiseInducedBias = np.zeros( + shape = (differenceImageWithDefocusIncrease.sizeAlongXAxis, differenceImageWithDefocusIncrease.sizeAlongYAxis) + ) + + # For now not working + # noiseInducedBias = noiseInducedBiasComputation( + # speckleImage, + # uniformImage, + # doubleGaussianFilter + # ) + + constrastFunctionSquared = ( + (differenceImageWithDefocusIncrease.standardDev**2 - noiseInducedBias) / + speckleImage.mean**2 + ) + + return np.sqrt(constrastFunctionSquared) + + +def rescaleImage(image: np.ndarray) -> np.ndarray: + """DOCS + """ + minValue = np.min(image) + if minValue < 0: + return image + abs(minValue) + else: + return image + + +def createImageFromPath(imagePath: str) -> np.ndarray: + """ + This function is used to read an image and output the data related to that image. If the image as negative values, + it will shift the image to all positive values. + + Parameters + ---------- + imagePath: string + Path of the image to read and get the pixel values from. + + Returns + ---------- + imageAsArray: np.ndarray + The image pixel values as a numpy 2 dimensionnal array. + """ + imageAsArray = tf.imread(imagePath) + + return rescaleImage(imageAsArray) + + +def createDifferenceImage(imageBeingSubstracted: ImageForHiLo, imageToSubstract: ImageForHiLo) -> ImageForHiLo: + """DOCS + """ + substractedImageData = rescaleImage(imageBeingSubstracted.data - imageToSubstract.data) + return ImageForHiLo(globalParams["windowValue"], imageArray = substractedImageData) + + +def noiseInducedBiasComputation( + speckleImage: ImageForHiLo, + uniformImage: ImageForHiLo, + bandpassFilter: FiltersForHiLo + ) -> np.ndarray: + """DOCS + """ + speckleImage.meanOfWholeImage() + uniformImage.meanOfWholeImage() + meanOfSpeckleImage, meanOfUniformImage = speckleImage.mean, uniformImage.mean + + filterIntegration = integrate2DArray( + bandpassFilter.kSpaceX[0], + bandpassFilter.kSpaceY.T[0], + np.abs(bandpassFilter.doubleGaussianFilter())**2 + ) + + sizeX, sizeY = speckleImage.sizeAlongXAxis, speckleImage.sizeAlongYAxis + + noiseBias = np.zeros(shape = (sizeX, sizeY)) + for x in range(sizeX): + for y in range(sizeY): + noiseBias[x, y] = ( + ( + ( + (globalParams["cameraGain"] * meanOfUniformImage[x, y]) + + (globalParams["cameraGain"] * meanOfSpeckleImage[x, y]) + + globalParams["readoutNoiseVariance"] + ) * filterIntegration + ) + ) + return noiseBias + + +def integrate2DArray(xDomain: np.ndarray, yDomain: np.ndarray, theArray: np.ndarray) -> float: + """DOCS + """ + firstDomainIntegration = simpson(theArray, yDomain) + secondDomainIntegration = simpson(firstDomainIntegration, xDomain) + return secondDomainIntegration + + +def cameraOTF(sizeAlongXAxis: int, sizeAlongYAxis: int) -> np.ndarray: + """DOCS + """ + pixelValues = np.zeros(shape = (sizeAlongXAxis, sizeAlongYAxis)) + for x in range(sizeAlongXAxis): + for y in range(sizeAlongYAxis): + newX = (2 * x - sizeAlongXAxis) * np.pi / sizeAlongXAxis + newY = (2 * y - sizeAlongYAxis) * np.pi / sizeAlongYAxis + if newX != 0 and newY != 0: + pixelValues[x, y] = ( + (math.sin(newX) * math.sin(newY)) / (newX * newY) + ) + elif newX == 0 and newY != 0: + pixelValues[x, y] = math.sin(newY) / newY + elif newX != 0 and newY == 0: + pixelValues[x, y] = math.sin(newX) / newX + elif newX == 0 and newY == 0: + pixelValues[x, y] = 1 + if pixelValues[x, y] < 0: + pixelValues[x, y] = 0 + return pixelValues + + +def imagingOTF(sizeAlongXAxis: int, sizeAlongYAxis: int, numAperture: float, wavelength: float) -> np.ndarray: + """DOCS + """ + bandwidth = 2 * numAperture / (wavelength * 1e-9) + scaleUnits = globalParams["magnification"] / (globalParams["pixelSize"] * 1e-6) / bandwidth + pixelValues = np.zeros(shape = (sizeAlongXAxis, sizeAlongYAxis)) + for x in range(sizeAlongXAxis): + for y in range(sizeAlongYAxis): + pixelValues[x, y] = scaleUnits * math.sqrt( + ( + ( + ( + (x + 0.5 - sizeAlongXAxis / 2)**2 + / (sizeAlongXAxis / 2 - 1)**2 + ) + + + ( + (y + 0.5 - sizeAlongYAxis / 2)**2 + / (sizeAlongYAxis / 2 - 1)**2 + ) + ) + ) + ) + if pixelValues[x, y] > 1: + pixelValues[x, y] = 1 + pixelValues[x, y] = 0.6366197723675814 * ( + ( + math.acos(pixelValues[x, y]) - pixelValues[x, y] + * math.sqrt(1 - pixelValues[x, y] * pixelValues[x, y]) + ) + ) + if pixelValues[x, y] < 0: + pixelValues[x, y] = 0 + return pixelValues + + +def estimateEta(bandpassFilter: FiltersForHiLo) -> float: + """DOCS + """ + camOTF = cameraOTF(bandpassFilter.sizeAlongKSpaceX, bandpassFilter.sizeAlongKSpaceY) + illuminationOTF = imagingOTF( + bandpassFilter.sizeAlongKSpaceX, + bandpassFilter.sizeAlongKSpaceY, + globalParams["illuminationAperture"], + globalParams["illuminationWavelength"] + ) + detectionOTF = imagingOTF( + bandpassFilter.sizeAlongKSpaceX, + bandpassFilter.sizeAlongKSpaceY, + globalParams["detectionAperture"], + globalParams["detectionWavelength"] + ) + bandpassFilterData = bandpassFilter.data + fudgeFactor = 1.2 + numerator = 0 + denominator = 0 + for x in range(bandpassFilter.sizeAlongKSpaceX): + for y in range(bandpassFilter.sizeAlongKSpaceY): + denominator += ( + (bandpassFilterData[x, y] * detectionOTF[x, y] * camOTF[x, y])**2 * abs(illuminationOTF[x, y]) + ) + numerator += illuminationOTF[x, y] + + eta = math.sqrt(numerator / denominator) * fudgeFactor + return eta + + +def computeLoImageOfHiLo( + uniformImage: ImageForHiLo, + speckleImage: ImageForHiLo, + bandpassFilter: FiltersForHiLo + ) -> ImageForHiLo: + """DOCS + """ + contrastFunction = contrastCalculation(uniformImage, speckleImage, bandpassFilter) + contrastTimesUniform = ImageForHiLo( + globalParams["windowValue"], + imageArray = contrastFunction * uniformImage.data + ) + lowpassFilter = FiltersForHiLo((sizeXForFilter, sizeYForFilter), "lowpass", globalParams["sigmaLP"]) + loImage = contrastTimesUniform.applyFilter(lowpassFilter) + return ImageForHiLo(globalParams["windowValue"], imageArray = loImage) + + +def computeHiImageOfHiLo(uniformImage: ImageForHiLo) -> ImageForHiLo: + """DOCS + """ + highpassFilter = FiltersForHiLo((sizeXForFilter, sizeYForFilter), "highpass", globalParams["sigmaLP"]) + hiImage = uniformImage.applyFilter(highpassFilter) + return ImageForHiLo(globalParams["windowValue"], imageArray = hiImage) + + +def checkIfSizeMatch(uniformImage: ImageForHiLo, speckleImage: ImageForHiLo) -> tuple: + """DOCS + """ + if ( + (uniformImage.sizeAlongXAxis != speckleImage.sizeAlongXAxis) or + (uniformImage.sizeAlongYAxis != uniformImage.sizeAlongYAxis) + ): + raise Exception("The size of the the speckle image and the uniform image do not match!") + else: + return uniformImage.sizeAlongXAxis, uniformImage.sizeAlongYAxis + + +def checkIfLengthMatch(allUniformPath: list, allSpecklePath: list) -> None: + """DOCS + """ + if len(allUniformPath) != len(allSpecklePath): + raise Exception("The number of uniform and speckle images is not the same!") + return + + +def createHiLoImage(uniformImage: ImageForHiLo, speckleImage: ImageForHiLo) -> ImageForHiLo: + """DOCS + """ + print("Creating HiLo image!") + # Step to compute eta + bandpassFilter = FiltersForHiLo((sizeXForFilter, sizeYForFilter), "doublegauss", globalParams["sigma"]) + eta = estimateEta(bandpassFilter) + globalParams["eta"] = eta + + # Step to compute Lo image of HiLo + loImage = computeLoImageOfHiLo(uniformImage, speckleImage, bandpassFilter) + + # Step to compute the Hi image of HiLo + hiImage = computeHiImageOfHiLo(uniformImage) + + # Build the HiLo final image + hiLoImage = ImageForHiLo( + globalParams["windowValue"], + imageArray = globalParams["eta"] * np.abs(loImage.data) + np.abs(hiImage.data) + ) + return hiLoImage + + +def sendMultiprocessingUnits(functionToSplit, paramsToIterateOver: list) -> np.ndarray: + """DOCS + """ + processes = mp.Pool() + resultingArray = processes.starmap(functionToSplit, paramsToIterateOver) + processes.close() + processes.join() + return resultingArray + + +def obtainUniformAndSpecklesFromDir(uniformDirectory: str, speckleDirectory: str) -> tuple: + """DOCS + """ + uniformFiles, speckleFiles = gu.readFilesInDirectory(uniformDirectory), gu.readFilesInDirectory(speckleDirectory) + return uniformFiles, speckleFiles + + +def buildUniformAndSpeckleObjectsMultipleTifs(allUniformPath: list, allSpecklePath: list) -> list: + """DOCS + """ + checkIfLengthMatch(allUniformPath, allSpecklePath) + allImageForHiLo = [ + ( + ImageForHiLo(imagePath = unif, samplingWindow = globalParams["windowValue"]), + ImageForHiLo(imagePath = speck, samplingWindow = globalParams["windowValue"]) + ) + for unif, speck in zip(allUniformPath, allSpecklePath) + ] + sizeXForFilter, sizeYForFilter = checkIfSizeMatch(allImageForHiLo[0][0], allImageForHiLo[0][1]) + return allImageForHiLo, sizeXForFilter, sizeYForFilter + + +def buildUniformAndSpeckleObjectsArray(unifData: list, speckData: list) -> list: + """DOCS + """ + allImageForHiLo = [ + ( + ImageForHiLo(imageArray = unif, samplingWindow = globalParams["windowValue"]), + ImageForHiLo(imageArray = speck, samplingWindow = globalParams["windowValue"]) + ) + for unif, speck in zip(unifData, speckData) + ] + sizeXForFilter, sizeYForFilter = checkIfSizeMatch(allImageForHiLo[0][0], allImageForHiLo[0][1]) + return allImageForHiLo, sizeXForFilter, sizeYForFilter + + +def buildHiLoImagesDirectory(resultDirectory: str, resultHiLoStackName: str, allHiLoImages: list) -> None: + """DOCS + """ + hiLoImagesData = [image.data for image in allHiLoImages] + directoryToInputResult = resultDirectory + "HiLoImages/" + gu.createDirectoryIfInexistant(directoryToInputResult) + tf.imwrite(directoryToInputResult + resultHiLoStackName, hiLoImagesData, dtype = "uint16") + return + + +def runWholeHiLoImageProcessingOnDir( + simParams: dict, + mainDir: str, + unifDirPath: str, + speckDirPath: str, + resultHiLoStackName: str + ) -> None: + """DOCS + """ + global globalParams + globalParams = simParams + + allUniformPath, allSpecklePath = obtainUniformAndSpecklesFromDir(unifDirPath, speckDirPath) + allImageForHiLo, xFilter, yFilter = buildUniformAndSpeckleObjectsMultipleTifs(allUniformPath, allSpecklePath) + + global sizeXForFilter + global sizeYForFilter + sizeXForFilter, sizeYForFilter = xFilter, yFilter + + hiLoImages = sendMultiprocessingUnits(createHiLoImage, allImageForHiLo) + buildHiLoImagesDirectory(mainDir, resultHiLoStackName, hiLoImages) + return hiLoImages + + +def runWholeHiLoImageProcessingOnPaths( + simParams: dict, + resultDir: str, + unifDataPath: list, + speckDataPath: list, + resultHiLoStackName: str + ) -> None: + """DOCS + """ + global globalParams + globalParams = simParams + + allImageForHiLo, xFilter, yFilter = buildUniformAndSpeckleObjectsMultipleTifs(unifDataPath, speckDataPath) + + global sizeXForFilter + global sizeYForFilter + sizeXForFilter, sizeYForFilter = xFilter, yFilter + + hiLoImages = sendMultiprocessingUnits(createHiLoImage, allImageForHiLo) + buildHiLoImagesDirectory(resultDir, resultHiLoStackName, hiLoImages) + return hiLoImages + + +if __name__ == "__main__": + pass diff --git a/HiLoProcessingUtilities/registrationModule.py b/HiLoProcessingUtilities/registrationModule.py new file mode 100644 index 0000000..27837d8 --- /dev/null +++ b/HiLoProcessingUtilities/registrationModule.py @@ -0,0 +1,112 @@ +""" +Module containing functions to do registration from one stack of chosen image to another +""" +import numpy as np +import subprocess as sub +import nrrd +import tifffile as tf +import os +import glob +from rich import print +import globalUtilities as gu + +def registrationToConsole( + whereToStoreResults: str, + movingReferentialStack: str, + fixedReferentialStack: str, + alignedFilename: str, + pathOfANTs: str + ) -> str: + """DOCS + """ + alignedStackName = whereToStoreResults + alignedFilename + print("Registration starting!") + os.system( + f"export ANTSPATH={pathOfANTs};" + + f"export PATH=$ANTSPATH:$PATH;" + + f"export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=18;" + + f"antsRegistration -d 3 --float 1 -v --output [,{alignedStackName}] --interpolation gaussian" + + f" --use-histogram-matching 0 -r [{fixedReferentialStack},{movingReferentialStack},1]" + + f" -t rigid[0.05]" + + f" -m MI[{fixedReferentialStack},{movingReferentialStack},1,32,Regular,0.25] -c [1000x500x250x125,1e-8,10]" + + f" --shrink-factors 12x8x4x2 --smoothing-sigmas 4x3x2x1vox" + + f" -t Affine[0.1] -m" + + f" MI[{fixedReferentialStack},{movingReferentialStack},1,32,Regular,0.25] --convergence" + + f" [1000x500x250x125,1e-8,10] --shrink-factors 12x8x4x2 --smoothing-sigmas 4x3x2x1vox" + + f" -t SyN[0.2,6,0] -m" + + f" CC[{fixedReferentialStack},{movingReferentialStack},1,2] -c [200x200x200x200x20,1e-8,10]" + + f" --shrink-factors 12x8x4x2x1 --smoothing-sigmas 4x3x2x1x0vox" + ) + removeTransformsFilesAfterRegistration(os.getcwd()) + return + + +def applyRegistrationToTwoImageStacks( + whereToStoreResults: str, + movingReferentialStackData: np.ndarray, + movingReferentialStackName: str, + fixedReferentialStackData: np.ndarray, + fixedReferentialStackName: str, + pathOfANTs: str + ) -> np.ndarray: + """DOCS + """ + registrationResults = whereToStoreResults + "registrationResults/" + gu.createDirectoryIfInexistant(registrationResults) + + alignedStackName = ( + movingReferentialStackName.split(".")[0] + + "WarpedOn" + + fixedReferentialStackName.split(".")[0] + + ".nrrd" + ) + saveStackAsNrrd( + registrationResults, + movingReferentialStackData, + movingReferentialStackName + ) + saveStackAsNrrd( + registrationResults, + fixedReferentialStackData, + fixedReferentialStackName + ) + registrationToConsole( + registrationResults, + registrationResults + movingReferentialStackName, + registrationResults + fixedReferentialStackName, + alignedStackName, + pathOfANTs + ) + movingWarpedOnFixed = readNrrdStackAndStaveAsTif(registrationResults, alignedStackName) + return movingWarpedOnFixed + + +def removeTransformsFilesAfterRegistration(directoryToRemove: str) -> None: + """DOCS + """ + filesInWorkingDir = glob.glob(directoryToRemove + "/*.*") + filesInWorkingDirToRemove = [ + file for file in filesInWorkingDir if file.split(".")[1] == "mat" or file.split(".")[1] == "nii" + ] + for file in filesInWorkingDirToRemove: + os.remove(file) + return + + +def saveStackAsNrrd(whereToSave: str, stackDataToSave: str, nameOfFileToSave: str) -> None: + """DOCS + """ + nrrd.write(whereToSave + nameOfFileToSave, np.array(stackDataToSave), index_order = "C") + return + + +def readNrrdStackAndStaveAsTif(whereToSave: str, nameOfFileToRead: str) -> np.ndarray: + """DOCS + """ + dataOfFile, headerOfFile = nrrd.read(whereToSave + nameOfFileToRead, index_order = "C") + tf.imwrite(whereToSave + nameOfFileToRead.split(".")[0] + ".tif", dataOfFile) + return dataOfFile + + +if __name__ == "__main__": + pass diff --git a/HiLoProcessingUtilities/wholeHiLoProcessFlow.py b/HiLoProcessingUtilities/wholeHiLoProcessFlow.py new file mode 100644 index 0000000..70e097e --- /dev/null +++ b/HiLoProcessingUtilities/wholeHiLoProcessFlow.py @@ -0,0 +1,41 @@ +import globalUtilities as gu +import processingHiLoImages as hilo +import motionCorrection as mc +import registrationModule as reg +from rich import print +import time +import tifffile as tf +import numpy as np + + +if __name__ == "__main__": + mainDirectory = "/Users/dcclab/Desktop/MainDir/" + timeAcqSparqDir = "/Users/dcclab/Desktop/MainDir/timeAcqSparqDir/" + resultDirectory = "/Users/dcclab/Desktop/MainDir/resultsOfProcessing/" + gu.createDirectoryIfInexistant(resultDirectory) + + allProcessedPath, allUnifPath, allSpeckPath = gu.extractDataFromSparqDirectory(timeAcqSparqDir) + globalConstants = gu.setParameterForHiLoComputation() + + # Compute functionnal HiLo images + functionnalHiLoImages = hilo.runWholeHiLoImageProcessingOnPaths( + globalConstants, + resultDirectory, + allUnifPath, + allSpeckPath, + f"functionnalHiLoImages.tif" + ) + + # Treatement of the time acquisition datas + imagePerZstackForTimeAcquisition = 26 + timeAcqInOrderOfPlanes, timeAcqDirectoryToCorrect = gu.createTifFilesOfPlanesFromTimeAcquisition( + resultDirectory, + functionnalHiLoImages, + imagePerZstackForTimeAcquisition + ) + # mc.wholeMotionCorrectionOnDirectory(timeAcqDirectoryToCorrect) # Motion correction on the time acq data + HiLoFunctionnalZstack = gu.temporalMeanOfTimeAcquisitionNoMc( + timeAcqDirectoryToCorrect, + resultDirectory, + "HiLoFunctionnalZstackTempMean.tif" + ) diff --git a/dataDirectory/speckleSingleFile/200plans-2um-15mW-500ms-speckle-001-Orca flash-Ch-1-TimeLapse_100-Time_99.393s-Frame_000000.tif b/dataDirectory/speckleSingleFile/200plans-2um-15mW-500ms-speckle-001-Orca flash-Ch-1-TimeLapse_100-Time_99.393s-Frame_000000.tif new file mode 100644 index 0000000..4d32514 Binary files /dev/null and b/dataDirectory/speckleSingleFile/200plans-2um-15mW-500ms-speckle-001-Orca flash-Ch-1-TimeLapse_100-Time_99.393s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_978-Time_122.224s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_978-Time_122.224s-Frame_000000.tif new file mode 100644 index 0000000..e6bb66f Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_978-Time_122.224s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_979-Time_122.349s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_979-Time_122.349s-Frame_000000.tif new file mode 100644 index 0000000..d00859a Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_979-Time_122.349s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_980-Time_122.474s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_980-Time_122.474s-Frame_000000.tif new file mode 100644 index 0000000..a0023f1 Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_980-Time_122.474s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_981-Time_122.599s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_981-Time_122.599s-Frame_000000.tif new file mode 100644 index 0000000..ab9bdaf Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_981-Time_122.599s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_982-Time_122.724s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_982-Time_122.724s-Frame_000000.tif new file mode 100644 index 0000000..1da5a6c Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_982-Time_122.724s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_983-Time_122.849s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_983-Time_122.849s-Frame_000000.tif new file mode 100644 index 0000000..ac38738 Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_983-Time_122.849s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_984-Time_122.974s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_984-Time_122.974s-Frame_000000.tif new file mode 100644 index 0000000..2474229 Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_984-Time_122.974s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_985-Time_123.099s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_985-Time_123.099s-Frame_000000.tif new file mode 100644 index 0000000..1b5a8ba Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_985-Time_123.099s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_986-Time_123.224s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_986-Time_123.224s-Frame_000000.tif new file mode 100644 index 0000000..515645a Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_986-Time_123.224s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_987-Time_123.349s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_987-Time_123.349s-Frame_000000.tif new file mode 100644 index 0000000..3296607 Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_987-Time_123.349s-Frame_000000.tif differ diff --git a/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_988-Time_123.474s-Frame_000000.tif b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_988-Time_123.474s-Frame_000000.tif new file mode 100644 index 0000000..a89cf53 Binary files /dev/null and b/dataDirectory/timeAcquisition/1000plans25plansparVol-10umStep-200mW-25ms-uniform-001-Orca flash-Ch-1-TimeLapse_988-Time_123.474s-Frame_000000.tif differ diff --git a/dataDirectory/uniformSingleFile/200plans-2um-15mW-500ms-uniform-001-Orca flash-Ch-1-TimeLapse_100-Time_99.148s-Frame_000000.tif b/dataDirectory/uniformSingleFile/200plans-2um-15mW-500ms-uniform-001-Orca flash-Ch-1-TimeLapse_100-Time_99.148s-Frame_000000.tif new file mode 100644 index 0000000..df4a51c Binary files /dev/null and b/dataDirectory/uniformSingleFile/200plans-2um-15mW-500ms-uniform-001-Orca flash-Ch-1-TimeLapse_100-Time_99.148s-Frame_000000.tif differ diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..7231fce --- /dev/null +++ b/requirements.txt @@ -0,0 +1,10 @@ +antspyx==0.3.8 +caiman==1.9.15 +Calimba==1.0 +matplotlib==3.5.1 +multiprocess==0.70.14 +numpy==1.21.6 +pynrrd==0.4.2 +rich==13.5.1 +scipy==1.8.0 +tifffile==2022.3.16