{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# imports\n", "import sys\n", "import pathlib\n", "import scipy\n", "import skimage\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "registration_dir_path = pathlib.Path().resolve().parent / \"matlab_translation\"\n", "if registration_dir_path not in sys.path:\n", " sys.path.append(str(registration_dir_path))\n", "\n", "__package__ = \"registration\"\n", "\n", "from .preprocessing_functions import get_ims, preprocessing\n", "from .global_registration import (\n", " _rescale_and_blur_image,\n", " _normalize_image,\n", " _reg_ims_com,\n", " _apply_transformation,\n", " _get_com,\n", " _bandpass_img,\n", ")\n", "from .calculate_dislocation import _get_bright_spot_centroid" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "image_dir = pathlib.Path(\"/Users/margareteminizer/Desktop/coda/sample_dataset/20x files/1x\")\n", "ref_im_name_stem = \"lungs_013\"\n", "moving_im_name_stem = \"lungs_011\"\n", "image_suffix = \".tif\"\n", "\n", "moving_im_0, moving_mask = get_ims(image_dir, moving_im_name_stem, image_suffix)\n", "ref_im_0, ref_mask = get_ims(image_dir, ref_im_name_stem, image_suffix)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "max_size = [1306, 1187]\n", "pad_all = 250\n", "\n", "ref_im, ref_im_g, ref_mask, _ = preprocessing(\n", " ref_im_0, ref_mask, max_size, pad_all\n", " )\n", "moving_im, moving_im_g, moving_mask, fillvals = preprocessing(\n", " moving_im_0, moving_mask, max_size, pad_all\n", ")" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "rescale = 6\n", "a_ref = _rescale_and_blur_image(ref_im_g, rescale)\n", "a_moving = _rescale_and_blur_image(moving_im_g, rescale)\n", "a_ref = _normalize_image(a_ref)\n", "a_moving = _normalize_image(a_moving)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "OUTPUT_DIR = pathlib.Path().resolve() / \"group_of_reg_plots\"\n", "OUTPUT_DIR.mkdir(exist_ok=True)\n", "\n", "def make_frame_plot(ref_im, moving_im, current_angle, best_angle, theta_range, theta_step, iteration, iangle, clip_off=50):\n", " # center of mass of the reference image\n", " ref_com = _get_com(ref_im)\n", " # apply the current rotation and COM translation\n", " im_moving_current_rot = _apply_transformation(moving_im, current_angle)\n", " com_trans_xy_current = ref_com - _get_com(im_moving_current_rot)\n", " moving_im_current = _apply_transformation(im_moving_current_rot, translation_xy=com_trans_xy_current)\n", " # apply the best rotation and COM translation\n", " im_moving_best_rot = _apply_transformation(moving_im, best_angle)\n", " com_trans_xy_best = ref_com - _get_com(im_moving_best_rot)\n", " moving_im_best = _apply_transformation(im_moving_best_rot, translation_xy=com_trans_xy_best)\n", " # make the Radon transforms\n", " theta = np.arange(-1.0 * theta_range, theta_range + theta_step, theta_step)\n", " radon_moving_current = _bandpass_img(skimage.transform.radon(moving_im_current, theta, circle=True), 1, 3)\n", " radon_moving_current = radon_moving_current[20:-20, :]\n", " # clip off the blank edges\n", " ref_im_c = ref_im[clip_off:-clip_off, clip_off:-clip_off]\n", " moving_im_current_c = moving_im_current[clip_off:-clip_off, clip_off:-clip_off]\n", " moving_im_best_c = moving_im_best[clip_off:-clip_off, clip_off:-clip_off]\n", " # normalize\n", " ref_im_norm = ref_im_c / np.percentile(ref_im_c, 99.0)\n", " moving_im_current_norm = moving_im_current_c / np.percentile(moving_im_current_c, 99.0)\n", " moving_im_best_norm = moving_im_best_c / np.percentile(moving_im_best_c, 99.0)\n", " # make overlays\n", " overlay_current = np.clip(np.stack((ref_im_norm, moving_im_current_norm, 0.5*(ref_im_norm + moving_im_current_norm)), axis=-1), 0., 1.)\n", " overlay_best = np.clip(np.stack((ref_im_norm, moving_im_best_norm, 0.5*(ref_im_norm + moving_im_best_norm)), axis=-1), 0., 1.)\n", " # make the overlay plots\n", " f, ax = plt.subplots(1, 2, figsize=(2*6., 6.*ref_im.shape[0]/ref_im.shape[1]))\n", " ax[0].imshow(overlay_best)\n", " ax[0].set_title(f\"Best rotation ({best_angle:.2f}°)\")\n", " ax[1].imshow(overlay_current)\n", " ax[1].set_title(f\"Current rotation ({current_angle:.2f}°)\")\n", " for axs in ax:\n", " axs.axis(\"off\")\n", " f.suptitle(f\"Iteration {iteration}\")\n", " f.tight_layout()\n", " #plt.show()\n", " plt.savefig(OUTPUT_DIR / f\"overlays_iter_{iteration:01d}_angle_{iangle:01d}.png\", bbox_inches=\"tight\")\n", " plt.close()\n", " # make the plot of the radon transform\n", " f, ax = plt.subplots(figsize=(6.*radon_moving_current.shape[1]/radon_moving_current.shape[0], 6.))\n", " ax.imshow(radon_moving_current, cmap=\"gray\")\n", " min_angle = current_angle - theta_range\n", " max_angle = current_angle + theta_range\n", " ax.set_title(f\"Current Radon Transform\\n({min_angle:.1f} - {max_angle:.1f} deg, {theta_step:.1f} deg steps)\")\n", " ax.axis(\"off\")\n", " # plt.show()\n", " plt.savefig(OUTPUT_DIR / f\"radons_iter_{iteration:01d}_angle_{iangle:01d}.png\", bbox_inches=\"tight\")\n", " plt.close()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "n_iters = 2\n", "\n", "angles = [0.0, 180.0, 90.0, 270.0]\n", "step = 45.0\n", "best_metric_value = 0.2\n", "best_rotation_angle = 0.0\n", "best_trans_xy = np.array([0.0, 0.0])\n", "iteration = 1\n", "while step >= 1:\n", " theta_range = step * 3 if step / 2 >= 1 else 90\n", " theta_step = step / 10.0 if step / 2 >= 1 else 0.5\n", " for iangle, angle in enumerate(angles):\n", " #### plot stuff below ####\n", " make_frame_plot(\n", " a_ref, a_moving, angle, best_rotation_angle, theta_range, theta_step, iteration, iangle\n", " )\n", " #### plot stuff above ####\n", " registered_im, rot_angle, trans_xy, rr_metric = _reg_ims_com(\n", " a_ref,\n", " a_moving,\n", " n_iters,\n", " angle,\n", " do_center_of_mass_translation=True,\n", " theta_range=theta_range,\n", " theta_step=theta_step,\n", " )\n", " if rr_metric != 0:\n", " aa = (a_ref > 0).astype(np.uint8) + (registered_im > 0).astype(np.uint8)\n", " rr_metric = np.sum(aa == 2) / np.sum(aa > 0)\n", " if rr_metric > best_metric_value:\n", " best_metric_value = rr_metric\n", " best_rotation_angle = rot_angle\n", " best_trans_xy = trans_xy\n", " # if best_metric_value > metric_threshold:\n", " # break\n", " iteration+=1\n", " # Refine angles around the best angle found\n", " angles = [best_rotation_angle + step * i for i in (0, -1, 1, -2, 2)]\n", " step /= 2 # Reduce step size for finer search" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "angles = [0.0, 180.0, 90.0, 270.0]\n", "step = 45.0\n", "best_metric_value = 0.2\n", "best_rotation_angle = 0.0\n", "best_trans_xy = np.array([0.0, 0.0])\n", "iteration = 1\n", "while step >= 1:\n", " theta_range = step * 3 if step / 2 >= 1 else 90\n", " theta_step = step / 10.0 if step / 2 >= 1 else 0.5\n", " for iangle, angle in enumerate(angles):\n", " #### plot stuff below ####\n", " make_frame_plot(\n", " a_ref, a_moving, angle, best_rotation_angle, theta_range, theta_step, iteration, iangle\n", " )\n", " #### plot stuff above ####\n", " registered_im, rot_angle, trans_xy, rr_metric = _reg_ims_com(\n", " a_ref,\n", " a_moving,\n", " n_iters,\n", " angle,\n", " do_center_of_mass_translation=True,\n", " theta_range=theta_range,\n", " theta_step=theta_step,\n", " )\n", " if rr_metric != 0:\n", " aa = (a_ref > 0).astype(np.uint8) + (registered_im > 0).astype(np.uint8)\n", " rr_metric = np.sum(aa == 2) / np.sum(aa > 0)\n", " if rr_metric > best_metric_value:\n", " best_metric_value = rr_metric\n", " best_rotation_angle = rot_angle\n", " best_trans_xy = trans_xy\n", " # if best_metric_value > metric_threshold:\n", " # break\n", " iteration+=1\n", " # Refine angles around the best angle found\n", " angles = [best_rotation_angle + step * i for i in (0, -1, 1, -2, 2)]\n", " step /= 2 # Reduce step size for finer search" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "angles = (\n", " [-2, 177, 87, 268, -1, 88, 269, 178]\n", " + list(range(-7, 8, 2))\n", " + list(range(179, 184))\n", " + list(range(89, 94))\n", " + list(range(270, 273))\n", ")\n", "best_metric_value = 0.2\n", "best_rotation_angle = 0.0\n", "best_trans_xy = np.array([0.0, 0.0])\n", "iteration = 1\n", "theta_range = 90\n", "theta_step = 0.5\n", "for iangle, angle in enumerate(angles):\n", " #### plot stuff below ####\n", " make_frame_plot(\n", " a_ref, a_moving, angle, best_rotation_angle, theta_range, theta_step, iteration, iangle\n", " )\n", " #### plot stuff above ####\n", " registered_im, rot_angle, trans_xy, rr_metric = _reg_ims_com(\n", " a_ref,\n", " a_moving,\n", " n_iters,\n", " angle,\n", " do_center_of_mass_translation=True,\n", " theta_range=theta_range,\n", " theta_step=theta_step,\n", " )\n", " if rr_metric != 0:\n", " aa = (a_ref > 0).astype(np.uint8) + (registered_im > 0).astype(np.uint8)\n", " rr_metric = np.sum(aa == 2) / np.sum(aa > 0)\n", " if rr_metric > best_metric_value:\n", " best_metric_value = rr_metric\n", " best_rotation_angle = rot_angle\n", " best_trans_xy = trans_xy\n", " # if best_metric_value > metric_threshold:\n", " # break" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "coda_scratch", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.15" } }, "nbformat": 4, "nbformat_minor": 2 }