From 85b2e794f0feef6e822d099c832465956965f34f Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Thu, 17 Apr 2025 12:51:10 +0200 Subject: [PATCH 01/23] generate deformed test image --- .gitignore | 5 + examples/.gitattributes | 1 + .../data/platy1_muscles_stardist_fixed.tif | 3 + .../data/platy1_muscles_stardist_moving.tif | 3 + examples/generate_test_data.py | 113 ++++++++++++++++++ 5 files changed, 125 insertions(+) create mode 100644 examples/.gitattributes create mode 100644 examples/data/platy1_muscles_stardist_fixed.tif create mode 100644 examples/data/platy1_muscles_stardist_moving.tif create mode 100644 examples/generate_test_data.py diff --git a/.gitignore b/.gitignore index 0a19790..54a3041 100644 --- a/.gitignore +++ b/.gitignore @@ -172,3 +172,8 @@ cython_debug/ # PyPI configuration file .pypirc + +# huge data +./examples/data/platy1_muscles_stardist.tif + +.DS_Store \ No newline at end of file diff --git a/examples/.gitattributes b/examples/.gitattributes new file mode 100644 index 0000000..d92417a --- /dev/null +++ b/examples/.gitattributes @@ -0,0 +1 @@ +*.tif filter=lfs diff=lfs merge=lfs -text diff --git a/examples/data/platy1_muscles_stardist_fixed.tif b/examples/data/platy1_muscles_stardist_fixed.tif new file mode 100644 index 0000000..d43e439 --- /dev/null +++ b/examples/data/platy1_muscles_stardist_fixed.tif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:cbbf90829a7830ec3cc57c048661ecddf864377b3c1c0fccbc15794b2f64aead +size 163013402 diff --git a/examples/data/platy1_muscles_stardist_moving.tif b/examples/data/platy1_muscles_stardist_moving.tif new file mode 100644 index 0000000..43600dc --- /dev/null +++ b/examples/data/platy1_muscles_stardist_moving.tif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f4157af6f678cdcee19d4cff13445d760aa9fdf540f9ededc5058c9c28d0b3c4 +size 114242416 diff --git a/examples/generate_test_data.py b/examples/generate_test_data.py new file mode 100644 index 0000000..c549973 --- /dev/null +++ b/examples/generate_test_data.py @@ -0,0 +1,113 @@ +import numpy as np +import tifffile as tif +import transforms3d as tf3d +from scipy.ndimage import affine_transform +import napari +from elf.wrapper.resized_volume import ResizedVolume + + +def downscale_seg(seg, factor): + new_shape = np.array(seg.shape) // 2 + downsampled_seg = ResizedVolume(seg, shape=new_shape)[:] + print("Downsampled shape", downsampled_seg.shape) + + return downsampled_seg + + +def rotate_seg(seg, angles): + # pad zeros so segmentation is not cut out off the image + shape = seg.shape + diagonal = int(np.ceil(np.linalg.norm(shape))) + print("Diagonal", diagonal) + + # Calculate how much padding is needed on each axis + pad_widths = [] + for dim in shape: + total_pad = diagonal - dim + before = total_pad // 2 + after = total_pad - before + pad_widths.append((before, after)) + + # Apply zero padding + padded = np.pad(seg, pad_width=pad_widths, mode='constant', constant_values=0) + print("New shape after padding:", padded.shape) + + rotation_matrix = tf3d.euler.euler2mat(*angles, axes='szyx') + print(rotation_matrix) + + # Compute center + center = np.array(padded.shape) / 2 + offset = center - rotation_matrix @ center + + # rotate the image around the center + rotated_seg = affine_transform( + padded, + matrix=rotation_matrix, + offset=offset, # offset to ensure the image fits within the new bounding box + order=0, # interpolation (use 0 for discrete/label data) + mode='constant', # fill mode + cval=0.0 # fill value (if constant mode) + ) + + print(f"New shape after rotation: {rotated_seg.shape}") + return rotated_seg + + +def remove_instances(seg, prob=0.05): + # Get all unique instance IDs, excluding background (assumed to be 0) + instance_ids = np.unique(seg) + instance_ids = instance_ids[instance_ids != 0] + + # Randomly select 10% of the instance IDs + num_to_remove = int(len(instance_ids) * prob) + print(f"Number of instances to remove: {num_to_remove}") + selected_ids = np.random.choice(instance_ids, size=num_to_remove, replace=False) + + # Create a mask for the selected IDs and set them to 0 + mask = np.isin(seg, selected_ids) + seg[mask] = 0 + + print(f"Number of instances left: {len(np.unique(seg))}") + return seg + + +def crop_to_bbox(seg): + """Crops a 3D volume to the minimal bounding box around all nonzero voxels.""" + # Find where the volume is nonzero (i.e., contains instances) + nonzero = np.argwhere(seg) + + # Compute bounding box from min to max index along each axis + z_min, y_min, x_min = nonzero.min(axis=0) + z_max, y_max, x_max = nonzero.max(axis=0) + 1 # +1 to include the max index + + # Crop the volume + cropped = seg[z_min:z_max, y_min:y_max, x_min:x_max] + return cropped + + +def main(): + seg = tif.imread("data/platy1_muscles_stardist.tif") + + # downsample image + factor = 4 + seg = downscale_seg(seg, factor) + tif.imwrite("./data/platy1_muscles_stardist_fixed.tif", seg) + # rotate image + angles = [np.deg2rad(155), np.deg2rad(30), np.deg2rad(65)] # z,y,x + rotated_seg = rotate_seg(seg, angles) + + # randomly remove instances + probability = 0.05 + seg_rm = remove_instances(rotated_seg, prob=probability) + seg_cropped = crop_to_bbox(seg_rm) + print("Cropped shape", seg_cropped.shape) + tif.imwrite("./data/platy1_muscles_stardist_moving.tif", seg_cropped) + + # visualize + v = napari.Viewer() + v.add_labels(seg_cropped, name="moving") + napari.run() + + +if __name__ == "__main__": + main() From 0fea233278fb1cb79a6ebb6bf429dd40bfe2631a Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Thu, 17 Apr 2025 14:01:20 +0200 Subject: [PATCH 02/23] save test deformations --- examples/.gitattributes | 1 + ...erate_test_data.py => deform_test_data.py} | 20 ++++++++++++++----- examples/reverse_transform.py | 0 3 files changed, 16 insertions(+), 5 deletions(-) rename examples/{generate_test_data.py => deform_test_data.py} (88%) create mode 100644 examples/reverse_transform.py diff --git a/examples/.gitattributes b/examples/.gitattributes index d92417a..c2e5313 100644 --- a/examples/.gitattributes +++ b/examples/.gitattributes @@ -1 +1,2 @@ *.tif filter=lfs diff=lfs merge=lfs -text +*.n5 filter=lfs diff=lfs merge=lfs -text diff --git a/examples/generate_test_data.py b/examples/deform_test_data.py similarity index 88% rename from examples/generate_test_data.py rename to examples/deform_test_data.py index c549973..b7fea54 100644 --- a/examples/generate_test_data.py +++ b/examples/deform_test_data.py @@ -1,10 +1,15 @@ import numpy as np import tifffile as tif import transforms3d as tf3d +import z5py from scipy.ndimage import affine_transform import napari from elf.wrapper.resized_volume import ResizedVolume +import sys +sys.path.append("..") +from matchmaker.vis import plot_overlay + def downscale_seg(seg, factor): new_shape = np.array(seg.shape) // 2 @@ -91,7 +96,9 @@ def main(): # downsample image factor = 4 seg = downscale_seg(seg, factor) - tif.imwrite("./data/platy1_muscles_stardist_fixed.tif", seg) + with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "w") as f: + f.create_dataset("seg", data=seg) + # rotate image angles = [np.deg2rad(155), np.deg2rad(30), np.deg2rad(65)] # z,y,x rotated_seg = rotate_seg(seg, angles) @@ -101,12 +108,15 @@ def main(): seg_rm = remove_instances(rotated_seg, prob=probability) seg_cropped = crop_to_bbox(seg_rm) print("Cropped shape", seg_cropped.shape) - tif.imwrite("./data/platy1_muscles_stardist_moving.tif", seg_cropped) + with z5py.File("./data/platy1_muscles_stardist_moving.n5", "w") as f: + f.create_dataset("seg", data=seg_cropped) + + plot_overlay(seg, seg_cropped) # visualize - v = napari.Viewer() - v.add_labels(seg_cropped, name="moving") - napari.run() + # v = napari.Viewer() + # v.add_labels(seg_cropped, name="moving") + # napari.run() if __name__ == "__main__": diff --git a/examples/reverse_transform.py b/examples/reverse_transform.py new file mode 100644 index 0000000..e69de29 From 4a40de817f019f90fb37335e4a4cd85365503e6d Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Wed, 23 Apr 2025 09:33:35 +0200 Subject: [PATCH 03/23] updated gitattributes file --- .gitattributes | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 .gitattributes diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..9a6d0d0 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,2 @@ +.n5 filter=lfs diff=lfs merge=lfs -text +.tif filter=lfs diff=lfs merge=lfs -text From c24fed7774c0158902052f0907be8b237355bd70 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Wed, 23 Apr 2025 09:40:39 +0200 Subject: [PATCH 04/23] new gitattributes --- .gitattributes | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitattributes b/.gitattributes index 9a6d0d0..5ee325b 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,2 +1,5 @@ .n5 filter=lfs diff=lfs merge=lfs -text .tif filter=lfs diff=lfs merge=lfs -text +.n5/** filter=lfs diff=lfs merge=lfs -text +examples/data filter=lfs diff=lfs merge=lfs -text +examples/data/** filter=lfs diff=lfs merge=lfs -text From 367ef06ac94c13f89788ca110c5cf8a1771b3761 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Wed, 23 Apr 2025 09:41:32 +0200 Subject: [PATCH 05/23] started pre alignment --- examples/deform_test_data.py | 84 ++++++--------------------------- examples/registration_test.py | 30 ++++++++++++ examples/reverse_deformation.py | 27 +++++++++++ examples/u.py | 57 ++++++++++++++++++++++ 4 files changed, 129 insertions(+), 69 deletions(-) create mode 100644 examples/registration_test.py create mode 100644 examples/reverse_deformation.py create mode 100644 examples/u.py diff --git a/examples/deform_test_data.py b/examples/deform_test_data.py index b7fea54..f8b94d5 100644 --- a/examples/deform_test_data.py +++ b/examples/deform_test_data.py @@ -1,14 +1,10 @@ import numpy as np import tifffile as tif -import transforms3d as tf3d import z5py -from scipy.ndimage import affine_transform -import napari from elf.wrapper.resized_volume import ResizedVolume -import sys -sys.path.append("..") -from matchmaker.vis import plot_overlay +from matchmaker.transform_utils import pad_img, get_rotation_matrix, rotate_img, crop_to_bbox +from matchmaker.vis import plot_three_slices, plot_overlay def downscale_seg(seg, factor): @@ -19,45 +15,6 @@ def downscale_seg(seg, factor): return downsampled_seg -def rotate_seg(seg, angles): - # pad zeros so segmentation is not cut out off the image - shape = seg.shape - diagonal = int(np.ceil(np.linalg.norm(shape))) - print("Diagonal", diagonal) - - # Calculate how much padding is needed on each axis - pad_widths = [] - for dim in shape: - total_pad = diagonal - dim - before = total_pad // 2 - after = total_pad - before - pad_widths.append((before, after)) - - # Apply zero padding - padded = np.pad(seg, pad_width=pad_widths, mode='constant', constant_values=0) - print("New shape after padding:", padded.shape) - - rotation_matrix = tf3d.euler.euler2mat(*angles, axes='szyx') - print(rotation_matrix) - - # Compute center - center = np.array(padded.shape) / 2 - offset = center - rotation_matrix @ center - - # rotate the image around the center - rotated_seg = affine_transform( - padded, - matrix=rotation_matrix, - offset=offset, # offset to ensure the image fits within the new bounding box - order=0, # interpolation (use 0 for discrete/label data) - mode='constant', # fill mode - cval=0.0 # fill value (if constant mode) - ) - - print(f"New shape after rotation: {rotated_seg.shape}") - return rotated_seg - - def remove_instances(seg, prob=0.05): # Get all unique instance IDs, excluding background (assumed to be 0) instance_ids = np.unique(seg) @@ -76,47 +33,36 @@ def remove_instances(seg, prob=0.05): return seg -def crop_to_bbox(seg): - """Crops a 3D volume to the minimal bounding box around all nonzero voxels.""" - # Find where the volume is nonzero (i.e., contains instances) - nonzero = np.argwhere(seg) - - # Compute bounding box from min to max index along each axis - z_min, y_min, x_min = nonzero.min(axis=0) - z_max, y_max, x_max = nonzero.max(axis=0) + 1 # +1 to include the max index - - # Crop the volume - cropped = seg[z_min:z_max, y_min:y_max, x_min:x_max] - return cropped - - def main(): seg = tif.imread("data/platy1_muscles_stardist.tif") # downsample image factor = 4 seg = downscale_seg(seg, factor) + seg_fixed = crop_to_bbox(seg) + print("Cropped shape", seg_fixed.shape) + # save downsampled image with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "w") as f: - f.create_dataset("seg", data=seg) + f.create_dataset("seg", data=seg_fixed) # rotate image + seg_padded = pad_img(seg) angles = [np.deg2rad(155), np.deg2rad(30), np.deg2rad(65)] # z,y,x - rotated_seg = rotate_seg(seg, angles) + rotation_matrix = get_rotation_matrix(angles, save_path="./data/rotation_matrix.txt") + rotated_seg = rotate_img(seg_padded, rotation_matrix) + seg_cropped = crop_to_bbox(rotated_seg) + print("Cropped shape", seg_cropped.shape) # randomly remove instances probability = 0.05 - seg_rm = remove_instances(rotated_seg, prob=probability) - seg_cropped = crop_to_bbox(seg_rm) - print("Cropped shape", seg_cropped.shape) + seg_moving = remove_instances(seg_cropped, prob=probability) with z5py.File("./data/platy1_muscles_stardist_moving.n5", "w") as f: - f.create_dataset("seg", data=seg_cropped) + f.create_dataset("seg", data=seg_moving) - plot_overlay(seg, seg_cropped) # visualize - # v = napari.Viewer() - # v.add_labels(seg_cropped, name="moving") - # napari.run() + plot_three_slices(seg_moving) + plot_overlay(seg_fixed, seg_moving) if __name__ == "__main__": diff --git a/examples/registration_test.py b/examples/registration_test.py new file mode 100644 index 0000000..80590fd --- /dev/null +++ b/examples/registration_test.py @@ -0,0 +1,30 @@ +import z5py +import numpy as np + +from matchmaker.prealignment import get_SVD_transform +from matchmaker.vis import plot_three_slices, plot_overlay +from matchmaker.transform_utils import rotate + + +def main(): + with z5py.File("./data/platy1_muscles_stardist_moving.n5", "r") as f: + seg_moving = f["seg"][:] + + with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: + seg_fixed = f["seg"][:] + + gc, Vt, pos_c = get_SVD_transform(seg_moving) + + seg_moving_rotated = rotate(seg_moving, np.linalg.inv(Vt)) # TODO: add gc as rotation center? + + plot_three_slices(seg_moving_rotated) + + plot_overlay( + seg_moving_rotated, + seg_fixed, + ) + # y sieht aus wie z: kacke + + +if __name__ == "__main__": + main() diff --git a/examples/reverse_deformation.py b/examples/reverse_deformation.py new file mode 100644 index 0000000..708781a --- /dev/null +++ b/examples/reverse_deformation.py @@ -0,0 +1,27 @@ +import numpy as np +import z5py +from matchmaker.vis import plot_three_slices, plot_overlay +from deform_test_data import pad_seg, crop_to_bbox, rotate_seg + + +def main(): + with z5py.File("./data/platy1_muscles_stardist_moving.n5", "r") as f: + seg = f["seg"][:] + + seg_padded = pad_seg(seg) + rotation_matrix = np.loadtxt("./data/rotation_matrix.txt") + rotated_seg = rotate_seg(seg_padded, np.linalg.inv(rotation_matrix)) + seg_cropped = crop_to_bbox(rotated_seg) + print("Cropped shape", seg_cropped.shape) + + plot_three_slices(seg_cropped) + + # compare with original image + with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: + seg_fixed = f["seg"][:] + + plot_overlay(seg_cropped, seg_fixed) + + +if __name__ == "__main__": + main() diff --git a/examples/u.py b/examples/u.py new file mode 100644 index 0000000..06510d0 --- /dev/null +++ b/examples/u.py @@ -0,0 +1,57 @@ + +from transforms3d.euler import euler2mat, mat2euler +from scipy.ndimage import affine_transform +import numpy as np + +def get_translation_transform(translation): + M = np.identity(4) + M[0:3, 3] = translation + return M + + +def get_affine_transform(translation, rotation, rotation_order="szyx"): + M = get_translation_transform(translation) + M[0:3, 0:3] = euler2mat(*rotation, rotation_order) + return M + + +def flip_dorsal(img, n5_fname): + """ + For dorsal-direction samples flip along z so they are same orientation + """ + if "dorsal" in n5_fname.lower(): + print("Flip dorsal image") + img_flipped = np.flip(img, axis=0) + return img_flipped + + else: + return img + +def move_along_z(margin, image_shape): + """ + After rotation the volume might not fit within the same Z range, + so move it up by some margin and make output image size bigger + """ + assert len(image_shape) == 3 + R_t = get_translation_transform([-margin, -margin, -margin]) + new_image_shape = list(image_shape) + new_image_shape = [sh + 2 * margin for sh in image_shape] + + return R_t, new_image_shape + + +R_t = get_translation_transform(-gc) +R_t_r = get_translation_transform(gc) + +R_rot = get_translation_transform([0, 0, 0]) +R_rot[0:3, 0:3] = np.linalg.inv(Vt) + +R_rot_90_z = get_affine_transform([0, 0, 0], [0, -np.pi / 2, 0], rotation_order="szyx") + +R_t_margin, new_image_shape = move_along_z(80, norm_img.shape) + +mat = R_t_r @ R_rot @ R_t @ R_t_margin +print(mat) +print(norm_img.shape) + +img_reg = affine_transform(img_isotropic, mat, order=0, output_shape=new_image_shape) \ No newline at end of file From 22d21ef704ff51a80b4227b5df55f5336486e8f0 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Wed, 23 Apr 2025 18:48:56 +0200 Subject: [PATCH 06/23] some linter corrections; started rotation to pre-align --- examples/registration_test.py | 25 ++-- matchmaker/align_rigid_elastix.py | 230 ++++++++++++++++-------------- matchmaker/data.py | 21 +-- matchmaker/elastix_utils.py | 73 +++++----- matchmaker/prealignment.py | 125 ++++++++-------- matchmaker/preprocessing.py | 31 ++++ matchmaker/transform_utils.py | 122 ++++++++++++++++ matchmaker/vis.py | 56 ++++---- 8 files changed, 448 insertions(+), 235 deletions(-) create mode 100644 matchmaker/preprocessing.py create mode 100644 matchmaker/transform_utils.py diff --git a/examples/registration_test.py b/examples/registration_test.py index 80590fd..b10e8b4 100644 --- a/examples/registration_test.py +++ b/examples/registration_test.py @@ -3,7 +3,9 @@ from matchmaker.prealignment import get_SVD_transform from matchmaker.vis import plot_three_slices, plot_overlay -from matchmaker.transform_utils import rotate +from matchmaker.transform_utils import rotate_with_padding, rotate_with_shape +from scipy.ndimage import affine_transform + def main(): @@ -15,15 +17,22 @@ def main(): gc, Vt, pos_c = get_SVD_transform(seg_moving) - seg_moving_rotated = rotate(seg_moving, np.linalg.inv(Vt)) # TODO: add gc as rotation center? + + + # seg_moving_rotated = rotate_with_padding(seg_moving, Vt.T) # TODO: add gc as rotation center? + # seg_moving_rotated2 = rotate_with_shape(seg_moving, Vt.T) + + import napari + v = napari.Viewer() + v.add_image(seg_moving_rotated2, name="fixed") + napari.run() - plot_three_slices(seg_moving_rotated) + # plot_three_slices(seg_moving_rotated) - plot_overlay( - seg_moving_rotated, - seg_fixed, - ) - # y sieht aus wie z: kacke + # plot_overlay( + # seg_fixed, + # seg_moving_rotated2, + # ) if __name__ == "__main__": diff --git a/matchmaker/align_rigid_elastix.py b/matchmaker/align_rigid_elastix.py index 320bf46..72c0e35 100644 --- a/matchmaker/align_rigid_elastix.py +++ b/matchmaker/align_rigid_elastix.py @@ -1,134 +1,154 @@ from pathlib import Path import numpy as np import argparse -from aicsimageio import readers - from platy_reg.n5_utils import read_volume, write_volume, get_attrs -from platy_reg.vis import plot_overlay - -from aicsimageio.aics_image import AICSImage +from matchmaker.vis import plot_overlay import logging import sys - -from platy_reg.preprocessing import percentile_norm -from platy_reg.elastix_utils import * - +from matchmaker import elastix_utils import itk -from skimage.filters import gaussian -def elastix_segm_rigid_alignment(fixed_img_np, fixed_resolution, moving_img_np, moving_resolution, log_dir): + +def elastix_segm_rigid_alignment( + fixed_img_np, fixed_resolution, moving_img_np, moving_resolution, log_dir +): """ Run rigid alignment of the ventral and dorsal datasets using elastix. """ - - logging.info(f"Do rigid transform of unnormalized images") + + logging.info("Do rigid transform of unnormalized images") logging.info(f"Fixed resolution {fixed_resolution}") logging.info(f"Moving resolution {moving_resolution}") - logging.info(f"Do rigid transform of unnormalized images") - + logging.info("Do rigid transform of unnormalized images") + fixed_img_semantic_np = (fixed_img_np > 0).astype(np.float32) moving_img_semantic_np = (moving_img_np > 0).astype(np.float32) - - fixed_img = itk_scalar_img(fixed_img_semantic_np, fixed_resolution) - moving_img = itk_scalar_img(moving_img_semantic_np, moving_resolution) - - logging.info(f"Fixed image") + + fixed_img = elastix_utils.itk_scalar_img(fixed_img_semantic_np, fixed_resolution) + moving_img = elastix_utils.itk_scalar_img(moving_img_semantic_np, moving_resolution) + + logging.info("Fixed image") logging.info(f"{fixed_img}") - logging.info(f"Moving image") + logging.info("Moving image") logging.info(f"{moving_img}") - + # plot_overlay(fixed_img_np, moving_img_np, log_dir / "intersample_segm_overlay_before_alignment.png") - - plot_overlay(itk_to_np_order(itk.GetArrayFromImage(fixed_img)), itk_to_np_order(itk.GetArrayFromImage(moving_img)), log_dir / "intersample_segm_overlay_before_alignment.png") - - parameter_map_paths = ["pipeline_steps/inter_sample_registration/ParameterMap_segm_rigid_registration_corr.txt"] - logging.info(f"Run rigid registration") - result_image, result_transform_parameters = run_registration(fixed_img, moving_img, parameter_map_paths, str(log_dir), log_name="elastix_log_rigid.log", set_threads=True) + + plot_overlay( + elastix_utils.itk_to_np_order(itk.GetArrayFromImage(fixed_img)), + elastix_utils.itk_to_np_order(itk.GetArrayFromImage(moving_img)), + log_dir / "intersample_segm_overlay_before_alignment.png", + ) + + parameter_map_paths = [ + "pipeline_steps/inter_sample_registration/ParameterMap_segm_rigid_registration_corr.txt" + ] + logging.info("Run rigid registration") + result_image, result_transform_parameters = elastix_utils.run_registration( + fixed_img, + moving_img, + parameter_map_paths, + str(log_dir), + log_name="elastix_log_rigid.log", + set_threads=True, + ) # serialize_parameter_object(result_transform_parameters, "ParameterMap_rigid_transform", log_dir) - - logging.info(f"Result image shape {result_image.shape}") - result_img_np = itk_to_np_order(itk.GetArrayFromImage(result_image)) - plot_overlay(itk_to_np_order(itk.GetArrayFromImage(fixed_img)), result_img_np, log_dir / f"intersample_segm_rigid_alignment_semantic.png") + logging.info(f"Result image shape {result_image.shape}") + result_img_np = elastix_utils.itk_to_np_order(itk.GetArrayFromImage(result_image)) + plot_overlay( + elastix_utils.itk_to_np_order(itk.GetArrayFromImage(fixed_img)), + result_img_np, + log_dir / "intersample_segm_rigid_alignment_semantic.png", + ) - logging.info(f"Apply transform to all channels") - result_img_np = apply_transform_chanwise(result_transform_parameters, moving_img_np, moving_resolution) + logging.info("Apply transform to all channels") + result_img_np = elastix_utils.apply_transform_chanwise( + result_transform_parameters, moving_img_np, moving_resolution + ) logging.info(f"Result image shape {result_img_np.shape}") - + return result_img_np def main(): - parser = argparse.ArgumentParser( - description="""Align rigidly moving segmentation volume to fixed volume. - """ + parser = argparse.ArgumentParser( + description="""Align rigidly moving segmentation volume to fixed volume.""" + ) + + parser.add_argument("fixed_n5", type=str, help="Path of the output n5") + parser.add_argument("fixed_key", type=str, help="Key of ventral dataset") + parser.add_argument("moving_n5", type=str, help="Path of the output n5") + parser.add_argument("moving_key", type=str, help="Key of ventral dataset") + parser.add_argument("output_n5", type=str, help="Path of the output n5") + parser.add_argument("output_key", type=str, help="Key to write aligned ventral dataset to in the output n5") + parser.add_argument("log_dir", type=str, help="Directory to store diagnostic plots etc") + parser.add_argument("log_path", type=str, help="Path to log file") + args = parser.parse_args() + + log_dir = Path(args.log_dir) + log_dir.mkdir(exist_ok=True) + + logging.basicConfig( + level=logging.INFO, + format="%(asctime)s [%(levelname)s] %(message)s", + handlers=[ + logging.FileHandler(args.log_path, mode="w"), + logging.StreamHandler(sys.stdout), + ], + datefmt="%Y-%m-%d %H:%M:%S", + ) + + logging.info("Read image file") + fixed_img_np = read_volume(args.fixed_n5, args.fixed_key) + + # pad_z = 50 + # fixed_img_np = np.pad(fixed_img_np, pad_width=((pad_z, pad_z), (0, 0), (0, 0))) + + # fixed_mask = gaussian(fixed_img_np, sigma=[6, 12, 12]) + # fixed_mask_bin = fixed_mask > np.quantile(fixed_mask, 0.7) + # fixed_img_np = fixed_img_np * fixed_mask_bin + + if args.fixed_n5 == args.moving_n5: + logging.info("Same n5 for fixed and moving image, not doing registration") + moving_img_np = fixed_img_np + + else: + moving_img_np = read_volume(args.moving_n5, args.moving_key) + # moving_mask = gaussian(moving_img_np, sigma=[6, 12, 12]) + # moving_mask_bin = moving_mask > np.quantile(moving_mask, 0.7) + # moving_img_np = moving_img_np * moving_mask_bin + + fixed_img_np = fixed_img_np.astype(np.float32) + moving_img_np = moving_img_np.astype(np.float32) + + logging.info("Start registration") + + # Dirty hack for the sample size mismatch + # For no apparent reason the size of the light samples seems to be ~20% larger than the EM + moving_resolution = get_attrs(args.moving_n5, args.moving_key)["resolution"] + moving_resolution = [r * 0.9 for r in moving_resolution] + + moving_img_np = elastix_segm_rigid_alignment( + fixed_img_np, + get_attrs(args.fixed_n5, args.fixed_key)["resolution"], + moving_img_np, + moving_resolution, + log_dir, ) - parser.add_argument("fixed_n5", type=str, help="Path of the output n5") - parser.add_argument("fixed_key", type=str, help="Key of ventral dataset") - parser.add_argument("moving_n5", type=str, help="Path of the output n5") - parser.add_argument("moving_key", type=str, help="Key of ventral dataset") - parser.add_argument("output_n5", type=str, help="Path of the output n5") - parser.add_argument("output_key", type=str, help="Key to write aligned ventral dataset to in the output n5") - parser.add_argument("log_dir", type=str, help="Directory to store diagnostic plots etc") - parser.add_argument("log_path", type=str, help="Path to log file") - args = parser.parse_args() - - log_dir = Path(args.log_dir) - log_dir.mkdir(exist_ok=True) - - logging.basicConfig(level=logging.INFO, - format="%(asctime)s [%(levelname)s] %(message)s", - handlers=[logging.FileHandler(args.log_path, mode="w"), - logging.StreamHandler(sys.stdout)], - datefmt='%Y-%m-%d %H:%M:%S') - - logging.info("Read image file") - fixed_img_np = read_volume(args.fixed_n5, args.fixed_key) - - # pad_z = 50 - # fixed_img_np = np.pad(fixed_img_np, pad_width=((pad_z, pad_z), (0, 0), (0, 0))) - - # fixed_mask = gaussian(fixed_img_np, sigma=[6, 12, 12]) - # fixed_mask_bin = fixed_mask > np.quantile(fixed_mask, 0.7) - # fixed_img_np = fixed_img_np * fixed_mask_bin - - - if args.fixed_n5 == args.moving_n5: - logging.info("Same n5 for fixed and moving image, not doing registration") - moving_img_np = fixed_img_np - - - else: - moving_img_np = read_volume(args.moving_n5, args.moving_key) - # moving_mask = gaussian(moving_img_np, sigma=[6, 12, 12]) - # moving_mask_bin = moving_mask > np.quantile(moving_mask, 0.7) - # moving_img_np = moving_img_np * moving_mask_bin - - fixed_img_np = fixed_img_np.astype(np.float32) - moving_img_np = moving_img_np.astype(np.float32) - - logging.info("Start registration") - - # Dirty hack for the sample size mismatch - # For no apparent reason the size of the light samples seems to be ~20% larger than the EM - moving_resolution = get_attrs(args.moving_n5, args.moving_key)["resolution"] - moving_resolution = [r * 0.9 for r in moving_resolution] - - moving_img_np = elastix_segm_rigid_alignment(fixed_img_np, - get_attrs(args.fixed_n5, args.fixed_key)["resolution"], - moving_img_np, - moving_resolution, - log_dir) - - logging.info("Write results") - attributes = dict(get_attrs(args.moving_n5, args.moving_key)) - attributes["resolution"] = get_attrs(args.fixed_n5, args.fixed_key)["resolution"] - - write_volume(args.output_n5, moving_img_np, args.output_key, chunks=(128, 512, 512), attrs=attributes) - - -if __name__=="__main__": + logging.info("Write results") + attributes = dict(get_attrs(args.moving_n5, args.moving_key)) + attributes["resolution"] = get_attrs(args.fixed_n5, args.fixed_key)["resolution"] + + write_volume( + args.output_n5, + moving_img_np, + args.output_key, + chunks=(128, 512, 512), + attrs=attributes, + ) + + +if __name__ == "__main__": main() - - diff --git a/matchmaker/data.py b/matchmaker/data.py index 8ce95a1..a9739f7 100644 --- a/matchmaker/data.py +++ b/matchmaker/data.py @@ -4,11 +4,12 @@ def create_point_cloud(segm): - props = regionprops_table(segm, properties=("label", 'centroid')) + props = regionprops_table(segm, properties=("label", "centroid")) props = pd.DataFrame(props) segm_labels = props["label"].to_numpy() centroid_columns = [col for col in props.columns if col.startswith("centroid")] pos = props[centroid_columns].to_numpy().astype(np.float32) + # TODO: center, maybe w/o labels return pos, segm_labels @@ -22,7 +23,7 @@ def write_pcd(pcd_df, pcd_path): pcd_df.to_csv(pcd_path, index=False) -def pcd_np_to_df(pos:np.array, segm_labels=None, reg_gt=None, matching=None): +def pcd_np_to_df(pos: np.array, segm_labels=None, reg_gt=None, matching=None): # print("pos len", len(pos)) pcd_df = pd.DataFrame() N = pos.shape[0] @@ -32,11 +33,15 @@ def pcd_np_to_df(pos:np.array, segm_labels=None, reg_gt=None, matching=None): pcd_df[f"coord_{col}"] = pos[:, col] if segm_labels is not None: - assert len(segm_labels) == N, print(f"Labels of shape {segm_labels.shape} don't correspond to coordinates of shape {pos.shape}") + assert len(segm_labels) == N, print( + f"Labels of shape {segm_labels.shape} don't correspond to coordinates of shape {pos.shape}" + ) pcd_df["segm_labels"] = segm_labels if reg_gt is not None: - assert len(reg_gt) == N, print(f"Labels of shape {reg_gt.shape} don't correspond to coordinates of shape {pos.shape}") + assert len(reg_gt) == N, print( + f"Labels of shape {reg_gt.shape} don't correspond to coordinates of shape {pos.shape}" + ) pcd_df["reg_gt"] = reg_gt pcd_df.set_index("order_idx", drop=False) @@ -67,12 +72,12 @@ def index_pairs_to_matching(): def write_index_pairs(pairs, pairs_path): - with open(pairs_path,"w+") as f: + with open(pairs_path, "w+") as f: for pair in pairs: - f.write(str(pair[0]) + "," + str(pair[1]) + '\n') + f.write(str(pair[0]) + "," + str(pair[1]) + "\n") def read_index_pairs(pairs_path): - with open(pairs_path,"r") as f: - pairs = [(int(l.split(",")[0]), int(l.split(",")[1])) for l in f.readlines()] + with open(pairs_path, "r") as f: + pairs = [(int(loc.split(",")[0]), int(loc.split(",")[1])) for loc in f.readlines()] return pairs diff --git a/matchmaker/elastix_utils.py b/matchmaker/elastix_utils.py index af791dd..60fab2d 100644 --- a/matchmaker/elastix_utils.py +++ b/matchmaker/elastix_utils.py @@ -1,19 +1,7 @@ from pathlib import Path import numpy as np -import argparse -from aicsimageio import readers - -from platy_reg.n5_utils import read_volume, write_volume, get_attrs -from platy_reg.vis import plot_overlay - -from aicsimageio.aics_image import AICSImage import logging -import sys - -import re - import itk -from platy_reg.preprocessing import percentile_norm def initial_alignment(ventral_img_np, dorsal_img_np): @@ -33,12 +21,12 @@ def itk_scalar_img(img: np.array, resolution, ch=0): ch: _description_. Defaults to 0. Returns: - Scalar ITK Image object for one of the channels. + Scalar ITK Image object for one of the channels. """ logging.info(f"Converting image of type {img.dtype} to ITK object") if img.ndim == 4: img = np.moveaxis(img[ch, :, :, :], 0, -1) - + else: img = np.moveaxis(img[:, :, :], 0, -1) @@ -57,7 +45,7 @@ def itk_to_np_order(img: np.array): img = np.swapaxes(img, 0, 3) img = np.swapaxes(img, 1, 2) return np.flip(img) - + def np_to_itk_order(img: np.array): # In numpy: CZYX @@ -68,18 +56,25 @@ def np_to_itk_order(img: np.array): img = np.swapaxes(img, 0, 3) img = np.swapaxes(img, 1, 2) return img - + def create_parameter_object(parameter_map_paths): parameter_object = itk.ParameterObject.New() for parameter_map_path in parameter_map_paths: parameter_object.AddParameterFile(parameter_map_path) print(parameter_object) - + return parameter_object -def run_registration(fixed_img, moving_img, parameter_map_paths, log_dir, log_name="elastix.log", set_threads=False): +def run_registration( + fixed_img, + moving_img, + parameter_map_paths, + log_dir, + log_name="elastix.log", + set_threads=False, +): logging.info("Start creating parameter object") parameter_object = create_parameter_object(parameter_map_paths) # Load Elastix Image Filter Object @@ -96,17 +91,26 @@ def run_registration(fixed_img, moving_img, parameter_map_paths, log_dir, log_na elastix_object.SetOutputDirectory(log_dir) elastix_object.SetLogFileName(log_name) logging.info("Set parameter map") - + logging.info("Start registration") elastix_object.UpdateLargestPossibleRegion() - + result_image = elastix_object.GetOutput() result_transform_parameters = elastix_object.GetTransformParameterObject() - + return result_image, result_transform_parameters -def run_pointset_registration(fixed_img, moving_img, parameter_map_paths, fixed_pointset, moving_pointset, log_dir, log_name="elastix.log", set_threads=False): +def run_pointset_registration( + fixed_img, + moving_img, + parameter_map_paths, + fixed_pointset, + moving_pointset, + log_dir, + log_name="elastix.log", + set_threads=False, +): logging.info("Start creating parameter object") parameter_object = create_parameter_object(parameter_map_paths) # Load Elastix Image Filter Object @@ -120,17 +124,17 @@ def run_pointset_registration(fixed_img, moving_img, parameter_map_paths, fixed_ elastix_object.SetLogToFile(True) elastix_object.SetOutputDirectory(log_dir) elastix_object.SetLogFileName(log_name) - + logging.info("Created registration object") logging.info(elastix_object) logging.info("Set parameter map") - + logging.info("Start registration") elastix_object.UpdateLargestPossibleRegion() - + result_image = elastix_object.GetOutput() result_transform_parameters = elastix_object.GetTransformParameterObject() - + return result_image, result_transform_parameters @@ -138,9 +142,11 @@ def serialize_parameter_object(parameter_object, prefix, write_dir): write_dir = Path(write_dir) for index in range(parameter_object.GetNumberOfParameterMaps()): parameter_map = parameter_object.GetParameterMap(index) - parameter_object.WriteParameterFile(parameter_map, write_dir / f"{prefix}.{index}.txt") - - + parameter_object.WriteParameterFile( + parameter_map, write_dir / f"{prefix}.{index}.txt" + ) + + def deserialize_parameter_object(prefix, cur_dir=Path("./")): parameter_files = sorted(list(cur_dir.glob(f"{prefix}*.txt"))) parameter_files = [str(fname) for fname in parameter_files] @@ -154,7 +160,7 @@ def create_transformix_object(transform_parameter_object): transformix_filter = itk.TransformixFilter[ImageType].New() transformix_filter.SetTransformParameterObject(transform_parameter_object) return transformix_filter - + def apply_transform(transformix_filter, moving_img): transformix_filter.SetMovingImage(moving_img) @@ -164,10 +170,11 @@ def apply_transform(transformix_filter, moving_img): output_img_np = itk_to_np_order(output_img_np) return output_img_np + def apply_transform_chanwise(transform_parameter_object, moving_img_np, resolution): - transformix_filter = create_transformix_object(transform_parameter_object) + transformix_filter = create_transformix_object(transform_parameter_object) result_img = [] - + if moving_img_np.ndim == 4: for chan in range(moving_img_np.shape[0]): moving_img = itk_scalar_img(moving_img_np, resolution, chan) @@ -177,5 +184,5 @@ def apply_transform_chanwise(transform_parameter_object, moving_img_np, resoluti else: moving_img = itk_scalar_img(moving_img_np, resolution, 0) result_img = apply_transform(transformix_filter, moving_img) - + return result_img diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index 40b8804..098d70c 100644 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -1,11 +1,11 @@ from pathlib import Path import numpy as np import matplotlib.pyplot as plt -from platy_reg.vis import plot_overlay, plot_three_slices import logging -from platy_reg.preprocessing import convert_to_point_cloud from scipy.ndimage import rotate from skimage.filters import gaussian +# from matchmaker.preprocessing import convert_to_point_cloud +from matchmaker.data import create_point_cloud def get_SVD_transform(img, plot_path=None, percentile_trsh=90): @@ -19,66 +19,73 @@ def get_SVD_transform(img, plot_path=None, percentile_trsh=90): Returns: Variance matrix and principal axes matrix. """ - trsh = np.percentile(img, percentile_trsh) - logging.info(f"Threshold used for converting to point cloud: {trsh}") - X = convert_to_point_cloud(img, trsh) - gc = X.mean(axis=0) + + pos, _ = create_point_cloud(img) + # how many pos? subsample? + # center coord with center of the image + + # else: + # trsh = np.percentile(img, percentile_trsh) + # logging.info(f"Threshold used for converting to point cloud: {trsh}") + # X = convert_to_point_cloud(img, trsh) + gc = pos.mean(axis=0) gc = np.array(img.shape) // 2 - X_c = X - gc - logging.info(f"Point cloud shape {X.shape}") + pos_c = pos - gc + logging.info(f"Point cloud shape {pos.shape}") logging.info(f"Point cloud center {gc}") - - if X.shape[1] == 3: - random_subset = np.random.choice(X.shape[0], 10000, replace=False) - X_subset = X_c[random_subset, 1:] - else: - X_subset = X_c[::100] - logging.info("Subset point cloud") - logging.info(f"Subset point cloud shape {X_subset.shape}") - + + # if X.shape[1] == 3: + # random_subset = np.random.choice(X.shape[0], 10000, replace=False) + # X_subset = X_c[random_subset, 1:] + # else: + # X_subset = X_c[::100] + # logging.info("Subset point cloud") + # logging.info(f"Subset point cloud shape {X_subset.shape}") + logging.info("Run SVD") - U, S, Vt = np.linalg.svd(X_subset) - # logging.info("U") - # logging.info(U) + U, S, Vt = np.linalg.svd(pos_c, full_matrices=False) + logging.info("U") + logging.info(U) logging.info("S") logging.info(str(S)) logging.info("Vt") logging.info(str(Vt)) - + logging.info("Rotate point cloud") - vr = X_subset @ Vt.T - + vr = pos_c @ Vt.T + plt.figure(figsize=(10, 5)) - plt.subplot(1,2,1) - plt.title('original vertices') - plt.scatter(X_subset[:, 0], X_subset[:, 1], alpha=0.2) - plt.subplot(1,2,2) - plt.title('rotated vertices') + plt.subplot(1, 2, 1) + plt.title("original vertices") + plt.scatter(pos_c[:, 0], pos_c[:, 1], alpha=0.2) + plt.subplot(1, 2, 2) + plt.title("rotated vertices") plt.scatter(vr[:, 0], vr[:, 1], alpha=0.1) if plot_path: - plt.savefig(plot_path, dpi=300) - - return gc, Vt - - + plt.savefig(plot_path, dpi=300) + + return gc, Vt, pos_c + + def orient_head(img, plot_path=None): - """Euristic to orient all samples "head up": calculate sum intensity profile along the Y axis, if max is closer to 0 then do nothing, else rotate 180 degrees. + """Euristic to orient all samples "head up": calculate sum intensity profile along the Y axis, + if max is closer to 0 then do nothing, else rotate 180 degrees. Args: img: DAPI volume - + Returns: True if rotation is needed. """ assert img.ndim == 3, f"Input image should have 3 dimensions, has {img.ndim}" int_profile = np.sum(img, axis=(0, 2)) - + plt.figure() plt.plot(int_profile) plt.xlabel("Coordinate") plt.ylabel("Sum intensity along Y axis") plt.savefig(plot_path, dpi=300) - + max_pos = int_profile.argmax() logging.info(f"Max position is {max_pos}, dimension shape is {img.shape[1]}") if max_pos < img.shape[1] // 2: @@ -87,7 +94,7 @@ def orient_head(img, plot_path=None): else: logging.info("Rotate 180 degree to align head position") return True - + def orient_sample(img, dapi_chan, plot_path, dorsal=False): """Preliminary orientation of the samples with body axis along Y, head closer to 0. @@ -100,40 +107,42 @@ def orient_sample(img, dapi_chan, plot_path, dorsal=False): Returns: oriented image """ - + # Rotate around Y if the volume is dorsal plot_path = Path(plot_path) - + if dorsal: - logging.info(f"Rotated dorsal sample to align Z") + logging.info("Rotated dorsal sample to align Z") img = img[:, ::-1, :, ::-1] - - + # Rotate in XY plane, because samples can be oriented randomly, not only along X or Y max_proj = np.max(img, axis=1)[dapi_chan, ...] plt.figure() plt.imshow(max_proj, cmap="Reds") plt.savefig(plot_path / "max_proj_input.png") - + logging.info(f"Smooth image with sigma={2}") img_smoothed = gaussian(img[dapi_chan, ::10, ::10, ::10], sigma=3) - gc, Vt = get_SVD_transform(img_smoothed, plot_path / "max_proj_point_cloud_random_angle.png", percentile_trsh=90) + gc, Vt = get_SVD_transform( + img_smoothed, + plot_path / "max_proj_point_cloud_random_angle.png", + percentile_trsh=90, + ) rot_angle = 90 - np.degrees(np.arctan2(Vt[0, 1], Vt[0, 0])) logging.info(f"Rotation angle to correct for random angle is {rot_angle}") img = rotate(img, rot_angle, axes=[-2, -1], order=3, mode="constant") logging.info(f"Rotated the input volume around Z by {rot_angle} degrees") - + plt.figure() plt.imshow(np.max(img, axis=1)[dapi_chan, ...], alpha=0.5, cmap="Blues") plt.savefig(plot_path / "max_proj_rotated_random_angle.png", dpi=300) - - + # Check again if the sample is along X or along Y # Make max projection and determine the direction of principal axes max_proj = np.max(img, axis=1)[dapi_chan, ...] img_smoothed = gaussian(img[dapi_chan, ::10, ::10, ::10], sigma=2) gc, Vt = get_SVD_transform(img_smoothed, plot_path / "max_proj_point_cloud.png") - + rot_angle = np.arccos(Vt[0, 0]) logging.info(f"Rotation angle is {rot_angle}") if np.abs(rot_angle - np.pi / 2) < np.pi / 4: @@ -141,24 +150,26 @@ def orient_sample(img, dapi_chan, plot_path, dorsal=False): plt.figure() plt.imshow(max_proj, cmap="Reds") img_rot = np.rot90(img, axes=(2, 3)) - - elif np.abs(rot_angle) < np.pi / 4: + + elif np.abs(rot_angle) < np.pi / 4: logging.info("Rotation is already correct") img_rot = img else: - logging.info(f"90 degree rotation can't be determined with first principal axis angle of {rot_angle * 180/ np.pi}") + logging.info( + f"90 degree rotation can't be determined with first principal axis angle of {rot_angle * 180/ np.pi}" + ) img_rot = img - + # Check if head is oriented correctly, else rotate 180 degrees - + if orient_head(img_rot[dapi_chan, ...], plot_path / "sum_intensity_profile.png"): img_reg = np.rot90(img_rot, k=2, axes=(2, 3)) else: img_reg = img_rot - + logging.info(f"Final image shape is {img_reg.shape}") plt.figure() plt.imshow(np.max(img_reg, axis=1)[dapi_chan, ...], alpha=0.5, cmap="Blues") plt.savefig(plot_path / "max_proj_rotated_final.png", dpi=300) - - return img_reg \ No newline at end of file + + return img_reg diff --git a/matchmaker/preprocessing.py b/matchmaker/preprocessing.py new file mode 100644 index 0000000..5f51629 --- /dev/null +++ b/matchmaker/preprocessing.py @@ -0,0 +1,31 @@ +import numpy as np + + +def zero_mean_unit_variance(img): + return (img - np.mean(img)) / np.std(img) + + +def percentile_norm(img, pmin, pmax, eps=1e-10, channelwise=False): + """ + Percentile normalization. Image is supposed to be C(Z)YX + + Args: + img: _description_ + pmin: in percents + pmax: in percents + eps: _description_. Defaults to 1e-10. + + Returns: + _description_ + """ + if channelwise: + norm_img = np.zeros_like(img) + for chan in range(norm_img.shape[0]): + pmin_val = np.percentile(img[chan, ...], pmin) + pmax_val = np.percentile(img[chan, ...], pmax) + norm_img[chan, ...] = (img[chan, ...] - pmin_val) / (pmax_val - pmin_val + eps) + return norm_img + else: + pmin = np.percentile(img, pmin) + pmax = np.percentile(img, pmax) + return (img - pmin) / (pmax - pmin + eps) diff --git a/matchmaker/transform_utils.py b/matchmaker/transform_utils.py new file mode 100644 index 0000000..1e7daae --- /dev/null +++ b/matchmaker/transform_utils.py @@ -0,0 +1,122 @@ +import numpy as np +import tifffile as tif +import transforms3d as tf3d +import z5py +from scipy.ndimage import affine_transform +from elf.wrapper.resized_volume import ResizedVolume + +from matchmaker.vis import plot_three_slices, plot_overlay + + +def downscale_seg(seg, factor): + new_shape = np.array(seg.shape) // factor + downsampled_seg = ResizedVolume(seg, shape=new_shape)[:] + print("Downsampled shape", downsampled_seg.shape) + + return downsampled_seg + + +def pad_img(img): + shape = img.shape + diagonal = int(np.ceil(np.linalg.norm(shape))) + print("Diagonal", diagonal) + + # Calculate how much padding is needed on each axis + pad_widths = [] + for dim in shape: + total_pad = diagonal - dim + before = total_pad // 2 + after = total_pad - before + pad_widths.append((before, after)) + + # Apply zero padding + padded = np.pad(img, pad_width=pad_widths, mode='constant', constant_values=0) + print("New shape after padding:", padded.shape) + + return padded + + +def get_rotation_matrix(angles, save_path=None): + rotation_matrix = tf3d.euler.euler2mat(*angles, axes='szyx') + print("Rotation matrix:\n", rotation_matrix) + if save_path is not None: + np.savetxt(save_path, rotation_matrix) + return rotation_matrix + + +def rotate_img(img, rot_matrix, output_shape=None): + # Compute center + center = np.array(img.shape) / 2 # TODO: rotate around center of mass? + offset = center - rot_matrix @ center + + # rotate the image around the center + rotated_img = affine_transform( + img, + matrix=rot_matrix, + output_shape=output_shape, # new shape after rotation + offset=offset, # offset to ensure the image fits within the new bounding box + order=0, # interpolation (use 0 for discrete/label data) + mode='constant', # fill mode + cval=0.0 # fill value (if constant mode) + ) + + print(f"Shape after rotation: {rotated_img.shape}") + return rotated_img + + +def crop_to_bbox(img): + """Crops a 3D volume to the minimal bounding box around all nonzero voxels.""" + # Find where the volume is nonzero (i.e., contains instances) + nonzero = np.argwhere(img) + + # Compute bounding box from min to max index along each axis + z_min, y_min, x_min = nonzero.min(axis=0) + z_max, y_max, x_max = nonzero.max(axis=0) + 1 # +1 to include the max index + + # Crop the volume + cropped = img[z_min:z_max, y_min:y_max, x_min:x_max] + return cropped + + +# NOTE rotate with padding +def rotate_with_padding(img, rot_matrix): + padded = pad_img(img) + rotated = rotate_img(padded, rot_matrix) + cropped = crop_to_bbox(rotated) + return cropped + + +# NOTE rotate with new shape +def get_rotated_shape(img, rotation_matrix): + # Step 1: Define the 8 corners of the original volume + dz, dy, dx = img.shape + center = np.array([dz, dy, dx]) / 2 + corners = np.array([ + [0, 0, 0], + [0, 0, dx], + [0, dy, 0], + [0, dy, dx], + [dz, 0, 0], + [dz, 0, dx], + [dz, dy, 0], + [dz, dy, dx] + ]) + centered_corners = corners - center # Center the corners around the origin + # Step 2: Apply the rotation matrix + rotated_corners = corners @ rotation_matrix # NOTE ohne T? + + # Step 3: Find min and max of the rotated corners + min_coords = rotated_corners.min(axis=0) + max_coords = rotated_corners.max(axis=0) + + # Step 4: Calculate the new shape + new_shape = np.ceil(max_coords - min_coords).astype(int) + + return new_shape + + +def rotate_with_shape(img, rot_matrix): + new_shape = get_rotated_shape(img, rot_matrix) + rotated = rotate_img(img, rot_matrix, output_shape=new_shape) + + return rotated diff --git a/matchmaker/vis.py b/matchmaker/vis.py index af78530..c2a954a 100644 --- a/matchmaker/vis.py +++ b/matchmaker/vis.py @@ -1,9 +1,18 @@ import matplotlib.pyplot as plt import numpy as np -from platy_reg.preprocessing import percentile_norm +from matchmaker.preprocessing import percentile_norm -def plot_three_slices(img, save_path, x_pos=None, y_pos=None, z_pos=None, cmap="Greys_r", max_pos=False, alpha=False): +def plot_three_slices( + img, + save_path=None, + x_pos=None, + y_pos=None, + z_pos=None, + cmap="Greys_r", + max_pos=False, + alpha=False, +): """Plot slices of a 3D image along each axis. Args: @@ -22,32 +31,31 @@ def plot_three_slices(img, save_path, x_pos=None, y_pos=None, z_pos=None, cmap=" y_pos = int(img.shape[1] // 2) if z_pos is None: z_pos = int(img.shape[0] // 2) - + if max_pos: z_pos, y_pos, x_pos = np.unravel_index(np.argmax(img), img.shape) - + if alpha: alpha = (img > 0).astype(np.float32) else: alpha = np.ones_like(img) plt.figure(figsize=(15, 5)) - plt.subplot(1,3,1) - plt.title(f'z slice at {z_pos}') + plt.subplot(1, 3, 1) + plt.title(f"z slice at {z_pos}") plt.imshow(img[z_pos, :, :], cmap=cmap, alpha=alpha[z_pos, :, :]) - plt.subplot(1,3,2) - plt.title(f'y slice at {y_pos}') + plt.subplot(1, 3, 2) + plt.title(f"y slice at {y_pos}") plt.imshow(img[:, y_pos, :], cmap=cmap, alpha=alpha[:, y_pos, :]) - plt.subplot(1,3,3) - plt.title(f'x slice at {x_pos}') + plt.subplot(1, 3, 3) + plt.title(f"x slice at {x_pos}") plt.imshow(img[:, :, x_pos], cmap=cmap, alpha=alpha[:, :, x_pos]) if save_path is None: plt.show() else: plt.savefig(save_path, dpi=300) plt.close() - - - + + def plot_overlay(img1, img2, save_path=None, x_pos=None, y_pos=None, z_pos=None): """Plot slices of two 3D images along each axis. @@ -61,29 +69,29 @@ def plot_overlay(img1, img2, save_path=None, x_pos=None, y_pos=None, z_pos=None) """ assert img1.ndim == 3 assert img2.ndim == 3 - + if x_pos is None: x_pos = min(int(img1.shape[2] // 2), int(img2.shape[2] // 2)) if y_pos is None: y_pos = min(int(img1.shape[1] // 2), int(img2.shape[1] // 2)) if z_pos is None: z_pos = min(int(img1.shape[0] // 2), int(img2.shape[0] // 2)) - + plt.figure(figsize=(30, 10), dpi=300) - plt.subplot(1,3,1) - plt.title(f'z slice at {z_pos}') + plt.subplot(1, 3, 1) + plt.title(f"z slice at {z_pos}") img1_alpha = (percentile_norm(img1, 0, 100) > 0) * 0.5 img2_alpha = (percentile_norm(img2, 0, 100) > 0) * 0.5 plt.imshow(img1[z_pos, :, :], cmap="Reds", alpha=img1_alpha[z_pos, :, :]) plt.imshow(img2[z_pos, :, :], cmap="Blues", alpha=img2_alpha[z_pos, :, :]) - - plt.subplot(1,3,2) - plt.title(f'y slice at {y_pos}') + + plt.subplot(1, 3, 2) + plt.title(f"y slice at {y_pos}") plt.imshow(img1[:, y_pos, :], cmap="Reds", alpha=img1_alpha[:, y_pos, :]) plt.imshow(img2[:, y_pos, :], cmap="Blues", alpha=img2_alpha[:, y_pos, :]) - - plt.subplot(1,3,3) - plt.title(f'x slice at {x_pos}') + + plt.subplot(1, 3, 3) + plt.title(f"x slice at {x_pos}") plt.imshow(img1[:, :, x_pos], cmap="Reds", alpha=img1_alpha[:, :, x_pos]) plt.imshow(img2[:, :, x_pos], cmap="Blues", alpha=img2_alpha[:, :, x_pos]) @@ -91,4 +99,4 @@ def plot_overlay(img1, img2, save_path=None, x_pos=None, y_pos=None, z_pos=None) plt.show() else: plt.savefig(save_path, dpi=300) - plt.close() \ No newline at end of file + plt.close() From d159f7c3205cc5aab86ebee113d3ba853f0c183d Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Thu, 24 Apr 2025 13:43:40 +0200 Subject: [PATCH 07/23] try rotate image with calculating new shape --- examples/registration_test.py | 5 +++-- matchmaker/transform_utils.py | 29 ++++++++++++++--------------- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/examples/registration_test.py b/examples/registration_test.py index b10e8b4..4bd44f9 100644 --- a/examples/registration_test.py +++ b/examples/registration_test.py @@ -20,11 +20,12 @@ def main(): # seg_moving_rotated = rotate_with_padding(seg_moving, Vt.T) # TODO: add gc as rotation center? - # seg_moving_rotated2 = rotate_with_shape(seg_moving, Vt.T) + seg_moving_rotated2 = rotate_with_shape(seg_moving, np.linalg.inv(Vt)) import napari v = napari.Viewer() - v.add_image(seg_moving_rotated2, name="fixed") + v.add_labels(seg_moving) + v.add_labels(seg_moving_rotated2) napari.run() # plot_three_slices(seg_moving_rotated) diff --git a/matchmaker/transform_utils.py b/matchmaker/transform_utils.py index 1e7daae..8c70fbb 100644 --- a/matchmaker/transform_utils.py +++ b/matchmaker/transform_utils.py @@ -1,12 +1,8 @@ import numpy as np -import tifffile as tif import transforms3d as tf3d -import z5py from scipy.ndimage import affine_transform from elf.wrapper.resized_volume import ResizedVolume -from matchmaker.vis import plot_three_slices, plot_overlay - def downscale_seg(seg, factor): new_shape = np.array(seg.shape) // factor @@ -44,15 +40,15 @@ def get_rotation_matrix(angles, save_path=None): return rotation_matrix -def rotate_img(img, rot_matrix, output_shape=None): +def rotate_img(img, rotation_matrix, output_shape=None): # Compute center - center = np.array(img.shape) / 2 # TODO: rotate around center of mass? - offset = center - rot_matrix @ center + center = np.array(img.shape) / 2 + offset = center - rotation_matrix @ center # rotate the image around the center rotated_img = affine_transform( img, - matrix=rot_matrix, + matrix=rotation_matrix, output_shape=output_shape, # new shape after rotation offset=offset, # offset to ensure the image fits within the new bounding box order=0, # interpolation (use 0 for discrete/label data) @@ -90,7 +86,7 @@ def rotate_with_padding(img, rot_matrix): def get_rotated_shape(img, rotation_matrix): # Step 1: Define the 8 corners of the original volume dz, dy, dx = img.shape - center = np.array([dz, dy, dx]) / 2 + corners = np.array([ [0, 0, 0], [0, 0, dx], @@ -101,13 +97,16 @@ def get_rotated_shape(img, rotation_matrix): [dz, dy, 0], [dz, dy, dx] ]) + center = np.array([dz, dy, dx]) / 2 # NOTE: or center with the center from SVD? centered_corners = corners - center # Center the corners around the origin # Step 2: Apply the rotation matrix - rotated_corners = corners @ rotation_matrix # NOTE ohne T? + rotated_corners = centered_corners @ rotation_matrix + + final_corners = rotated_corners + center # Translate back to original position # Step 3: Find min and max of the rotated corners - min_coords = rotated_corners.min(axis=0) - max_coords = rotated_corners.max(axis=0) + min_coords = final_corners.min(axis=0) + max_coords = final_corners.max(axis=0) # Step 4: Calculate the new shape new_shape = np.ceil(max_coords - min_coords).astype(int) @@ -115,8 +114,8 @@ def get_rotated_shape(img, rotation_matrix): return new_shape -def rotate_with_shape(img, rot_matrix): - new_shape = get_rotated_shape(img, rot_matrix) - rotated = rotate_img(img, rot_matrix, output_shape=new_shape) +def rotate_with_shape(img, rotation_matrix): + new_shape = get_rotated_shape(img, rotation_matrix) + rotated = rotate_img(img, rotation_matrix, output_shape=new_shape) return rotated From 76030a465344dd5c741a4b58750312dc3bfffc72 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Thu, 24 Apr 2025 13:45:12 +0200 Subject: [PATCH 08/23] start applying transform --- examples/registration_test.py | 19 +++++++------------ matchmaker/prealignment.py | 2 +- 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/examples/registration_test.py b/examples/registration_test.py index 4bd44f9..2acef21 100644 --- a/examples/registration_test.py +++ b/examples/registration_test.py @@ -4,8 +4,6 @@ from matchmaker.prealignment import get_SVD_transform from matchmaker.vis import plot_three_slices, plot_overlay from matchmaker.transform_utils import rotate_with_padding, rotate_with_shape -from scipy.ndimage import affine_transform - def main(): @@ -15,25 +13,22 @@ def main(): with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: seg_fixed = f["seg"][:] - gc, Vt, pos_c = get_SVD_transform(seg_moving) - - + gc, Vt = get_SVD_transform(seg_moving) - # seg_moving_rotated = rotate_with_padding(seg_moving, Vt.T) # TODO: add gc as rotation center? - seg_moving_rotated2 = rotate_with_shape(seg_moving, np.linalg.inv(Vt)) + seg_moving_rotated = rotate_with_shape(seg_moving, np.linalg.inv(Vt)) import napari v = napari.Viewer() v.add_labels(seg_moving) - v.add_labels(seg_moving_rotated2) + v.add_labels(seg_moving_rotated) napari.run() # plot_three_slices(seg_moving_rotated) - # plot_overlay( - # seg_fixed, - # seg_moving_rotated2, - # ) + plot_overlay( + seg_fixed, + seg_moving_rotated, + ) if __name__ == "__main__": diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index 098d70c..683252b 100644 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -64,7 +64,7 @@ def get_SVD_transform(img, plot_path=None, percentile_trsh=90): if plot_path: plt.savefig(plot_path, dpi=300) - return gc, Vt, pos_c + return gc, Vt def orient_head(img, plot_path=None): From af6c69a9cd7bb0db0e1e7b36ae07f82b85d5991a Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Fri, 25 Apr 2025 14:27:15 +0200 Subject: [PATCH 09/23] added code to compute and apply pre-alignment of the volumes --- examples/registration_test.py | 43 +++++++++++++++++----- matchmaker/transform_utils.py | 69 +++++++++++++++++++++++------------ 2 files changed, 79 insertions(+), 33 deletions(-) diff --git a/examples/registration_test.py b/examples/registration_test.py index 2acef21..8579b0e 100644 --- a/examples/registration_test.py +++ b/examples/registration_test.py @@ -2,8 +2,8 @@ import numpy as np from matchmaker.prealignment import get_SVD_transform +from matchmaker.transform_utils import get_transformation_matrix, rotate_img, crop_to_bbox from matchmaker.vis import plot_three_slices, plot_overlay -from matchmaker.transform_utils import rotate_with_padding, rotate_with_shape def main(): @@ -13,21 +13,46 @@ def main(): with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: seg_fixed = f["seg"][:] + # 1) align moving image with PCs gc, Vt = get_SVD_transform(seg_moving) + T, new_shape = get_transformation_matrix(seg_moving, gc, Vt) + seg_moving_rotated = rotate_img(seg_moving, T, output_shape=new_shape) + seg_moving_rotated = crop_to_bbox(seg_moving_rotated) - seg_moving_rotated = rotate_with_shape(seg_moving, np.linalg.inv(Vt)) + plot_three_slices( + seg_moving_rotated, + z_pos=236, + y_pos=181, + x_pos=153, + save_path="./data/plots/moving_rotated.png", + ) + + # 2) align fixed image with PCs + gc, Vt = get_SVD_transform(seg_fixed) + T, new_shape = get_transformation_matrix(seg_fixed, gc, Vt) + seg_fixed_rotated = rotate_img(seg_fixed, T, output_shape=new_shape) + seg_fixed_rotated = crop_to_bbox(seg_fixed_rotated) + + plot_three_slices( + seg_fixed_rotated, + z_pos=236, + y_pos=181, + x_pos=153, + save_path="./data/plots/fixed_rotated.png", + ) - import napari - v = napari.Viewer() - v.add_labels(seg_moving) - v.add_labels(seg_moving_rotated) - napari.run() + # 3) test backtransform + seg_fixed_inv = rotate_img(seg_fixed_rotated, np.linalg.inv(T), output_shape=seg_fixed.shape) + # NOTE: cropped & backtransform? new_shape center is not the same anymore + # how to align the rotated images then? - # plot_three_slices(seg_moving_rotated) + plot_three_slices(seg_fixed_inv, save_path="./data/plots/fixed_rotated_inv.png") + # 4) plot overlay of prealigned volumes plot_overlay( - seg_fixed, + seg_fixed_rotated, seg_moving_rotated, + save_path="./data/plots/overlay_prealignment.png", ) diff --git a/matchmaker/transform_utils.py b/matchmaker/transform_utils.py index 8c70fbb..3db4346 100644 --- a/matchmaker/transform_utils.py +++ b/matchmaker/transform_utils.py @@ -32,20 +32,15 @@ def pad_img(img): return padded -def get_rotation_matrix(angles, save_path=None): - rotation_matrix = tf3d.euler.euler2mat(*angles, axes='szyx') - print("Rotation matrix:\n", rotation_matrix) - if save_path is not None: - np.savetxt(save_path, rotation_matrix) - return rotation_matrix +# def get_rotation_matrix(angles, save_path=None): +# rotation_matrix = tf3d.euler.euler2mat(*angles, axes='szyx') +# print("Rotation matrix:\n", rotation_matrix) +# if save_path is not None: +# np.savetxt(save_path, rotation_matrix) +# return rotation_matrix -def rotate_img(img, rotation_matrix, output_shape=None): - # Compute center - center = np.array(img.shape) / 2 - offset = center - rotation_matrix @ center - - # rotate the image around the center +def rotate_img(img, rotation_matrix, output_shape=None, offset=None): rotated_img = affine_transform( img, matrix=rotation_matrix, @@ -97,25 +92,51 @@ def get_rotated_shape(img, rotation_matrix): [dz, dy, 0], [dz, dy, dx] ]) - center = np.array([dz, dy, dx]) / 2 # NOTE: or center with the center from SVD? - centered_corners = corners - center # Center the corners around the origin - # Step 2: Apply the rotation matrix - rotated_corners = centered_corners @ rotation_matrix - final_corners = rotated_corners + center # Translate back to original position + # Step 2: Apply the rotation matrix + rotated_corners = corners @ rotation_matrix[0:3, 0:3] # Apply rotation # Step 3: Find min and max of the rotated corners - min_coords = final_corners.min(axis=0) - max_coords = final_corners.max(axis=0) + min_coords = rotated_corners.min(axis=0) + max_coords = rotated_corners.max(axis=0) # Step 4: Calculate the new shape - new_shape = np.ceil(max_coords - min_coords).astype(int) + new_shape = (np.ceil(max_coords - min_coords).astype(int)) return new_shape -def rotate_with_shape(img, rotation_matrix): - new_shape = get_rotated_shape(img, rotation_matrix) - rotated = rotate_img(img, rotation_matrix, output_shape=new_shape) +# def get_affine_transform(translation, rotation, rotation_order="szyx"): +# M = get_translation_matrix(translation) +# M[0:3, 0:3] = tf3d.euler.euler2mat(*rotation, rotation_order) +# return M + + +def get_translation_matrix(translation): + M = np.identity(4) + M[0:3, 3] = translation + return M + + +def get_rotation_matrix(R): + M = np.identity(4) + M[0:3, 0:3] = R + return M + + +### now again: what do I need? + +def get_transformation_matrix(img, gc, Vt): + # 1. center image on origin + center_to_origin = get_translation_matrix(gc) + # 2. rotate image + rot = get_rotation_matrix(Vt.T) + # 3. get new shape + new_shape = get_rotated_shape(img, rot) + # 4. center image on new shape + new_shape_center = np.array(new_shape) // 2 + center_to_new_shape = get_translation_matrix(-new_shape_center) + # 5. combine all transforms: get transformation matrix + R = center_to_origin @ rot @ center_to_new_shape - return rotated + return R, new_shape From 0c14dfdebf92ea036a53fba8a23373fa853e6ce8 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Mon, 28 Apr 2025 13:24:56 +0200 Subject: [PATCH 10/23] calculated pre-alignment with pca of toy example --- .gitignore | 2 ++ examples/data/rotation_matrix.txt | 3 +++ examples/registration_test.py | 13 ++----------- 3 files changed, 7 insertions(+), 11 deletions(-) create mode 100644 examples/data/rotation_matrix.txt diff --git a/.gitignore b/.gitignore index 54a3041..b2e0298 100644 --- a/.gitignore +++ b/.gitignore @@ -175,5 +175,7 @@ cython_debug/ # huge data ./examples/data/platy1_muscles_stardist.tif +*.n5/ +*.png .DS_Store \ No newline at end of file diff --git a/examples/data/rotation_matrix.txt b/examples/data/rotation_matrix.txt new file mode 100644 index 0000000..207754e --- /dev/null +++ b/examples/data/rotation_matrix.txt @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:dcdd4056e70e24809604b56e3c8efb544a04faeb8b2432b6a8feeb71e484d6b9 +size 231 diff --git a/examples/registration_test.py b/examples/registration_test.py index 8579b0e..ed80dd2 100644 --- a/examples/registration_test.py +++ b/examples/registration_test.py @@ -2,7 +2,7 @@ import numpy as np from matchmaker.prealignment import get_SVD_transform -from matchmaker.transform_utils import get_transformation_matrix, rotate_img, crop_to_bbox +from matchmaker.transform_utils import get_transformation_matrix, rotate_img from matchmaker.vis import plot_three_slices, plot_overlay @@ -17,13 +17,9 @@ def main(): gc, Vt = get_SVD_transform(seg_moving) T, new_shape = get_transformation_matrix(seg_moving, gc, Vt) seg_moving_rotated = rotate_img(seg_moving, T, output_shape=new_shape) - seg_moving_rotated = crop_to_bbox(seg_moving_rotated) plot_three_slices( seg_moving_rotated, - z_pos=236, - y_pos=181, - x_pos=153, save_path="./data/plots/moving_rotated.png", ) @@ -31,20 +27,15 @@ def main(): gc, Vt = get_SVD_transform(seg_fixed) T, new_shape = get_transformation_matrix(seg_fixed, gc, Vt) seg_fixed_rotated = rotate_img(seg_fixed, T, output_shape=new_shape) - seg_fixed_rotated = crop_to_bbox(seg_fixed_rotated) + # seg_fixed_rotated = crop_to_bbox(seg_fixed_rotated) plot_three_slices( seg_fixed_rotated, - z_pos=236, - y_pos=181, - x_pos=153, save_path="./data/plots/fixed_rotated.png", ) # 3) test backtransform seg_fixed_inv = rotate_img(seg_fixed_rotated, np.linalg.inv(T), output_shape=seg_fixed.shape) - # NOTE: cropped & backtransform? new_shape center is not the same anymore - # how to align the rotated images then? plot_three_slices(seg_fixed_inv, save_path="./data/plots/fixed_rotated_inv.png") From d50c75ecaca4502530af8b18fe20eb12124c36a7 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Mon, 28 Apr 2025 14:51:50 +0200 Subject: [PATCH 11/23] update deformation of test data --- environment.yml | 1 + examples/.gitattributes | 2 - examples/data/transformation_matrix.txt | 3 ++ examples/deform_test_data.py | 34 ++++++------- examples/registration_test.py | 19 +++++-- examples/reverse_deformation.py | 20 ++++---- examples/reverse_transform.py | 0 examples/u.py | 57 --------------------- matchmaker/prealignment.py | 41 +++++---------- matchmaker/transform_utils.py | 67 +++++++++---------------- 10 files changed, 86 insertions(+), 158 deletions(-) delete mode 100644 examples/.gitattributes create mode 100644 examples/data/transformation_matrix.txt delete mode 100644 examples/reverse_transform.py delete mode 100644 examples/u.py diff --git a/environment.yml b/environment.yml index b5abc49..645abb6 100644 --- a/environment.yml +++ b/environment.yml @@ -16,6 +16,7 @@ dependencies: - optuna - proxsuite - cvxpy + - transforms3d - ipykernel - seaborn diff --git a/examples/.gitattributes b/examples/.gitattributes deleted file mode 100644 index c2e5313..0000000 --- a/examples/.gitattributes +++ /dev/null @@ -1,2 +0,0 @@ -*.tif filter=lfs diff=lfs merge=lfs -text -*.n5 filter=lfs diff=lfs merge=lfs -text diff --git a/examples/data/transformation_matrix.txt b/examples/data/transformation_matrix.txt new file mode 100644 index 0000000..a4ae13a --- /dev/null +++ b/examples/data/transformation_matrix.txt @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d5e9bd5d72c68d4b8fa4d68b37382c836e24bdf87e0893222cdb564c8bf9ab3e +size 406 diff --git a/examples/deform_test_data.py b/examples/deform_test_data.py index f8b94d5..7173aa7 100644 --- a/examples/deform_test_data.py +++ b/examples/deform_test_data.py @@ -1,20 +1,12 @@ import numpy as np import tifffile as tif import z5py -from elf.wrapper.resized_volume import ResizedVolume +import transforms3d as tf3d -from matchmaker.transform_utils import pad_img, get_rotation_matrix, rotate_img, crop_to_bbox +from matchmaker.transform_utils import downscale_seg, get_transformation_matrix, rotate_img, crop_to_bbox from matchmaker.vis import plot_three_slices, plot_overlay -def downscale_seg(seg, factor): - new_shape = np.array(seg.shape) // 2 - downsampled_seg = ResizedVolume(seg, shape=new_shape)[:] - print("Downsampled shape", downsampled_seg.shape) - - return downsampled_seg - - def remove_instances(seg, prob=0.05): # Get all unique instance IDs, excluding background (assumed to be 0) instance_ids = np.unique(seg) @@ -41,27 +33,35 @@ def main(): seg = downscale_seg(seg, factor) seg_fixed = crop_to_bbox(seg) print("Cropped shape", seg_fixed.shape) + # save downsampled image with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "w") as f: f.create_dataset("seg", data=seg_fixed) # rotate image - seg_padded = pad_img(seg) - angles = [np.deg2rad(155), np.deg2rad(30), np.deg2rad(65)] # z,y,x - rotation_matrix = get_rotation_matrix(angles, save_path="./data/rotation_matrix.txt") - rotated_seg = rotate_img(seg_padded, rotation_matrix) - seg_cropped = crop_to_bbox(rotated_seg) - print("Cropped shape", seg_cropped.shape) + center = np.array(seg_fixed.shape) // 2 + rotation = tf3d.euler.euler2mat(*[np.deg2rad(155), np.deg2rad(30), np.deg2rad(65)], axes="szyx") # z,y,x + + T, new_shape = get_transformation_matrix(seg_fixed, center, rotation, save_path="./data/transformation_matrix.txt") + seg_moving = rotate_img(seg_fixed, T, output_shape=new_shape) + + # seg_padded = pad_img(seg) + # angles = [np.deg2rad(155), np.deg2rad(30), np.deg2rad(65)] # z,y,x + # rotation_matrix = get_rotation_matrix(angles, save_path="./data/rotation_matrix.txt") + # rotated_seg = rotate_img(seg_padded, rotation_matrix) + # seg_cropped = crop_to_bbox(rotated_seg) + # print("Cropped shape", seg_cropped.shape) # randomly remove instances probability = 0.05 - seg_moving = remove_instances(seg_cropped, prob=probability) + seg_moving = remove_instances(seg_moving, prob=probability) with z5py.File("./data/platy1_muscles_stardist_moving.n5", "w") as f: f.create_dataset("seg", data=seg_moving) # visualize plot_three_slices(seg_moving) + plot_three_slices(seg_fixed) plot_overlay(seg_fixed, seg_moving) diff --git a/examples/registration_test.py b/examples/registration_test.py index ed80dd2..154242c 100644 --- a/examples/registration_test.py +++ b/examples/registration_test.py @@ -1,7 +1,7 @@ import z5py import numpy as np -from matchmaker.prealignment import get_SVD_transform +from matchmaker.prealignment import get_SVD_transform, orient_head from matchmaker.transform_utils import get_transformation_matrix, rotate_img from matchmaker.vis import plot_three_slices, plot_overlay @@ -10,6 +10,8 @@ def main(): with z5py.File("./data/platy1_muscles_stardist_moving.n5", "r") as f: seg_moving = f["seg"][:] + plot_three_slices(seg_moving, save_path="./data/plots/moving.png") + with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: seg_fixed = f["seg"][:] @@ -18,6 +20,13 @@ def main(): T, new_shape = get_transformation_matrix(seg_moving, gc, Vt) seg_moving_rotated = rotate_img(seg_moving, T, output_shape=new_shape) + # Check if head is oriented correctly, else rotate 180 degrees + rotate_head = orient_head(seg_moving_rotated, save_path="./data/plots/moving_orient_head.png") + # NOTE: already enough? + if rotate_head: + print("Rotate head 180 degrees ...") + seg_moving_rotated = np.rot90(seg_moving_rotated, k=2) + plot_three_slices( seg_moving_rotated, save_path="./data/plots/moving_rotated.png", @@ -27,7 +36,11 @@ def main(): gc, Vt = get_SVD_transform(seg_fixed) T, new_shape = get_transformation_matrix(seg_fixed, gc, Vt) seg_fixed_rotated = rotate_img(seg_fixed, T, output_shape=new_shape) - # seg_fixed_rotated = crop_to_bbox(seg_fixed_rotated) + + rotate_head = orient_head(seg_fixed_rotated) + if rotate_head: + print("Rotate head 180 degrees ...") + seg_fixed_rotated = np.rot90(seg_fixed_rotated, k=2) plot_three_slices( seg_fixed_rotated, @@ -35,7 +48,7 @@ def main(): ) # 3) test backtransform - seg_fixed_inv = rotate_img(seg_fixed_rotated, np.linalg.inv(T), output_shape=seg_fixed.shape) + seg_fixed_inv = rotate_img(seg_fixed_rotated, np.linalg.inv(T), output_shape=seg_fixed.shape) plot_three_slices(seg_fixed_inv, save_path="./data/plots/fixed_rotated_inv.png") diff --git a/examples/reverse_deformation.py b/examples/reverse_deformation.py index 708781a..5f163bb 100644 --- a/examples/reverse_deformation.py +++ b/examples/reverse_deformation.py @@ -1,26 +1,28 @@ import numpy as np import z5py from matchmaker.vis import plot_three_slices, plot_overlay -from deform_test_data import pad_seg, crop_to_bbox, rotate_seg +from matchmaker.transform_utils import rotate_img def main(): with z5py.File("./data/platy1_muscles_stardist_moving.n5", "r") as f: - seg = f["seg"][:] + seg_moving = f["seg"][:] - seg_padded = pad_seg(seg) - rotation_matrix = np.loadtxt("./data/rotation_matrix.txt") - rotated_seg = rotate_seg(seg_padded, np.linalg.inv(rotation_matrix)) - seg_cropped = crop_to_bbox(rotated_seg) - print("Cropped shape", seg_cropped.shape) + # compare with original image + with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: + seg_fixed = f["seg"][:] + + T = np.loadtxt("./data/transformation_matrix.txt") + seg_moving_reverse = rotate_img(seg_moving, np.linalg.inv(T), output_shape=seg_fixed.shape) - plot_three_slices(seg_cropped) + plot_three_slices(seg_moving_reverse) # compare with original image with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: seg_fixed = f["seg"][:] - plot_overlay(seg_cropped, seg_fixed) + plot_three_slices(seg_fixed) + plot_overlay(seg_fixed, seg_moving_reverse) if __name__ == "__main__": diff --git a/examples/reverse_transform.py b/examples/reverse_transform.py deleted file mode 100644 index e69de29..0000000 diff --git a/examples/u.py b/examples/u.py deleted file mode 100644 index 06510d0..0000000 --- a/examples/u.py +++ /dev/null @@ -1,57 +0,0 @@ - -from transforms3d.euler import euler2mat, mat2euler -from scipy.ndimage import affine_transform -import numpy as np - -def get_translation_transform(translation): - M = np.identity(4) - M[0:3, 3] = translation - return M - - -def get_affine_transform(translation, rotation, rotation_order="szyx"): - M = get_translation_transform(translation) - M[0:3, 0:3] = euler2mat(*rotation, rotation_order) - return M - - -def flip_dorsal(img, n5_fname): - """ - For dorsal-direction samples flip along z so they are same orientation - """ - if "dorsal" in n5_fname.lower(): - print("Flip dorsal image") - img_flipped = np.flip(img, axis=0) - return img_flipped - - else: - return img - -def move_along_z(margin, image_shape): - """ - After rotation the volume might not fit within the same Z range, - so move it up by some margin and make output image size bigger - """ - assert len(image_shape) == 3 - R_t = get_translation_transform([-margin, -margin, -margin]) - new_image_shape = list(image_shape) - new_image_shape = [sh + 2 * margin for sh in image_shape] - - return R_t, new_image_shape - - -R_t = get_translation_transform(-gc) -R_t_r = get_translation_transform(gc) - -R_rot = get_translation_transform([0, 0, 0]) -R_rot[0:3, 0:3] = np.linalg.inv(Vt) - -R_rot_90_z = get_affine_transform([0, 0, 0], [0, -np.pi / 2, 0], rotation_order="szyx") - -R_t_margin, new_image_shape = move_along_z(80, norm_img.shape) - -mat = R_t_r @ R_rot @ R_t @ R_t_margin -print(mat) -print(norm_img.shape) - -img_reg = affine_transform(img_isotropic, mat, order=0, output_shape=new_image_shape) \ No newline at end of file diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index 683252b..6b99226 100644 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -8,7 +8,7 @@ from matchmaker.data import create_point_cloud -def get_SVD_transform(img, plot_path=None, percentile_trsh=90): +def get_SVD_transform(img, save_path=None, percentile_trsh=90): """Convert image to point cloud by thresholding, then run SVD on resulting point cloud. Args: @@ -21,27 +21,12 @@ def get_SVD_transform(img, plot_path=None, percentile_trsh=90): """ pos, _ = create_point_cloud(img) - # how many pos? subsample? - # center coord with center of the image - - # else: - # trsh = np.percentile(img, percentile_trsh) - # logging.info(f"Threshold used for converting to point cloud: {trsh}") - # X = convert_to_point_cloud(img, trsh) gc = pos.mean(axis=0) gc = np.array(img.shape) // 2 pos_c = pos - gc logging.info(f"Point cloud shape {pos.shape}") logging.info(f"Point cloud center {gc}") - # if X.shape[1] == 3: - # random_subset = np.random.choice(X.shape[0], 10000, replace=False) - # X_subset = X_c[random_subset, 1:] - # else: - # X_subset = X_c[::100] - # logging.info("Subset point cloud") - # logging.info(f"Subset point cloud shape {X_subset.shape}") - logging.info("Run SVD") U, S, Vt = np.linalg.svd(pos_c, full_matrices=False) logging.info("U") @@ -61,13 +46,13 @@ def get_SVD_transform(img, plot_path=None, percentile_trsh=90): plt.subplot(1, 2, 2) plt.title("rotated vertices") plt.scatter(vr[:, 0], vr[:, 1], alpha=0.1) - if plot_path: - plt.savefig(plot_path, dpi=300) + if save_path: + plt.savefig(save_path, dpi=300) return gc, Vt -def orient_head(img, plot_path=None): +def orient_head(img, save_path=None): """Euristic to orient all samples "head up": calculate sum intensity profile along the Y axis, if max is closer to 0 then do nothing, else rotate 180 degrees. @@ -84,7 +69,7 @@ def orient_head(img, plot_path=None): plt.plot(int_profile) plt.xlabel("Coordinate") plt.ylabel("Sum intensity along Y axis") - plt.savefig(plot_path, dpi=300) + plt.savefig(save_path, dpi=300) max_pos = int_profile.argmax() logging.info(f"Max position is {max_pos}, dimension shape is {img.shape[1]}") @@ -96,7 +81,7 @@ def orient_head(img, plot_path=None): return True -def orient_sample(img, dapi_chan, plot_path, dorsal=False): +def orient_sample(img, dapi_chan, save_path, dorsal=False): """Preliminary orientation of the samples with body axis along Y, head closer to 0. Args: @@ -109,7 +94,7 @@ def orient_sample(img, dapi_chan, plot_path, dorsal=False): """ # Rotate around Y if the volume is dorsal - plot_path = Path(plot_path) + save_path = Path(save_path) if dorsal: logging.info("Rotated dorsal sample to align Z") @@ -119,13 +104,13 @@ def orient_sample(img, dapi_chan, plot_path, dorsal=False): max_proj = np.max(img, axis=1)[dapi_chan, ...] plt.figure() plt.imshow(max_proj, cmap="Reds") - plt.savefig(plot_path / "max_proj_input.png") + plt.savefig(save_path / "max_proj_input.png") logging.info(f"Smooth image with sigma={2}") img_smoothed = gaussian(img[dapi_chan, ::10, ::10, ::10], sigma=3) gc, Vt = get_SVD_transform( img_smoothed, - plot_path / "max_proj_point_cloud_random_angle.png", + save_path / "max_proj_point_cloud_random_angle.png", percentile_trsh=90, ) rot_angle = 90 - np.degrees(np.arctan2(Vt[0, 1], Vt[0, 0])) @@ -135,13 +120,13 @@ def orient_sample(img, dapi_chan, plot_path, dorsal=False): plt.figure() plt.imshow(np.max(img, axis=1)[dapi_chan, ...], alpha=0.5, cmap="Blues") - plt.savefig(plot_path / "max_proj_rotated_random_angle.png", dpi=300) + plt.savefig(save_path / "max_proj_rotated_random_angle.png", dpi=300) # Check again if the sample is along X or along Y # Make max projection and determine the direction of principal axes max_proj = np.max(img, axis=1)[dapi_chan, ...] img_smoothed = gaussian(img[dapi_chan, ::10, ::10, ::10], sigma=2) - gc, Vt = get_SVD_transform(img_smoothed, plot_path / "max_proj_point_cloud.png") + gc, Vt = get_SVD_transform(img_smoothed, save_path / "max_proj_point_cloud.png") rot_angle = np.arccos(Vt[0, 0]) logging.info(f"Rotation angle is {rot_angle}") @@ -162,7 +147,7 @@ def orient_sample(img, dapi_chan, plot_path, dorsal=False): # Check if head is oriented correctly, else rotate 180 degrees - if orient_head(img_rot[dapi_chan, ...], plot_path / "sum_intensity_profile.png"): + if orient_head(img_rot[dapi_chan, ...], save_path / "sum_intensity_profile.png"): img_reg = np.rot90(img_rot, k=2, axes=(2, 3)) else: img_reg = img_rot @@ -170,6 +155,6 @@ def orient_sample(img, dapi_chan, plot_path, dorsal=False): logging.info(f"Final image shape is {img_reg.shape}") plt.figure() plt.imshow(np.max(img_reg, axis=1)[dapi_chan, ...], alpha=0.5, cmap="Blues") - plt.savefig(plot_path / "max_proj_rotated_final.png", dpi=300) + plt.savefig(save_path / "max_proj_rotated_final.png", dpi=300) return img_reg diff --git a/matchmaker/transform_utils.py b/matchmaker/transform_utils.py index 3db4346..b07b750 100644 --- a/matchmaker/transform_utils.py +++ b/matchmaker/transform_utils.py @@ -1,6 +1,6 @@ import numpy as np -import transforms3d as tf3d from scipy.ndimage import affine_transform +import transforms3d as tf3d from elf.wrapper.resized_volume import ResizedVolume @@ -32,29 +32,6 @@ def pad_img(img): return padded -# def get_rotation_matrix(angles, save_path=None): -# rotation_matrix = tf3d.euler.euler2mat(*angles, axes='szyx') -# print("Rotation matrix:\n", rotation_matrix) -# if save_path is not None: -# np.savetxt(save_path, rotation_matrix) -# return rotation_matrix - - -def rotate_img(img, rotation_matrix, output_shape=None, offset=None): - rotated_img = affine_transform( - img, - matrix=rotation_matrix, - output_shape=output_shape, # new shape after rotation - offset=offset, # offset to ensure the image fits within the new bounding box - order=0, # interpolation (use 0 for discrete/label data) - mode='constant', # fill mode - cval=0.0 # fill value (if constant mode) - ) - - print(f"Shape after rotation: {rotated_img.shape}") - return rotated_img - - def crop_to_bbox(img): """Crops a 3D volume to the minimal bounding box around all nonzero voxels.""" # Find where the volume is nonzero (i.e., contains instances) @@ -69,15 +46,6 @@ def crop_to_bbox(img): return cropped -# NOTE rotate with padding -def rotate_with_padding(img, rot_matrix): - padded = pad_img(img) - rotated = rotate_img(padded, rot_matrix) - cropped = crop_to_bbox(rotated) - return cropped - - -# NOTE rotate with new shape def get_rotated_shape(img, rotation_matrix): # Step 1: Define the 8 corners of the original volume dz, dy, dx = img.shape @@ -106,27 +74,24 @@ def get_rotated_shape(img, rotation_matrix): return new_shape -# def get_affine_transform(translation, rotation, rotation_order="szyx"): -# M = get_translation_matrix(translation) -# M[0:3, 0:3] = tf3d.euler.euler2mat(*rotation, rotation_order) -# return M - - def get_translation_matrix(translation): M = np.identity(4) M[0:3, 3] = translation return M +def get_rotation(angles): + R = tf3d.euler.euler2mat(*angles, axes='szyx') + return R + + def get_rotation_matrix(R): M = np.identity(4) M[0:3, 0:3] = R return M -### now again: what do I need? - -def get_transformation_matrix(img, gc, Vt): +def get_transformation_matrix(img, gc, Vt, save_path=None): # 1. center image on origin center_to_origin = get_translation_matrix(gc) # 2. rotate image @@ -139,4 +104,22 @@ def get_transformation_matrix(img, gc, Vt): # 5. combine all transforms: get transformation matrix R = center_to_origin @ rot @ center_to_new_shape + if save_path is not None: + np.savetxt(save_path, R) + return R, new_shape + + +def rotate_img(img, rotation_matrix, output_shape=None, offset=None): + rotated_img = affine_transform( + img, + matrix=rotation_matrix, + output_shape=output_shape, # new shape after rotation + offset=offset, # offset to ensure the image fits within the new bounding box + order=0, # interpolation (use 0 for discrete/label data) + mode='constant', # fill mode + cval=0.0 # fill value (if constant mode) + ) + + print(f"Shape after rotation: {rotated_img.shape}") + return rotated_img From 2fc7c6f6c7f3325c4ee72c55dcaa6171082a1238 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Thu, 5 Jun 2025 13:47:28 +0200 Subject: [PATCH 12/23] update gitignore --- .gitignore | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index b2e0298..d73e375 100644 --- a/.gitignore +++ b/.gitignore @@ -178,4 +178,6 @@ cython_debug/ *.n5/ *.png -.DS_Store \ No newline at end of file +.DS_Store + +./examples/data/test/* \ No newline at end of file From 5cb17c2a152d2708181ae3880fc9ebce5430d9eb Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Thu, 5 Jun 2025 13:48:44 +0200 Subject: [PATCH 13/23] added head orientation alignment to prealignment and transformation matrix --- environment.yml | 1 + examples/deform_test_data.py | 34 +++++--- examples/registration_test.py | 71 ++++++++--------- matchmaker/align_rigid_elastix.py | 2 +- matchmaker/file_test.py | 9 +++ matchmaker/mobie_export.py | 79 ++++++++++++++++++ matchmaker/n5_utils.py | 73 +++++++++++++++++ matchmaker/prealignment.py | 128 +++++++++++++----------------- matchmaker/transform_utils.py | 6 +- 9 files changed, 275 insertions(+), 128 deletions(-) create mode 100644 matchmaker/file_test.py create mode 100644 matchmaker/mobie_export.py create mode 100644 matchmaker/n5_utils.py mode change 100644 => 100755 matchmaker/prealignment.py diff --git a/environment.yml b/environment.yml index 645abb6..cc8a81f 100644 --- a/environment.yml +++ b/environment.yml @@ -30,5 +30,6 @@ dependencies: - pyscipopt - click - pyyaml + - itk diff --git a/examples/deform_test_data.py b/examples/deform_test_data.py index 7173aa7..005c64f 100644 --- a/examples/deform_test_data.py +++ b/examples/deform_test_data.py @@ -3,7 +3,12 @@ import z5py import transforms3d as tf3d -from matchmaker.transform_utils import downscale_seg, get_transformation_matrix, rotate_img, crop_to_bbox +from matchmaker.transform_utils import ( + downscale_seg, + get_transformation_matrix, + rotate_img, + crop_to_bbox, +) from matchmaker.vis import plot_three_slices, plot_overlay @@ -38,20 +43,20 @@ def main(): with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "w") as f: f.create_dataset("seg", data=seg_fixed) + # save also as tiff + tif.imwrite("./data/platy1_muscles_stardist_fixed.tif", seg_fixed) + # rotate image center = np.array(seg_fixed.shape) // 2 - rotation = tf3d.euler.euler2mat(*[np.deg2rad(155), np.deg2rad(30), np.deg2rad(65)], axes="szyx") # z,y,x + rotation = tf3d.euler.euler2mat( + *[np.deg2rad(155), np.deg2rad(30), np.deg2rad(65)], axes="szyx" + ) - T, new_shape = get_transformation_matrix(seg_fixed, center, rotation, save_path="./data/transformation_matrix.txt") + T, new_shape = get_transformation_matrix( + seg_fixed, center, rotation, save_path="./data/transformation_matrix.txt" + ) seg_moving = rotate_img(seg_fixed, T, output_shape=new_shape) - # seg_padded = pad_img(seg) - # angles = [np.deg2rad(155), np.deg2rad(30), np.deg2rad(65)] # z,y,x - # rotation_matrix = get_rotation_matrix(angles, save_path="./data/rotation_matrix.txt") - # rotated_seg = rotate_img(seg_padded, rotation_matrix) - # seg_cropped = crop_to_bbox(rotated_seg) - # print("Cropped shape", seg_cropped.shape) - # randomly remove instances probability = 0.05 seg_moving = remove_instances(seg_moving, prob=probability) @@ -59,10 +64,13 @@ def main(): with z5py.File("./data/platy1_muscles_stardist_moving.n5", "w") as f: f.create_dataset("seg", data=seg_moving) + # save also as tiff + tif.imwrite("./data/platy1_muscles_stardist_moving.tif", seg_moving) + # visualize - plot_three_slices(seg_moving) - plot_three_slices(seg_fixed) - plot_overlay(seg_fixed, seg_moving) + plot_three_slices(seg_moving, save_path="./data/plots/seg_moving.png") + plot_three_slices(seg_fixed, save_path="./data/plots/seg_fixed.png") + plot_overlay(seg_fixed, seg_moving, save_path="./data/plots/seg_overlay.png") if __name__ == "__main__": diff --git a/examples/registration_test.py b/examples/registration_test.py index 154242c..2f634e0 100644 --- a/examples/registration_test.py +++ b/examples/registration_test.py @@ -1,62 +1,55 @@ +import os import z5py import numpy as np -from matchmaker.prealignment import get_SVD_transform, orient_head +from matchmaker.prealignment import prealign_sample from matchmaker.transform_utils import get_transformation_matrix, rotate_img from matchmaker.vis import plot_three_slices, plot_overlay def main(): - with z5py.File("./data/platy1_muscles_stardist_moving.n5", "r") as f: + output_dir = "./data/test" + if not os.path.exists(f"{output_dir}/plots"): + os.makedirs(f"{output_dir}/plots") + + # prealign moving image + moving_input = "./data/platy1_muscles_stardist_moving.n5" + with z5py.File(moving_input, "r") as f: seg_moving = f["seg"][:] - plot_three_slices(seg_moving, save_path="./data/plots/moving.png") + plot_three_slices(seg_moving, save_path=f"{output_dir}/plots/moving.png") + + seg_moving_prealigned = prealign_sample(seg_moving, file_name="moving", save_path=output_dir) + + with z5py.File(f"{output_dir}/moving_prealigned.n5", "w") as f: + f.create_dataset("seg", data=seg_moving_prealigned, compression="gzip") - with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: + # prealign fixed image + fixed_input = "./data/platy1_muscles_stardist_fixed.n5" + with z5py.File(fixed_input, "r") as f: seg_fixed = f["seg"][:] - # 1) align moving image with PCs - gc, Vt = get_SVD_transform(seg_moving) - T, new_shape = get_transformation_matrix(seg_moving, gc, Vt) - seg_moving_rotated = rotate_img(seg_moving, T, output_shape=new_shape) - - # Check if head is oriented correctly, else rotate 180 degrees - rotate_head = orient_head(seg_moving_rotated, save_path="./data/plots/moving_orient_head.png") - # NOTE: already enough? - if rotate_head: - print("Rotate head 180 degrees ...") - seg_moving_rotated = np.rot90(seg_moving_rotated, k=2) - - plot_three_slices( - seg_moving_rotated, - save_path="./data/plots/moving_rotated.png", - ) + plot_three_slices(seg_fixed, save_path=f"{output_dir}/plots/fixed.png") - # 2) align fixed image with PCs - gc, Vt = get_SVD_transform(seg_fixed) - T, new_shape = get_transformation_matrix(seg_fixed, gc, Vt) - seg_fixed_rotated = rotate_img(seg_fixed, T, output_shape=new_shape) + # prealign moving image + seg_fixed_prealigned = prealign_sample(seg_fixed, file_name="fixed", save_path=output_dir) - rotate_head = orient_head(seg_fixed_rotated) - if rotate_head: - print("Rotate head 180 degrees ...") - seg_fixed_rotated = np.rot90(seg_fixed_rotated, k=2) + with z5py.File(f"{output_dir}/fixed_prealigned.n5", "w") as f: + f.create_dataset("seg", data=seg_fixed_prealigned, compression="gzip") - plot_three_slices( - seg_fixed_rotated, - save_path="./data/plots/fixed_rotated.png", + # test backtransform + T = np.loadtxt(f"{output_dir}/fixed_T_prealignment.txt") + seg_fixed_inv = rotate_img( + seg_fixed_prealigned, np.linalg.inv(T), output_shape=seg_fixed.shape ) - # 3) test backtransform - seg_fixed_inv = rotate_img(seg_fixed_rotated, np.linalg.inv(T), output_shape=seg_fixed.shape) - - plot_three_slices(seg_fixed_inv, save_path="./data/plots/fixed_rotated_inv.png") + plot_three_slices(seg_fixed_inv, save_path=f"{output_dir}/plots/fixed_prealigned_inv.png") - # 4) plot overlay of prealigned volumes + # plot overlay of prealigned volumes plot_overlay( - seg_fixed_rotated, - seg_moving_rotated, - save_path="./data/plots/overlay_prealignment.png", + seg_fixed_prealigned, + seg_moving_prealigned, + save_path=f"{output_dir}/plots/overlay_prealignment.png", ) diff --git a/matchmaker/align_rigid_elastix.py b/matchmaker/align_rigid_elastix.py index 72c0e35..2bacee8 100644 --- a/matchmaker/align_rigid_elastix.py +++ b/matchmaker/align_rigid_elastix.py @@ -1,7 +1,7 @@ from pathlib import Path import numpy as np import argparse -from platy_reg.n5_utils import read_volume, write_volume, get_attrs +from matchmaker.n5_utils import read_volume, write_volume, get_attrs from matchmaker.vis import plot_overlay import logging import sys diff --git a/matchmaker/file_test.py b/matchmaker/file_test.py new file mode 100644 index 0000000..8ac0f47 --- /dev/null +++ b/matchmaker/file_test.py @@ -0,0 +1,9 @@ +import z5py +import napari + +with z5py.File("/Users/marei/git-repositories/matchmaker/examples/data/test/mobie_project/platy1_muscles_stardist/images/ome-zarr/original.ome.zarr", "r") as f: + seg = f["s0"][:] + +v = napari.Viewer() +v.add_image(seg) +napari.run() \ No newline at end of file diff --git a/matchmaker/mobie_export.py b/matchmaker/mobie_export.py new file mode 100644 index 0000000..e6cd503 --- /dev/null +++ b/matchmaker/mobie_export.py @@ -0,0 +1,79 @@ +# export results as mobie project +# menu 1: dataset +# menu 2: fixed_segmentation; then original, pre-aligned, ... +# menu 3: moving_segmentation; then original, pre-aligned, ... +import os +import mobie +from matchmaker.n5_utils import get_attrs +# from mobie import add_segmentation + + +# def add_segmentation_to_mobie(input_path, input_key, seg_name): + +# resolution = get_attrs(input_path, input_key)["resolution"] +# unit = "micrometer" +# chunks = (64, 64, 64) +# scale_factors = 4 * [[2, 2, 2]] + +# chunks = (128, 128, 128) +# max_jobs = 8 + +# add_segmentation( +# input_path=input_path, +# input_key=input_key, +# root=MOBIE_FOLDER, +# dataset_name=DS_NAME, +# segmentation_name=seg_name, +# resolution=resolution, +# scale_factors=scale_factors, +# chunks=chunks, +# file_format="ome.zarr", +# max_jobs=max_jobs, +# ) + +# def process_segmentation(input_path, input_key, seg_name, scale): +# metadata = read_dataset_metadata(os.path.join(MOBIE_FOLDER, DS_NAME)) +# if seg_name in metadata["sources"]: +# update_segmentation(input_path, input_key, seg_name, scale) +# else: +# add_segmentation_to_mobie(input_path, input_key, seg_name, scale) + + +def create_mobie_project(input_path, input_key, output_dir): + if not os.path.exists(f"{output_dir}/mobie_project"): + os.makedirs(f"{output_dir}/mobie_project") + + # Set parameters for MOBIE + mobie_folder = f"{output_dir}/mobie_project" + dataset_name = "platy1_muscles_stardist" + breakpoint() + # resolution = get_attrs(input_path, input_key)["resolution"] + chunks = (64, 64, 64) + scale_factors = 4 * [[2, 2, 2]] + + mobie.add_segmentation( + input_path=input_path, + input_key=input_key, + root=mobie_folder, + dataset_name=dataset_name, + segmentation_name="original", + resolution=[1,1,1], + scale_factors=scale_factors, + chunks=chunks, + menu_name="fixed", + file_format="ome.zarr", + is_default_dataset=True + ) + + +def main(): + input_path = "/Users/marei/git-repositories/matchmaker/examples/data/platy1_muscles_stardist_fixed.n5" + input_key = "seg" + output_dir = "/Users/marei/git-repositories/matchmaker/examples/data/test" + + create_mobie_project(input_path, input_key, output_dir) + print(f"MoBIE project created at {output_dir}/mobie_project") + + +if __name__ == "__main__": + main() diff --git a/matchmaker/n5_utils.py b/matchmaker/n5_utils.py new file mode 100644 index 0000000..a899bc2 --- /dev/null +++ b/matchmaker/n5_utils.py @@ -0,0 +1,73 @@ +import z5py +from pathlib import PurePath +import numpy as np + + +def print_key_tree(f: z5py.File): + print(f"Key structure of z5 file {f.filename}") + f.visititems(lambda name, obj: print(name)) + + +def read_volume( + f: z5py.File, key: str, roi: np.lib.index_tricks.IndexExpression = np.s_[:] +): + print(type(f)) + if isinstance(f, (str, PurePath)): + f = z5py.File(f, "r") + + try: + ds = f[key] + except KeyError: + print(f"No key {key} in file {f.filename}") + print_key_tree(f) + return None + + ds.n_threads = 8 + print(f"Reading roi {roi} of volume {key} from {f.filename}") + vol = ds[roi] + print(f"Read volume with shape {vol.shape}, data type {vol.dtype}") + + return vol + + +def get_attrs(f: z5py.File, key: str): + print(type(f)) + if isinstance(f, (str, PurePath)): + f = z5py.File(f, "a") + + try: + ds = f[key] + except KeyError: + print(f"No key {key} in file {f.filename}") + print_key_tree(f) + return None + + return ds.attrs + + +def write_volume(f, arr: np.array, key, chunks=(1, 512, 512), attrs=None): + shape = arr.shape + compression = "gzip" + dtype = arr.dtype + + print(type(f)) + if isinstance(f, (str, PurePath)): + f = z5py.File(f, "a") + + if key not in f.keys(): + print(f"Created dataset {key}") + ds = f.create_dataset( + key, shape=shape, compression=compression, chunks=chunks, dtype=dtype + ) + else: + print(f"Overwriting {key}") + ds = f[key] + + ds.n_threads = 8 + print(f"Writing array to {key}") + ds[:] = arr + + if attrs is not None: + print("Assigning attributes of the dataset") + for key in attrs.keys(): + ds.attrs[key] = attrs[key] diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py old mode 100644 new mode 100755 index 6b99226..3c9066e --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -1,11 +1,12 @@ -from pathlib import Path import numpy as np import matplotlib.pyplot as plt import logging -from scipy.ndimage import rotate -from skimage.filters import gaussian -# from matchmaker.preprocessing import convert_to_point_cloud +import z5py +import os from matchmaker.data import create_point_cloud +from matchmaker.transform_utils import get_rotated_shape, get_transformation_matrix, rotate_img +from matchmaker.vis import plot_three_slices +import click def get_SVD_transform(img, save_path=None, percentile_trsh=90): @@ -69,7 +70,10 @@ def orient_head(img, save_path=None): plt.plot(int_profile) plt.xlabel("Coordinate") plt.ylabel("Sum intensity along Y axis") - plt.savefig(save_path, dpi=300) + if save_path is not None: + plt.savefig(save_path, dpi=300) + else: + plt.show() max_pos = int_profile.argmax() logging.info(f"Max position is {max_pos}, dimension shape is {img.shape[1]}") @@ -81,80 +85,60 @@ def orient_head(img, save_path=None): return True -def orient_sample(img, dapi_chan, save_path, dorsal=False): - """Preliminary orientation of the samples with body axis along Y, head closer to 0. +def prealign_sample(img, file_name, save_path): - Args: - img: input volume - dapi_chan: channel to use for registration - dorsal: if True, rotate around Y to align with ventral in Z direction + # align segmentation with PCs + gc, Vt = get_SVD_transform(img) + T, new_shape = get_transformation_matrix(img, gc, Vt, save_path=f"{save_path}/{file_name}_T_prealignment.txt") + img_rotated = rotate_img(img, T, output_shape=new_shape) - Returns: - oriented image - """ + # Check if head is oriented correctly, else rotate 180 degrees + rotate_head = orient_head( + img_rotated, f"{save_path}/plots/{file_name}_intensity_profile.png" + ) + if rotate_head: + print("Rotate head 180 degrees ...") + # img_rotated = np.rot90(img_rotated, k=2) # NOTE: use rotation matrix instead + R_3x3 = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 1]]) + img_center = 0.5 * (np.array(img_rotated.shape)-1) + offset = img_center - R_3x3 @ img_center + R = np.eye(4) + R[:3, :3] = R_3x3 + R[:3, 3] = offset + img_rotated = rotate_img(img_rotated, R, output_shape=img_rotated.shape) + + # update transformation matrix + T = T @ R + np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) + + plot_three_slices( + img_rotated, + save_path=f"{save_path}/plots/{file_name}_prealigned.png" + ) - # Rotate around Y if the volume is dorsal - save_path = Path(save_path) + return img_rotated - if dorsal: - logging.info("Rotated dorsal sample to align Z") - img = img[:, ::-1, :, ::-1] - # Rotate in XY plane, because samples can be oriented randomly, not only along X or Y - max_proj = np.max(img, axis=1)[dapi_chan, ...] - plt.figure() - plt.imshow(max_proj, cmap="Reds") - plt.savefig(save_path / "max_proj_input.png") - - logging.info(f"Smooth image with sigma={2}") - img_smoothed = gaussian(img[dapi_chan, ::10, ::10, ::10], sigma=3) - gc, Vt = get_SVD_transform( - img_smoothed, - save_path / "max_proj_point_cloud_random_angle.png", - percentile_trsh=90, - ) - rot_angle = 90 - np.degrees(np.arctan2(Vt[0, 1], Vt[0, 0])) - logging.info(f"Rotation angle to correct for random angle is {rot_angle}") - img = rotate(img, rot_angle, axes=[-2, -1], order=3, mode="constant") - logging.info(f"Rotated the input volume around Z by {rot_angle} degrees") +@click.command() +@click.option("-i", "--input-path", required=True, help="Input file") +@click.option("-k", "--input-key", required=True, help="Input key") +@click.option("-o", "--output-dir", required=True, help="Output directory") +def main(input_path, input_key, output_dir): + with z5py.File(input_path, "r") as f: + img = f[input_key][:] + file_name = os.path.splitext(os.path.basename(input_path))[0] - plt.figure() - plt.imshow(np.max(img, axis=1)[dapi_chan, ...], alpha=0.5, cmap="Blues") - plt.savefig(save_path / "max_proj_rotated_random_angle.png", dpi=300) - - # Check again if the sample is along X or along Y - # Make max projection and determine the direction of principal axes - max_proj = np.max(img, axis=1)[dapi_chan, ...] - img_smoothed = gaussian(img[dapi_chan, ::10, ::10, ::10], sigma=2) - gc, Vt = get_SVD_transform(img_smoothed, save_path / "max_proj_point_cloud.png") - - rot_angle = np.arccos(Vt[0, 0]) - logging.info(f"Rotation angle is {rot_angle}") - if np.abs(rot_angle - np.pi / 2) < np.pi / 4: - logging.info("Rotate 90 degrees") - plt.figure() - plt.imshow(max_proj, cmap="Reds") - img_rot = np.rot90(img, axes=(2, 3)) - - elif np.abs(rot_angle) < np.pi / 4: - logging.info("Rotation is already correct") - img_rot = img - else: - logging.info( - f"90 degree rotation can't be determined with first principal axis angle of {rot_angle * 180/ np.pi}" - ) - img_rot = img + if not os.path.exists(f"{output_dir}/plots"): + os.makedirs(f"{output_dir}/plots") - # Check if head is oriented correctly, else rotate 180 degrees + img_prealigned = prealign_sample(img, file_name, save_path=output_dir) - if orient_head(img_rot[dapi_chan, ...], save_path / "sum_intensity_profile.png"): - img_reg = np.rot90(img_rot, k=2, axes=(2, 3)) - else: - img_reg = img_rot + with z5py.File(f"{output_dir}/{file_name}_prealigned.n5", "w") as f: + f.create_dataset(input_key, data=img_prealigned, compression="gzip") + + # create MoBIE project + # create_mobie_project() - logging.info(f"Final image shape is {img_reg.shape}") - plt.figure() - plt.imshow(np.max(img_reg, axis=1)[dapi_chan, ...], alpha=0.5, cmap="Blues") - plt.savefig(save_path / "max_proj_rotated_final.png", dpi=300) - return img_reg +if __name__ == "__main__": + main() diff --git a/matchmaker/transform_utils.py b/matchmaker/transform_utils.py index b07b750..794ca1c 100644 --- a/matchmaker/transform_utils.py +++ b/matchmaker/transform_utils.py @@ -102,12 +102,12 @@ def get_transformation_matrix(img, gc, Vt, save_path=None): new_shape_center = np.array(new_shape) // 2 center_to_new_shape = get_translation_matrix(-new_shape_center) # 5. combine all transforms: get transformation matrix - R = center_to_origin @ rot @ center_to_new_shape + T = center_to_origin @ rot @ center_to_new_shape if save_path is not None: - np.savetxt(save_path, R) + np.savetxt(save_path, T) - return R, new_shape + return T, new_shape def rotate_img(img, rotation_matrix, output_shape=None, offset=None): From dc7a506d14e78c2c8e464a828ab011e0bbb53e86 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Thu, 5 Jun 2025 13:50:28 +0200 Subject: [PATCH 14/23] started to export results to mobie --- examples/registration_test.py | 2 +- matchmaker/mobie_export.py | 3 +-- matchmaker/prealignment.py | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/examples/registration_test.py b/examples/registration_test.py index 2f634e0..c3334a1 100644 --- a/examples/registration_test.py +++ b/examples/registration_test.py @@ -3,7 +3,7 @@ import numpy as np from matchmaker.prealignment import prealign_sample -from matchmaker.transform_utils import get_transformation_matrix, rotate_img +from matchmaker.transform_utils import rotate_img from matchmaker.vis import plot_three_slices, plot_overlay diff --git a/matchmaker/mobie_export.py b/matchmaker/mobie_export.py index e6cd503..21aaacb 100644 --- a/matchmaker/mobie_export.py +++ b/matchmaker/mobie_export.py @@ -4,7 +4,6 @@ # menu 3: moving_segmentation; then original, pre-aligned, ... import os import mobie -from matchmaker.n5_utils import get_attrs # from mobie import add_segmentation @@ -57,7 +56,7 @@ def create_mobie_project(input_path, input_key, output_dir): root=mobie_folder, dataset_name=dataset_name, segmentation_name="original", - resolution=[1,1,1], + resolution=[1, 1, 1], scale_factors=scale_factors, chunks=chunks, menu_name="fixed", diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index 3c9066e..f0107ce 100755 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -4,7 +4,7 @@ import z5py import os from matchmaker.data import create_point_cloud -from matchmaker.transform_utils import get_rotated_shape, get_transformation_matrix, rotate_img +from matchmaker.transform_utils import get_transformation_matrix, rotate_img from matchmaker.vis import plot_three_slices import click From 62c93ef24a2f8d9f1739d441d7b8b766f1e0617f Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Fri, 1 Aug 2025 14:41:20 +0200 Subject: [PATCH 15/23] set up prealignment + mobie export --- README.md | 15 ++--- examples/registration_test.py | 4 +- examples/reverse_deformation.py | 3 + matchmaker/apply.py | 1 + matchmaker/compute_registration.py | 103 +++++++++++++++++++++++++++++ matchmaker/mobie_export.py | 13 ++-- matchmaker/prealignment.py | 52 ++++++++++++--- 7 files changed, 167 insertions(+), 24 deletions(-) create mode 100644 matchmaker/compute_registration.py diff --git a/README.md b/README.md index 2083036..700c282 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# matchmaker +# 💞 Matchmaker Tool for segmentation-based deformable registration and object matching @@ -7,10 +7,8 @@ Tool for segmentation-based deformable registration and object matching ### Input - **Fixed image**: 3D instance segmentation in `n5` + resolution - - **Moving image**: 3D instance segmentation `n5` + resolution -- **Registration parameters yaml**: parameters of the `matchmaker` pipeline Expected image shape: ZYX @@ -23,11 +21,12 @@ Expected image shape: ZYX ### Registration steps -0. Create point cloud -1. PCA pre-alignment ---> **rigid transform matrix** -2. Manual input: should the image be flipped? -3. Rigid pre-alignment with Elastix OR with rigid CPD ---> **rigid transform matrix** -4. Coherent point drift +**1. PCA pre-alignment**: rigid alignment of fixed and moving image using the PCs \ + `prealignment.py --input_path ... --input_key ... --output_dir --mobie_export --type ...` \ + output: prealigned images + **rigid transformation matrix** +2. Rigid pre-alignment with Elastix OR with rigid CPD ---> **rigid transform matrix** + +**3. Coherent point drift** 5. Matching points with mixed integer programming 6. Deformable registration with Elastix with distance between keypoints in loss and rigidity penalty ---> **B-spline coefficients** diff --git a/examples/registration_test.py b/examples/registration_test.py index c3334a1..fc31f18 100644 --- a/examples/registration_test.py +++ b/examples/registration_test.py @@ -12,6 +12,7 @@ def main(): if not os.path.exists(f"{output_dir}/plots"): os.makedirs(f"{output_dir}/plots") + ####################### # prealign moving image moving_input = "./data/platy1_muscles_stardist_moving.n5" with z5py.File(moving_input, "r") as f: @@ -24,6 +25,7 @@ def main(): with z5py.File(f"{output_dir}/moving_prealigned.n5", "w") as f: f.create_dataset("seg", data=seg_moving_prealigned, compression="gzip") + ####################### # prealign fixed image fixed_input = "./data/platy1_muscles_stardist_fixed.n5" with z5py.File(fixed_input, "r") as f: @@ -31,12 +33,12 @@ def main(): plot_three_slices(seg_fixed, save_path=f"{output_dir}/plots/fixed.png") - # prealign moving image seg_fixed_prealigned = prealign_sample(seg_fixed, file_name="fixed", save_path=output_dir) with z5py.File(f"{output_dir}/fixed_prealigned.n5", "w") as f: f.create_dataset("seg", data=seg_fixed_prealigned, compression="gzip") + ####################### # test backtransform T = np.loadtxt(f"{output_dir}/fixed_T_prealignment.txt") seg_fixed_inv = rotate_img( diff --git a/examples/reverse_deformation.py b/examples/reverse_deformation.py index 5f163bb..9849da8 100644 --- a/examples/reverse_deformation.py +++ b/examples/reverse_deformation.py @@ -5,6 +5,9 @@ def main(): + ''' + Reverse the deformation of a sample by applying the inverse transformation matrix. + ''' with z5py.File("./data/platy1_muscles_stardist_moving.n5", "r") as f: seg_moving = f["seg"][:] diff --git a/matchmaker/apply.py b/matchmaker/apply.py index e69de29..f4a22ab 100644 --- a/matchmaker/apply.py +++ b/matchmaker/apply.py @@ -0,0 +1 @@ +# apply transformations \ No newline at end of file diff --git a/matchmaker/compute_registration.py b/matchmaker/compute_registration.py new file mode 100644 index 0000000..5614d5c --- /dev/null +++ b/matchmaker/compute_registration.py @@ -0,0 +1,103 @@ +import os +import z5py +import numpy as np + +from matchmaker.prealignment import prealignment_per_image +from matchmaker.mobie_export import export_to_mobie +from matchmaker.transform_utils import rotate_img +from matchmaker.vis import plot_three_slices, plot_overlay + + +def create_mobie_project( + fixed_input_path, + fixed_input_key, + moving_input_path, + moving_input_key, + output_dir, + ): + + # export fixed image to MoBIE + fixed_file_name = os.path.splitext(os.path.basename(fixed_input_path))[0] + export_to_mobie( + fixed_input_path, + fixed_input_key, + output_dir, + segmentation_name=f"{fixed_file_name}_original", + menu_name="fixed", + ) + + # export moving image to MoBIE + fixed_file_name = os.path.splitext(os.path.basename(moving_input_path))[0] + export_to_mobie( + moving_input_path, + moving_input_key, + output_dir, + segmentation_name=f"{fixed_file_name}_original", + menu_name="moving", + ) + + +def run_prealignment( + fixed_input_path, + fixed_input_key, + moving_input_path, + moving_input_key, + output_dir, + mobie_export=True, + ): + + fixed_prealigned = prealignment_per_image( + fixed_input_path, + fixed_input_key, + output_dir, + mobie_export, + image_type="fixed" + ) + + moving_prealigned = prealignment_per_image( + moving_input_path, + moving_input_key, + output_dir, + mobie_export, + image_type="moving" + ) + + plot_overlay( + fixed_prealigned, + moving_prealigned, + save_path=f"{output_dir}/plots/overlay_prealignment.png", + ) + + return fixed_prealigned, moving_prealigned + + +def main(): + fixed_input_path = "../examples/data/platy1_muscles_stardist_fixed.n5" + fixed_input_key = "seg" + moving_input_path = "../examples/data/platy1_muscles_stardist_moving.n5" + moving_input_key = "seg" + + output_dir = "../examples/data/test" + if not os.path.exists(f"{output_dir}/plots"): + os.makedirs(f"{output_dir}/plots") + + create_mobie_project( + fixed_input_path, + fixed_input_key, + moving_input_path, + moving_input_key, + output_dir, + ) + + fixed_prealigned, moving_prealigned = run_prealignment( + fixed_input_path, + fixed_input_key, + moving_input_path, + moving_input_key, + output_dir, + mobie_export=True, + ) + + +if __name__ == "__main__": + main() diff --git a/matchmaker/mobie_export.py b/matchmaker/mobie_export.py index 21aaacb..c4b6a78 100644 --- a/matchmaker/mobie_export.py +++ b/matchmaker/mobie_export.py @@ -38,14 +38,13 @@ # add_segmentation_to_mobie(input_path, input_key, seg_name, scale) -def create_mobie_project(input_path, input_key, output_dir): +def export_to_mobie(input_path, input_key, output_dir, segmentation_name, menu_name): if not os.path.exists(f"{output_dir}/mobie_project"): os.makedirs(f"{output_dir}/mobie_project") # Set parameters for MOBIE mobie_folder = f"{output_dir}/mobie_project" dataset_name = "platy1_muscles_stardist" - breakpoint() # resolution = get_attrs(input_path, input_key)["resolution"] chunks = (64, 64, 64) scale_factors = 4 * [[2, 2, 2]] @@ -55,22 +54,24 @@ def create_mobie_project(input_path, input_key, output_dir): input_key=input_key, root=mobie_folder, dataset_name=dataset_name, - segmentation_name="original", + segmentation_name=segmentation_name, resolution=[1, 1, 1], scale_factors=scale_factors, chunks=chunks, - menu_name="fixed", + menu_name=menu_name, file_format="ome.zarr", is_default_dataset=True ) def main(): - input_path = "/Users/marei/git-repositories/matchmaker/examples/data/platy1_muscles_stardist_fixed.n5" + input_path = "/Users/marei/git-repositories/matchmaker/examples/CLI_test/platy1_muscles_stardist_fixed_prealigned.n5" input_key = "seg" output_dir = "/Users/marei/git-repositories/matchmaker/examples/data/test" - create_mobie_project(input_path, input_key, output_dir) + file_name = os.path.splitext(os.path.basename(input_path))[0] + + export_to_mobie(input_path, input_key, output_dir, segmentation_name=f"{file_name}_prealigned", menu_name="fixed") print(f"MoBIE project created at {output_dir}/mobie_project") diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index f0107ce..6822125 100755 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -4,6 +4,7 @@ import z5py import os from matchmaker.data import create_point_cloud +from matchmaker.mobie_export import export_to_mobie from matchmaker.transform_utils import get_transformation_matrix, rotate_img from matchmaker.vis import plot_three_slices import click @@ -86,8 +87,17 @@ def orient_head(img, save_path=None): def prealign_sample(img, file_name, save_path): + """ + Pre-align a sample segmentation with its principal components. + + Args: + img: Segmentation volume + file_name: Name of the sample + save_path: Folder to save results - # align segmentation with PCs + Returns: + Pre-aligned segmentation volume. + """ gc, Vt = get_SVD_transform(img) T, new_shape = get_transformation_matrix(img, gc, Vt, save_path=f"{save_path}/{file_name}_T_prealignment.txt") img_rotated = rotate_img(img, T, output_shape=new_shape) @@ -98,7 +108,6 @@ def prealign_sample(img, file_name, save_path): ) if rotate_head: print("Rotate head 180 degrees ...") - # img_rotated = np.rot90(img_rotated, k=2) # NOTE: use rotation matrix instead R_3x3 = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 1]]) img_center = 0.5 * (np.array(img_rotated.shape)-1) offset = img_center - R_3x3 @ img_center @@ -119,25 +128,50 @@ def prealign_sample(img, file_name, save_path): return img_rotated -@click.command() -@click.option("-i", "--input-path", required=True, help="Input file") -@click.option("-k", "--input-key", required=True, help="Input key") -@click.option("-o", "--output-dir", required=True, help="Output directory") -def main(input_path, input_key, output_dir): - with z5py.File(input_path, "r") as f: +def prealignment_per_image(input_path, input_key, output_dir, mobie_export, image_type): + # load input image + with z5py.File(input_path, 'r') as f: img = f[input_key][:] file_name = os.path.splitext(os.path.basename(input_path))[0] if not os.path.exists(f"{output_dir}/plots"): os.makedirs(f"{output_dir}/plots") + # plot three slices of input image + plot_three_slices(img, save_path=f"{output_dir}/plots/{image_type}.png") + + # pre-align image img_prealigned = prealign_sample(img, file_name, save_path=output_dir) + # save pre-aligned image with z5py.File(f"{output_dir}/{file_name}_prealigned.n5", "w") as f: f.create_dataset(input_key, data=img_prealigned, compression="gzip") + # export pre-aligned image to MoBIE # create MoBIE project - # create_mobie_project() + if mobie_export: + export_to_mobie( + input_path=f"{output_dir}/{file_name}_prealigned.n5", + input_key=input_key, + output_dir=output_dir, + segmentation_name=f"{file_name}_prealigned", + menu_name=image_type + ) + + return img_prealigned + + +@click.command() +@click.option("-i", "--input_path", required=True, help="Input file") +@click.option("-k", "--input_key", required=True, help="Input key") +@click.option("-o", "--output_dir", required=True, help="Output directory") +@click.option("-m", "--mobie_export", required=False, is_flag=True, help="MoBIE export") +@click.option("-t", "--image_type", required=False, default="moving", help="Image Type: fixed or moving") +def main(input_path, input_key, output_dir, mobie_export, image_type): + if mobie_export and (image_type is None): + raise click.UsageError("--image_type is required when --mobie_export is set.") + + img_prealigned = prealignment_per_image(input_path, input_key, output_dir, mobie_export, image_type) if __name__ == "__main__": From d9e13ff6c42b0108e7fd89e4fc545c41527f1dcd Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Mon, 4 Aug 2025 18:09:44 +0200 Subject: [PATCH 16/23] update .gitignore --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index d73e375..e8ca84b 100644 --- a/.gitignore +++ b/.gitignore @@ -180,4 +180,4 @@ cython_debug/ .DS_Store -./examples/data/test/* \ No newline at end of file +./examples/data/* \ No newline at end of file From 668817131ad03f56b01cd5426e3e7f8e80fea414 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Mon, 4 Aug 2025 18:10:40 +0200 Subject: [PATCH 17/23] also include rigid alignment computation and upload to mobie --- examples/clean_up.py | 35 ++++++++ examples/deform_test_data.py | 25 ++++-- examples/reverse_deformation.py | 19 ++-- examples/view_data.py | 11 +++ matchmaker/align_rigid_elastix.py | 117 +++++++++++++------------ matchmaker/compute_registration.py | 122 +++++++------------------- matchmaker/elastix_utils.py | 2 - matchmaker/file_test.py | 9 -- matchmaker/mobie_export.py | 29 +++++++ matchmaker/n5_utils.py | 19 ++-- matchmaker/prealignment.py | 134 +++++++++++++++++++++++------ 11 files changed, 316 insertions(+), 206 deletions(-) create mode 100644 examples/clean_up.py create mode 100644 examples/view_data.py delete mode 100644 matchmaker/file_test.py diff --git a/examples/clean_up.py b/examples/clean_up.py new file mode 100644 index 0000000..75ba3b6 --- /dev/null +++ b/examples/clean_up.py @@ -0,0 +1,35 @@ +import json + +# clean up dataset.json +input_path = "./data/test/mobie_project/platy1_muscles_stardist/dataset.json" + +# Define the key you want to delete +key_to_delete = "platy1_muscles_stardist_moving_prealigned_rigid_aligned" + +# Load the JSON file +with open(input_path, "r") as file: + dataset_dict = json.load(file) + + +# Recursively delete all entries with the specified key +def delete_key_recursively(obj, key): + if isinstance(obj, dict): + # Remove the key if it exists in the current dictionary + obj.pop(key, None) + # Recurse into the dictionary values + for value in obj.values(): + delete_key_recursively(value, key) + elif isinstance(obj, list): + # Recurse into each item in the list + for item in obj: + delete_key_recursively(item, key) + + +# Call the recursive function on the loaded JSON data +delete_key_recursively(dataset_dict, key_to_delete) + +# Save the modified JSON back to a file +with open(input_path, "w") as file: + json.dump(dataset_dict, file, indent=4) + +print(f"Entries with key '{key_to_delete}' have been deleted.") diff --git a/examples/deform_test_data.py b/examples/deform_test_data.py index 005c64f..262ea2c 100644 --- a/examples/deform_test_data.py +++ b/examples/deform_test_data.py @@ -1,6 +1,5 @@ import numpy as np import tifffile as tif -import z5py import transforms3d as tf3d from matchmaker.transform_utils import ( @@ -10,6 +9,7 @@ crop_to_bbox, ) from matchmaker.vis import plot_three_slices, plot_overlay +from matchmaker.n5_utils import write_volume def remove_instances(seg, prob=0.05): @@ -39,9 +39,15 @@ def main(): seg_fixed = crop_to_bbox(seg) print("Cropped shape", seg_fixed.shape) - # save downsampled image - with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "w") as f: - f.create_dataset("seg", data=seg_fixed) + # save downsampled, fixed image + attributes = {"resolution": [1, 1, 1]} + write_volume( + f="./data/platy1_muscles_stardist_fixed.n5", + arr=seg_fixed, + key="seg", + chunks=(128, 512, 512), + attrs=attributes, + ) # save also as tiff tif.imwrite("./data/platy1_muscles_stardist_fixed.tif", seg_fixed) @@ -61,8 +67,15 @@ def main(): probability = 0.05 seg_moving = remove_instances(seg_moving, prob=probability) - with z5py.File("./data/platy1_muscles_stardist_moving.n5", "w") as f: - f.create_dataset("seg", data=seg_moving) + # save moving image + attributes = {"resolution": [1, 1, 1]} + write_volume( + f="./data/platy1_muscles_stardist_moving.n5", + arr=seg_moving, + key="seg", + chunks=(128, 512, 512), + attrs=attributes, + ) # save also as tiff tif.imwrite("./data/platy1_muscles_stardist_moving.tif", seg_moving) diff --git a/examples/reverse_deformation.py b/examples/reverse_deformation.py index 9849da8..38dc87b 100644 --- a/examples/reverse_deformation.py +++ b/examples/reverse_deformation.py @@ -1,29 +1,28 @@ import numpy as np -import z5py from matchmaker.vis import plot_three_slices, plot_overlay from matchmaker.transform_utils import rotate_img +from matchmaker.n5_utils import read_volume def main(): ''' Reverse the deformation of a sample by applying the inverse transformation matrix. ''' - with z5py.File("./data/platy1_muscles_stardist_moving.n5", "r") as f: - seg_moving = f["seg"][:] + seg_moving = read_volume( + f="./data/platy1_muscles_stardist_moving.n5", + key="seg", + ) # compare with original image - with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: - seg_fixed = f["seg"][:] + seg_fixed = read_volume( + f="./data/platy1_muscles_stardist_fixed.n5", + key="seg", + ) T = np.loadtxt("./data/transformation_matrix.txt") seg_moving_reverse = rotate_img(seg_moving, np.linalg.inv(T), output_shape=seg_fixed.shape) plot_three_slices(seg_moving_reverse) - - # compare with original image - with z5py.File("./data/platy1_muscles_stardist_fixed.n5", "r") as f: - seg_fixed = f["seg"][:] - plot_three_slices(seg_fixed) plot_overlay(seg_fixed, seg_moving_reverse) diff --git a/examples/view_data.py b/examples/view_data.py new file mode 100644 index 0000000..2a77a6b --- /dev/null +++ b/examples/view_data.py @@ -0,0 +1,11 @@ +from matchmaker.n5_utils import read_volume +import napari + +seg_fixed = read_volume("./data/test/platy1_muscles_stardist_fixed_prealigned.n5", key="seg") +seg_moving = read_volume("./data/test/platy1_muscles_stardist_moving_prealigned.n5", key="seg") +# seg_moving = seg_moving[:, :, ::-1] + +v = napari.Viewer() +v.add_labels(seg_fixed) +v.add_labels(seg_moving) +napari.run() \ No newline at end of file diff --git a/matchmaker/align_rigid_elastix.py b/matchmaker/align_rigid_elastix.py index 2bacee8..33b54f1 100644 --- a/matchmaker/align_rigid_elastix.py +++ b/matchmaker/align_rigid_elastix.py @@ -6,11 +6,14 @@ import logging import sys from matchmaker import elastix_utils +from matchmaker.mobie_export import export_to_mobie import itk +import click +import os def elastix_segm_rigid_alignment( - fixed_img_np, fixed_resolution, moving_img_np, moving_resolution, log_dir + fixed_img_np, fixed_resolution, moving_img_np, moving_resolution, output_dir ): """ Run rigid alignment of the ventral and dorsal datasets using elastix. @@ -19,7 +22,6 @@ def elastix_segm_rigid_alignment( logging.info("Do rigid transform of unnormalized images") logging.info(f"Fixed resolution {fixed_resolution}") logging.info(f"Moving resolution {moving_resolution}") - logging.info("Do rigid transform of unnormalized images") fixed_img_semantic_np = (fixed_img_np > 0).astype(np.float32) moving_img_semantic_np = (moving_img_np > 0).astype(np.float32) @@ -32,36 +34,33 @@ def elastix_segm_rigid_alignment( logging.info("Moving image") logging.info(f"{moving_img}") - # plot_overlay(fixed_img_np, moving_img_np, log_dir / "intersample_segm_overlay_before_alignment.png") - plot_overlay( elastix_utils.itk_to_np_order(itk.GetArrayFromImage(fixed_img)), elastix_utils.itk_to_np_order(itk.GetArrayFromImage(moving_img)), - log_dir / "intersample_segm_overlay_before_alignment.png", + f"{output_dir}/plots/intersample_segm_overlay_before_alignment.png", ) parameter_map_paths = [ - "pipeline_steps/inter_sample_registration/ParameterMap_segm_rigid_registration_corr.txt" + "../ParameterMap_segm_rigid_registration_corr.txt" ] logging.info("Run rigid registration") result_image, result_transform_parameters = elastix_utils.run_registration( fixed_img, moving_img, parameter_map_paths, - str(log_dir), + output_dir, log_name="elastix_log_rigid.log", set_threads=True, ) - # serialize_parameter_object(result_transform_parameters, "ParameterMap_rigid_transform", log_dir) logging.info(f"Result image shape {result_image.shape}") result_img_np = elastix_utils.itk_to_np_order(itk.GetArrayFromImage(result_image)) plot_overlay( elastix_utils.itk_to_np_order(itk.GetArrayFromImage(fixed_img)), result_img_np, - log_dir / "intersample_segm_rigid_alignment_semantic.png", + f"{output_dir}/plots/intersample_segm_rigid_alignment_semantic.png", ) - + # NOTE: difference between result_img_np before and after applying transform? logging.info("Apply transform to all channels") result_img_np = elastix_utils.apply_transform_chanwise( result_transform_parameters, moving_img_np, moving_resolution @@ -71,84 +70,90 @@ def elastix_segm_rigid_alignment( return result_img_np -def main(): - parser = argparse.ArgumentParser( - description="""Align rigidly moving segmentation volume to fixed volume.""" - ) - - parser.add_argument("fixed_n5", type=str, help="Path of the output n5") - parser.add_argument("fixed_key", type=str, help="Key of ventral dataset") - parser.add_argument("moving_n5", type=str, help="Path of the output n5") - parser.add_argument("moving_key", type=str, help="Key of ventral dataset") - parser.add_argument("output_n5", type=str, help="Path of the output n5") - parser.add_argument("output_key", type=str, help="Key to write aligned ventral dataset to in the output n5") - parser.add_argument("log_dir", type=str, help="Directory to store diagnostic plots etc") - parser.add_argument("log_path", type=str, help="Path to log file") - args = parser.parse_args() - - log_dir = Path(args.log_dir) - log_dir.mkdir(exist_ok=True) +def run_rigid_alignment( + fixed_path, + fixed_key, + moving_path, + moving_key, + output_dir, + mobie_export +): + if not os.path.exists(f"{output_dir}/plots"): + os.makedirs(f"{output_dir}/plots") logging.basicConfig( level=logging.INFO, format="%(asctime)s [%(levelname)s] %(message)s", handlers=[ - logging.FileHandler(args.log_path, mode="w"), + logging.FileHandler(f"{output_dir}/rigid_alignment.log", mode="w"), logging.StreamHandler(sys.stdout), ], datefmt="%Y-%m-%d %H:%M:%S", ) logging.info("Read image file") - fixed_img_np = read_volume(args.fixed_n5, args.fixed_key) - - # pad_z = 50 - # fixed_img_np = np.pad(fixed_img_np, pad_width=((pad_z, pad_z), (0, 0), (0, 0))) - - # fixed_mask = gaussian(fixed_img_np, sigma=[6, 12, 12]) - # fixed_mask_bin = fixed_mask > np.quantile(fixed_mask, 0.7) - # fixed_img_np = fixed_img_np * fixed_mask_bin + fixed_img_np = read_volume(fixed_path, fixed_key) + print(fixed_img_np.dtype) - if args.fixed_n5 == args.moving_n5: + if fixed_path == moving_path: logging.info("Same n5 for fixed and moving image, not doing registration") moving_img_np = fixed_img_np else: - moving_img_np = read_volume(args.moving_n5, args.moving_key) - # moving_mask = gaussian(moving_img_np, sigma=[6, 12, 12]) - # moving_mask_bin = moving_mask > np.quantile(moving_mask, 0.7) - # moving_img_np = moving_img_np * moving_mask_bin + moving_img_np = read_volume(moving_path, moving_key) fixed_img_np = fixed_img_np.astype(np.float32) moving_img_np = moving_img_np.astype(np.float32) logging.info("Start registration") - # Dirty hack for the sample size mismatch - # For no apparent reason the size of the light samples seems to be ~20% larger than the EM - moving_resolution = get_attrs(args.moving_n5, args.moving_key)["resolution"] - moving_resolution = [r * 0.9 for r in moving_resolution] - moving_img_np = elastix_segm_rigid_alignment( - fixed_img_np, - get_attrs(args.fixed_n5, args.fixed_key)["resolution"], - moving_img_np, - moving_resolution, - log_dir, + fixed_img_np=fixed_img_np, + fixed_resolution=get_attrs(fixed_path, fixed_key)["resolution"], + moving_img_np=moving_img_np, + moving_resolution=get_attrs(moving_path, moving_key)["resolution"], + output_dir=output_dir, ) - + moving_img_np = moving_img_np.astype(np.uint16) logging.info("Write results") - attributes = dict(get_attrs(args.moving_n5, args.moving_key)) - attributes["resolution"] = get_attrs(args.fixed_n5, args.fixed_key)["resolution"] + attributes = dict(get_attrs(moving_path, moving_key)) + + file_name = os.path.splitext(os.path.basename(moving_path))[0] + file_name = file_name.removesuffix("_prealigned") + print(file_name) write_volume( - args.output_n5, + f"{output_dir}/{file_name}_rigid_aligned.n5", moving_img_np, - args.output_key, + moving_key, chunks=(128, 512, 512), attrs=attributes, ) + # export rigid alignment to mobie + if mobie_export: + export_to_mobie( + input_path=f"{output_dir}/{file_name}_rigid_aligned.n5", + input_key=moving_key, + output_dir=output_dir, + segmentation_name=f"{file_name}_rigid_aligned", + menu_name="moving" + ) + + +@click.command() +@click.option("-fi", "--fixed_path", required=True, help="Fixed input .n5 file") +@click.option("-fk", "--fixed_key", required=True, help="Fixed input key") +@click.option("-mi", "--moving_path", required=True, help="Moving input .n5 file") +@click.option("-mk", "--moving_key", required=True, help="Moving input key") +@click.option("-o", "--output_dir", required=True, help="Output directory") +@click.option("-m", "--mobie_export", required=False, is_flag=True, help="MoBIE export") +def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export): + + run_rigid_alignment(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export) + if __name__ == "__main__": main() + + # python align_rigid_elastix.py ../examples/data/test/platy1_muscles_stardist_fixed_prealigned.n5 seg ../examples/data/test/platy1_muscles_stardist_moving_prealigned.n5 seg ../examples/data/test/platy1_muscles_stardist_moving_rigid_aligned.n5 seg diff --git a/matchmaker/compute_registration.py b/matchmaker/compute_registration.py index 5614d5c..373975b 100644 --- a/matchmaker/compute_registration.py +++ b/matchmaker/compute_registration.py @@ -1,101 +1,43 @@ import os -import z5py import numpy as np - -from matchmaker.prealignment import prealignment_per_image -from matchmaker.mobie_export import export_to_mobie -from matchmaker.transform_utils import rotate_img -from matchmaker.vis import plot_three_slices, plot_overlay - - -def create_mobie_project( - fixed_input_path, - fixed_input_key, - moving_input_path, - moving_input_key, - output_dir, - ): - - # export fixed image to MoBIE - fixed_file_name = os.path.splitext(os.path.basename(fixed_input_path))[0] - export_to_mobie( - fixed_input_path, - fixed_input_key, - output_dir, - segmentation_name=f"{fixed_file_name}_original", - menu_name="fixed", - ) - - # export moving image to MoBIE - fixed_file_name = os.path.splitext(os.path.basename(moving_input_path))[0] - export_to_mobie( - moving_input_path, - moving_input_key, - output_dir, - segmentation_name=f"{fixed_file_name}_original", - menu_name="moving", - ) - - -def run_prealignment( - fixed_input_path, - fixed_input_key, - moving_input_path, - moving_input_key, - output_dir, - mobie_export=True, - ): - - fixed_prealigned = prealignment_per_image( - fixed_input_path, - fixed_input_key, - output_dir, - mobie_export, - image_type="fixed" - ) - - moving_prealigned = prealignment_per_image( - moving_input_path, - moving_input_key, - output_dir, - mobie_export, - image_type="moving" - ) - - plot_overlay( - fixed_prealigned, - moving_prealigned, - save_path=f"{output_dir}/plots/overlay_prealignment.png", - ) - - return fixed_prealigned, moving_prealigned - - -def main(): - fixed_input_path = "../examples/data/platy1_muscles_stardist_fixed.n5" - fixed_input_key = "seg" - moving_input_path = "../examples/data/platy1_muscles_stardist_moving.n5" - moving_input_key = "seg" - - output_dir = "../examples/data/test" +import click +from matchmaker.prealignment import run_prealignment +from matchmaker.mobie_export import create_mobie_project, export_to_mobie + + +@click.command() +@click.option("-fi", "--fixed_path", required=True, help="Fixed input .n5 file") +@click.option("-fk", "--fixed_key", required=True, help="Fixed input key") +@click.option("-mi", "--moving_path", required=True, help="Moving input .n5 file") +@click.option("-mk", "--moving_key", required=True, help="Moving input key") +@click.option("-o", "--output_dir", required=True, help="Output directory") +@click.option("-m", "--mobie_export", required=False, is_flag=True, help="MoBIE export") +def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export): + # fixed_input_path = "../examples/data/platy1_muscles_stardist_fixed.n5" + # fixed_input_key = "seg" + # moving_input_path = "../examples/data/platy1_muscles_stardist_moving.n5" + # moving_input_key = "seg" + + # output_dir = "../examples/data/test" if not os.path.exists(f"{output_dir}/plots"): os.makedirs(f"{output_dir}/plots") - create_mobie_project( - fixed_input_path, - fixed_input_key, - moving_input_path, - moving_input_key, - output_dir, - ) + if mobie_export: + create_mobie_project( + fixed_path, + fixed_key, + moving_path, + moving_key, + output_dir, + ) fixed_prealigned, moving_prealigned = run_prealignment( - fixed_input_path, - fixed_input_key, - moving_input_path, - moving_input_key, + fixed_path, + fixed_key, + moving_path, + moving_key, output_dir, - mobie_export=True, + mobie_export, ) diff --git a/matchmaker/elastix_utils.py b/matchmaker/elastix_utils.py index 60fab2d..738b67f 100644 --- a/matchmaker/elastix_utils.py +++ b/matchmaker/elastix_utils.py @@ -81,8 +81,6 @@ def run_registration( elastix_object = itk.ElastixRegistrationMethod.New(fixed_img, moving_img) logging.info(elastix_object) logging.info("Created registration object") - # elastix_object.SetFixedImage(fixed_image) - # elastix_object.SetMovingImage(moving_image) elastix_object.SetParameterObject(parameter_object) if set_threads: elastix_object.SetNumberOfThreads(32) diff --git a/matchmaker/file_test.py b/matchmaker/file_test.py deleted file mode 100644 index 8ac0f47..0000000 --- a/matchmaker/file_test.py +++ /dev/null @@ -1,9 +0,0 @@ -import z5py -import napari - -with z5py.File("/Users/marei/git-repositories/matchmaker/examples/data/test/mobie_project/platy1_muscles_stardist/images/ome-zarr/original.ome.zarr", "r") as f: - seg = f["s0"][:] - -v = napari.Viewer() -v.add_image(seg) -napari.run() \ No newline at end of file diff --git a/matchmaker/mobie_export.py b/matchmaker/mobie_export.py index c4b6a78..735f67d 100644 --- a/matchmaker/mobie_export.py +++ b/matchmaker/mobie_export.py @@ -38,6 +38,35 @@ # add_segmentation_to_mobie(input_path, input_key, seg_name, scale) +def create_mobie_project( + fixed_input_path, + fixed_input_key, + moving_input_path, + moving_input_key, + output_dir, +): + + # export fixed image to MoBIE + fixed_file_name = os.path.splitext(os.path.basename(fixed_input_path))[0] + export_to_mobie( + fixed_input_path, + fixed_input_key, + output_dir, + segmentation_name=f"{fixed_file_name}_original", + menu_name="fixed", + ) + + # export moving image to MoBIE + moving_file_name = os.path.splitext(os.path.basename(moving_input_path))[0] + export_to_mobie( + moving_input_path, + moving_input_key, + output_dir, + segmentation_name=f"{moving_file_name}_original", + menu_name="moving", + ) + + def export_to_mobie(input_path, input_key, output_dir, segmentation_name, menu_name): if not os.path.exists(f"{output_dir}/mobie_project"): os.makedirs(f"{output_dir}/mobie_project") diff --git a/matchmaker/n5_utils.py b/matchmaker/n5_utils.py index a899bc2..a822b2d 100644 --- a/matchmaker/n5_utils.py +++ b/matchmaker/n5_utils.py @@ -11,7 +11,7 @@ def print_key_tree(f: z5py.File): def read_volume( f: z5py.File, key: str, roi: np.lib.index_tricks.IndexExpression = np.s_[:] ): - print(type(f)) + # print(type(f)) if isinstance(f, (str, PurePath)): f = z5py.File(f, "r") @@ -23,15 +23,15 @@ def read_volume( return None ds.n_threads = 8 - print(f"Reading roi {roi} of volume {key} from {f.filename}") + # print(f"Reading roi {roi} of volume {key} from {f.filename}") vol = ds[roi] - print(f"Read volume with shape {vol.shape}, data type {vol.dtype}") + # print(f"Read volume with shape {vol.shape}, data type {vol.dtype}") return vol def get_attrs(f: z5py.File, key: str): - print(type(f)) + # print(type(f)) if isinstance(f, (str, PurePath)): f = z5py.File(f, "a") @@ -50,24 +50,25 @@ def write_volume(f, arr: np.array, key, chunks=(1, 512, 512), attrs=None): compression = "gzip" dtype = arr.dtype - print(type(f)) if isinstance(f, (str, PurePath)): f = z5py.File(f, "a") if key not in f.keys(): - print(f"Created dataset {key}") + # print(f"Created dataset {key}") ds = f.create_dataset( key, shape=shape, compression=compression, chunks=chunks, dtype=dtype ) else: - print(f"Overwriting {key}") + # print(f"Overwriting {key}") ds = f[key] ds.n_threads = 8 - print(f"Writing array to {key}") + # print(f"Writing array to {key}") ds[:] = arr + print(f"Dataset {key} written to {f.filename}") + if attrs is not None: - print("Assigning attributes of the dataset") + # print("Assigning attributes of the dataset") for key in attrs.keys(): ds.attrs[key] = attrs[key] diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index 6822125..128d1f0 100755 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -1,22 +1,22 @@ import numpy as np import matplotlib.pyplot as plt import logging -import z5py import os +import click + from matchmaker.data import create_point_cloud from matchmaker.mobie_export import export_to_mobie from matchmaker.transform_utils import get_transformation_matrix, rotate_img -from matchmaker.vis import plot_three_slices -import click +from matchmaker.n5_utils import read_volume, get_attrs, write_volume +from matchmaker.vis import plot_three_slices, plot_overlay -def get_SVD_transform(img, save_path=None, percentile_trsh=90): +def get_SVD_transform(img, save_path=None): """Convert image to point cloud by thresholding, then run SVD on resulting point cloud. Args: img: _description_ plot_path: _description_. Defaults to None. - percentile_trsh: _description_. Defaults to 90. Returns: Variance matrix and principal axes matrix. @@ -54,7 +54,7 @@ def get_SVD_transform(img, save_path=None, percentile_trsh=90): return gc, Vt -def orient_head(img, save_path=None): +def orient_yaxis(img, save_path=None): """Euristic to orient all samples "head up": calculate sum intensity profile along the Y axis, if max is closer to 0 then do nothing, else rotate 180 degrees. @@ -65,11 +65,11 @@ def orient_head(img, save_path=None): True if rotation is needed. """ assert img.ndim == 3, f"Input image should have 3 dimensions, has {img.ndim}" - int_profile = np.sum(img, axis=(0, 2)) + int_profile = np.sum(img, axis=(0, 2)) # profile along y (2nd dimension) plt.figure() plt.plot(int_profile) - plt.xlabel("Coordinate") + plt.xlabel("Y Coordinate") plt.ylabel("Sum intensity along Y axis") if save_path is not None: plt.savefig(save_path, dpi=300) @@ -86,6 +86,37 @@ def orient_head(img, save_path=None): return True +def orient_xaxis(img, save_path=None): + """Heuristic to orient all samples 'head left': calculate sum intensity profile along the X axis. + If max is closer to 0, then do nothing, else flip along x (mirror). + + Returns: + True if mirroring is needed. + """ + assert img.ndim == 3, f"Input image should have 3 dimensions, has {img.ndim}" + int_profile = np.sum(img, axis=(0, 1)) # profile along x (3rd dimension) + + plt.figure() + plt.plot(int_profile) + plt.xlabel("X coordinate") + plt.ylabel("Sum intensity along X axis") + if save_path is not None: + plt.savefig(save_path, dpi=300) + else: + plt.show() + + max_pos = int_profile.argmax() + logging.info(f"[orient_left] Max position is {max_pos}, shape[2] = {img.shape[2]}") + if max_pos < img.shape[2] // 2: + logging.info("Correct left-right orientation") + return False + else: + logging.info("Mirror along x-axis to align left-right") + return True + +# NOTE: also needed for z axis? + + def prealign_sample(img, file_name, save_path): """ Pre-align a sample segmentation with its principal components. @@ -103,11 +134,11 @@ def prealign_sample(img, file_name, save_path): img_rotated = rotate_img(img, T, output_shape=new_shape) # Check if head is oriented correctly, else rotate 180 degrees - rotate_head = orient_head( - img_rotated, f"{save_path}/plots/{file_name}_intensity_profile.png" + rotate_y = orient_yaxis( + img_rotated, f"{save_path}/plots/{file_name}_intensity_profile_y.png" ) - if rotate_head: - print("Rotate head 180 degrees ...") + if rotate_y: + print("Rotate head 180 degrees/ mirror on y axis ...") R_3x3 = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 1]]) img_center = 0.5 * (np.array(img_rotated.shape)-1) offset = img_center - R_3x3 @ img_center @@ -120,6 +151,24 @@ def prealign_sample(img, file_name, save_path): T = T @ R np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) + # Check if left-right is oriented correctly, else mirror on x axis + rotate_x = orient_xaxis( + img_rotated, f"{save_path}/plots/{file_name}_intensity_profile_x.png" + ) + if rotate_x: + print("Mirror on x axis ...") + R_3x3 = np.array([[1, 0, 0], [0, 1, 0], [0, 0, -1]]) + img_center = 0.5 * (np.array(img_rotated.shape)-1) + offset = img_center - R_3x3 @ img_center + R = np.eye(4) + R[:3, :3] = R_3x3 + R[:3, 3] = offset + img_rotated = rotate_img(img_rotated, R, output_shape=img_rotated.shape) + + # update transformation matrix + T = T @ R + np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) + plot_three_slices( img_rotated, save_path=f"{save_path}/plots/{file_name}_prealigned.png" @@ -130,8 +179,7 @@ def prealign_sample(img, file_name, save_path): def prealignment_per_image(input_path, input_key, output_dir, mobie_export, image_type): # load input image - with z5py.File(input_path, 'r') as f: - img = f[input_key][:] + img = read_volume(input_path, input_key) file_name = os.path.splitext(os.path.basename(input_path))[0] if not os.path.exists(f"{output_dir}/plots"): @@ -144,11 +192,15 @@ def prealignment_per_image(input_path, input_key, output_dir, mobie_export, imag img_prealigned = prealign_sample(img, file_name, save_path=output_dir) # save pre-aligned image - with z5py.File(f"{output_dir}/{file_name}_prealigned.n5", "w") as f: - f.create_dataset(input_key, data=img_prealigned, compression="gzip") + attributes = dict(get_attrs(input_path, input_key)) + write_volume( + f=f"{output_dir}/{file_name}_prealigned.n5", + arr=img_prealigned, + key=input_key, + attrs=attributes + ) # export pre-aligned image to MoBIE - # create MoBIE project if mobie_export: export_to_mobie( input_path=f"{output_dir}/{file_name}_prealigned.n5", @@ -161,17 +213,51 @@ def prealignment_per_image(input_path, input_key, output_dir, mobie_export, imag return img_prealigned +def run_prealignment( + fixed_input_path, + fixed_input_key, + moving_input_path, + moving_input_key, + output_dir, + mobie_export=True, +): + + fixed_prealigned = prealignment_per_image( + fixed_input_path, + fixed_input_key, + output_dir, + mobie_export, + image_type="fixed" + ) + + moving_prealigned = prealignment_per_image( + moving_input_path, + moving_input_key, + output_dir, + mobie_export, + image_type="moving", + ) + print(fixed_prealigned.dtype) + + plot_overlay( + fixed_prealigned, + moving_prealigned, + save_path=f"{output_dir}/plots/overlay_prealignment.png", + ) + + # return fixed_prealigned, moving_prealigned + + @click.command() -@click.option("-i", "--input_path", required=True, help="Input file") -@click.option("-k", "--input_key", required=True, help="Input key") +@click.option("-fi", "--fixed_path", required=True, help="Fixed input .n5 file") +@click.option("-fk", "--fixed_key", required=True, help="Fixed input key") +@click.option("-mi", "--moving_path", required=True, help="Moving input .n5 file") +@click.option("-mk", "--moving_key", required=True, help="Moving input key") @click.option("-o", "--output_dir", required=True, help="Output directory") @click.option("-m", "--mobie_export", required=False, is_flag=True, help="MoBIE export") -@click.option("-t", "--image_type", required=False, default="moving", help="Image Type: fixed or moving") -def main(input_path, input_key, output_dir, mobie_export, image_type): - if mobie_export and (image_type is None): - raise click.UsageError("--image_type is required when --mobie_export is set.") +def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export): - img_prealigned = prealignment_per_image(input_path, input_key, output_dir, mobie_export, image_type) + run_prealignment(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export) if __name__ == "__main__": From 82134ecfd85a4f09a70d461365dc796a7fafb4fe Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Wed, 6 Aug 2025 10:53:17 +0200 Subject: [PATCH 18/23] prealingment before updating overall structure --- matchmaker/prealignment.py | 177 +++++++++++++++++++++---------------- 1 file changed, 99 insertions(+), 78 deletions(-) diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index 128d1f0..f0cea95 100755 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -3,6 +3,8 @@ import logging import os import click +import sys +import napari from matchmaker.data import create_point_cloud from matchmaker.mobie_export import export_to_mobie @@ -29,10 +31,10 @@ def get_SVD_transform(img, save_path=None): logging.info(f"Point cloud shape {pos.shape}") logging.info(f"Point cloud center {gc}") - logging.info("Run SVD") + logging.info("Run SVD ...") U, S, Vt = np.linalg.svd(pos_c, full_matrices=False) logging.info("U") - logging.info(U) + logging.info(str(U)) logging.info("S") logging.info(str(S)) logging.info("Vt") @@ -41,6 +43,7 @@ def get_SVD_transform(img, save_path=None): logging.info("Rotate point cloud") vr = pos_c @ Vt.T + # TODO: update axis labels plt.figure(figsize=(10, 5)) plt.subplot(1, 2, 1) plt.title("original vertices") @@ -54,67 +57,53 @@ def get_SVD_transform(img, save_path=None): return gc, Vt -def orient_yaxis(img, save_path=None): - """Euristic to orient all samples "head up": calculate sum intensity profile along the Y axis, - if max is closer to 0 then do nothing, else rotate 180 degrees. - - Args: - img: DAPI volume - - Returns: - True if rotation is needed. - """ +def orient_axis(img, axis, save_path=None): + logging.info(f"Orient along axis {axis}") assert img.ndim == 3, f"Input image should have 3 dimensions, has {img.ndim}" - int_profile = np.sum(img, axis=(0, 2)) # profile along y (2nd dimension) + if axis == 0: + int_profile = np.sum(img, axis=(1, 2)) # profile along z + R_3x3 = np.array([[1, 0, 0], [0, -1, 0], [0, 0, -1]]) + if axis == 1: + int_profile = np.sum(img, axis=(0, 2)) # profile along y + R_3x3 = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]]) + if axis == 2: + int_profile = np.sum(img, axis=(0, 1)) # profile along x + R_3x3 = np.array([[1, 0, 0], [0, -1, 0], [0, 0, 1]]) plt.figure() plt.plot(int_profile) plt.xlabel("Y Coordinate") - plt.ylabel("Sum intensity along Y axis") + plt.ylabel(f"Sum intensity along axis = {axis}") if save_path is not None: plt.savefig(save_path, dpi=300) else: plt.show() max_pos = int_profile.argmax() - logging.info(f"Max position is {max_pos}, dimension shape is {img.shape[1]}") + logging.info(f"Max position is {max_pos}, dimension shape is {img.shape[axis]}") if max_pos < img.shape[1] // 2: - logging.info("Correct head orientation") + logging.info("Correct orientation") return False else: - logging.info("Rotate 180 degree to align head position") + logging.info("Rotate 180 degree to align") return True -def orient_xaxis(img, save_path=None): - """Heuristic to orient all samples 'head left': calculate sum intensity profile along the X axis. - If max is closer to 0, then do nothing, else flip along x (mirror). - - Returns: - True if mirroring is needed. - """ - assert img.ndim == 3, f"Input image should have 3 dimensions, has {img.ndim}" - int_profile = np.sum(img, axis=(0, 1)) # profile along x (3rd dimension) +# def orient_sample(img, axis, save_path=None): +# R_3x3 = orient_axis(img, axis=axis) #, save_path=save_path) +# if R_3x3 is not None: +# img_center = 0.5 * (np.array(img_rotated.shape)-1) +# offset = img_center - R_3x3 @ img_center +# R = np.eye(4) +# R[:3, :3] = R_3x3 +# R[:3, 3] = offset +# img_rotated = rotate_img(img_rotated, R, output_shape=img_rotated.shape) - plt.figure() - plt.plot(int_profile) - plt.xlabel("X coordinate") - plt.ylabel("Sum intensity along X axis") - if save_path is not None: - plt.savefig(save_path, dpi=300) - else: - plt.show() - - max_pos = int_profile.argmax() - logging.info(f"[orient_left] Max position is {max_pos}, shape[2] = {img.shape[2]}") - if max_pos < img.shape[2] // 2: - logging.info("Correct left-right orientation") - return False - else: - logging.info("Mirror along x-axis to align left-right") - return True - -# NOTE: also needed for z axis? +# # update transformation matrix +# T = T @ R +# np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) + +# return img_rotated def prealign_sample(img, file_name, save_path): @@ -133,31 +122,10 @@ def prealign_sample(img, file_name, save_path): T, new_shape = get_transformation_matrix(img, gc, Vt, save_path=f"{save_path}/{file_name}_T_prealignment.txt") img_rotated = rotate_img(img, T, output_shape=new_shape) - # Check if head is oriented correctly, else rotate 180 degrees - rotate_y = orient_yaxis( - img_rotated, f"{save_path}/plots/{file_name}_intensity_profile_y.png" - ) - if rotate_y: - print("Rotate head 180 degrees/ mirror on y axis ...") - R_3x3 = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 1]]) - img_center = 0.5 * (np.array(img_rotated.shape)-1) - offset = img_center - R_3x3 @ img_center - R = np.eye(4) - R[:3, :3] = R_3x3 - R[:3, 3] = offset - img_rotated = rotate_img(img_rotated, R, output_shape=img_rotated.shape) - - # update transformation matrix - T = T @ R - np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) - - # Check if left-right is oriented correctly, else mirror on x axis - rotate_x = orient_xaxis( - img_rotated, f"{save_path}/plots/{file_name}_intensity_profile_x.png" - ) - if rotate_x: - print("Mirror on x axis ...") + if np.linalg.det(Vt.T) < 0: + print("WARNING: V includes a reflection (mirroring)") R_3x3 = np.array([[1, 0, 0], [0, 1, 0], [0, 0, -1]]) + logging.info("Mirror back...") img_center = 0.5 * (np.array(img_rotated.shape)-1) offset = img_center - R_3x3 @ img_center R = np.eye(4) @@ -169,6 +137,22 @@ def prealign_sample(img, file_name, save_path): T = T @ R np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) + # check orientation in every three dimensions/ if heads up + # for axis in [0, 1, 2]: + # R_3x3 = orient_axis(img_rotated, axis=axis, save_path=f"{save_path}/plots/{file_name}_intensity_profile_{axis}.png") + # if R_3x3 is not None: + # logging.info(f"Rotate axis {axis} 180 degrees to align...") + # img_center = 0.5 * (np.array(img_rotated.shape)-1) + # offset = img_center - R_3x3 @ img_center + # R = np.eye(4) + # R[:3, :3] = R_3x3 + # R[:3, 3] = offset + # img_rotated = rotate_img(img_rotated, R, output_shape=img_rotated.shape) + + # # update transformation matrix + # T = T @ R + # np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) + plot_three_slices( img_rotated, save_path=f"{save_path}/plots/{file_name}_prealigned.png" @@ -178,20 +162,17 @@ def prealign_sample(img, file_name, save_path): def prealignment_per_image(input_path, input_key, output_dir, mobie_export, image_type): - # load input image + logging.info("Read image file") img = read_volume(input_path, input_key) file_name = os.path.splitext(os.path.basename(input_path))[0] - if not os.path.exists(f"{output_dir}/plots"): - os.makedirs(f"{output_dir}/plots") - # plot three slices of input image plot_three_slices(img, save_path=f"{output_dir}/plots/{image_type}.png") - # pre-align image + logging.info("Prealign image") img_prealigned = prealign_sample(img, file_name, save_path=output_dir) - # save pre-aligned image + logging.info("Save prealigned image") attributes = dict(get_attrs(input_path, input_key)) write_volume( f=f"{output_dir}/{file_name}_prealigned.n5", @@ -202,6 +183,7 @@ def prealignment_per_image(input_path, input_key, output_dir, mobie_export, imag # export pre-aligned image to MoBIE if mobie_export: + logging.info("Export prealigned image to MoBIE") export_to_mobie( input_path=f"{output_dir}/{file_name}_prealigned.n5", input_key=input_key, @@ -221,7 +203,19 @@ def run_prealignment( output_dir, mobie_export=True, ): + if not os.path.exists(f"{output_dir}/plots"): + os.makedirs(f"{output_dir}/plots") + logging.basicConfig( + level=logging.INFO, + format="%(asctime)s [%(levelname)s] %(message)s", + handlers=[ + logging.FileHandler(f"{output_dir}/prealignment.log", mode="w"), + logging.StreamHandler(sys.stdout), + ], + datefmt="%Y-%m-%d %H:%M:%S", + ) + logging.info("Start prealignment of fixed image ...") fixed_prealigned = prealignment_per_image( fixed_input_path, fixed_input_key, @@ -230,6 +224,7 @@ def run_prealignment( image_type="fixed" ) + logging.info("Start prealignment of moving image ...") moving_prealigned = prealignment_per_image( moving_input_path, moving_input_key, @@ -237,7 +232,35 @@ def run_prealignment( mobie_export, image_type="moving", ) - print(fixed_prealigned.dtype) + # check orientation + for axis in range(3): + rotate_axis_fixed = orient_axis(moving_prealigned, axis=axis, save_path=f"{output_dir}/plots/moving_intensity_profile_{axis}.png") + rotate_axis_moving = orient_axis(fixed_prealigned, axis=axis, save_path=f"{output_dir}/plots/fixed_intensity_profile_{axis}.png") + + if not rotate_axis_fixed and not rotate_axis_moving: + print(f"Correct orientation in axis {axis}") + else: + print(f"Rotate axis {axis} 180 degrees to align...") + R_3x3 = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 1]]) + logging.info(f"Rotate axis {axis} 180 degrees to align...") + img_center = 0.5 * (np.array(moving_prealigned.shape)-1) + offset = img_center - R_3x3 @ img_center + R = np.eye(4) + R[:3, :3] = R_3x3 + R[:3, 3] = offset + moving_prealigned = rotate_img(moving_prealigned, R, output_shape=moving_prealigned.shape) + + # update transformation matrix + # T = T @ R + # np.savetxt(f"{output_dir}/moving_T_prealignment.txt", T) + + logging.info("Prealignment done.") + + v = napari.Viewer() + v.add_labels(fixed_prealigned) + v.add_labels(moving_prealigned) + napari.run() + breakpoint() plot_overlay( fixed_prealigned, @@ -245,8 +268,6 @@ def run_prealignment( save_path=f"{output_dir}/plots/overlay_prealignment.png", ) - # return fixed_prealigned, moving_prealigned - @click.command() @click.option("-fi", "--fixed_path", required=True, help="Fixed input .n5 file") From c5d8bd15bcddf4248cc39b97500b2259851e3275 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Wed, 6 Aug 2025 14:03:33 +0200 Subject: [PATCH 19/23] update gitignore --- examples/.gitignore | 1 + 1 file changed, 1 insertion(+) create mode 100644 examples/.gitignore diff --git a/examples/.gitignore b/examples/.gitignore new file mode 100644 index 0000000..460aa0e --- /dev/null +++ b/examples/.gitignore @@ -0,0 +1 @@ +./data \ No newline at end of file From 6980898b30ac46a74a91e3558be0f8b9682026a8 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Wed, 6 Aug 2025 14:04:52 +0200 Subject: [PATCH 20/23] put together prealignment, rigid alignment and mobie export --- matchmaker/align_rigid_elastix.py | 5 +- matchmaker/compute_registration.py | 27 +++- matchmaker/mobie_export.py | 61 ++++---- matchmaker/prealignment.py | 239 ++++++++++++++--------------- 4 files changed, 176 insertions(+), 156 deletions(-) diff --git a/matchmaker/align_rigid_elastix.py b/matchmaker/align_rigid_elastix.py index 33b54f1..9294e84 100644 --- a/matchmaker/align_rigid_elastix.py +++ b/matchmaker/align_rigid_elastix.py @@ -1,6 +1,4 @@ -from pathlib import Path import numpy as np -import argparse from matchmaker.n5_utils import read_volume, write_volume, get_attrs from matchmaker.vis import plot_overlay import logging @@ -156,4 +154,5 @@ def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_expor if __name__ == "__main__": main() - # python align_rigid_elastix.py ../examples/data/test/platy1_muscles_stardist_fixed_prealigned.n5 seg ../examples/data/test/platy1_muscles_stardist_moving_prealigned.n5 seg ../examples/data/test/platy1_muscles_stardist_moving_rigid_aligned.n5 seg +# python align_rigid_elastix.py -fi ../examples/data/test/platy1_muscles_stardist_fixed_prealigned.n5 -fk seg +# -mi ../examples/data/test/platy1_muscles_stardist_moving_prealigned.n5 -mk seg -o ../examples/data/test -m diff --git a/matchmaker/compute_registration.py b/matchmaker/compute_registration.py index 373975b..902e1aa 100644 --- a/matchmaker/compute_registration.py +++ b/matchmaker/compute_registration.py @@ -1,8 +1,8 @@ import os -import numpy as np import click +from matchmaker.mobie_export import create_mobie_project from matchmaker.prealignment import run_prealignment -from matchmaker.mobie_export import create_mobie_project, export_to_mobie +from matchmaker.align_rigid_elastix import run_rigid_alignment @click.command() @@ -31,7 +31,7 @@ def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_expor output_dir, ) - fixed_prealigned, moving_prealigned = run_prealignment( + run_prealignment( fixed_path, fixed_key, moving_path, @@ -40,6 +40,27 @@ def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_expor mobie_export, ) + fixed_prealigned_path = f"{output_dir}/{os.path.splitext(os.path.basename(fixed_path))[0]}_prealigned.n5" + moving_prealigned_path = f"{output_dir}/{os.path.splitext(os.path.basename(moving_path))[0]}_prealigned.n5" + + run_rigid_alignment( + fixed_prealigned_path, + fixed_key, + moving_prealigned_path, + moving_key, + output_dir, + mobie_export + ) + + # fixed_rigid_aligned_path = f"{output_dir}/{os.path.splitext(os.path.basename(fixed_path))[0]}_rigid_aligned.n5" + # moving_rigid_aligned_path = f"{output_dir}/{os.path.splitext(os.path.basename(moving_path))[0]}_rigid_aligned.n5" + + # run_cpd(...) + if __name__ == "__main__": main() + + +# python compute_registration.py -fi ../examples/data/platy1_muscles_stardist_fixed.n5 -fk seg +# -mi ../examples/data/platy1_muscles_stardist_moving.n5 -mk seg -o ../examples/data/test -m \ No newline at end of file diff --git a/matchmaker/mobie_export.py b/matchmaker/mobie_export.py index 735f67d..5c3c87d 100644 --- a/matchmaker/mobie_export.py +++ b/matchmaker/mobie_export.py @@ -3,46 +3,47 @@ # menu 2: fixed_segmentation; then original, pre-aligned, ... # menu 3: moving_segmentation; then original, pre-aligned, ... import os +import json import mobie # from mobie import add_segmentation -# def add_segmentation_to_mobie(input_path, input_key, seg_name): +def update_default_view(dataset_json_path, new_segmentation_name): + """ + Update the name and sources for the segmentation in the 'default' view + in a MoBIE dataset.json file. -# resolution = get_attrs(input_path, input_key)["resolution"] -# unit = "micrometer" -# chunks = (64, 64, 64) -# scale_factors = 4 * [[2, 2, 2]] + Args: + dataset_json_path (str): Path to the dataset.json file. + new_segmentation_name (str): New segmentation name to set in the default view. + """ + if not os.path.exists(dataset_json_path): + raise FileNotFoundError(f"Could not find: {dataset_json_path}") -# chunks = (128, 128, 128) -# max_jobs = 8 + with open(dataset_json_path, "r") as f: + data = json.load(f) -# add_segmentation( -# input_path=input_path, -# input_key=input_key, -# root=MOBIE_FOLDER, -# dataset_name=DS_NAME, -# segmentation_name=seg_name, -# resolution=resolution, -# scale_factors=scale_factors, -# chunks=chunks, -# file_format="ome.zarr", -# max_jobs=max_jobs, -# ) + views = data.get("views", {}) + default_view = views.get("default", {}) -# def process_segmentation(input_path, input_key, seg_name, scale): -# metadata = read_dataset_metadata(os.path.join(MOBIE_FOLDER, DS_NAME)) -# if seg_name in metadata["sources"]: -# update_segmentation(input_path, input_key, seg_name, scale) -# else: -# add_segmentation_to_mobie(input_path, input_key, seg_name, scale) + source_displays = default_view.get("sourceDisplays", []) + for display in source_displays: + if "segmentationDisplay" in display: + display["segmentationDisplay"]["name"] = new_segmentation_name + display["segmentationDisplay"]["sources"] = [new_segmentation_name] + + # Save the updated dataset.json + with open(dataset_json_path, "w") as f: + json.dump(data, f, indent=2) + + print(f"Updated default view to use segmentation: '{new_segmentation_name}'") def create_mobie_project( fixed_input_path, - fixed_input_key, + fixed_key, moving_input_path, - moving_input_key, + moving_key, output_dir, ): @@ -50,7 +51,7 @@ def create_mobie_project( fixed_file_name = os.path.splitext(os.path.basename(fixed_input_path))[0] export_to_mobie( fixed_input_path, - fixed_input_key, + fixed_key, output_dir, segmentation_name=f"{fixed_file_name}_original", menu_name="fixed", @@ -60,7 +61,7 @@ def create_mobie_project( moving_file_name = os.path.splitext(os.path.basename(moving_input_path))[0] export_to_mobie( moving_input_path, - moving_input_key, + moving_key, output_dir, segmentation_name=f"{moving_file_name}_original", menu_name="moving", @@ -73,7 +74,7 @@ def export_to_mobie(input_path, input_key, output_dir, segmentation_name, menu_n # Set parameters for MOBIE mobie_folder = f"{output_dir}/mobie_project" - dataset_name = "platy1_muscles_stardist" + dataset_name = "platy1_muscles_stardist" # TODO: specify as input # resolution = get_attrs(input_path, input_key)["resolution"] chunks = (64, 64, 64) scale_factors = 4 * [[2, 2, 2]] diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index f0cea95..e119f0e 100755 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -4,10 +4,9 @@ import os import click import sys -import napari from matchmaker.data import create_point_cloud -from matchmaker.mobie_export import export_to_mobie +from matchmaker.mobie_export import export_to_mobie, update_default_view from matchmaker.transform_utils import get_transformation_matrix, rotate_img from matchmaker.n5_utils import read_volume, get_attrs, write_volume from matchmaker.vis import plot_three_slices, plot_overlay @@ -58,17 +57,14 @@ def get_SVD_transform(img, save_path=None): def orient_axis(img, axis, save_path=None): - logging.info(f"Orient along axis {axis}") + # logging.info(f"Orient along axis {axis}") assert img.ndim == 3, f"Input image should have 3 dimensions, has {img.ndim}" if axis == 0: int_profile = np.sum(img, axis=(1, 2)) # profile along z - R_3x3 = np.array([[1, 0, 0], [0, -1, 0], [0, 0, -1]]) if axis == 1: int_profile = np.sum(img, axis=(0, 2)) # profile along y - R_3x3 = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]]) if axis == 2: int_profile = np.sum(img, axis=(0, 1)) # profile along x - R_3x3 = np.array([[1, 0, 0], [0, -1, 0], [0, 0, 1]]) plt.figure() plt.plot(int_profile) @@ -81,7 +77,7 @@ def orient_axis(img, axis, save_path=None): max_pos = int_profile.argmax() logging.info(f"Max position is {max_pos}, dimension shape is {img.shape[axis]}") - if max_pos < img.shape[1] // 2: + if max_pos < img.shape[axis] // 2: logging.info("Correct orientation") return False else: @@ -89,23 +85,6 @@ def orient_axis(img, axis, save_path=None): return True -# def orient_sample(img, axis, save_path=None): -# R_3x3 = orient_axis(img, axis=axis) #, save_path=save_path) -# if R_3x3 is not None: -# img_center = 0.5 * (np.array(img_rotated.shape)-1) -# offset = img_center - R_3x3 @ img_center -# R = np.eye(4) -# R[:3, :3] = R_3x3 -# R[:3, 3] = offset -# img_rotated = rotate_img(img_rotated, R, output_shape=img_rotated.shape) - -# # update transformation matrix -# T = T @ R -# np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) - -# return img_rotated - - def prealign_sample(img, file_name, save_path): """ Pre-align a sample segmentation with its principal components. @@ -123,7 +102,7 @@ def prealign_sample(img, file_name, save_path): img_rotated = rotate_img(img, T, output_shape=new_shape) if np.linalg.det(Vt.T) < 0: - print("WARNING: V includes a reflection (mirroring)") + logging.warning("V includes a reflection (mirroring)") R_3x3 = np.array([[1, 0, 0], [0, 1, 0], [0, 0, -1]]) logging.info("Mirror back...") img_center = 0.5 * (np.array(img_rotated.shape)-1) @@ -137,69 +116,14 @@ def prealign_sample(img, file_name, save_path): T = T @ R np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) - # check orientation in every three dimensions/ if heads up - # for axis in [0, 1, 2]: - # R_3x3 = orient_axis(img_rotated, axis=axis, save_path=f"{save_path}/plots/{file_name}_intensity_profile_{axis}.png") - # if R_3x3 is not None: - # logging.info(f"Rotate axis {axis} 180 degrees to align...") - # img_center = 0.5 * (np.array(img_rotated.shape)-1) - # offset = img_center - R_3x3 @ img_center - # R = np.eye(4) - # R[:3, :3] = R_3x3 - # R[:3, 3] = offset - # img_rotated = rotate_img(img_rotated, R, output_shape=img_rotated.shape) - - # # update transformation matrix - # T = T @ R - # np.savetxt(f"{save_path}/{file_name}_T_prealignment.txt", T) - - plot_three_slices( - img_rotated, - save_path=f"{save_path}/plots/{file_name}_prealigned.png" - ) - - return img_rotated - - -def prealignment_per_image(input_path, input_key, output_dir, mobie_export, image_type): - logging.info("Read image file") - img = read_volume(input_path, input_key) - file_name = os.path.splitext(os.path.basename(input_path))[0] - - # plot three slices of input image - plot_three_slices(img, save_path=f"{output_dir}/plots/{image_type}.png") - - logging.info("Prealign image") - img_prealigned = prealign_sample(img, file_name, save_path=output_dir) - - logging.info("Save prealigned image") - attributes = dict(get_attrs(input_path, input_key)) - write_volume( - f=f"{output_dir}/{file_name}_prealigned.n5", - arr=img_prealigned, - key=input_key, - attrs=attributes - ) - - # export pre-aligned image to MoBIE - if mobie_export: - logging.info("Export prealigned image to MoBIE") - export_to_mobie( - input_path=f"{output_dir}/{file_name}_prealigned.n5", - input_key=input_key, - output_dir=output_dir, - segmentation_name=f"{file_name}_prealigned", - menu_name=image_type - ) - - return img_prealigned + return img_rotated, T def run_prealignment( - fixed_input_path, - fixed_input_key, - moving_input_path, - moving_input_key, + fixed_path, + fixed_key, + moving_path, + moving_key, output_dir, mobie_export=True, ): @@ -216,51 +140,102 @@ def run_prealignment( datefmt="%Y-%m-%d %H:%M:%S", ) logging.info("Start prealignment of fixed image ...") - fixed_prealigned = prealignment_per_image( - fixed_input_path, - fixed_input_key, + fixed_img = read_volume(fixed_path, fixed_key) + fixed_file_name = os.path.splitext(os.path.basename(fixed_path))[0] + + plot_three_slices( + fixed_img, + save_path=f"{output_dir}/plots/{fixed_file_name}_fixed.png" + ) + + fixed_prealigned, _ = prealign_sample( + fixed_img, + fixed_file_name, output_dir, - mobie_export, - image_type="fixed" ) logging.info("Start prealignment of moving image ...") - moving_prealigned = prealignment_per_image( - moving_input_path, - moving_input_key, + moving_img = read_volume(moving_path, moving_key) + moving_file_name = os.path.splitext(os.path.basename(moving_path))[0] + + plot_three_slices( + moving_img, + save_path=f"{output_dir}/plots/{moving_file_name}_moving.png" + ) + + moving_prealigned, T_moving = prealign_sample( + moving_img, + moving_file_name, output_dir, - mobie_export, - image_type="moving", ) - # check orientation + + # check orientation (if moving fits to fixed) + R_3x3 = np.eye(3) + change_orientation = False for axis in range(3): - rotate_axis_fixed = orient_axis(moving_prealigned, axis=axis, save_path=f"{output_dir}/plots/moving_intensity_profile_{axis}.png") - rotate_axis_moving = orient_axis(fixed_prealigned, axis=axis, save_path=f"{output_dir}/plots/fixed_intensity_profile_{axis}.png") - - if not rotate_axis_fixed and not rotate_axis_moving: - print(f"Correct orientation in axis {axis}") - else: + rotate_axis_fixed = orient_axis( + fixed_prealigned, + axis=axis, + save_path=f"{output_dir}/plots/fixed_intensity_profile_{axis}.png", + ) + rotate_axis_moving = orient_axis( + moving_prealigned, + axis=axis, + save_path=f"{output_dir}/plots/moving_intensity_profile_{axis}.png", + ) + + if rotate_axis_fixed or rotate_axis_moving: print(f"Rotate axis {axis} 180 degrees to align...") - R_3x3 = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 1]]) - logging.info(f"Rotate axis {axis} 180 degrees to align...") - img_center = 0.5 * (np.array(moving_prealigned.shape)-1) - offset = img_center - R_3x3 @ img_center - R = np.eye(4) - R[:3, :3] = R_3x3 - R[:3, 3] = offset - moving_prealigned = rotate_img(moving_prealigned, R, output_shape=moving_prealigned.shape) - - # update transformation matrix - # T = T @ R - # np.savetxt(f"{output_dir}/moving_T_prealignment.txt", T) - + change_orientation = True + R_3x3[axis, axis] = -1 + else: + print(f"Correct orientation in axis {axis}") + + if not change_orientation: + logging.info("Correct orientation") + else: + logging.info("Rotate moving image to align orientation...") + logging.info(str(R_3x3)) + img_center = 0.5 * (np.array(moving_prealigned.shape)-1) + offset = img_center - R_3x3 @ img_center + R = np.eye(4) + R[:3, :3] = R_3x3 + R[:3, 3] = offset + moving_prealigned = rotate_img(moving_prealigned, R, output_shape=moving_prealigned.shape) + + # update transformation matrix + T_moving = T_moving @ R + np.savetxt(f"{output_dir}/moving_T_prealignment.txt", T_moving) + logging.info("Prealignment done.") - v = napari.Viewer() - v.add_labels(fixed_prealigned) - v.add_labels(moving_prealigned) - napari.run() - breakpoint() + logging.info("Save prealigned fixed image") + attributes = dict(get_attrs(fixed_path, fixed_key)) + write_volume( + f=f"{output_dir}/{fixed_file_name}_prealigned.n5", + arr=fixed_prealigned, + key=fixed_key, + attrs=attributes + ) + + logging.info("Save prealigned moving image") + attributes = dict(get_attrs(moving_path, moving_key)) + write_volume( + f=f"{output_dir}/{moving_file_name}_prealigned.n5", + arr=moving_prealigned, + key=moving_key, + attrs=attributes + ) + + plot_three_slices( + fixed_prealigned, + save_path=f"{output_dir}/plots/{fixed_file_name}_prealigned.png" + ) + + plot_three_slices( + moving_prealigned, + save_path=f"{output_dir}/plots/{moving_file_name}_prealigned.png" + ) plot_overlay( fixed_prealigned, @@ -268,6 +243,30 @@ def run_prealignment( save_path=f"{output_dir}/plots/overlay_prealignment.png", ) + if mobie_export: + logging.info("Export prealigned images to MoBIE") + export_to_mobie( + input_path=f"{output_dir}/{fixed_file_name}_prealigned.n5", + input_key=fixed_key, + output_dir=output_dir, + segmentation_name=f"{fixed_file_name}_prealigned", + menu_name="fixed" + ) + update_default_view( + dataset_json_path=f"{output_dir}/mobie_project/platy1_muscles_stardist/dataset.json", + # TODO: make name independent + new_segmentation_name=f"{fixed_file_name}_prealigned" + ) + + export_to_mobie( + input_path=f"{output_dir}/{moving_file_name}_prealigned.n5", + input_key=moving_key, + output_dir=output_dir, + segmentation_name=f"{moving_file_name}_prealigned", + menu_name="moving" + ) + logging.info("MoBIE export done.") + @click.command() @click.option("-fi", "--fixed_path", required=True, help="Fixed input .n5 file") From ba92aeb28aa6ad6ccd6b227e7e75e40fd07986eb Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Thu, 7 Aug 2025 09:47:37 +0200 Subject: [PATCH 21/23] update readme, docstrings and logging --- .gitignore | 5 +- README.md | 20 ++++--- examples/.gitignore | 1 - matchmaker/align_rigid_elastix.py | 83 ++++++++++++++++++---------- matchmaker/compute_registration.py | 38 ++++++++++--- matchmaker/elastix_utils.py | 2 +- matchmaker/mobie_export.py | 55 ++++++++++++++----- matchmaker/n5_utils.py | 8 --- matchmaker/prealignment.py | 86 ++++++++++++++++++++---------- matchmaker/transform_utils.py | 22 ++++++++ matchmaker/vis.py | 6 ++- 11 files changed, 229 insertions(+), 97 deletions(-) delete mode 100644 examples/.gitignore diff --git a/.gitignore b/.gitignore index e8ca84b..a00d509 100644 --- a/.gitignore +++ b/.gitignore @@ -174,10 +174,11 @@ cython_debug/ .pypirc # huge data -./examples/data/platy1_muscles_stardist.tif +*.tif *.n5/ *.png +*.ome.zarr/ .DS_Store -./examples/data/* \ No newline at end of file +./examples/data/test/ \ No newline at end of file diff --git a/README.md b/README.md index 700c282..aa447e6 100644 --- a/README.md +++ b/README.md @@ -21,14 +21,22 @@ Expected image shape: ZYX ### Registration steps -**1. PCA pre-alignment**: rigid alignment of fixed and moving image using the PCs \ - `prealignment.py --input_path ... --input_key ... --output_dir --mobie_export --type ...` \ - output: prealigned images + **rigid transformation matrix** -2. Rigid pre-alignment with Elastix OR with rigid CPD ---> **rigid transform matrix** +**1. PCA pre-alignment**: alignment of fixed and moving image to the PCs \ + `prealignment.py --fixed_path ... --fixed_key ... --moving_path ... --moving_key ... --output_dir ... --mobie_export --dataset_name ...` \ + output: prealigned images + **transformation matrix** + +**2. Rigid pre-alignment with Elastix** \ +`apply_rigid_elastix.py --fixed_path ... --fixed_key ... --moving_path ... --moving_key ... --output_dir ... --mobie_export --dataset_name ...` \ +output: rigid alinged moving image + **rigid transformation matrix** **3. Coherent point drift** -5. Matching points with mixed integer programming -6. Deformable registration with Elastix with distance between keypoints in loss and rigidity penalty ---> **B-spline coefficients** + +**4. Matching points with mixed integer programming** + +**5. Deformable registration with Elastix** \ +with distance between keypoints in loss and rigidity penalty ---> **B-spline coefficients** + + ## Apply registration to other images diff --git a/examples/.gitignore b/examples/.gitignore deleted file mode 100644 index 460aa0e..0000000 --- a/examples/.gitignore +++ /dev/null @@ -1 +0,0 @@ -./data \ No newline at end of file diff --git a/matchmaker/align_rigid_elastix.py b/matchmaker/align_rigid_elastix.py index 9294e84..69f03fb 100644 --- a/matchmaker/align_rigid_elastix.py +++ b/matchmaker/align_rigid_elastix.py @@ -32,16 +32,10 @@ def elastix_segm_rigid_alignment( logging.info("Moving image") logging.info(f"{moving_img}") - plot_overlay( - elastix_utils.itk_to_np_order(itk.GetArrayFromImage(fixed_img)), - elastix_utils.itk_to_np_order(itk.GetArrayFromImage(moving_img)), - f"{output_dir}/plots/intersample_segm_overlay_before_alignment.png", - ) - parameter_map_paths = [ "../ParameterMap_segm_rigid_registration_corr.txt" ] - logging.info("Run rigid registration") + logging.info("Run rigid registration with elastix") result_image, result_transform_parameters = elastix_utils.run_registration( fixed_img, moving_img, @@ -56,7 +50,7 @@ def elastix_segm_rigid_alignment( plot_overlay( elastix_utils.itk_to_np_order(itk.GetArrayFromImage(fixed_img)), result_img_np, - f"{output_dir}/plots/intersample_segm_rigid_alignment_semantic.png", + f"{output_dir}/plots/overlay_after_rigid_alignment.png", ) # NOTE: difference between result_img_np before and after applying transform? logging.info("Apply transform to all channels") @@ -74,24 +68,37 @@ def run_rigid_alignment( moving_path, moving_key, output_dir, - mobie_export + mobie_export, + dataset_name, ): + """ + Perform rigid alignment of a moving image to a fixed image using Elastix. + + This function reads the fixed and moving images from the specified paths, + performs a rigid alignment using Elastix, and saves the aligned moving image + to the output directory. If the MoBIE export flag is set, the aligned image + is also exported to a MoBIE project. + + Args: + fixed_path (str): Path to the fixed image .n5 file. + fixed_key (str): Key to the fixed image data in the .n5 file. + moving_path (str): Path to the moving image .n5 file. + moving_key (str): Key to the moving image data in the .n5 file. + output_dir (str): Directory where the aligned image should be saved. + mobie_export (bool): Flag indicating whether to export to a MoBIE project. + dataset_name (str): Name of the dataset for the MoBIE export. + + Returns: + np.ndarray: The rigidly aligned moving image. + """ + if not os.path.exists(f"{output_dir}/plots"): os.makedirs(f"{output_dir}/plots") - logging.basicConfig( - level=logging.INFO, - format="%(asctime)s [%(levelname)s] %(message)s", - handlers=[ - logging.FileHandler(f"{output_dir}/rigid_alignment.log", mode="w"), - logging.StreamHandler(sys.stdout), - ], - datefmt="%Y-%m-%d %H:%M:%S", - ) + logging.info("Start rigid alignment") logging.info("Read image file") fixed_img_np = read_volume(fixed_path, fixed_key) - print(fixed_img_np.dtype) if fixed_path == moving_path: logging.info("Same n5 for fixed and moving image, not doing registration") @@ -103,7 +110,7 @@ def run_rigid_alignment( fixed_img_np = fixed_img_np.astype(np.float32) moving_img_np = moving_img_np.astype(np.float32) - logging.info("Start registration") + logging.info("Compute rigid alignment ...") moving_img_np = elastix_segm_rigid_alignment( fixed_img_np=fixed_img_np, @@ -113,27 +120,27 @@ def run_rigid_alignment( output_dir=output_dir, ) moving_img_np = moving_img_np.astype(np.uint16) - logging.info("Write results") - attributes = dict(get_attrs(moving_path, moving_key)) + logging.info("Save rigid aligned moving image") + attributes = dict(get_attrs(moving_path, moving_key)) file_name = os.path.splitext(os.path.basename(moving_path))[0] file_name = file_name.removesuffix("_prealigned") - print(file_name) write_volume( - f"{output_dir}/{file_name}_rigid_aligned.n5", - moving_img_np, - moving_key, - chunks=(128, 512, 512), + f=f"{output_dir}/{file_name}_rigid_aligned.n5", + arr=moving_img_np, + key=moving_key, attrs=attributes, ) # export rigid alignment to mobie if mobie_export: + logging.info("Export rigid aligned moving image to MoBIE") export_to_mobie( input_path=f"{output_dir}/{file_name}_rigid_aligned.n5", input_key=moving_key, output_dir=output_dir, + dataset_name=dataset_name, segmentation_name=f"{file_name}_rigid_aligned", menu_name="moving" ) @@ -148,11 +155,29 @@ def run_rigid_alignment( @click.option("-m", "--mobie_export", required=False, is_flag=True, help="MoBIE export") def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export): - run_rigid_alignment(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export) + logging.basicConfig( + level=logging.INFO, + format="%(asctime)s [%(levelname)s] %(message)s", + handlers=[ + logging.FileHandler(f"{output_dir}/rigid_alignment.log", mode="w"), + logging.StreamHandler(sys.stdout), + ], + datefmt="%Y-%m-%d %H:%M:%S", + ) + + run_rigid_alignment( + fixed_path, + fixed_key, + moving_path, + moving_key, + output_dir, + mobie_export, + dataset_name="platy1_muscles_stardist", + ) if __name__ == "__main__": main() -# python align_rigid_elastix.py -fi ../examples/data/test/platy1_muscles_stardist_fixed_prealigned.n5 -fk seg +# python align_rigid_elastix.py -fi ../examples/data/test/platy1_muscles_stardist_fixed_prealigned.n5 -fk seg # -mi ../examples/data/test/platy1_muscles_stardist_moving_prealigned.n5 -mk seg -o ../examples/data/test -m diff --git a/matchmaker/compute_registration.py b/matchmaker/compute_registration.py index 902e1aa..ae52bc8 100644 --- a/matchmaker/compute_registration.py +++ b/matchmaker/compute_registration.py @@ -1,5 +1,7 @@ import os +import sys import click +import logging from matchmaker.mobie_export import create_mobie_project from matchmaker.prealignment import run_prealignment from matchmaker.align_rigid_elastix import run_rigid_alignment @@ -13,15 +15,36 @@ @click.option("-o", "--output_dir", required=True, help="Output directory") @click.option("-m", "--mobie_export", required=False, is_flag=True, help="MoBIE export") def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export): - # fixed_input_path = "../examples/data/platy1_muscles_stardist_fixed.n5" - # fixed_input_key = "seg" - # moving_input_path = "../examples/data/platy1_muscles_stardist_moving.n5" - # moving_input_key = "seg" + """ + Main function to perform registration of moving image to fixed image. + + This function orchestrates the sequence of steps required to register a moving image to a fixed image. + It handles the creation of necessary directories, configures logging, and invokes prealignment and + rigid alignment functions. Optionally, it can create a MoBIE project for visualization. + + Args: + fixed_path (str): Path to the fixed input .n5 file. + fixed_key (str): Key to the fixed image data in the .n5 file. + moving_path (str): Path to the moving input .n5 file. + moving_key (str): Key to the moving image data in the .n5 file. + output_dir (str): Directory where the results should be saved. + mobie_export (bool): Flag indicating whether to export results to a MoBIE project. + """ - # output_dir = "../examples/data/test" if not os.path.exists(f"{output_dir}/plots"): os.makedirs(f"{output_dir}/plots") + logging.basicConfig( + level=logging.INFO, + format="%(asctime)s [%(levelname)s] %(message)s", + handlers=[ + logging.FileHandler(f"{output_dir}/registration.log", mode="w"), + logging.StreamHandler(sys.stdout), + ], + datefmt="%Y-%m-%d %H:%M:%S", + ) + + dataset_name = "platy1_muscles_stardist" if mobie_export: create_mobie_project( fixed_path, @@ -29,6 +52,7 @@ def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_expor moving_path, moving_key, output_dir, + dataset_name ) run_prealignment( @@ -38,6 +62,7 @@ def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_expor moving_key, output_dir, mobie_export, + dataset_name, ) fixed_prealigned_path = f"{output_dir}/{os.path.splitext(os.path.basename(fixed_path))[0]}_prealigned.n5" @@ -49,7 +74,8 @@ def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_expor moving_prealigned_path, moving_key, output_dir, - mobie_export + mobie_export, + dataset_name, ) # fixed_rigid_aligned_path = f"{output_dir}/{os.path.splitext(os.path.basename(fixed_path))[0]}_rigid_aligned.n5" diff --git a/matchmaker/elastix_utils.py b/matchmaker/elastix_utils.py index 738b67f..55c2b4c 100644 --- a/matchmaker/elastix_utils.py +++ b/matchmaker/elastix_utils.py @@ -141,7 +141,7 @@ def serialize_parameter_object(parameter_object, prefix, write_dir): for index in range(parameter_object.GetNumberOfParameterMaps()): parameter_map = parameter_object.GetParameterMap(index) parameter_object.WriteParameterFile( - parameter_map, write_dir / f"{prefix}.{index}.txt" + parameter_map, write_dir / f"{prefix}_{index}.txt" ) diff --git a/matchmaker/mobie_export.py b/matchmaker/mobie_export.py index 5c3c87d..f8ba736 100644 --- a/matchmaker/mobie_export.py +++ b/matchmaker/mobie_export.py @@ -1,11 +1,8 @@ -# export results as mobie project -# menu 1: dataset -# menu 2: fixed_segmentation; then original, pre-aligned, ... -# menu 3: moving_segmentation; then original, pre-aligned, ... import os import json +import logging import mobie -# from mobie import add_segmentation +from matchmaker.n5_utils import get_attrs def update_default_view(dataset_json_path, new_segmentation_name): @@ -36,7 +33,7 @@ def update_default_view(dataset_json_path, new_segmentation_name): with open(dataset_json_path, "w") as f: json.dump(data, f, indent=2) - print(f"Updated default view to use segmentation: '{new_segmentation_name}'") + logging.info(f"Updated default view to use segmentation: '{new_segmentation_name}'") def create_mobie_project( @@ -45,7 +42,18 @@ def create_mobie_project( moving_input_path, moving_key, output_dir, + dataset_name ): + """ + Create initial MoBIE project with the fixed and moving images. + + Args: + fixed_input_path (str): Path to the fixed image n5 file. + fixed_key (str): Key to the fixed image data in the n5 file. + moving_input_path (str): Path to the moving image n5 file. + moving_key (str): Key to the moving image data in the n5 file. + output_dir (str): Directory where the MoBIE project should be saved. + """ # export fixed image to MoBIE fixed_file_name = os.path.splitext(os.path.basename(fixed_input_path))[0] @@ -53,6 +61,7 @@ def create_mobie_project( fixed_input_path, fixed_key, output_dir, + dataset_name, segmentation_name=f"{fixed_file_name}_original", menu_name="fixed", ) @@ -63,19 +72,31 @@ def create_mobie_project( moving_input_path, moving_key, output_dir, + dataset_name, segmentation_name=f"{moving_file_name}_original", menu_name="moving", ) + logging.info("Created initial MoBIE project with fixed and moving images.") + +def export_to_mobie(input_path, input_key, output_dir, dataset_name, segmentation_name, menu_name): + """ + Export segmentation from n5 file to MoBIE project. -def export_to_mobie(input_path, input_key, output_dir, segmentation_name, menu_name): + Args: + input_path (str): Path to the n5 file containing the segmentation. + input_key (str): Key to the segmentation data in the n5 file. + output_dir (str): Directory where the MoBIE project should be saved. + dataset_name (str): Name of the MoBIE dataset. + segmentation_name (str): Name of the segmentation in the MoBIE project. + menu_name (str): Name of the menu in the MoBIE project. + """ if not os.path.exists(f"{output_dir}/mobie_project"): os.makedirs(f"{output_dir}/mobie_project") # Set parameters for MOBIE mobie_folder = f"{output_dir}/mobie_project" - dataset_name = "platy1_muscles_stardist" # TODO: specify as input - # resolution = get_attrs(input_path, input_key)["resolution"] + resolution = get_attrs(input_path, input_key)["resolution"] chunks = (64, 64, 64) scale_factors = 4 * [[2, 2, 2]] @@ -85,23 +106,31 @@ def export_to_mobie(input_path, input_key, output_dir, segmentation_name, menu_n root=mobie_folder, dataset_name=dataset_name, segmentation_name=segmentation_name, - resolution=[1, 1, 1], + resolution=resolution, scale_factors=scale_factors, chunks=chunks, menu_name=menu_name, file_format="ome.zarr", is_default_dataset=True ) + logging.info(f"Added segmentation: {segmentation_name}") def main(): - input_path = "/Users/marei/git-repositories/matchmaker/examples/CLI_test/platy1_muscles_stardist_fixed_prealigned.n5" + input_path = "../examples/CLI_test/platy1_muscles_stardist_fixed_prealigned.n5" input_key = "seg" - output_dir = "/Users/marei/git-repositories/matchmaker/examples/data/test" + output_dir = "../examples/data/test" file_name = os.path.splitext(os.path.basename(input_path))[0] - export_to_mobie(input_path, input_key, output_dir, segmentation_name=f"{file_name}_prealigned", menu_name="fixed") + export_to_mobie( + input_path, + input_key, + output_dir, + dataset_name="platy1_muscles_stardist", + segmentation_name=f"{file_name}_prealigned", + menu_name="fixed", + ) print(f"MoBIE project created at {output_dir}/mobie_project") diff --git a/matchmaker/n5_utils.py b/matchmaker/n5_utils.py index a822b2d..aea3fe4 100644 --- a/matchmaker/n5_utils.py +++ b/matchmaker/n5_utils.py @@ -11,7 +11,6 @@ def print_key_tree(f: z5py.File): def read_volume( f: z5py.File, key: str, roi: np.lib.index_tricks.IndexExpression = np.s_[:] ): - # print(type(f)) if isinstance(f, (str, PurePath)): f = z5py.File(f, "r") @@ -23,15 +22,12 @@ def read_volume( return None ds.n_threads = 8 - # print(f"Reading roi {roi} of volume {key} from {f.filename}") vol = ds[roi] - # print(f"Read volume with shape {vol.shape}, data type {vol.dtype}") return vol def get_attrs(f: z5py.File, key: str): - # print(type(f)) if isinstance(f, (str, PurePath)): f = z5py.File(f, "a") @@ -54,21 +50,17 @@ def write_volume(f, arr: np.array, key, chunks=(1, 512, 512), attrs=None): f = z5py.File(f, "a") if key not in f.keys(): - # print(f"Created dataset {key}") ds = f.create_dataset( key, shape=shape, compression=compression, chunks=chunks, dtype=dtype ) else: - # print(f"Overwriting {key}") ds = f[key] ds.n_threads = 8 - # print(f"Writing array to {key}") ds[:] = arr print(f"Dataset {key} written to {f.filename}") if attrs is not None: - # print("Assigning attributes of the dataset") for key in attrs.keys(): ds.attrs[key] = attrs[key] diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index e119f0e..2095f4e 100755 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -42,13 +42,12 @@ def get_SVD_transform(img, save_path=None): logging.info("Rotate point cloud") vr = pos_c @ Vt.T - # TODO: update axis labels plt.figure(figsize=(10, 5)) plt.subplot(1, 2, 1) - plt.title("original vertices") + plt.title("Original Vertices") plt.scatter(pos_c[:, 0], pos_c[:, 1], alpha=0.2) plt.subplot(1, 2, 2) - plt.title("rotated vertices") + plt.title("Rotated Vertices") plt.scatter(vr[:, 0], vr[:, 1], alpha=0.1) if save_path: plt.savefig(save_path, dpi=300) @@ -57,7 +56,18 @@ def get_SVD_transform(img, save_path=None): def orient_axis(img, axis, save_path=None): - # logging.info(f"Orient along axis {axis}") + """ + Plot the sum intensity along the given axis and return True if the maximum is closer to the + upper boundary than the lower boundary, False otherwise. + + Args: + img: 3D image + axis: Axis to sum along + save_path: Path to save the plot to. If None, show the plot instead. + + Returns: + True if the maximum is closer to the upper boundary than the lower boundary, False otherwise. + """ assert img.ndim == 3, f"Input image should have 3 dimensions, has {img.ndim}" if axis == 0: int_profile = np.sum(img, axis=(1, 2)) # profile along z @@ -68,7 +78,7 @@ def orient_axis(img, axis, save_path=None): plt.figure() plt.plot(int_profile) - plt.xlabel("Y Coordinate") + plt.xlabel(f"Axis {axis} Coordinate") plt.ylabel(f"Sum intensity along axis = {axis}") if save_path is not None: plt.savefig(save_path, dpi=300) @@ -78,10 +88,8 @@ def orient_axis(img, axis, save_path=None): max_pos = int_profile.argmax() logging.info(f"Max position is {max_pos}, dimension shape is {img.shape[axis]}") if max_pos < img.shape[axis] // 2: - logging.info("Correct orientation") return False else: - logging.info("Rotate 180 degree to align") return True @@ -125,20 +133,13 @@ def run_prealignment( moving_path, moving_key, output_dir, - mobie_export=True, + mobie_export, + dataset_name ): if not os.path.exists(f"{output_dir}/plots"): os.makedirs(f"{output_dir}/plots") - logging.basicConfig( - level=logging.INFO, - format="%(asctime)s [%(levelname)s] %(message)s", - handlers=[ - logging.FileHandler(f"{output_dir}/prealignment.log", mode="w"), - logging.StreamHandler(sys.stdout), - ], - datefmt="%Y-%m-%d %H:%M:%S", - ) + logging.info("Start prealignment") logging.info("Start prealignment of fixed image ...") fixed_img = read_volume(fixed_path, fixed_key) fixed_file_name = os.path.splitext(os.path.basename(fixed_path))[0] @@ -163,6 +164,12 @@ def run_prealignment( save_path=f"{output_dir}/plots/{moving_file_name}_moving.png" ) + plot_overlay( + fixed_img, + moving_img, + save_path=f"{output_dir}/plots/overlay_before.png", + ) + moving_prealigned, T_moving = prealign_sample( moving_img, moving_file_name, @@ -170,18 +177,19 @@ def run_prealignment( ) # check orientation (if moving fits to fixed) + logging.info("Check axis orientation ...") R_3x3 = np.eye(3) change_orientation = False for axis in range(3): rotate_axis_fixed = orient_axis( fixed_prealigned, axis=axis, - save_path=f"{output_dir}/plots/fixed_intensity_profile_{axis}.png", + save_path=f"{output_dir}/plots/fixed_prealigned_intensity_profile_{axis}.png", ) rotate_axis_moving = orient_axis( moving_prealigned, axis=axis, - save_path=f"{output_dir}/plots/moving_intensity_profile_{axis}.png", + save_path=f"{output_dir}/plots/moving_prealigned_intensity_profile_{axis}.png", ) if rotate_axis_fixed or rotate_axis_moving: @@ -189,10 +197,10 @@ def run_prealignment( change_orientation = True R_3x3[axis, axis] = -1 else: - print(f"Correct orientation in axis {axis}") + print(f"Correct orientation in axis {axis}.") if not change_orientation: - logging.info("Correct orientation") + logging.info("Correct orientation.") else: logging.info("Rotate moving image to align orientation...") logging.info(str(R_3x3)) @@ -205,11 +213,11 @@ def run_prealignment( # update transformation matrix T_moving = T_moving @ R - np.savetxt(f"{output_dir}/moving_T_prealignment.txt", T_moving) + np.savetxt(f"{output_dir}/{moving_file_name}_T_prealignment.txt", T_moving) logging.info("Prealignment done.") - logging.info("Save prealigned fixed image") + logging.info("Save prealigned fixed image ...") attributes = dict(get_attrs(fixed_path, fixed_key)) write_volume( f=f"{output_dir}/{fixed_file_name}_prealigned.n5", @@ -218,7 +226,7 @@ def run_prealignment( attrs=attributes ) - logging.info("Save prealigned moving image") + logging.info("Save prealigned moving image ...") attributes = dict(get_attrs(moving_path, moving_key)) write_volume( f=f"{output_dir}/{moving_file_name}_prealigned.n5", @@ -240,21 +248,22 @@ def run_prealignment( plot_overlay( fixed_prealigned, moving_prealigned, - save_path=f"{output_dir}/plots/overlay_prealignment.png", + save_path=f"{output_dir}/plots/overlay_after_prealignment.png", ) if mobie_export: - logging.info("Export prealigned images to MoBIE") + logging.info("Export prealigned images to MoBIE ...") export_to_mobie( input_path=f"{output_dir}/{fixed_file_name}_prealigned.n5", input_key=fixed_key, output_dir=output_dir, + dataset_name=dataset_name, segmentation_name=f"{fixed_file_name}_prealigned", menu_name="fixed" ) + logging.info("Change default dataset to prealigned fixed image.") update_default_view( - dataset_json_path=f"{output_dir}/mobie_project/platy1_muscles_stardist/dataset.json", - # TODO: make name independent + dataset_json_path=f"{output_dir}/mobie_project/{dataset_name}/dataset.json", new_segmentation_name=f"{fixed_file_name}_prealigned" ) @@ -262,6 +271,7 @@ def run_prealignment( input_path=f"{output_dir}/{moving_file_name}_prealigned.n5", input_key=moving_key, output_dir=output_dir, + dataset_name=dataset_name, segmentation_name=f"{moving_file_name}_prealigned", menu_name="moving" ) @@ -277,7 +287,25 @@ def run_prealignment( @click.option("-m", "--mobie_export", required=False, is_flag=True, help="MoBIE export") def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export): - run_prealignment(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export) + logging.basicConfig( + level=logging.INFO, + format="%(asctime)s [%(levelname)s] %(message)s", + handlers=[ + logging.FileHandler(f"{output_dir}/prealignment.log", mode="w"), + logging.StreamHandler(sys.stdout), + ], + datefmt="%Y-%m-%d %H:%M:%S", + ) + + run_prealignment( + fixed_path, + fixed_key, + moving_path, + moving_key, + output_dir, + mobie_export, + dataset_name="platy1_muscles_stardist", + ) if __name__ == "__main__": diff --git a/matchmaker/transform_utils.py b/matchmaker/transform_utils.py index 794ca1c..f05878f 100644 --- a/matchmaker/transform_utils.py +++ b/matchmaker/transform_utils.py @@ -111,6 +111,28 @@ def get_transformation_matrix(img, gc, Vt, save_path=None): def rotate_img(img, rotation_matrix, output_shape=None, offset=None): + """ + Rotate an image using a given rotation matrix. + + Parameters + ---------- + img : array + The 3D image to be rotated. + rotation_matrix : array + A 3x3 or 4x4 rotation matrix. + output_shape : tuple, optional + The desired output shape of the rotated image. If not given, the output shape + will be determined from the rotation matrix. + offset : tuple, optional + The offset to apply to the rotated image to ensure it fits within the new + bounding box. If not given, the offset will be determined from the rotation + matrix. + + Returns + ------- + rotated_img : array + The rotated image with the desired output shape (if given). + """ rotated_img = affine_transform( img, matrix=rotation_matrix, diff --git a/matchmaker/vis.py b/matchmaker/vis.py index c2a954a..d5b6a3f 100644 --- a/matchmaker/vis.py +++ b/matchmaker/vis.py @@ -13,7 +13,8 @@ def plot_three_slices( max_pos=False, alpha=False, ): - """Plot slices of a 3D image along each axis. + """ + Plot slices of a 3D image along each axis. Args: img: _description_ @@ -57,7 +58,8 @@ def plot_three_slices( def plot_overlay(img1, img2, save_path=None, x_pos=None, y_pos=None, z_pos=None): - """Plot slices of two 3D images along each axis. + """ + Plot slices of two 3D images along each axis. Args: img1: _description_target_shape = (25, 22, 29) From eda094070f57005f111ed798763c2d8f4afda7d0 Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Thu, 7 Aug 2025 09:50:18 +0200 Subject: [PATCH 22/23] update gitignore --- .gitignore | 7 ++++--- examples/data/platy1_muscles_stardist_fixed.tif | 4 ++-- examples/data/platy1_muscles_stardist_moving.tif | 4 ++-- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/.gitignore b/.gitignore index a00d509..30dc145 100644 --- a/.gitignore +++ b/.gitignore @@ -173,12 +173,13 @@ cython_debug/ # PyPI configuration file .pypirc +# test results +examples/data/test/ + # huge data *.tif *.n5/ *.png *.ome.zarr/ -.DS_Store - -./examples/data/test/ \ No newline at end of file +.DS_Store \ No newline at end of file diff --git a/examples/data/platy1_muscles_stardist_fixed.tif b/examples/data/platy1_muscles_stardist_fixed.tif index d43e439..02e6673 100644 --- a/examples/data/platy1_muscles_stardist_fixed.tif +++ b/examples/data/platy1_muscles_stardist_fixed.tif @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:cbbf90829a7830ec3cc57c048661ecddf864377b3c1c0fccbc15794b2f64aead -size 163013402 +oid sha256:bf565c9e84d198c58dc92d432a7b4a4ea2b7519e62b4d8f81a8d1ee85ad83ee9 +size 13164762 diff --git a/examples/data/platy1_muscles_stardist_moving.tif b/examples/data/platy1_muscles_stardist_moving.tif index 43600dc..a611bd2 100644 --- a/examples/data/platy1_muscles_stardist_moving.tif +++ b/examples/data/platy1_muscles_stardist_moving.tif @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:f4157af6f678cdcee19d4cff13445d760aa9fdf540f9ededc5058c9c28d0b3c4 -size 114242416 +oid sha256:ea6991f9a3d168a8dfd9ec0e8b83d38f73de778f8d1bc83a206d251f0ad54038 +size 60386230 From 189ca38e9ca80870b327943e9cf89c4b2a2565cc Mon Sep 17 00:00:00 2001 From: Marei Freitag Date: Wed, 24 Sep 2025 10:32:50 +0200 Subject: [PATCH 23/23] update Readme --- README.md | 24 ++++++- ...{clean_up.py => mobie_dataset_clean_up.py} | 0 matchmaker/prealignment.py | 70 +++++++++++++++++++ 3 files changed, 92 insertions(+), 2 deletions(-) rename examples/{clean_up.py => mobie_dataset_clean_up.py} (100%) diff --git a/README.md b/README.md index aa447e6..9cfd649 100644 --- a/README.md +++ b/README.md @@ -17,17 +17,37 @@ Expected image shape: ZYX - **QC plots** - **Files with all transforms** - **Table of correspondence between instances in moving and fixed instance segmentations** +- **Logging file: `registration.log`** +- **Optional: Mobie project saved at `{output_dir}/mobie_project/`** ### Registration steps **1. PCA pre-alignment**: alignment of fixed and moving image to the PCs \ `prealignment.py --fixed_path ... --fixed_key ... --moving_path ... --moving_key ... --output_dir ... --mobie_export --dataset_name ...` \ - output: prealigned images + **transformation matrix** + +Outputs: +- prealigned images: `{file_name}_prealigned.n5` +- transformation matrixes + - `{file_name}_fixed_T_prealignment.txt` + - `{file_name}_moving_T_prealignment.txt` (maybe final one is `moving_T_prealignment.txt`, couldn't figure this out) +- plots: + - slice per dimension before pre-alignment: `{file_name}_fixed.png`, `{file_name}_moving.png` + - slice per dimension after pre-alignment: `{file_name}_fixed_prealigned.png`, `{file_name}_moving_prealigned.png` + - overlay of slice per dimension after pre-alignment: `overlay_prealignment.png` + - intensity profiles per axis and volume: `fixed_intensity_profile_{axis}.png`, `moving_intensity_profile_{axis}.png` **2. Rigid pre-alignment with Elastix** \ `apply_rigid_elastix.py --fixed_path ... --fixed_key ... --moving_path ... --moving_key ... --output_dir ... --mobie_export --dataset_name ...` \ -output: rigid alinged moving image + **rigid transformation matrix** + +Outputs: +- rigid alinged moving image: `{file_name}_rigid_aligned.n5` +- rigid transformation matrix (Elastix outputs): `result.0.mhd`, `result.0.raw`, `TransformParameters.0.txt` +- logging file: `elastix_log_rigid.log` +- plots: + - `intersample_segm_overlay_before_alignment.png` + - `intersample_segm_rigid_alignment_semantic.png` + **3. Coherent point drift** diff --git a/examples/clean_up.py b/examples/mobie_dataset_clean_up.py similarity index 100% rename from examples/clean_up.py rename to examples/mobie_dataset_clean_up.py diff --git a/matchmaker/prealignment.py b/matchmaker/prealignment.py index 2095f4e..81ddef9 100755 --- a/matchmaker/prealignment.py +++ b/matchmaker/prealignment.py @@ -136,6 +136,59 @@ def run_prealignment( mobie_export, dataset_name ): + """ + Run prealignment of a fixed and moving 3D image volume. + + This function reads two volumetric datasets (a fixed and a moving image), + performs prealignment to roughly register them into a common space, and + saves diagnostic plots, transformation matrices, and prealigned volumes. + It also checks axis orientation consistency between the two images and + applies corrective rotations if necessary. Optionally, the results can be + exported into a MoBIE project for interactive visualization. + + Steps performed: + 1. Load fixed and moving volumes. + 2. Plot reference slices and overlays before alignment. + 3. Apply prealignment to both volumes. + 4. Check and correct axis orientations if required. + 5. Save prealigned volumes and transformation matrices. + 6. Generate plots before and after prealignment. + 7. Optionally export results to a MoBIE project. + + Parameters + ---------- + fixed_path : str + Path to the fixed volume file (e.g. N5, OME-Zarr). + fixed_key : str + Dataset key inside the fixed volume file. + moving_path : str + Path to the moving volume file (e.g. N5, OME-Zarr). + moving_key : str + Dataset key inside the moving volume file. + output_dir : str + Directory where outputs (plots, volumes, transformations) will be saved. + mobie_export : bool + If True, export prealigned images to a MoBIE project. + dataset_name : str + Name of the MoBIE dataset (used only if `mobie_export=True`). + + Outputs + ------- + - Plots of slices and overlays before and after prealignment, saved in + ``{output_dir}/plots/``. + - Prealigned fixed and moving volumes saved as N5 containers in + ``{output_dir}/``. + - Transformation matrix for the moving image saved as a text file. + - (Optional) Exported MoBIE project with updated views. + + Notes + ----- + - The function assumes the input volumes are large 3D datasets. + - Axis orientation is checked via intensity profile analysis; axes may be + flipped by 180° if misaligned. + - The MoBIE export modifies the `dataset.json` to set the prealigned fixed + volume as the default view. + """ if not os.path.exists(f"{output_dir}/plots"): os.makedirs(f"{output_dir}/plots") @@ -286,7 +339,24 @@ def run_prealignment( @click.option("-o", "--output_dir", required=True, help="Output directory") @click.option("-m", "--mobie_export", required=False, is_flag=True, help="MoBIE export") def main(fixed_path, fixed_key, moving_path, moving_key, output_dir, mobie_export): + """ + Perform prealignment of moving image to fixed image. + + This function orchestrates the sequence of steps required to prealign a moving image to a fixed image. + It handles the creation of necessary directories, configures logging, and invokes prealignment functions. + Optionally, it can create a MoBIE project for visualization. + Args: + fixed_path (str): Path to the fixed input .n5 file. + fixed_key (str): Key to the fixed image data in the .n5 file. + moving_path (str): Path to the moving input .n5 file. + moving_key (str): Key to the moving image data in the .n5 file. + output_dir (str): Directory where the results should be saved. + mobie_export (bool): Flag indicating whether to export results to a MoBIE project. + + Returns: + None + """ logging.basicConfig( level=logging.INFO, format="%(asctime)s [%(levelname)s] %(message)s",