diff --git a/README.md b/README.md index 523e8de..0bbd911 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,8 @@ Use pip to install the libraries ``` pip install git+https://github.com/virtualpregnancy/repro-image.git ``` -To use pip to install the libraries when you have cloned them onto your machine (note this is best to do within a virtual environment) +To use pip to install the libraries when you have cloned them onto your machine and are planning to develop them +(note this is best to do within a virtual environment): ``` pip install -e /path/to/reproimage ``` @@ -15,4 +16,5 @@ If you are running a system with multiple python versions, it can be safest to r ``` python -m pip ``` -To ensure that your pip command is explicitly tied to the correct python interpreter. +To ensure that your pip command is explicitly tied to the correct python interpreter. This is implicitly true if you are +already running python from a virtual environment that has been activated. diff --git a/src/reproimage/Graph.py b/src/reproimage/Graph.py index dd71058..4935fb3 100644 --- a/src/reproimage/Graph.py +++ b/src/reproimage/Graph.py @@ -24,6 +24,8 @@ def junction_node_subgraph(graph): :return: a series reduced tree with no degree 2 nodes """ jgraph = copy.deepcopy(graph) + print(nx.is_frozen(jgraph)) + jgraph = nx.Graph(jgraph) # Iterate through the edges and remove attributes for u, v, attrs in jgraph.edges(data=True): for attr_key in list(attrs.keys()): @@ -424,4 +426,24 @@ def get_num_junction_nodes(graph): for node in graph.nodes: if deg[node] > 2: degree_node_set.append(node) - return len(degree_node_set), degree_node_set \ No newline at end of file + return len(degree_node_set), degree_node_set + +def get_points_from_exnode(exnode_obj): + dimension = 3 + + points = np.zeros([exnode_obj.num_nodes+1, dimension]) + print(f"Processing nodes") + for node in exnode_obj.sections[0].nodes: + points[node.number, :] = node.values + return points + +def edges_from_exelem(exelem_obj): + edges = list() + print(f"Processing elements") + for element in exelem_obj.elements: + edges.append(element.nodes) + return edges + +def direct_UAG_from_inlet(UAG, inlet): + directed_graph = nx.bfs_tree(UAG, inlet) + return directed_graph \ No newline at end of file diff --git a/src/reproimage/image_operations.py b/src/reproimage/image_operations.py index e87ed5f..a3b0be3 100644 --- a/src/reproimage/image_operations.py +++ b/src/reproimage/image_operations.py @@ -1,7 +1,9 @@ -from reproimage.utils import cubify, calculate_fractal_dimension +from reproimage.utils import cubify, calculate_fractal_dimension, get_image_metadata_from_log_file, Suppressor import SimpleITK as sitk from scipy.stats import variation import numpy as np +import re, os, sys +from math import ceil def efficient_largest_ccmp_filter(image: sitk.Image): """ @@ -113,4 +115,263 @@ def n_largest_components_filter(img, n_components = 0): largest_comp = sitk.Threshold(ccmp_sort, upper=float(n_components), outsideValue=float(0)) largest_comp = sitk.Cast(largest_comp, sitk.sitkUInt8) - return largest_comp \ No newline at end of file + return largest_comp + +def image_stack_to_volume(reconstruction_directory, output_directory, sample_identifier, isotropic_downsample_dim=8, + grid_size=[8, 8, 8], specified_mosaic_piece_size=[], raw_im_stack_size=[], raw_im_stack_spacing=[]): + if not os.path.isdir(output_directory): + print(f"Error code 2: Output directory provided is not a real directory ({output_directory})") + exit(2) + + # script parameters + grid_x, grid_y, grid_z = grid_size + + root_dir = reconstruction_directory + f_name = output_directory + f_identifier = sample_identifier + + grid_map = sitk.Image([grid_z, grid_x, grid_y], sitk.sitkUInt16) + + if root_dir == '': + print("Script terminated as no directory was provided.") + exit(2) + + # calculate image parameters for the output image + try: + log_file = os.path.join(root_dir, + [f for f in os.listdir(root_dir) if f.endswith('rec.log') and not f.startswith('._')][ + 0]) + raw_im_stack_size, raw_im_stack_spacing, fileformat = get_image_metadata_from_log_file(log_file) + n_images = raw_im_stack_size[2] + except IndexError: + print( + "No log file available, script will run with no spatial information, and with the best possible approximation of stack size using the files contained within your chosen directory") + n_images = len(os.listdir(root_dir)) + if raw_im_stack_spacing == []: + raw_im_stack_spacing = [1, 1, 1] + # image size parameters based off first image in the stack + im_address_for_params = os.path.join(root_dir, os.listdir(root_dir)[-1]) + with Suppressor(): + im_for_params = sitk.ReadImage(im_address_for_params) + if raw_im_stack_size == []: + raw_im_stack_size = list(im_for_params.GetSize()) + raw_im_stack_size.append(n_images) + fileformat = '.' + os.path.split(im_address_for_params)[-1].split('.')[-1] + + # create list of filenames within target directory + filelist = [f for f in os.listdir(root_dir) if f.endswith(fileformat)] # tif code + filelist.sort(key=lambda f: int(re.sub('\D', '', f)), reverse=True) + filelist = list(filter(lambda f: not (f.endswith('spr' + fileformat)), filelist)) # tif code + filelist = list(filter(lambda f: not ('pp' in f), filelist)) # tif code + filelist = list(filter(lambda f: not ('prev' in f), filelist)) # tif code + + if len(filelist) == 0: + print("There are no suitable files in this directory") + exit(1) + + resamplinglayer_size = isotropic_downsample_dim + number_of_stacks_in_mosaic_layer = ceil(n_images / resamplinglayer_size) + stack_remainder = n_images % resamplinglayer_size + + output_spacing = tuple([dim_spacing * isotropic_downsample_dim for dim_spacing in raw_im_stack_spacing]) + + mosaic_piece_size = [0.0, 0.0, 0.0] + output_size_ds = [0, 0, 0] + + output_size_ds[0] = int(raw_im_stack_size[0] * raw_im_stack_spacing[0] / output_spacing[0] + .5) + output_size_ds[1] = int(raw_im_stack_size[1] * raw_im_stack_spacing[1] / output_spacing[1] + .5) + output_size_ds[2] = int(raw_im_stack_size[2] * raw_im_stack_spacing[2] / output_spacing[2] + .5) + + if specified_mosaic_piece_size != []: + output_size_ds[0] = output_size_ds[0] + ( + specified_mosaic_piece_size[0] - output_size_ds[0] % specified_mosaic_piece_size[0]) + output_size_ds[1] = output_size_ds[1] + ( + specified_mosaic_piece_size[1] - output_size_ds[1] % specified_mosaic_piece_size[1]) + output_size_ds[2] = output_size_ds[2] + ( + specified_mosaic_piece_size[2] - output_size_ds[2] % specified_mosaic_piece_size[2]) + + mosaic_piece_size[0] = int(output_size_ds[0] / grid_x) + mosaic_piece_size[1] = int(output_size_ds[1] / grid_y) + mosaic_piece_size[2] = int(output_size_ds[2] / grid_z) + im_size = os.stat(os.path.join(root_dir, filelist[1])).st_size # size of one full size .bmp image in bytes + + # resampling summary info + print(f'There are {n_images} \'.{fileformat}\' files in the reconstruction directory') + print(f'The entire image stack is {round(im_size * n_images / 1024 ** 3, 2)} GB') + print( + f'The stack is composed of {n_images} images that are {raw_im_stack_size[0]} by {raw_im_stack_size[1]} pixels, with spacing of {raw_im_stack_spacing} µm') + + # create output image + result_img = sitk.Image(mosaic_piece_size, 1, 1) # SIZE, pixelID, numberOfComponentsPerPixel + result_img.SetSpacing(output_spacing) + result_img.SetOrigin([0.0, 0.0, 0.0]) + result_img.SetDirection((1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) # identity matrix + + ds_size = [output_size_ds[0], output_size_ds[1], mosaic_piece_size[2]] + + # initialise resample image filter + resample = sitk.ResampleImageFilter() + resample.SetOutputDirection((1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) # Identity matrix + resample.SetOutputSpacing(output_spacing) + resample.SetOutputOrigin([0.0, 0.0, 0.0]) + resample.SetTransform(sitk.Transform()) + resample.SetInterpolator(sitk.sitkNearestNeighbor) + + extract = sitk.ExtractImageFilter() + Reader = sitk.ImageSeriesReader() + #######TEMP + debug = False + origins = [] + lin_idx = 0 + global_z_layer = 0 # work variable to keep track of what z slice the resampling algorithm is up to, specifically for when the resampled image stacks are added into the output image dataframe + + # resample the original images stack by stack + print("Resampling image data") + print( + f"The raw reconstruction image stack will be sliced into stacks of {resamplinglayer_size} images for the resampling process") + print("Processing . . . ") + + for z_grid_index in range(0, grid_z): + + downsample_img = sitk.Image(ds_size, 1, 1) # SIZE, pixelID, numberOfComponentsPerPixel + downsample_img.SetSpacing(output_spacing) + downsample_img.SetOrigin([0.0, 0.0, 0.0]) + downsample_img.SetDirection((1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) # identity matrix + + images_left_in_stack = len(filelist) + if images_left_in_stack > mosaic_piece_size[2] * isotropic_downsample_dim: + images_remaining = mosaic_piece_size[2] * isotropic_downsample_dim # -skip_end# n_images + else: + images_remaining = images_left_in_stack + number_of_stacks_in_mosaic_layer = mosaic_piece_size[2] * isotropic_downsample_dim / resamplinglayer_size + stack_number = 1 + + while images_remaining > 0: # this while loop loads in the entirety of the mosaic layer and either resamples it + # or does not, depending on the resmapling factor chosesn in the function execution + if images_remaining < resamplinglayer_size: + resamplinglayer_size = images_remaining + images_remaining -= resamplinglayer_size + if debug: + print('Images remaining', images_remaining) + # pop names off the sorted list of image filenames, used to parameterize the reader object + names = [] # clear the names list + for n in range(0, resamplinglayer_size): # each images i na stack + names.append(os.path.join(root_dir, filelist.pop(-1)).replace("\\", "/")) + if debug: + print('Names', names) + Reader.SetFileNames(names) + with Suppressor(): + temp_image = Reader.Execute() # read in stack of images + + if temp_image.GetPixelIDValue() != sitk.sitkUInt8: + if temp_image.GetPixelIDValue() == 13: + temp_image = sitk.VectorIndexSelectionCast(temp_image, 0, 1) + + if temp_image.GetPixelIDValue() != sitk.sitkUInt8: + temp_image = sitk.RescaleIntensity(temp_image, 0, 255) + temp_image = sitk.Cast(temp_image, sitk.sitkUInt8) + + temp_image.SetSpacing( + raw_im_stack_spacing) # bitmaps have no inherent spatial information, so it needs to be set using data read from the log file\ + + if isotropic_downsample_dim > 1: # this section encapsulates the bulk of what is needed to do the image resample if not mosaicing + res_stack_size = ds_size[:2] + res_stack_size.append(int(temp_image.GetSize()[2] * raw_im_stack_spacing[2] / output_spacing[2] + .5)) + resample.SetSize(res_stack_size) + temp_downsize = resample.Execute(temp_image) + + # write downsampled stack into output image + pasteFilter = sitk.PasteImageFilter() + pasteFilter.SetSourceSize(temp_downsize.GetSize()) + pasteFilter.SetSourceIndex([0, 0, 0]) + pasteFilter.SetDestinationIndex( + [0, 0, stack_number - 1]) # use of the global z layer should only happen when grid size is 1 + + # write resampled sub-stack into the output image + downsample_img = pasteFilter.Execute(downsample_img, temp_downsize) + else: + # just write downsampled stack into output image + pasteFilter = sitk.PasteImageFilter() + pasteFilter.SetSourceSize(temp_image.GetSize()) + pasteFilter.SetSourceIndex([0, 0, 0]) + pasteFilter.SetDestinationIndex([0, 0, stack_number - 1]) + downsample_img = pasteFilter.Execute(downsample_img, temp_image) + + global_z_layer += int(resamplinglayer_size / isotropic_downsample_dim) # update current z-slice level + + print( + f"\rProcessed Image stack {stack_number - 1}/{number_of_stacks_in_mosaic_layer}, in layer: {z_grid_index + 1}/{grid_z}, {lin_idx + 1}/{grid_x * grid_y * grid_z} pieces have been written to file.", + end='') + sys.stdout.flush() + + stack_number += 1 # + + start_x = 0 + for i in range(0, grid_x): + start_y = 0 + for j in range(0, grid_y): + # lin_idx programmed to incremented by grid_size**2 as we move down each mosaic level down the z-index + lin_idx = z_grid_index * (grid_x * grid_y) + i * (grid_y) ** 1 + j * (grid_y) ** 0 + grid_map[z_grid_index, i, j] = lin_idx + # Select same subregion using ExtractImageFilter + extract.SetSize([mosaic_piece_size[0], mosaic_piece_size[1], mosaic_piece_size[2]]) + extract.SetIndex([start_x, start_y, 0]) + # print([start_x, start_y, global_layer]) + origins.append((lin_idx, [start_x * output_spacing[0], start_y * output_spacing[1], + z_grid_index * number_of_stacks_in_mosaic_layer * output_spacing[2]])) + + # write resampled sub-stack into the output image + result_img = extract.Execute(downsample_img) + + # write resampled image to file, TODO: improve readability of origin handling + result_img.SetMetaData("xyzt_units", "3") + result_img.SetMetaData("origin", str([start_x, start_y, + z_grid_index * number_of_stacks_in_mosaic_layer * isotropic_downsample_dim])) + result_img.SetOrigin([start_x * output_spacing[0], start_y * output_spacing[1], + z_grid_index * number_of_stacks_in_mosaic_layer * output_spacing[2]]) + result_img.SetSpacing(output_spacing) + # statistics_image_filter.Execute(result_img) + # print(statistics_image_filter.GetMaximum()) + write = True + if write: + try: + append = f_identifier + im_f_name = os.path.join(f_name, str(lin_idx) + '_' + append + '.nii') + sitk.WriteImage(result_img, im_f_name) + print( + f"\rProcessed Image stack {stack_number - 1}/{number_of_stacks_in_mosaic_layer}, in layer: {z_grid_index + 1}/{grid_z}, {lin_idx + 1}/{grid_x * grid_y * grid_z} pieces have been written to file.", + end='') + sys.stdout.flush() + except RuntimeError: + append = str(lin_idx) + '_' + print("Filename provided is invalid") + print(f"Saving as {lin_idx}_resampled_stack.nii") + im_f_name = append + 'resampled_stack.nii' + sitk.WriteImage(result_img, im_f_name) + + start_y = start_y + mosaic_piece_size[1] + start_x = start_x + mosaic_piece_size[0] + if grid_size != [1, 1, 1]: + with open(os.path.join(output_directory, f_identifier + '_mosaic_mapping.txt'), 'w') as f: + for mosaic_piece in origins: + f.write(f'{mosaic_piece[0]}, {mosaic_piece[1]}\n') + f.write(f'mosaic piece spacing, {output_spacing}\n') + f.write(f'total mosaic size, {output_size_ds}\n') + f.write(f'mosiac piece size, {mosaic_piece_size}\n') + f.write(f'Grid size, {grid_size}\n') + f.write(f'identifier,{f_identifier}\n') + f.write(f'original reconstruction location, {root_dir}\n') + f.write(f'downsampling factor, {isotropic_downsample_dim}\n') + else: + with open(os.path.join(output_directory, f_identifier + '_resampling_summary.txt'), 'w') as f: + f.write(f'resampled image spacing, {output_spacing}\n') + f.write(f'resampled image size, {output_size_ds}\n') + f.write(f'Grid size, {grid_size}\n') + f.write(f'identifier,{f_identifier}\n') + f.write(f'original reconstruction location, {root_dir}\n') + f.write(f'downsampling factor, {isotropic_downsample_dim}\n') + + +def downsample(reconstruction_directory, output_directory, sample_identifier, isotropic_downsample_dim): + return image_stack_to_volume(reconstruction_directory, output_directory, sample_identifier, isotropic_downsample_dim, + grid_size=[1, 1, 1]) \ No newline at end of file diff --git a/src/reproimage/utils.py b/src/reproimage/utils.py index b11f7e5..60619f7 100644 --- a/src/reproimage/utils.py +++ b/src/reproimage/utils.py @@ -1,5 +1,5 @@ import inspect - +import sys import numpy as np from pathlib import Path from scipy.signal import find_peaks @@ -88,7 +88,6 @@ def mean_wave(x_values, y_values, verbose=False): trough_loc = trough_loc[-1] trough_indices.append(all_trough_indices[trough_loc]) - interpolated_waves = [] if verbose: plt.figure(figsize=(10, 6)) @@ -128,7 +127,7 @@ def mean_wave(x_values, y_values, verbose=False): # Convert the list of interpolated waves to a NumPy array for calculations interpolated_waves_np = np.vstack(interpolated_waves) - + # subsampling # Calculate the initial average and standard deviation average_wave = np.mean(interpolated_waves_np, axis=0) std_wave = np.std(interpolated_waves_np, axis=0) @@ -142,8 +141,8 @@ def mean_wave(x_values, y_values, verbose=False): for wave_index, waveform in enumerate(interpolated_waves): plt.plot(x_common, waveform, label=f"Waveform {wave_index}") plt.plot(x_common, average_wave, label='Average wave', linestyle='-.') - plt.plot(x_common, average_wave+std_wave, label='Average wave + SD', linestyle='--') - plt.plot(x_common, average_wave-std_wave, label='Average wave - SD', linestyle='--') + plt.plot(x_common, average_wave+0.2*amplitude_of_ave, label='Average wave + ampl/5', linestyle='--') + plt.plot(x_common, average_wave-0.2*amplitude_of_ave, label='Average wave - ampl/5', linestyle='--') plt.legend() plt.grid(True) plt.show() @@ -154,6 +153,7 @@ def mean_wave(x_values, y_values, verbose=False): filtered_waves = [] excluded_waves = [] count_excluded = 0 + for wave in interpolated_waves_np: # Calculate the percentage of points that meet the OR condition within_range = (wave >= (average_wave - 0.2*amplitude_of_ave)) & ( @@ -163,19 +163,20 @@ def mean_wave(x_values, y_values, verbose=False): if percentage_within_range >= threshold_percentage: filtered_waves.append(wave) else: - count_excluded =+ 1 + count_excluded += 1 excluded_waves.append(wave) if verbose: print("Waves filtered, num excluded", count_excluded) if verbose: - print(f"{len(interpolated_waves)} waveforms included in the calculationg for the average waveform," - f"using a cutoff proportion of {threshold_percentage} % for points within one standard deviation of the " - f"raw native waveform") + print(f"{len(filtered_waves)} waveforms included in the Calculation for the average waveform," + f"using a cutoff proportion of {threshold_percentage} % for points within 20% of the mean waveform with " + f"average waveform amplitude as the distance raw native waveform") # Recalculate the average and standard deviation with the filtered waves filtered_waves_np = np.vstack(filtered_waves) new_average_wave = np.mean(filtered_waves_np, axis=0) new_std_wave = np.std(filtered_waves_np, axis=0) + if verbose: plt.figure(figsize=(10, 6)) plt.title("Average waveform") @@ -343,4 +344,54 @@ def Generate_timeDelta(self, delta_key, delta_value): "Your units must be one of 'day', 'days', 'month', 'months', 'year', 'years'") def to_numpy(self): - return self.dataframe.to_numpy() \ No newline at end of file + return self.dataframe.to_numpy() + +def get_image_metadata_from_log_file(file): + """ + This function takes in the file pointer for a reconstruction log file, reads the file n pixels in x,y,z directions, + and reads the spatial data for the reconstructed images + """ + with open(file) as f: + im_size = [0, 0, 0] + spacing = [] + for x in f: + if "Result Image Width" in x: + im_size[0] = int(x.split('=')[-1]) + elif "Result Image Height" in x: + im_size[1] = int(x.split('=')[-1]) + elif "Sections Count" in x: + im_size[2] = int(x.split('=')[-1]) + elif "Image Pixel Size" in x: + spacing = tuple([float( + x.split("=")[-1])]) * 3 # known bug, this line picks up first camera pixel, then actual pixel size + elif "Result File Type" in x: + formatString = x.split("=")[-1] + if "TIFF" in formatString or "TIF" in formatString: + fileformat = ".tif" + elif "BMP" in formatString: + fileformat = ".bmp" + else: + print(f"Error: no database match for {formatString}") + + return im_size, spacing, fileformat + +class Suppressor(object): + + def __enter__(self): + self.stdout = sys.stdout + sys.stdout = self + + def __exit__(self, type, value, traceback): + sys.stdout = self.stdout + if type is not None: + pass + + def write(self, _): + pass # ignore all writes + + def flush(self): + pass # required by some code that calls flush + +def invert_dictionary(dictionary): + inverted_map = {v: k for k, v in dictionary.items()} + return inverted_map