From d269fff6c6fbbe60b56ee9bfed0266c3fede6733 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Thu, 11 Sep 2025 11:05:05 +1200 Subject: [PATCH 01/25] variance based filtering for waveforms --- src/reproimage/utils.py | 43 ++++++++++++++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 9 deletions(-) diff --git a/src/reproimage/utils.py b/src/reproimage/utils.py index b11f7e5..6bec0aa 100644 --- a/src/reproimage/utils.py +++ b/src/reproimage/utils.py @@ -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,8 @@ 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) - + interpolated_waves_np = jackknife_variance_filter(interpolated_waves_np) # filter based on variance using jackknife + # 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 +142,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 +154,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 +164,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 +345,27 @@ 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 jackknife_variance_filter(data): + converged = False + while not converged: + pop_var = np.var(data,axis=0) + subgroup_vars = [] + for i in range(data.shape[0]): + sample = np.delete(data,i, axis=0) + subgroup_vars.append(np.var(sample, axis=0)) + + subgroup_vars = np.asarray(subgroup_vars) + subgroup_vars = np.abs(pop_var-subgroup_vars) + subgroup_vars = np.mean(subgroup_vars, axis=1) + + feature = int(np.argmax(subgroup_vars)) + + var_without_feature = np.mean(np.var(np.delete(data, feature, axis=0), axis=0)) + if subgroup_vars[feature] >= var_without_feature * 1.5: + data = np.delete(data, feature, axis=0) + else: + converged = True + return data + From 2eeda203b789c0e987bd6438d163dfa25c812de6 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Fri, 19 Sep 2025 16:14:47 +1200 Subject: [PATCH 02/25] utility functions for mosaic creation --- src/reproimage/utils.py | 42 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/src/reproimage/utils.py b/src/reproimage/utils.py index 6bec0aa..cf06306 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 @@ -369,3 +369,43 @@ def jackknife_variance_filter(data): converged = True return data +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 + # Do normal exception handling \ No newline at end of file From f821385c5b1df1274a205046c7506ee2a429581b Mon Sep 17 00:00:00 2001 From: tjac799 Date: Fri, 19 Sep 2025 16:16:30 +1200 Subject: [PATCH 03/25] mosaic creation and downsampling for large image stacks --- src/reproimage/image_operations.py | 265 ++++++++++++++++++++++++++++- 1 file changed, 263 insertions(+), 2 deletions(-) diff --git a/src/reproimage/image_operations.py b/src/reproimage/image_operations.py index e87ed5f..0108fb8 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([raw_im_stack_spacing[0] * isotropic_downsample_dim]) * 3 + 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 From ebd31ab26e579d05c9db88ca8670cfd264ce59a5 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Mon, 29 Sep 2025 14:20:48 +1300 Subject: [PATCH 04/25] Readme changes regarding developer installation --- README.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) 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. From eaae91fd3746ead5f107aaaf55987198178410e4 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Mon, 29 Sep 2025 14:34:35 +1300 Subject: [PATCH 05/25] utility functions for mosaic creation --- src/reproimage/utils.py | 45 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 43 insertions(+), 2 deletions(-) diff --git a/src/reproimage/utils.py b/src/reproimage/utils.py index b11f7e5..d4e293d 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 @@ -343,4 +343,45 @@ 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 + # Do normal exception handling \ No newline at end of file From b6afee66493fe20c6390aa0b707e3c5192e2db99 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Fri, 19 Sep 2025 16:16:30 +1200 Subject: [PATCH 06/25] mosaic creation and downsampling for large image stacks --- src/reproimage/image_operations.py | 265 ++++++++++++++++++++++++++++- 1 file changed, 263 insertions(+), 2 deletions(-) diff --git a/src/reproimage/image_operations.py b/src/reproimage/image_operations.py index e87ed5f..0108fb8 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([raw_im_stack_spacing[0] * isotropic_downsample_dim]) * 3 + 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 From 8aa446639d79c33842a6f56c1e4d6f31edeed102 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Mon, 29 Sep 2025 14:20:48 +1300 Subject: [PATCH 07/25] Readme changes regarding developer installation --- README.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) 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. From 739d1068085816ad1ef7076785ec9c88546c6958 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Tue, 30 Sep 2025 16:47:09 +1300 Subject: [PATCH 08/25] functions for creation of sub-trees --- src/reproimage/Graph.py | 24 +++++++++++++++++++++++- src/reproimage/utils.py | 5 ++++- 2 files changed, 27 insertions(+), 2 deletions(-) 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/utils.py b/src/reproimage/utils.py index 1a487ca..7953b26 100644 --- a/src/reproimage/utils.py +++ b/src/reproimage/utils.py @@ -385,4 +385,7 @@ def __exit__(self, type, value, traceback): sys.stdout = self.stdout if type is not None: pass - # Do normal exception handling \ No newline at end of file + +def invert_dictionary(dictionary): + inverted_map = {v: k for k, v in dictionary.items()} + return inverted_map From 83853f3e6cc254ec531eff63ce119651ef9a0a9a Mon Sep 17 00:00:00 2001 From: tjac799 Date: Tue, 30 Sep 2025 17:09:19 +1300 Subject: [PATCH 09/25] modification to suppresor class so it works --- src/reproimage/utils.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/reproimage/utils.py b/src/reproimage/utils.py index 7953b26..60619f7 100644 --- a/src/reproimage/utils.py +++ b/src/reproimage/utils.py @@ -386,6 +386,12 @@ def __exit__(self, type, value, traceback): 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 From c973abdd1c00cc0d2829cdb63ec3d63847aaf6e9 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Thu, 27 Nov 2025 15:55:13 +1300 Subject: [PATCH 10/25] generalises handling of anisotropy for anisotropic im stacks with no log file --- src/reproimage/image_operations.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/reproimage/image_operations.py b/src/reproimage/image_operations.py index 0108fb8..6ec2ea5 100644 --- a/src/reproimage/image_operations.py +++ b/src/reproimage/image_operations.py @@ -117,7 +117,6 @@ def n_largest_components_filter(img, n_components = 0): 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): @@ -174,7 +173,7 @@ def image_stack_to_volume(reconstruction_directory, output_directory, sample_ide number_of_stacks_in_mosaic_layer = ceil(n_images / resamplinglayer_size) stack_remainder = n_images % resamplinglayer_size - output_spacing = tuple([raw_im_stack_spacing[0] * isotropic_downsample_dim]) * 3 + output_spacing = tuple([dim_spacing * isotropic_downsample_dim for dim_spacing in raw_im_stack_spacing]) * 3 mosaic_piece_size = [0.0, 0.0, 0.0] output_size_ds = [0, 0, 0] From d26dbe09926c0128e238cdb0fd47e08a75338010 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Thu, 27 Nov 2025 16:08:37 +1300 Subject: [PATCH 11/25] fix --- src/reproimage/image_operations.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/reproimage/image_operations.py b/src/reproimage/image_operations.py index 6ec2ea5..a3b0be3 100644 --- a/src/reproimage/image_operations.py +++ b/src/reproimage/image_operations.py @@ -173,7 +173,8 @@ def image_stack_to_volume(reconstruction_directory, output_directory, sample_ide 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]) * 3 + 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] From 5353033f3d7f713452037aa55180a46560890abe Mon Sep 17 00:00:00 2001 From: tjac799 Date: Wed, 14 May 2025 12:49:04 +1200 Subject: [PATCH 12/25] specifying pyvista version in pyproject.toml --- pyproject.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 346b2f7..eb5b76c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,5 +9,4 @@ description = "A package for processing of placental images for the development readme = "README.md" requires-python = ">=3.9" classifiers = ["Programming Language :: Python :: 3", "License :: OSI Approved :: APACHE License", "Operating System :: OS Independent",] -dependencies = ["SimpleITK", "networkx", "scipy", "scikit-image", "numba", "pyvista"] - +dependencies = ["SimpleITK", "networkx", "scipy", "scikit-image", "numba", "pyvista>=0.43"] \ No newline at end of file From 6d091f15825179fdec2067fb5a071d05f9364b94 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Fri, 11 Jul 2025 14:39:03 +1200 Subject: [PATCH 13/25] dependency update --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index eb5b76c..1fd8c44 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,4 +9,4 @@ description = "A package for processing of placental images for the development readme = "README.md" requires-python = ">=3.9" classifiers = ["Programming Language :: Python :: 3", "License :: OSI Approved :: APACHE License", "Operating System :: OS Independent",] -dependencies = ["SimpleITK", "networkx", "scipy", "scikit-image", "numba", "pyvista>=0.43"] \ No newline at end of file +dependencies = ["SimpleITK", "networkx", "scipy", "scikit-image", "numba", "pyvista>=0.43", "pandas"] \ No newline at end of file From d9a9911d59db4a2c7681045158380d279f4a9bb5 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Fri, 11 Jul 2025 14:40:44 +1200 Subject: [PATCH 14/25] update query wrapper to accept local scope of calling frame --- src/reproimage/utils.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/reproimage/utils.py b/src/reproimage/utils.py index 99d8c2a..2fc3fb0 100644 --- a/src/reproimage/utils.py +++ b/src/reproimage/utils.py @@ -1,3 +1,5 @@ +import inspect + import numpy as np from pathlib import Path from scipy.signal import find_peaks @@ -81,9 +83,10 @@ def mean_wave(x_values, y_values, verbose=False): ## Want to find the last local minima that has occured before a main peak trough_indices = [] for peak in peak_indices: - trough_loc = np.where(all_trough_indices Date: Tue, 5 Aug 2025 14:22:56 +1200 Subject: [PATCH 15/25] init py file causing issues with uninstalling --- src/__init__.py | 1 - 1 file changed, 1 deletion(-) delete mode 100644 src/__init__.py diff --git a/src/__init__.py b/src/__init__.py deleted file mode 100644 index 8d1c8b6..0000000 --- a/src/__init__.py +++ /dev/null @@ -1 +0,0 @@ - From a9f45a31af837d0add9bd1cf2f84433f8f7c39f3 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Tue, 5 Aug 2025 14:23:11 +1200 Subject: [PATCH 16/25] whitespace adjustment --- src/reproimage/utils.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/reproimage/utils.py b/src/reproimage/utils.py index 2fc3fb0..b11f7e5 100644 --- a/src/reproimage/utils.py +++ b/src/reproimage/utils.py @@ -168,9 +168,6 @@ def mean_wave(x_values, y_values, verbose=False): 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 " From 34823def63e6be217af4eaa2599fd1c9d3c19de7 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Tue, 5 Aug 2025 14:26:35 +1200 Subject: [PATCH 17/25] visualising a single graph with no field --- src/reproimage/Visualisation.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/reproimage/Visualisation.py b/src/reproimage/Visualisation.py index 5046ed9..5d729ed 100644 --- a/src/reproimage/Visualisation.py +++ b/src/reproimage/Visualisation.py @@ -1,6 +1,6 @@ import numpy as np import pyvista as pv - +import networkx as nx def remap_node_field_for_vis(graph, field): """ graph - nx.Graph @@ -34,7 +34,7 @@ def generate_visualisation_arrays(coords, edges): connections_with_padding = np.vstack((padding, remapped_connections.T)).T return new_coords, connections_with_padding -def visualise_graph_and_field(graph, coords, field, field_name='radii', title='', need_remap = True): +def visualise_graph_and_field(graph: nx.graph, coords: np.array, field, field_name='radii', title='', need_remap = True): if title == '': title = f'{field_name} visualisation' if need_remap: @@ -50,4 +50,13 @@ def visualise_graph_and_field(graph, coords, field, field_name='radii', title='' plotter.add_mesh(pod_tube, render_lines_as_tubes=True, show_scalar_bar=False) plotter.add_scalar_bar(field_name, position_x=0.25) plotter.camera_position = 'xz' + plotter.show() + +def visualise_graph(graph: nx.graph, coords: np.array, radius: float = 1.0): + vis_coords, vis_connections = generate_visualisation_arrays(coords, + np.array(graph.edges())) + plotter = pv.Plotter() + pod = pv.PolyData(vis_coords, lines=vis_connections, n_lines=vis_connections.shape[0]) + pod_tube = pod.tube(radius=radius) + plotter.add_mesh(pod_tube, render_lines_as_tubes=True) plotter.show() \ No newline at end of file From cfa3ecb5e9265221dd39e1802c06e3744bcd1743 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Mon, 29 Sep 2025 14:34:35 +1300 Subject: [PATCH 18/25] utility functions for mosaic creation --- src/reproimage/utils.py | 45 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 43 insertions(+), 2 deletions(-) diff --git a/src/reproimage/utils.py b/src/reproimage/utils.py index b11f7e5..d4e293d 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 @@ -343,4 +343,45 @@ 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 + # Do normal exception handling \ No newline at end of file From 96a97b7b61ea932b72742c02e503e11e945f8b65 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Fri, 19 Sep 2025 16:16:30 +1200 Subject: [PATCH 19/25] mosaic creation and downsampling for large image stacks --- src/reproimage/image_operations.py | 265 ++++++++++++++++++++++++++++- 1 file changed, 263 insertions(+), 2 deletions(-) diff --git a/src/reproimage/image_operations.py b/src/reproimage/image_operations.py index e87ed5f..0108fb8 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([raw_im_stack_spacing[0] * isotropic_downsample_dim]) * 3 + 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 From 54c7065a5e31e62df3640224cfab2fae9b6593f3 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Mon, 29 Sep 2025 14:20:48 +1300 Subject: [PATCH 20/25] Readme changes regarding developer installation --- README.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) 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. From 514172a98fa2d06d46a20000cdafdf964b6820c6 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Mon, 29 Sep 2025 14:20:48 +1300 Subject: [PATCH 21/25] Readme changes regarding developer installation From dd907ca6d0f3065f7868ee7648cbbf2b5f4ba8ed Mon Sep 17 00:00:00 2001 From: tjac799 Date: Tue, 30 Sep 2025 16:47:09 +1300 Subject: [PATCH 22/25] functions for creation of sub-trees --- src/reproimage/Graph.py | 24 +++++++++++++++++++++++- src/reproimage/utils.py | 5 ++++- 2 files changed, 27 insertions(+), 2 deletions(-) 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/utils.py b/src/reproimage/utils.py index d4e293d..f568b2c 100644 --- a/src/reproimage/utils.py +++ b/src/reproimage/utils.py @@ -384,4 +384,7 @@ def __exit__(self, type, value, traceback): sys.stdout = self.stdout if type is not None: pass - # Do normal exception handling \ No newline at end of file + +def invert_dictionary(dictionary): + inverted_map = {v: k for k, v in dictionary.items()} + return inverted_map From db50f0a9f6dfda25095ae640b4ae557c3749d999 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Tue, 30 Sep 2025 17:09:19 +1300 Subject: [PATCH 23/25] modification to suppresor class so it works --- src/reproimage/utils.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/reproimage/utils.py b/src/reproimage/utils.py index f568b2c..f8d4cad 100644 --- a/src/reproimage/utils.py +++ b/src/reproimage/utils.py @@ -385,6 +385,12 @@ def __exit__(self, type, value, traceback): 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 From 87b241187712b2d83b717d80d829934015068a78 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Thu, 27 Nov 2025 15:55:13 +1300 Subject: [PATCH 24/25] generalises handling of anisotropy for anisotropic im stacks with no log file --- src/reproimage/image_operations.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/reproimage/image_operations.py b/src/reproimage/image_operations.py index 0108fb8..6ec2ea5 100644 --- a/src/reproimage/image_operations.py +++ b/src/reproimage/image_operations.py @@ -117,7 +117,6 @@ def n_largest_components_filter(img, n_components = 0): 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): @@ -174,7 +173,7 @@ def image_stack_to_volume(reconstruction_directory, output_directory, sample_ide number_of_stacks_in_mosaic_layer = ceil(n_images / resamplinglayer_size) stack_remainder = n_images % resamplinglayer_size - output_spacing = tuple([raw_im_stack_spacing[0] * isotropic_downsample_dim]) * 3 + output_spacing = tuple([dim_spacing * isotropic_downsample_dim for dim_spacing in raw_im_stack_spacing]) * 3 mosaic_piece_size = [0.0, 0.0, 0.0] output_size_ds = [0, 0, 0] From 9664aec4fcb40e612d8c110218c0c5d80d9a0627 Mon Sep 17 00:00:00 2001 From: tjac799 Date: Thu, 27 Nov 2025 16:08:37 +1300 Subject: [PATCH 25/25] fix --- src/reproimage/image_operations.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/reproimage/image_operations.py b/src/reproimage/image_operations.py index 6ec2ea5..a3b0be3 100644 --- a/src/reproimage/image_operations.py +++ b/src/reproimage/image_operations.py @@ -173,7 +173,8 @@ def image_stack_to_volume(reconstruction_directory, output_directory, sample_ide 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]) * 3 + 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]