Exp_data_reader.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import sys\n",
    "# sys.path.insert(0, '/Users/austin/Documents/GitHub/pyTEMlib')\n",
    "\n",
    "# sys.path.insert(0, '/Users/austin/Documents/GitHub/SciFiReaders')\n",
    "# import SciFiReaders\n",
    "import pyTEMlib\n",
    "from pyTEMlib import file_tools as ft\n",
    "from pyTEMlib import image_tools as it\n",
    "%matplotlib ipympl\n",
    "\n",
    "\n",
    "from scipy.spatial import KDTree\n",
    "import matplotlib.cm as cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pyTEMlib.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Graphene example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Graphene exp data\n",
    "path = '/Users/austin/Dropbox/graphene_images/Hao_Wang_colab/'\n",
    "files = os.listdir(path)\n",
    "files = [f for f in files if f.endswith('.emd')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dates = [file.split(' ')[2].split('_')[0] for file in files]\n",
    "dates = list(set(dates))\n",
    "\n",
    "date = dates[0]\n",
    "sorted_files = [file for file in files if '.emd' in file]\n",
    "\n",
    "# file indexes 2,8,14 on dates[4] are Moires"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = sorted_files[14]\n",
    "\n",
    "dset = ft.open_file(path + file)\n",
    "image = dset['Channel_000']\n",
    "viw = image.plot(cmap='gray')\n",
    "\n",
    "pixel_size = image.original_metadata['BinaryResult']['PixelSize']['width']\n",
    "pixel_size = np.array(pixel_size)\n",
    "print('Pixel size: ', pixel_size, )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im_array = np.array(image)\n",
    "\n",
    "np.savez('Moire.npz', im_array=im_array, pixel_size=pixel_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "power_spectrum = pyTEMlib.image_tools.power_spectrum(image, smoothing=1)\n",
    "\n",
    "power_spectrum.view_metadata()\n",
    "view = power_spectrum.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find spots in power spectrum\n",
    "# import blob finder\n",
    "from skimage.feature import blob_log\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find blobs\n",
    "blobs = blob_log(power_spectrum, min_sigma=1, max_sigma=10, num_sigma=10, threshold=0.5)\n",
    "\n",
    "# plot blobs\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(power_spectrum, cmap='gray')\n",
    "for blob in blobs:\n",
    "    y, x, r = blob\n",
    "    c = plt.Circle((x, y), r, color='red', linewidth=2, fill=False)\n",
    "    ax.add_patch(c)\n",
    "\n",
    "ax.set_xlim(power_spectrum.shape[1]//3, 2*power_spectrum.shape[1]//3)\n",
    "ax.set_ylim(power_spectrum.shape[0]//3, 2*power_spectrum.shape[0]//3)\n",
    "ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get inner blobs:\n",
    "inner_blobs = blob_log(power_spectrum, min_sigma=1, max_sigma=10, num_sigma=10, threshold=0.52)\n",
    "# filter by distance from center\n",
    "# center is blob closest to COM of the blobs\n",
    "center = np.array([power_spectrum.shape[0]//2, power_spectrum.shape[1]//2])\n",
    "distances = [np.linalg.norm(np.array([blob[0], blob[1]]) - center) for blob in inner_blobs]\n",
    "inner_blobs = [inner_blobs[i] for i in range(len(inner_blobs)) if distances[i] < 30]\n",
    "inner_blobs = [inner_blobs[i] for i in range(len(inner_blobs)) if distances[i] > 18]\n",
    "\n",
    "# get outer blobs\n",
    "outer_blobs = blob_log(power_spectrum, min_sigma=1, max_sigma=10, num_sigma=10, threshold=0.5)\n",
    "# filter by distance from center\n",
    "# center is blob closest to COM of the blobs\n",
    "center = np.array([power_spectrum.shape[0]//2, power_spectrum.shape[1]//2])\n",
    "distances = [np.linalg.norm(np.array([blob[0], blob[1]]) - center) for blob in outer_blobs]\n",
    "outer_blobs = [outer_blobs[i] for i in range(len(outer_blobs)) if distances[i] > 30]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_blobs = inner_blobs + outer_blobs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(power_spectrum, cmap='gray')\n",
    "for blob in inner_blobs:\n",
    "    y, x, r = blob\n",
    "    c = plt.Circle((x, y), r, color='red', linewidth=2, fill=False)\n",
    "    plt.gca().add_patch(c)\n",
    "\n",
    "for blob in outer_blobs:\n",
    "    y, x, r = blob\n",
    "    c = plt.Circle((x, y), r, color='blue', linewidth=2, fill=False)\n",
    "    plt.gca().add_patch(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get angles for all the blobs, for grouping them\n",
    "angles = [np.arctan2(blob[0] - center[0], blob[1] - center[1]) for blob in all_blobs]\n",
    "angles = np.array(angles) * 180 / np.pi + 180\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the remainder after dividing by 60\n",
    "angles = angles % 60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(angles, bins=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spots_A = [all_blobs[i] for i in range(len(all_blobs)) if angles[i] < 25]\n",
    "spots_B = [all_blobs[i] for i in range(len(all_blobs)) if angles[i] > 25]\n",
    "\n",
    "spots_A = np.array(spots_A)\n",
    "spots_B = np.array(spots_B)\n",
    "\n",
    "distances_A = [np.linalg.norm(np.array([blob[0], blob[1]]) - center) for blob in spots_A]\n",
    "distances_B = [np.linalg.norm(np.array([blob[0], blob[1]]) - center) for blob in spots_B]\n",
    "\n",
    "# this is the bit to make it checkerboard\n",
    "outer_A = [spots_A[i] for i in range(len(spots_A)) if distances_A[i] > 30]\n",
    "outer_B = [spots_B[i] for i in range(len(spots_B)) if distances_B[i] > 30]\n",
    "\n",
    "inner_A = [spots_B[i] for i in range(len(spots_B)) if distances_B[i] < 30]\n",
    "inner_B = [spots_A[i] for i in range(len(spots_A)) if distances_A[i] < 30]\n",
    "\n",
    "\n",
    "spots_A = np.array(inner_A + outer_A)\n",
    "spots_B = np.array(inner_B + outer_B)\n",
    "\n",
    "\n",
    "spots_A = spots_A[:, 0:2]\n",
    "spots_B = spots_B[:, 0:2]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(power_spectrum, cmap='gray')\n",
    "plt.scatter([spot[1] for spot in spots_A], [spot[0] for spot in spots_A], color='red')\n",
    "plt.scatter([spot[1] for spot in spots_B], [spot[0] for spot in spots_B], color='blue')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reciprocal_scale = np.array([ft.get_slope(power_spectrum.u.values), ft.get_slope(power_spectrum.v.values)])\n",
    "spots_A = spots_A * reciprocal_scale + [power_spectrum.u.values[0], power_spectrum.v.values[0]]\n",
    "spots_B = spots_B * reciprocal_scale + [power_spectrum.u.values[0], power_spectrum.v.values[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_A = pyTEMlib.image_tools.adaptive_fourier_filter(image, spots=spots_A, low_pass=1, reflection_radius=1)\n",
    "view = filtered_A.plot(cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_power_spectrum = pyTEMlib.image_tools.power_spectrum(filtered_A, smoothing=0)\n",
    "view = filtered_power_spectrum.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_B = pyTEMlib.image_tools.adaptive_fourier_filter(image, spots=spots_B, low_pass=1, reflection_radius=1)\n",
    "view = filtered_B.plot(cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in each of these, get the atoms and bond distances\n",
    "# get the atoms\n",
    "atoms_A = blob_log(filtered_A, min_sigma=4, max_sigma=10, num_sigma=10, threshold=0.9)\n",
    "atoms_A = np.array(atoms_A[:, 0:2])\n",
    "\n",
    "atoms_B = blob_log(filtered_B, min_sigma=4, max_sigma=10, num_sigma=10, threshold=0.9)\n",
    "atoms_B = np.array(atoms_B[:, 0:2])\n",
    "\n",
    "fig, ax = plt.subplots(1,2, sharex=True, sharey=True, figsize = (15,7))   \n",
    "ax[0].imshow(image, cmap='gray')\n",
    "ax[1].imshow(image, cmap='gray')\n",
    "ax[0].scatter([atom[1] for atom in atoms_A], [atom[0] for atom in atoms_A], color='red',s=10)\n",
    "ax[1].scatter([atom[1] for atom in atoms_B], [atom[0] for atom in atoms_B], color='blue',s=10)\n",
    "\n",
    "for a in ax:\n",
    "    a.axis('off')\n",
    "fig.tight_layout()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import KDTree\n",
    "from scipy.spatial import KDTree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the distance between the centroids\n",
    "dist_A = []\n",
    "tree = KDTree(atoms_A)\n",
    "distances, indices = tree.query(atoms_A, k=3)\n",
    "nearest_distances = distances[:, 1:] * float(pixel_size) * 1e10 # angstroms\n",
    "dist_A.append(nearest_distances.flatten())\n",
    "dist_A = np.array(dist_A).flatten()\n",
    "\n",
    "dist_B = []\n",
    "tree = KDTree(atoms_B)\n",
    "distances, indices = tree.query(atoms_B, k=3)\n",
    "nearest_distances = distances[:, 1:] * float(pixel_size) * 1e10 # angstroms\n",
    "dist_B.append(nearest_distances.flatten())\n",
    "dist_B = np.array(dist_B).flatten()\n",
    "\n",
    "dist_total = np.concatenate((dist_A, dist_B))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Median bond distance: ', np.median(dist_total))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(dist_total, bins=40, color='gray', alpha=0.5)\n",
    "plt.xlim(0.8,2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WSSe Exp data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# path = '/Users/austin/Desktop/Projects/WS2_twist/WS2_1deg/20230315/WS2_1deg/20230315 1045 STEM HAADF 9.20 Mx.emd'\n",
    "path = '/Users/austin/Desktop/Projects/WS2_twist/WS2_15deg/20230303/20230303 1703 STEM HAADF 4.60 Mx.emd'\n",
    "dset = ft.open_file(path)\n",
    "image = dset['Channel_000']\n",
    "\n",
    "# normalize the image\n",
    "image = image - image.min()\n",
    "image = image/image.max()\n",
    "im_array = np.array(image)\n",
    "\n",
    "\n",
    "view = image.plot()\n",
    "\n",
    "pixel_size = dset['Channel_000'].original_metadata['BinaryResult']['PixelSize']['width'] # m/pixel\n",
    "pixel_size = float(pixel_size) * 1e10 # Angstrom/pixel\n",
    "\n",
    "print('Pixel size: ', pixel_size, ' Angstroms/pixel')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# scipy ndimage zoom to get the right pixel size\n",
    "from scipy.ndimage import zoom\n",
    "\n",
    "zoom_factor = 2\n",
    "zoomed_im = zoom(im_array, zoom_factor, order=3)\n",
    "pixel_size = pixel_size / zoom_factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(zoomed_im, cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rigid_registered_dataset = it.rigid_registration(image)\n",
    "view = rigid_registered_dataset.plot()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im_array = np.array(rigid_registered_dataset).sum(axis=0)\n",
    "im_array.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# im_array = np.array(rigid_registered_dataset).sum(axis=0)\n",
    "\n",
    "# np.savez('WSSe_haadf_2.npz', im_array=zoomed_im, pixel_size=pixel_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drift = rigid_registered_dataset.metadata['drift']\n",
    "polynom_degree = 2 # 1 is linear fit, 2 is parabolic fit, ...\n",
    "\n",
    "x = np.linspace(0,drift.shape[0]-1,drift.shape[0])\n",
    "\n",
    "line_fit_x = np.polyfit(x, drift[:,0], polynom_degree)\n",
    "poly_x = np.poly1d(line_fit_x)\n",
    "line_fit_y = np.polyfit(x, drift[:,1], polynom_degree)\n",
    "poly_y = np.poly1d(line_fit_y)\n",
    "\n",
    "plt.figure()\n",
    "plt.axhline(color = 'gray')\n",
    "plt.plot(x, drift[:,0], label = 'drift x')\n",
    "plt.plot(x, drift[:,1], label = 'drift y')\n",
    "plt.plot(x, poly_x(x),  label = 'fit_drift_x')\n",
    "plt.plot(x, poly_y(x),  label = 'fit_drift_y')\n",
    "\n",
    "plt.legend();\n",
    "ax_pixels = plt.gca()\n",
    "ax_pixels.step(1, 1)\n",
    "\n",
    "scaleX = (rigid_registered_dataset.x[1]-rigid_registered_dataset.x[0])*1000.  #in pm\n",
    "\n",
    "ax_pm = ax_pixels.twinx()\n",
    "x_1, x_2 = ax_pixels.get_ylim()\n",
    "\n",
    "ax_pm.set_ylim(x_1*scaleX, x_2*scaleX)\n",
    "\n",
    "ax_pixels.set_ylabel('drift [pixels]')\n",
    "ax_pm.set_ylabel('drift [pm]')\n",
    "ax_pixels.set_xlabel('image number');\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# non_rigid_registered = it.demon_registration(rigid_registered_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## get some more images of WSSe\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WSSe Exp data \n",
    "path = '/Users/austin/Desktop/Projects/WS2_twist/WS2_1deg/20230315/WS2_1deg/20230315 1045 STEM HAADF 9.20 Mx.emd'\n",
    "dset = ft.open_file(path)\n",
    "image = dset['Channel_000']\n",
    "\n",
    "# normalize the image\n",
    "image = image - image.min()\n",
    "image = image/image.max()\n",
    "im_array = np.array(image)\n",
    "\n",
    "view = image.plot()\n",
    "\n",
    "pixel_size = dset['Channel_000'].original_metadata['BinaryResult']['PixelSize']['width'] # m/pixel\n",
    "pixel_size = float(pixel_size) * 1e10 # Angstrom/pixel\n",
    "\n",
    "print('Pixel size: ', pixel_size, ' Angstroms/pixel')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rigid_registered_dataset = it.rigid_registration(image)\n",
    "\n",
    "sum_im = rigid_registered_dataset.sum(axis=0)\n",
    "\n",
    "# np.savez('WSSe_haadf.npz', im_array=im_array, pixel_size=pixel_size)\n",
    "\n",
    "view = rigid_registered_dataset.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_im.data_type = 'image'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "power_spectrum = pyTEMlib.image_tools.power_spectrum(sum_im, smoothing=1)\n",
    "\n",
    "power_spectrum.view_metadata()\n",
    "view = power_spectrum.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.feature import blob_log\n",
    "im_array = np.array(sum_im)\n",
    "\n",
    "# find blobs\n",
    "blobs = blob_log(power_spectrum, min_sigma=1, max_sigma=10, num_sigma=10, threshold=0.5)\n",
    "\n",
    "# plot blobs\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(power_spectrum, cmap='gray')\n",
    "for blob in blobs:\n",
    "    y, x, r = blob\n",
    "    c = plt.Circle((x, y), r, color='red', linewidth=2, fill=False)\n",
    "    ax.add_patch(c)\n",
    "\n",
    "ax.set_xlim(power_spectrum.shape[1]//3, 2*power_spectrum.shape[1]//3)\n",
    "ax.set_ylim(power_spectrum.shape[0]//3, 2*power_spectrum.shape[0]//3)\n",
    "ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get inner blobs:\n",
    "inner_blobs = blob_log(power_spectrum, min_sigma=1, max_sigma=10, num_sigma=10, threshold=0.5)\n",
    "# filter by distance from center\n",
    "# center is blob closest to COM of the blobs\n",
    "center = np.array([power_spectrum.shape[0]//2, power_spectrum.shape[1]//2])\n",
    "distances = [np.linalg.norm(np.array([blob[0], blob[1]]) - center) for blob in inner_blobs]\n",
    "inner_blobs = [inner_blobs[i] for i in range(len(inner_blobs)) if distances[i] < 40]\n",
    "inner_blobs = [inner_blobs[i] for i in range(len(inner_blobs)) if distances[i] > 0]\n",
    "\n",
    "# get outer blobs\n",
    "outer_blobs = blob_log(power_spectrum, min_sigma=1, max_sigma=10, num_sigma=10, threshold=0.5)\n",
    "# filter by distance from center\n",
    "# center is blob closest to COM of the blobs\n",
    "center = np.array([power_spectrum.shape[0]//2, power_spectrum.shape[1]//2])\n",
    "distances = [np.linalg.norm(np.array([blob[0], blob[1]]) - center) for blob in outer_blobs]\n",
    "outer_blobs = [outer_blobs[i] for i in range(len(outer_blobs)) if distances[i] > 40 and distances[i] < 80]\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_blobs = inner_blobs + outer_blobs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(power_spectrum, cmap='gray')\n",
    "for blob in inner_blobs:\n",
    "    y, x, r = blob\n",
    "    c = plt.Circle((x, y), r, color='red', linewidth=2, fill=False)\n",
    "    plt.gca().add_patch(c)\n",
    "\n",
    "\n",
    "for blob in outer_blobs:\n",
    "    y, x, r = blob\n",
    "    c = plt.Circle((x, y), r, color='blue', linewidth=2, fill=False)\n",
    "    plt.gca().add_patch(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the nearest neighbor for each blob\n",
    "tree = KDTree(all_blobs)\n",
    "distances, indices = tree.query(all_blobs, k=2)  # k=2 to include the point itself as the closest\n",
    "pairs = {(min(i, indices[i][1]), max(i, indices[i][1])) for i in range(len(all_blobs))}\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "colors = cm.rainbow(np.linspace(0, 1, len(pairs)))\n",
    "for idx, pair in enumerate(pairs):\n",
    "    plt.scatter([all_blobs[pair[0]][1]], [all_blobs[pair[0]][0]], color=colors[idx])\n",
    "    plt.scatter([all_blobs[pair[1]][1]], [all_blobs[pair[1]][0]], color=colors[idx])\n",
    "plt.scatter(center[1], center[0], color='black', marker='x')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smaller_angle_coords = []\n",
    "larger_angle_coords = []\n",
    "\n",
    "for pair in pairs:\n",
    "    angle1 = np.arctan2(all_blobs[pair[0]][0] - center[0], all_blobs[pair[0]][1] - center[1]) * 180 / np.pi\n",
    "    angle2 = np.arctan2(all_blobs[pair[1]][0] - center[0], all_blobs[pair[1]][1] - center[1]) * 180 / np.pi\n",
    "\n",
    "    if angle1 < 0:\n",
    "        angle1 += 360\n",
    "    if angle2 < 0:\n",
    "        angle2 += 360\n",
    "\n",
    "    # Compare angles to determine which is smaller in the circular sense\n",
    "    angle_diff = (angle2 - angle1 + 360) % 360\n",
    "\n",
    "    avg_dist_from_center = (np.linalg.norm(np.array(all_blobs[pair[0]])[:2] - center) + np.linalg.norm(np.array(all_blobs[pair[1]])[:2] - center)) / 2\n",
    "\n",
    "    if avg_dist_from_center < 40:\n",
    "        if angle_diff < 180:\n",
    "            smaller_angle_coords.append(all_blobs[pair[0]])\n",
    "            larger_angle_coords.append(all_blobs[pair[1]])\n",
    "        else:\n",
    "            smaller_angle_coords.append(all_blobs[pair[1]])\n",
    "            larger_angle_coords.append(all_blobs[pair[0]])\n",
    "\n",
    "    else:\n",
    "        if angle_diff > 180:\n",
    "            smaller_angle_coords.append(all_blobs[pair[0]])\n",
    "            larger_angle_coords.append(all_blobs[pair[1]])\n",
    "        else:\n",
    "            smaller_angle_coords.append(all_blobs[pair[1]])\n",
    "            larger_angle_coords.append(all_blobs[pair[0]])\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(power_spectrum, cmap='gray')\n",
    "\n",
    "for blob in smaller_angle_coords:\n",
    "    y, x, r = blob\n",
    "    c = plt.Circle((x, y), r, color='red', linewidth=2, fill=False)\n",
    "    plt.gca().add_patch(c)\n",
    "\n",
    "for blob in larger_angle_coords:\n",
    "    y, x, r = blob\n",
    "    c = plt.Circle((x, y), r, color='blue', linewidth=2, fill=False)\n",
    "    plt.gca().add_patch(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spots_A = np.array(smaller_angle_coords)[:,:2]\n",
    "spots_B = np.array(larger_angle_coords)[:,:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reciprocal_scale = np.array([ft.get_slope(power_spectrum.u.values), ft.get_slope(power_spectrum.v.values)])\n",
    "spots_A = spots_A * reciprocal_scale + [power_spectrum.u.values[0], power_spectrum.v.values[0]]\n",
    "spots_B = spots_B * reciprocal_scale + [power_spectrum.u.values[0], power_spectrum.v.values[0]]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_A = pyTEMlib.image_tools.adaptive_fourier_filter(sum_im, spots=spots_A, low_pass=1, reflection_radius=1)\n",
    "view = filtered_A.plot(cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_power_spectrum = pyTEMlib.image_tools.power_spectrum(filtered_A, smoothing=0)\n",
    "view = filtered_power_spectrum.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_B = pyTEMlib.image_tools.adaptive_fourier_filter(sum_im, spots=spots_B, low_pass=1, reflection_radius=1)\n",
    "view = filtered_B.plot(cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_power_spectrum = pyTEMlib.image_tools.power_spectrum(filtered_B, smoothing=0)\n",
    "view = filtered_power_spectrum.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in each of these, get the atoms and bond distances\n",
    "# get the atoms\n",
    "atoms_A = blob_log(filtered_A, min_sigma=4, max_sigma=10, num_sigma=10, threshold=0.01)\n",
    "atoms_A = np.array(atoms_A[:, 0:2])\n",
    "\n",
    "atoms_B = blob_log(filtered_B, min_sigma=4, max_sigma=10, num_sigma=10, threshold=0.08)\n",
    "atoms_B = np.array(atoms_B[:, 0:2])\n",
    "\n",
    "fig, ax = plt.subplots(1,2, sharex=True, sharey=True, figsize = (15,7))   \n",
    "ax[0].imshow(filtered_A, cmap='gray')\n",
    "ax[1].imshow(filtered_B, cmap='gray')\n",
    "ax[0].scatter([atom[1] for atom in atoms_A], [atom[0] for atom in atoms_A], color='red',s=10)\n",
    "ax[1].scatter([atom[1] for atom in atoms_B], [atom[0] for atom in atoms_B], color='blue',s=10)\n",
    "\n",
    "for a in ax:\n",
    "    a.axis('off')\n",
    "fig.tight_layout()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.draw import disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im = np.array(sum_im)\n",
    "# Radius of the circle\n",
    "r = 3\n",
    "\n",
    "# Create an empty list to store the summed intensities\n",
    "intensities = []\n",
    "\n",
    "# Set up plot for visualizing\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(im, cmap='gray')\n",
    "\n",
    "# Iterate through each point\n",
    "for point in atoms_A:\n",
    "    row, col = point\n",
    "    # Create a circular mask centered at the point\n",
    "    rr, cc = disk((row, col), r, shape=im.shape)\n",
    "    \n",
    "    # Sum the intensity within the circle\n",
    "    intensity = np.sum(im[rr, cc])\n",
    "    intensities.append(intensity)\n",
    "    \n",
    "    # Visualize the circle on the image\n",
    "    circle = plt.Circle((col, row), r, color='red', fill=False)\n",
    "    ax.add_patch(circle)\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(intensities, bins=20, alpha=0.5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im = np.array(sum_im)\n",
    "# Radius of the circle\n",
    "r = 3\n",
    "\n",
    "# Create an empty list to store the summed intensities\n",
    "intensities = []\n",
    "\n",
    "# Set up plot for visualizing\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(im, cmap='gray')\n",
    "\n",
    "# Iterate through each point\n",
    "for point in atoms_B:\n",
    "    row, col = point\n",
    "    # Create a circular mask centered at the point\n",
    "    rr, cc = disk((row, col), r, shape=im.shape)\n",
    "    \n",
    "    # Sum the intensity within the circle\n",
    "    intensity = np.sum(im[rr, cc])\n",
    "    intensities.append(intensity)\n",
    "    \n",
    "    # Visualize the circle on the image\n",
    "    circle = plt.Circle((col, row), r, color='b', fill=False)\n",
    "    ax.add_patch(circle)\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(intensities, bins=20, alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im = np.array(filtered_A)\n",
    "# Radius of the circle\n",
    "r = 3\n",
    "\n",
    "# Create an empty list to store the summed intensities\n",
    "intensities = []\n",
    "\n",
    "# Set up plot for visualizing\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(im, cmap='gray')\n",
    "\n",
    "# Iterate through each point\n",
    "for point in atoms_A:\n",
    "    row, col = point\n",
    "    # Create a circular mask centered at the point\n",
    "    rr, cc = disk((row, col), r, shape=im.shape)\n",
    "    \n",
    "    # Sum the intensity within the circle\n",
    "    intensity = np.sum(im[rr, cc])\n",
    "    intensities.append(intensity)\n",
    "    \n",
    "    # Visualize the circle on the image\n",
    "    circle = plt.Circle((col, row), r, color='r', fill=False)\n",
    "    ax.add_patch(circle)\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(intensities, bins=50, alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im = np.array(filtered_B)\n",
    "# Radius of the circle\n",
    "r = 3\n",
    "\n",
    "# Create an empty list to store the summed intensities\n",
    "intensities = []\n",
    "\n",
    "# Set up plot for visualizing\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(im, cmap='gray')\n",
    "\n",
    "# Iterate through each point\n",
    "for point in atoms_B:\n",
    "    row, col = point\n",
    "    # Create a circular mask centered at the point\n",
    "    rr, cc = disk((row, col), r, shape=im.shape)\n",
    "    \n",
    "    # Sum the intensity within the circle\n",
    "    intensity = np.sum(im[rr, cc])\n",
    "    intensities.append(intensity)\n",
    "    \n",
    "    # Visualize the circle on the image\n",
    "    circle = plt.Circle((col, row), r, color='r', fill=False)\n",
    "    ax.add_patch(circle)\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(intensities, bins=50, alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pyTEMlib",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}