notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3DCT Refactor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "# markers_3d: correlation points in 3D image (FM)\n",
    "# markers_2d: correlation points in 2D image (FIB/SEM)\n",
    "# spots_3d: points of interest in 3D image (FM)\n",
    "# rotation_center: rotation center of the 3D image (FM)\n",
    "    ## The default value is set as the center of a cube with an edge length equal to the longest edge of the image volume\n",
    "# results_file: path to save the results (.txt)\n",
    "# imageProps: dictionary with image properties (voxel size, image size, etc.)\n",
    "\n",
    "# installer napari: req \n",
    "# #conda install -c conda-forge libstdcxx-ng\n",
    "# pip install -e . \n",
    "\n",
    "\n",
    "from tdct import correlation\n",
    "from pprint import pprint\n",
    "import os\n",
    "import csv\n",
    "import tifffile as tff\n",
    "import numpy as np\n",
    "\n",
    "PATH = \"/home/patrick/development/data/CORRELATION/3dct/3D_correlation_test_dataset\"\n",
    "PATH = \"/home/patrick/github/3DCT/3D_correlation_test_dataset\"\n",
    "\n",
    "\n",
    "fib_coord_filename = os.path.join(PATH, \"IB_coordinates.txt\")\n",
    "fm_coord_filename = os.path.join(PATH, \"LM_coordinates_withPOI.txt\")\n",
    "\n",
    "fib_image_filename = os.path.join(PATH, \"IB_image.tif\")\n",
    "fm_image_filename = os.path.join(PATH, \"LM_image_stack_reslized.tif\")\n",
    "\n",
    "def parse_coordinate_file(filename: str, delimiter: str = \"\\t\") -> list:\n",
    "    coords: list = []\n",
    "    with open(filename) as csv_file:\n",
    "        for row in csv.reader(csv_file, delimiter=delimiter):\n",
    "            coords.append([field for field in row])\n",
    "    return coords\n",
    "\n",
    "fib_coordinates = parse_coordinate_file(fib_coord_filename)\n",
    "fm_coordinates = parse_coordinate_file(fm_coord_filename)\n",
    "\n",
    "fib_coordinates = np.array(fib_coordinates, dtype=np.float32)\n",
    "fm_coordinates = np.array(fm_coordinates, dtype=np.float32)\n",
    "\n",
    "print(f\"FIB Coordinates:\")\n",
    "pprint(fib_coordinates)\n",
    "print(f\"FM Coordinates:\")\n",
    "pprint(fm_coordinates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract the number of coordinates for each type (fib, fm, poi)\n",
    "\n",
    "nrows_fib = len(fib_coordinates)\n",
    "nrows_fm = len(fm_coordinates)\n",
    "nrows_poi = nrows_fm - nrows_fib\n",
    "\n",
    "print(f\"Number of FIB coordinates: {nrows_fib}\")\n",
    "print(f\"Number of FM coordinates: {nrows_fm}\")\n",
    "print(f\"Number of POI coordinates: {nrows_poi}\")\n",
    "\n",
    "# points of interest are excess coordinates in the FM coordinate file\n",
    "fib_coordinates_corr = fib_coordinates\n",
    "fm_coordinates_corr = fm_coordinates[:nrows_fib]\n",
    "poi_coordinates = fm_coordinates[nrows_fib:]\n",
    "\n",
    "print(fib_coordinates_corr)\n",
    "print(fm_coordinates_corr)\n",
    "print(poi_coordinates)\n",
    "\n",
    "# points are selected in the image coordinates (0, 0 top left)\n",
    "# z level is in pixels as well\n",
    "# no need to convert to physical coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rotation center is the center of the image volume\n",
    "\n",
    "fm_image = tff.imread(fm_image_filename)\n",
    "fib_image = tff.imread(fib_image_filename)\n",
    "\n",
    "print(fm_image.shape)\n",
    "print(fib_image.shape)\n",
    "\n",
    "halfmax_dim = max(fm_image.shape) * 0.5\n",
    "rotation_center = (halfmax_dim, halfmax_dim, halfmax_dim)\n",
    "print(rotation_center)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract fib pixel size from tiff metadata\n",
    "# open tiff file\n",
    "with tff.TiffReader(fib_image_filename) as tif:\n",
    "    for page in tif.pages:\n",
    "        for tag in page.tags.values():\n",
    "            if tag.name == 'FEI_HELIOS':\n",
    "                md = tag.value\n",
    "\n",
    "pprint(md[\"Scan\"][\"PixelWidth\"])    \n",
    "pixel_size = md[\"Scan\"][\"PixelWidth\"]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# image properties: used for converting coordinates to physical units (microscope coordinates)\n",
    "\n",
    "print(f\"FM Shape: {fm_image.shape}\")\n",
    "print(f\"Initial FIB Shape: {fib_image.shape}\")\n",
    "\n",
    "if fib_image.ndim == 3:\n",
    "    fib_shape = fib_image.shape[:-1]\n",
    "    print(f\"2D FIB Shape: {fib_shape}\")\n",
    "\n",
    "# standard resolutions: (old), needs to be modernised\n",
    "standard_res = [(442, 512), (884, 1024), (1768, 2048), (3536, 4096), (7072, 8192)]\n",
    "\n",
    "# find closest resolution\n",
    "res = min(standard_res, key=lambda x: abs(x[0] - fib_shape[0]))\n",
    "\n",
    "if fib_shape != res:\n",
    "    print(f\"Clipped to closest standard resolution: {res}\")\n",
    "\n",
    "# fib image shape minus metadata, fib_pixelsize (microns), fm_image_shape\n",
    "image_props = [res, pixel_size*1e6, fm_image.shape]\n",
    "print(f\"Image Properties: {image_props}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fib_coordinates_corr[:, :2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "\n",
    "timestamp = datetime.datetime.now().strftime('%Y-%m-%d_%H-%M-%S')\n",
    "results_file = os.path.join(PATH, f\"{timestamp}_correlation.txt\")\n",
    "\n",
    "cor_ret = ret = correlation.main(markers_3d=fm_coordinates_corr,\n",
    "                                markers_2d=fib_coordinates_corr,\n",
    "                                spots_3d=poi_coordinates,\n",
    "                                rotation_center=rotation_center,\n",
    "                                results_file=results_file,\n",
    "                                imageProps=image_props\n",
    "                                )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# transf, transf3d, spots2d, deltad2d, cm_3d_markers, mod_trans\n",
    "\n",
    "from tdct.app import save_correlation_data, parse_correlation_result\n",
    "\n",
    "# input data\n",
    "input_data = {\n",
    "    \"fib_coordinates\": fib_coordinates_corr.tolist(),\n",
    "    \"fm_coordinates\": fm_coordinates_corr.tolist(),\n",
    "    \"poi_coordinates\": poi_coordinates.tolist(),\n",
    "    \"image_properties\": {\n",
    "        \"fib_image_filename\": fib_image_filename,\n",
    "        \"fib_image_shape\": list(image_props[0]),\n",
    "        \"fib_pixel_size_um\": float(image_props[1]),\n",
    "        \"fm_image_filename\": fm_image_filename,\n",
    "        \"fm_image_shape\": list(fm_image.shape),\n",
    "    },\n",
    "    \"rotation_center\": list(rotation_center),\n",
    "    \"rotation_center_custom\": list(rotation_center)\n",
    "}\n",
    "print(input_data[\"image_properties\"][\"fib_image_shape\"])\n",
    "# output data\n",
    "correlation_data = parse_correlation_result(cor_ret, input_data)\n",
    "\n",
    "# full correlation data\n",
    "full_correlation_data = {\n",
    "    \"metadata\": {\n",
    "        \"timestamp\": datetime.datetime.now().strftime('%Y-%m-%d_%H-%M-%S'),\n",
    "        \"data_path\": PATH,\n",
    "        \"csv_path\": \"\", \n",
    "        \"project_path\": \"\",\n",
    "    },\n",
    "    \"correlation\": correlation_data,\n",
    "\n",
    "}\n",
    "save_correlation_data(full_correlation_data, PATH)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cor_ret[0] # transformation matrix\n",
    "# cor_ret[1] # reprojected 3D points\n",
    "# cor_ret[2] # spots2d\n",
    "# cor_ret[3] # difference reprojected 3D points and 2D points\n",
    "# cor_ret[4] # cm_3d_markers\n",
    "# cor_ret[5] # translation around rotation center (custom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(cor_ret[1].T)\n",
    "\n",
    "print(fib_coordinates_corr[:, :2] + cor_ret[3].T)\n",
    "\n",
    "# convert from q list to matrix\n",
    "# q = np.array(full_correlation_data[\"rotation_quaternion\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error = full_correlation_data[\"correlation\"][\"output\"][\"error\"]\n",
    "\n",
    "delta_2d = np.array(error[\"delta_2d\"]).T\n",
    "\n",
    "print(delta_2d.shape)\n",
    "# loop through all points and and print (x, y)\n",
    "for dx, dy in delta_2d:\n",
    "    \n",
    "    # fmt as .2f str\n",
    "    dx, dy = f\"{dx:.2f}\", f\"{dy:.2f}\"\n",
    "\n",
    "    print(f\"({dx}, {dy})\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pixel_size = tif.pages[0].tags['FEI_HELIOS'].value['Scan']['PixelWidth']\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import tifffile as tff\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "import numpy as np\n",
    "import os\n",
    "FIB_IMAGE_PATH = \"3D_correlation_test_dataset/IB_image.tif\"\n",
    "\n",
    "\n",
    "# things to consider: \n",
    "# extract pixel size\n",
    "# convert to grayscale (if rgb)\n",
    "# trim metadata bar\n",
    "import logging\n",
    "\n",
    "logging.basicConfig(level=logging.INFO)\n",
    "from tdct.app import load_and_parse_fib_image\n",
    "\n",
    "fib_image, pixel_size = load_and_parse_fib_image(FIB_IMAGE_PATH)\n",
    "\n",
    "plt.imshow(fib_image, cmap='gray')\n",
    "plt.show()\n",
    "\n",
    "print(f\"Pixel Size: {pixel_size*1e6}\")\n",
    "\n",
    "# fig, ax = plt.subplots(1, 3, figsize=(15, 3))\n",
    "\n",
    "# plt.suptitle(f\"Image: {os.path.basename(FIB_IMAGE_PATH)}, Pixel Size: {pixel_size*1e9:.2f} nm\")\n",
    "\n",
    "# ax[0].imshow(fib_image, cmap='gray')\n",
    "# ax[0].set_title(f\"Original: {fib_image.shape}\")\n",
    "# ax[1].imshow(img, cmap='gray')\n",
    "# ax[1].set_title(f\"Base: {img.shape}\")\n",
    "# ax[2].imshow(trimmed_img, cmap='gray')\n",
    "# ax[2].set_title(f\"Trimmed: {trimmed_img.shape}\")\n",
    "# plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Refactor Testing\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import os\n",
    "import tifffile as tff\n",
    "import matplotlib.pyplot as plt\n",
    "import datetime\n",
    "import numpy as np\n",
    "from tdct import correlation\n",
    "from pprint import pprint\n",
    "from tdct.app import load_and_parse_fib_image, parse_coordinates\n",
    "\n",
    "\n",
    "PATH = \"/home/patrick/github/3DCT/3D_correlation_test_dataset\"\n",
    "\n",
    "\n",
    "fib_coord_filename = os.path.join(PATH, \"IB_coordinates.txt\")\n",
    "fm_coord_filename = os.path.join(PATH, \"LM_coordinates_withPOI.txt\")\n",
    "\n",
    "fib_image_filename = os.path.join(PATH, \"IB_image.tif\")\n",
    "fm_image_filename = os.path.join(PATH, \"LM_image_stack_reslized.tif\")\n",
    "\n",
    "\n",
    "fib_coordinates, fm_coordinates = parse_coordinates(fib_coord_filename, fm_coord_filename)\n",
    "\n",
    "# extract the number of coordinates for each type (fib, fm, poi)\n",
    "\n",
    "nrows_fib = len(fib_coordinates)\n",
    "nrows_fm = len(fm_coordinates)\n",
    "nrows_poi = nrows_fm - nrows_fib\n",
    "\n",
    "print(f\"Number of FIB coordinates: {nrows_fib}\")\n",
    "print(f\"Number of FM coordinates: {nrows_fm}\")\n",
    "print(f\"Number of POI coordinates: {nrows_poi}\")\n",
    "\n",
    "# points of interest are excess coordinates in the FM coordinate file\n",
    "fib_coordinates = fib_coordinates\n",
    "poi_coordinates = fm_coordinates[nrows_fib:]\n",
    "fm_coordinates = fm_coordinates[:nrows_fib]\n",
    "\n",
    "\n",
    "\n",
    "print(\"-\"*50)\n",
    "print(f\"FIB Coordinates:\")\n",
    "pprint(fib_coordinates)\n",
    "print(f\"FM Coordinates:\")\n",
    "pprint(fm_coordinates)\n",
    "print(f\"POI Coordinates:\")\n",
    "pprint(poi_coordinates)\n",
    "\n",
    "fib_image, pixel_size = load_and_parse_fib_image(fib_image_filename)\n",
    "fm_image = tff.imread(fm_image_filename)\n",
    "\n",
    "print(\"-\"*50)\n",
    "print(f\"FIB Image Shape: {fib_image.shape}\")\n",
    "print(f\"Pixel Size: {pixel_size*1e6}um\")\n",
    "print(f\"FM Shape: {fm_image.shape}\")\n",
    "\n",
    "print(\"-\"*50)\n",
    "\n",
    "# rotation center\n",
    "halfmax_dim = max(fm_image.shape) * 0.5\n",
    "rotation_center = (halfmax_dim, halfmax_dim, halfmax_dim)\n",
    "print(f\"Rotation Center: {rotation_center}\")\n",
    "\n",
    "\n",
    "# set image properties\n",
    "image_props = [fib_image.shape, pixel_size*1e6, fm_image.shape]\n",
    "print(f\"Image Properties: {image_props}\")\n",
    "\n",
    "print(\"-\"*50)\n",
    "\n",
    "timestamp = datetime.datetime.now().strftime('%Y-%m-%d_%H-%M-%S')\n",
    "results_file = os.path.join(PATH, f\"{timestamp}_correlation.txt\")\n",
    "print(f\"Results File: {results_file}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cor_ret = correlation.main(markers_3d=fm_coordinates,\n",
    "                                markers_2d=fib_coordinates,\n",
    "                                spots_3d=poi_coordinates,\n",
    "                                rotation_center=rotation_center,\n",
    "                                results_file=results_file,\n",
    "                                imageProps=image_props\n",
    "                                )\n",
    "\n",
    "cor_ret_2 = correlation.correlate(markers_3d=fm_coordinates,\n",
    "                                markers_2d=fib_coordinates,\n",
    "                                poi_3d=poi_coordinates,\n",
    "                                rotation_center=rotation_center,\n",
    "                                imageProps=image_props\n",
    "                                )\n",
    "\n",
    "# correlation.save_results(cor_ret_2, results_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_data = {\n",
    "\"fib_coordinates\": fib_coordinates.tolist(),\n",
    "\"fm_coordinates\": fm_coordinates.tolist(),\n",
    "\"poi_coordinates\": poi_coordinates.tolist(),\n",
    "\"image_properties\": {\n",
    "    \"fib_image_filename\": fib_image_filename,\n",
    "    \"fib_image_shape\": list(image_props[0]),\n",
    "    \"fib_pixel_size_um\": float(image_props[1]),\n",
    "    \"fm_image_filename\": fm_image_filename,\n",
    "    \"fm_image_shape\": list(image_props[2]),\n",
    "},\n",
    "\"rotation_center\": list(rotation_center),\n",
    "\"rotation_center_custom\": list(rotation_center),\n",
    "}\n",
    "\n",
    "from tdct.app import parse_correlation_result\n",
    "\n",
    "# output data\n",
    "correlation_data = parse_correlation_result(\n",
    "    cor_ret, input_data=input_data\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tdct.io import parse_correlation_result_v2\n",
    "corr2 = parse_correlation_result_v2(cor_ret_2, input_data=input_data)\n",
    "\n",
    "pprint(corr2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pprint(correlation_data[\"output\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pre_cor_ret = cor_ret\n",
    "\n",
    "print(\"-\"*50)\n",
    "\n",
    "poi2d = cor_ret[2].T\n",
    "poi3d = cor_ret[1].T\n",
    "d2d = cor_ret[3].T\n",
    "mod_trans = cor_ret[5]\n",
    "\n",
    "cor_ret_2 = cor_ret_2[\"output\"]\n",
    "poi2d_2 = cor_ret_2[\"reprojected_2d_poi\"].T\n",
    "poi3d_2 = cor_ret_2[\"reprojected_3d_coordinates\"].T\n",
    "d2d_2 = cor_ret_2[\"reprojection_error\"].T\n",
    "mod_trans_2 = cor_ret_2[\"modified_translation\"]\n",
    "\n",
    "\n",
    "# assert the arrays are equal\n",
    "# np.testing.assert_array_equal(poi2d, poi2d_2)\n",
    "np.testing.assert_array_equal(poi3d, poi3d_2)\n",
    "assert np.allclose(d2d, d2d_2)\n",
    "assert np.allclose(mod_trans, mod_trans_2)\n",
    "\n",
    "print(\"Arrays are equal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility Refactor\n",
    "\n",
    "getzgauss\n",
    "\n",
    "getxyzopt\n",
    "\n",
    "interpolation\n",
    "\n",
    "multi-channel support"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# getzgauss\n",
    "\n",
    "import tifffile as tff\n",
    "import numpy as np\n",
    "import logging\n",
    "\n",
    "from tdct.beadPos import getzGauss\n",
    "from tdct.util import get_z_gauss, multi_channel_get_z_guass\n",
    "image = tff.imread(fm_image_filename)\n",
    "\n",
    "logging.basicConfig(level=logging.INFO)\n",
    "# get z position of the bead\n",
    "\n",
    "for i, coord in enumerate(fm_coordinates):\n",
    "    x, y, z_init = coord\n",
    "\n",
    "    z = getzGauss(x=x, y=y, img=image)  # threshVal, cutout\n",
    "    _, z2, _ = get_z_gauss(x=x, y=y, image=image, show=False)\n",
    "    _, mcz, _ = multi_channel_get_z_guass(x=x, y=y, image=image, show=False)\n",
    "\n",
    "    logging.info(f\"Initial Z: {z_init}, Calculated Z: {z}, Calculated Z2: {z2}, MCZ: {mcz}\")\n",
    "    assert np.isclose(z, z2, atol=1e-2)\n",
    "    assert np.isclose(z, mcz, atol=1e-2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compArray = np.array([\n",
    "        [[0, 0, 0, 0, 0],\n",
    "            [0, 0, 0, 0, 0],\n",
    "            [0, 0, 0, 0, 0],\n",
    "            [0, 0, 0, 0, 0],\n",
    "            [0, 0, 0, 0, 0]],\n",
    "        [[17, 17, 17, 17, 17],\n",
    "            [17, 17, 17, 17, 17],\n",
    "            [17, 17, 17, 17, 17],\n",
    "            [17, 17, 17, 17, 17],\n",
    "            [17, 17, 17, 17, 17]],\n",
    "        [[34, 34, 34, 34, 34],\n",
    "            [34, 34, 34, 34, 34],\n",
    "            [34, 34, 34, 34, 34],\n",
    "            [34, 34, 34, 34, 34],\n",
    "            [34, 34, 34, 34, 34]],\n",
    "        [[51, 51, 51, 51, 51],\n",
    "            [51, 51, 51, 51, 51],\n",
    "            [51, 51, 51, 51, 51],\n",
    "            [51, 51, 51, 51, 51],\n",
    "            [51, 51, 51, 51, 51]],\n",
    "        [[67, 67, 67, 67, 67],\n",
    "            [67, 67, 67, 67, 67],\n",
    "            [67, 67, 67, 67, 67],\n",
    "            [67, 67, 67, 67, 67],\n",
    "            [67, 67, 67, 67, 67]],\n",
    "        [[84, 84, 84, 84, 84],\n",
    "            [84, 84, 84, 84, 84],\n",
    "            [84, 84, 84, 84, 84],\n",
    "            [84, 84, 84, 84, 84],\n",
    "            [84, 84, 84, 84, 84]],\n",
    "        [[101, 101, 101, 101, 101],\n",
    "            [101, 101, 101, 101, 101],\n",
    "            [101, 101, 101, 101, 101],\n",
    "            [101, 101, 101, 101, 101],\n",
    "            [101, 101, 101, 101, 101]],\n",
    "        [[117, 117, 117, 117, 117],\n",
    "            [117, 117, 117, 117, 117],\n",
    "            [117, 117, 117, 117, 117],\n",
    "            [117, 117, 117, 117, 117],\n",
    "            [117, 117, 117, 117, 117]],\n",
    "        [[134, 134, 134, 134, 134],\n",
    "            [134, 134, 134, 134, 134],\n",
    "            [134, 134, 134, 134, 134],\n",
    "            [134, 134, 134, 134, 134],\n",
    "            [134, 134, 134, 134, 134]],\n",
    "        [[0, 0, 0, 0, 0],\n",
    "            [0, 0, 0, 0, 0],\n",
    "            [0, 0, 0, 0, 0],\n",
    "            [0, 0, 0, 0, 0],\n",
    "            [0, 0, 0, 0, 0]]], dtype=np.uint8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# interpolation\n",
    "from tdct import stackProcessing\n",
    "from tdct.util import interpolate_z_stack, showgraph_\n",
    "calcArray = np.ones([4,5,5],dtype=np.uint8)\n",
    "calcArray[1] += 50\n",
    "calcArray[2] += 100\n",
    "calcArray[3] += 150\n",
    "\n",
    "retArray = stackProcessing.interpol(calcArray, 300., 100., \"spline\", showgraph=False)\n",
    "new_array = interpolate_z_stack(calcArray, 300., 100., method=\"spline\") \n",
    "\n",
    "print(f\"Shapes: {retArray.shape}, {compArray.shape}, {new_array.shape}\")\n",
    "print(f\"Calc array same: {np.array_equal(retArray, compArray)}\")\n",
    "print(f\"New array equal: {np.array_equal(new_array, compArray)}\")\n",
    "\n",
    "\n",
    "# check why the arrays are not equal\n",
    "\n",
    "# for ret, comp in zip(retArray, compArray):\n",
    "#     print(ret)\n",
    "#     print(comp)\n",
    "#     print(np.array_equal(ret, comp))\n",
    "\n",
    "# TODO: linear first array is not correct, but spline is...??\n",
    "# TODO: add test for multi-channel squeezing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "showgraph_(calcArray, 300., 100.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(new_array.shape)\n",
    "print(interpolated.shape)\n",
    "print(fm_image.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z0 = 500e-9\n",
    "z1 = 77.4e-9\n",
    "print(z0 / z1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "# multi-channel images\n",
    "PATH = \"/home/patrick/github/3DCT/3D_correlation_test_dataset/test-image.ome.tiff\"\n",
    "import logging\n",
    "import tifffile as tff\n",
    "import numpy as np\n",
    "from tdct.util import interpolate_z_stack, multi_channel_interpolation\n",
    "\n",
    "logging.basicConfig(level=logging.INFO)\n",
    "\n",
    "px_in = 2.5e-7\n",
    "px_out = 3.89922e-8\n",
    "\n",
    "image = tff.imread(PATH)\n",
    "print(image.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "print(\"starting spline interpolation\")\n",
    "spline_inter = multi_channel_interpolation(image, px_in, px_out, method=\"spline\")\n",
    "\n",
    "print(\"starting fast cubic interpolation\")\n",
    "cubic_inter = multi_channel_interpolation(image, px_in, px_out, method=\"fast-cubic\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import napari\n",
    "\n",
    "viewer = napari.view_image(spline_inter, name=\"Spline Interpolated Image\", scale=(px_out, px_out, px_out))\n",
    "viewer.add_image(cubic_inter, scale=(px_out, px_out, px_out), name=\"Cubic Interpolated Image\")\n",
    "viewer.add_image(image, scale=(px_in, px_out, px_out), name=\"Original Image\")\n",
    "\n",
    "viewer.scale_bar.visible = True\n",
    "viewer.scale_bar.unit = \"m\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(cubic_inter.shape)\n",
    "print(spline_inter.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tdct.util import multi_channel_get_z_guass, get_z_gauss\n",
    "\n",
    "# get z position of the bead\n",
    "x=497\n",
    "y=253\n",
    "zvals = multi_channel_get_z_guass(image, x=x, y=y) \n",
    "zval = get_z_gauss(x=x, y=y, image=image[1], show=False)\n",
    "print(zvals)\n",
    "print(zval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = napari.view_image(image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zvals = multi_channel_get_z_guass(image[1], x=x, y=y) \n",
    "print(zvals)\n",
    "\n",
    "print(zvals[0].dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read Images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from tdct.app import load_and_parse_fib_image, load_image_and_metadata\n",
    "\n",
    "PATH = \"/home/patrick/github/3DCT/3D_correlation_test_dataset/fib_002.jpg\"\n",
    "\n",
    "# image, md = load_image_and_metadata(PATH)\n",
    "\n",
    "# print(md)\n",
    "# print (image.shape)\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import PIL.Image as Image\n",
    "\n",
    "\n",
    "image = Image.open(PATH)\n",
    "print(image.size)\n",
    "\n",
    "image = np.array(image)\n",
    "\n",
    "\n",
    "plt.imshow(image)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "fov = 50e-6\n",
    "width = 1536\n",
    "\n",
    "pixel_size = fov / width\n",
    "print(pixel_size)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save as tiff\n",
    "tff.imwrite(\"fib_002.tif\", image)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_and_parse_fib_image(\"/home/patrick/github/3DCT/3D_correlation_test_dataset/fib_002.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "3.2552083333333335e-08 * 1e6\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## XYZ Targeting\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import os\n",
    "import glob\n",
    "\n",
    "import matplotlib.pyplot as plt \n",
    "from tdct.util import multi_channel_get_z_guass, multi_channel_zyx_targeting\n",
    "from tdct.io import load_and_parse_fm_image\n",
    "\n",
    "PATH = \"/home/patrick/development/data/CORRELATION/3dct/3D_correlation_test_dataset\"\n",
    "filenames = glob.glob(os.path.join(PATH, \"*ome.tiff\"))\n",
    "\n",
    "print(filenames)\n",
    "\n",
    "image, md = load_and_parse_fm_image(filenames[0])\n",
    "print(image.shape)\n",
    "print(md)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image1 = image[1]\n",
    "xinit, yinit = 510, 255\n",
    "\n",
    "zval, zidx, zsigma = multi_channel_get_z_guass(image1, x=xinit, y=yinit)\n",
    "print(zval, zidx, zsigma)\n",
    "\n",
    "plt.imshow(image1[round(zidx)], cmap='gray')\n",
    "plt.plot(xinit, yinit, \"y+\", ms=50)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "xinit, yinit = 515, 235\n",
    "ch, (x,y,z) = multi_channel_zyx_targeting(image, xinit, yinit)\n",
    "\n",
    "print(ch, x, y, z)\n",
    "\n",
    "plt.imshow(image[ch, round(z)], cmap='gray')\n",
    "plt.plot(xinit, yinit, \"c+\", ms=50)\n",
    "plt.plot(x, y, \"y+\", ms=50)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "3dct-env",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}