import numpy as np import math from numba import jit from numba.core.errors import NumbaDeprecationWarning, NumbaPendingDeprecationWarning import warnings import pywt from skimage.filters import gaussian from skimage.transform import rescale # warnings.simplefilter('ignore', category=NumbaDeprecationWarning) # warnings.simplefilter('ignore', category=NumbaPendingDeprecationWarning) def check_dims(image): """ Checks image dimentions and adds a 3rd axis in case it is 2D :param image: image to add axis if needed :return: 3D image """ if len(image.shape) < 2 or len(image.shape) > 3: raise ValueError() if len(image.shape) < 3: image = np.expand_dims(image, axis=0) return image # @jit def min_max_normalization(image, low, high): """ Min/Max Normalization of the image using low and high percentiles and following the formula: I_new = (I - min) * (Max_new - Min_new) / (Max - Min) + Min_new :param image: stack or image to be normalized with low and high percentile :param low: lower percentile for normalization :param high: higher percentile for normalization :return: normalized image using the defined percentiles """ # image = check_dims(image) if low is not None: p_low = low else: p_low = 1.0 if high is not None: p_high = high else: p_high = 99.0 minI = np.percentile(image, p_low) maxI = np.percentile(image, p_high) image_norm = (image - minI) / (maxI - minI + 1e-20) return image_norm # @jit def absolute_laplacian(image, do_gauss): """ Wrapper(ish) function Computes the Absolute Laplacian of each slice of a 3D stack and returns an array with values :param image: image or stack where to compute the absolute laplacians :return: array shape=(len(z)) containing the absolute laplacian for each slice """ # n_px_slice = (image.shape[1] - 2) * (image.shape[2] - 2) abs_laplacian_array = np.zeros((image.shape[0]), dtype='float32') for z in range(image.shape[0]): image_slice = image[z].copy() if do_gauss: image_slice = gaussian(image_slice, sigma=1) abs_laplacian_array[z] = absolute_laplacian_slice(image_slice) return abs_laplacian_array def absolute_laplacian_scaled(image, scale_factor): """ TODO add detailed description :param image: :param scale_factor: :return: """ assert scale_factor <= 1, print("Scale factor needs to be smaller or equal to 1") abs_laplacian_array = np.zeros((image.shape[0]), dtype='float32') for z in range(image.shape[0]): image_slice = image[z].copy() if scale_factor != 1: image_slice = rescale(image_slice, scale_factor) abs_laplacian_array[z] = absolute_laplacian_slice(image_slice) return abs_laplacian_array @jit def absolute_laplacian_slice(image): """ Computes the Absolute Laplacian of a single image and returns value :param image: nd array of single image where to compute the absolute laplacian :return: absolute laplacian value """ n_px = (image.shape[0] - 2) * (image.shape[1] - 2) laplacian = 0 for y in range(1, image.shape[0] - 1): # skip first and last pixel in y for x in range(1, image.shape[1] - 1): # skip first and last pixel in x laplacian += abs(2 * image[y][x] - image[y][x - 1] - image[y][x + 1]) +\ abs(2 * image[y][x] - image[y - 1][x] - image[y + 1][x]) return np.float32(laplacian / n_px) @jit def quant_first_derivative(image): """ Computes the 1st derivative of each pixel and returns results as an nd array with shape(z, y-2, x-2) :param image: image where to compute the second derivatives in an image :return: nd array image where the pixels are the 1st derivative at each pixel """ # image = check_dims(image) image_derivative = np.zeros_like(image) for z in range(image.shape[0]): for y in range(1, image.shape[1] - 1): for x in range(1, image.shape[2] - 1): image_derivative[z][y][z] = abs(image[z][y][x] - image[z][y][x - 1]) + abs(image[z][y][x] - image[z][y - 1][x]) # next step reduces size of output to exclude zeros on borders which are not derivatives image_derivative = image_derivative[:, 1:(image_derivative.shape[1] - 1), 1:(image_derivative.shape[2] - 1)] return image_derivative def quant_second_derivative(image, do_gauss): """ Computes the 2nd derivative of each pixel and returns results as an nd array with shape(z, y-2, x-2) :param image: image where to compute the second derivatives in an image :return: nd array with smaller X and Y dims image where the pixels are the 2nd derivative at each pixel """ derivative_array = np.zeros((image.shape[0], image.shape[1]-2, image.shape[2]-2), dtype='float32') for z in range(image.shape[0]): image_slice = image[z].copy() if do_gauss: image_slice = gaussian(image_slice, sigma=1) derivative_array[z] = quant_second_derivative_slice(image_slice) return derivative_array @jit def quant_second_derivative_slice(image): """ Computes the 2nd derivative of each pixel and returns results as an nd array with shape(z, y-2, x-2) :param image: image where to compute the second derivatives in an image :return: nd array image where the pixels are the 2nd derivative at each pixel """ image_derivative = np.zeros_like(image) for y in range(1, image.shape[0] - 1): for x in range(1, image.shape[1] - 1): image_derivative[y][x] = (abs( 2 * image[y][x] - image[y][x - 1] - image[y][x + 1] ) + abs( 2 * image[y][x] - image[y - 1][x] - image[y + 1][x]) ) / 2 # next step reduces size of output to exclude zeros on borders which are not derivatives image_derivative = image_derivative[1:(image_derivative.shape[0] - 1), 1:(image_derivative.shape[1] - 1)] return image_derivative # @jit def total_variation(image, do_gauss): """ Computes Total Variation metric for the whole stack :param image: nd array of a 3D image where to compute the Total Variation :return: Total Variation score of image stack """ tv_array = np.zeros((image.shape[0]), dtype='float32') for z in range(image.shape[0]): image_slice = image[z].copy() if do_gauss: image_slice = gaussian(image_slice, sigma=1) tv_array[z] = total_variation_slice(image_slice) return tv_array @jit def total_variation_slice(image): """ Computes the Total Variation of each slice seperately and returns as an array :param image: nd array from 3D image :return: array with Total Variation scores for each slice """ n_px_slice = (image.shape[0] - 2) * (image.shape[1] - 2) tv = 0.0 for y in range(1, image.shape[0] - 1): for x in range(1, image.shape[1] - 1): tv += math.sqrt(pow(image[y][x + 1] - image[y][x - 1], 2) + pow(image[y + 1][x] - image[y - 1][x], 2)) return np.float32(tv / n_px_slice) def wavelet_contrast_index(image, p_low, p_high): """ Wrapper function for contrast index quantification, checking image dimensions Following the quantification for contrast using wavelet decomposition/analysis seen in Albright et al. 2023 PNAS. And defined as contrast_index = log(w_{95} / w_{50}) :param image: image to measure contrast in :param p_low: lowest percentile to normalize for contrast index :param p_high: top percentile to measure contrast, should not be sensitive to noise in image :return: contrast index """ contrast_indexes = np.zeros((image.shape[0]), dtype='float32') for z in range(image.shape[0]): image_slice = image[z].copy() contrast_indexes[z] = wavelet_contrast_index_slice(image_slice, p_low, p_high) return contrast_indexes def wavelet_contrast_index_slice(image, p_low, p_high): """ Following the quantification for contrast using wavelet decomposition/analysis seen in Albright et al. 2023 PNAS. And defined as contrast_index = log(w_{95} / w_{50}) :param image: image to measure contrast in :param p_low: lowest percentile to normalize for contrast index :param p_high: top percentile to measure contrast, should not be sensitive to noise in image :return: contrast index """ wavelet_coefficients = pywt.wavedec2(image, 'haar', level=4) coefficient_array, coefficient_slices = pywt.coeffs_to_array(wavelet_coefficients[:-1]) coefficients_absolute = np.absolute(coefficient_array) coeff_p_low = np.percentile(coefficients_absolute, p_low) coeff_p_high = np.percentile(coefficients_absolute, p_high) contrast_index = math.log(coeff_p_high / coeff_p_low) return contrast_index def percentile_contrast_index(image, p_low, p_high): """ TODO add description :param image: :param p_low: :param p_high: :return: """ assert p_low < p_high, print("Low percentile needs to be smaller than High percentile") assert 0 <= p_low <= 100.0, print("Lower percentile needs to be between [0-100]") assert 0 <= p_high <= 100.0, print("Higher percentile needs to be between [0-100]") contrast_indexes = np.zeros((image.shape[0]), dtype='float32') for z in range(image.shape[0]): image_slice = image[z] perc_low = np.percentile(image_slice, p_low) perc_high = np.percentile(image_slice, p_high) try: contrast_indexes[z] = math.log(perc_high / perc_low) except: contrast_indexes[z] = np.nan return contrast_indexes def PSNR(gt, pred, range_=4095.0): """ TODO add description Based on code from PPN2V repo for PSNR quantification :param gt: :param pred: :param range_: :return: """ assert gt.shape[0] == pred.shape[0], print("Stacks don't have the same number of z planes") psnr_indices = np.zeros((gt.shape[0]), dtype='float32') for z in range(gt.shape[0]): mse = np.mean((gt[z] - pred[z])**2) psnr = 20 * np.log10((range_)/np.sqrt(mse)) psnr_indices[z] = psnr return psnr_indices def MSE(gt, pred): """ TODO add description """ assert gt.shape == pred.shape, print("Stacks don't have the same dimensions") mse_values = np.zeros((gt.shape[0]), dtype='float32') for z in range(gt.shape[0]): mse_values[z] = np.mean((gt[z] - pred[z])**2) return mse_values def MAE(gt, pred): """ TODO add description """ assert gt.shape == pred.shape, print("Stacks don't have the same dimensions") mae_values = np.zeros((gt.shape[0]), dtype='float32') for z in range(gt.shape[0]): mae_values[z] = np.mean((gt[z] - pred[z])) return mae_values