analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of a TEM image of a grain boundary\n",
    "* <a href=\"#download\">Downloading data</a>\n",
    "* <a href=\"#open\">Exploring the data</a>\n",
    "* <a href=\"#analysis\">Analysis of the dataset</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"download\"></a>1. Downloading the data if necessary and checking the hash\n",
    "This we do using bash as it's a bit easier than with python and could also be done in a shell outside the jupyter notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# create the data folder in the current folder if it doesn't exist\n",
    "if [ ! -d data ]; then\n",
    "   mkdir data\n",
    "fi\n",
    "\n",
    "if [ ! -f data/dataset.emd ]; then\n",
    "    echo \"Data not yet present, downloading\"\n",
    "    wget https://owncloud.gwdg.de/index.php/s/utJfj0388mp8W1S/download -O data/dataset.emd\n",
    "else\n",
    "    echo \"Data already present\"\n",
    "fi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's also important to check the hash of a downloaded file to verify it is the file we want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "hash=$(sha256sum data/dataset.emd | awk -F\" \" '{print $1}')\n",
    "if [ $hash == 777a5f480c5b10288bb9c83c12be440408e7dd25620adc56b5fad27bbfc65d05 ]; then\n",
    "    echo \"Hash is ok\"\n",
    "else\n",
    "    echo \"Wrong hash! This may be the wrong file, or it may have been tampered with.\"\n",
    "fi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"open\"></a> 2. Opening and exploring the data\n",
    "We will work with an `.emd` file, a semi-proprietary HDF5 based file format created by Thermo Fisher as output format for their Velox software.\n",
    "The details of this file format are not important; we can use hyperspy as a standard interface with many different microscopy data formats in the same way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hyperspy.api as hs\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file actually contains a list of different \"signals\".\n",
    "Each signal has data shape of 39x512x512, meaning 39 frames of 512x512 images.\n",
    "BF, HAADF, DF4 and DF2 are individual detectors in the microscope, representing different signals collected during the experiment.\n",
    "The DCFI datasets are actually datasets derived from the other datasets by image processing in Velox, specifically rigid registration to align the images and correct for drift."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = hs.load(\"data/dataset.emd\")\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[1].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that there is drift during the image acquisition.\n",
    "We could attempt to align the images ourselves but the DCFI option in Velox does a decent job so we continue with this.\n",
    "The last frame in the DCFI dataset is usually what we want - it represents the average of all the aligned frames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = data[3].inav[-1]\n",
    "image.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close(\"all\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can explore the acquisition conditions, microscope parameters and other metadata found in the original file through the `metadata` and `original_metadata` attribute.\n",
    "For scale information (related to the axes) we look at the `axes_manager` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# metadata in the format of hyperspy. Only a subset of all the metadata parsed from the file, but will have the same structure for all different types of files.\n",
    "image.metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# raw metadata as parsed from the file \n",
    "image.original_metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image.axes_manager"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can index into these structures just like python dictionaries to query information we care about.\n",
    "For example, convergence angle is a parameter we often care about as this determines the shape of the electron probe and hence the image formation physics.\n",
    "We might need this for simulations, but in this case we just show it's possible to get the information.\n",
    "\n",
    "Of course the data needs to be in the file to be able to query it.\n",
    "If the file format does not store critical acquisition parameters, it's important to keep track of it manually through, for example, electronic lab notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data is loaded as strings, we must cast to float\n",
    "convergence_angle = float(image.original_metadata[\"Optics\"][\"BeamConvergence\"])*1000  # data is in radians, we usually work with mrad\n",
    "c2_aperture_size = float(image.original_metadata[\"Optics\"][\"Apertures\"][\"Aperture-1\"][\"Diameter\"])*1e6  # data is in meters, we like to work with micrometers\n",
    "print(\"Convergence angle (mrad):\", convergence_angle)\n",
    "print(\"Aperture radius (micron):\", c2_aperture_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"analysis\"></a> 3. Perform some atomic column based analysis\n",
    "Just as a basic example (and to work in some machine learning) let's try to identify the atomic columns in the image and segment them based on nearest neighbor information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import atomap.api as am\n",
    "from atomap.atom_finding_refining import _remove_too_close_atoms\n",
    "import hyperspy.api as hs\n",
    "from skimage.filters import gaussian\n",
    "from skimage.exposure import rescale_intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we do some image processing and find the coordinates in the image with peak finding.\n",
    "Then we refine using center of mass (implemented in atomap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we only take a subset of the image near the grain boundary\n",
    "half_x = image.data.shape[1]//2\n",
    "half_y = image.data.shape[0]//2\n",
    "half_wx = 100\n",
    "half_wy = 200\n",
    "#sub_image = image.data[half_y-half_wy:half_y+half_wy, half_x-half_wx:half_y+half_wx]\n",
    "sub_image = image.data[:, half_x-half_wx:half_y+half_wx]\n",
    "# because we mainly care about the positioning of the hex rings, we smooth and find the centers of the hex rings instead of the individual columns\n",
    "sub_image_hex = gaussian(sub_image, 0.8)\n",
    "sub_image_hex = rescale_intensity(-sub_image_hex, out_range=(0, 1))\n",
    "\n",
    "im_peaks = sub_image_hex\n",
    "atom_positions = am.get_atom_positions(hs.signals.Signal2D(im_peaks), separation=1)\n",
    "atom_positions = _remove_too_close_atoms(atom_positions, 5)\n",
    "sublattice = am.Sublattice(atom_positions, image=im_peaks)\n",
    "sublattice.find_nearest_neighbors()\n",
    "sublattice.refine_atom_positions_using_center_of_mass()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(atom_positions.shape)\n",
    "fig, ax = plt.subplots(figsize=(5, 10))\n",
    "ax.imshow(sub_image, cmap=\"Greys_r\", vmax = np.percentile(sub_image, 98), vmin = np.percentile(sub_image, 2))\n",
    "ax.scatter(sublattice.atom_positions[:,0], sublattice.atom_positions[:,1], s=5, c=\"red\")\n",
    "ax.axis(\"off\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To know how large we should make our search distance we can look at the histogram of the distance of the n'th nearest neighbor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.spatial import cKDTree\n",
    "# check the distribution of nearest neighbor spacings\n",
    "columns = sublattice.atom_positions\n",
    "tree_c = cKDTree(columns)\n",
    "# which nearest neighbor\n",
    "NN = 1\n",
    "nn_distance = tree_c.query(columns, k=NN+1)[0][:,NN]\n",
    "print(f\"Median {NN}th nearest neighbor distance (pixels): \", np.median(nn_distance))\n",
    "fig, ax = plt.subplots()\n",
    "_ = ax.hist(nn_distance, bins=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close(\"all\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Features we might be interested in:\n",
    "* 0: the average distance of first layer of nearest neighbors\n",
    "* 1: the standard deviation in the this distance\n",
    "* 2: the average angle towards consecutive nearest neighbors\n",
    "* 3: the standard deviation of these angles\n",
    "* 4: the \"rotation\" of the voronoi cell by looking at the smallest angle with the horizontal\n",
    "* 5: the number of nearest neighbors within a certain distance\n",
    "* 6: the closest nearest neighbor\n",
    "* 7: the furthest nearest neighbor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_nn_features(coordinates, max_distance=12, k=10):\n",
    "    # A distance of 12 seems reasonable to separate first and second layer atoms\n",
    "    nearest_neighbor_indices = tree_c.query(coordinates, k=k, distance_upper_bound=max_distance)[1][:, 1:]\n",
    "    features = np.zeros((coordinates.shape[0], 8), dtype=np.float64)\n",
    "    for point_index, nns in enumerate(nearest_neighbor_indices):\n",
    "        # get rid of all the indices that don't correspond to a point\n",
    "        nns = nns[nns != coordinates.shape[0]]\n",
    "        # get the coordinates of those points\n",
    "        nn_coords = coordinates[nns]\n",
    "        # original coordinates\n",
    "        point_coords = coordinates[point_index]\n",
    "        # the difference vectors\n",
    "        difference_vectors = nn_coords - point_coords\n",
    "        # length information\n",
    "        lengths = np.linalg.norm(difference_vectors, axis=1)\n",
    "        average_length = np.mean(lengths)\n",
    "        std_length = np.std(lengths)\n",
    "        # angle information\n",
    "        angles = np.arctan2(difference_vectors[:,1], difference_vectors[:,0]) % (2*np.pi)\n",
    "        sorted_angles = np.array(np.sort(angles), angles.min()+2*np.pi)   # add the first angle again + 2*pi because we care about increments\n",
    "        angle_increments = np.diff(sorted_angles)\n",
    "        average_increment = np.mean(angle_increments)\n",
    "        std_increment = np.std(angle_increments)\n",
    "        min_angle = np.min(angles)\n",
    "        # number of vertices\n",
    "        vertices = nns.shape[0]\n",
    "        # fill in the values into the array\n",
    "        features[point_index, 0] = average_length\n",
    "        features[point_index, 1] = std_length\n",
    "        features[point_index, 2] = average_increment\n",
    "        features[point_index, 3] = std_increment\n",
    "        features[point_index, 4] = min_angle\n",
    "        features[point_index, 5] = vertices\n",
    "        features[point_index, 6] = lengths.min()\n",
    "        features[point_index, 7] = lengths.max()\n",
    "    return features\n",
    "        \n",
    "features = get_nn_features(columns, max_distance=11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot a histogram of the number of nearest neighbors found in the search distance. Most will have 6 (hexagons). The points at the edge will have <5 usually so we discard them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.hist(features[:,5], 50)\n",
    "ax.set_yscale(\"log\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we only care about points not on the edge where we have at least 5 nearest neighbors\n",
    "off_edge = features[:,5] >= 5\n",
    "filtered_features = features[off_edge]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we can plot various features on the image to see their distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the raw features to see the effect\n",
    "FI = 5\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 10))\n",
    "\n",
    "ax.imshow(sub_image, cmap=\"Greys_r\", vmax = np.percentile(sub_image, 98), vmin = np.percentile(sub_image, 2))\n",
    "\n",
    "filtered_coordinates = sublattice.atom_positions[off_edge]\n",
    "points = ax.scatter(filtered_coordinates[:,0], filtered_coordinates[:,1], s=20, c=filtered_features[:,FI], cmap=\"viridis\")\n",
    "fig.colorbar(points, ax=ax)\n",
    "\n",
    "ax.set_xlim(0,sub_image.shape[1])\n",
    "ax.set_ylim(0,sub_image.shape[0])\n",
    "ax.axis(\"off\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now want to cluster based on these features. Since they have different ranges we rescale them to a similar range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
    "scaler = MinMaxScaler()\n",
    "scaled_features = scaler.fit_transform(filtered_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With PCA we check how significant each of our features are to the principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "dimension_reduction = PCA(8)\n",
    "xy = dimension_reduction.fit_transform(scaled_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "gg = ax.imshow(dimension_reduction.components_, cmap=\"coolwarm\")\n",
    "fig.colorbar(gg)\n",
    "ax.set_xlabel(\"Feature\")\n",
    "ax.set_ylabel(\"Principal component\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then attempt a clustering using various clustering algorithms implemented in scikit learn. Meanshift seemed to give a decent answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.cluster import DBSCAN, KMeans, SpectralClustering, MeanShift\n",
    "\n",
    "cluster_model_1 = SpectralClustering(5)  # optimal kmeans\n",
    "#cluster_model_1 = MeanShift()\n",
    "#datafit = PolynomialFeatures(2).fit_transform(scaled_features)\n",
    "# it seems that the most significant features for good clustering are the standard deviation in length, the angle of the cell, and the number of nearest neighbors\n",
    "datafit = scaled_features[:, [1, 4, 5]]\n",
    "cluster_model_1.fit(datafit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.set_aspect(\"equal\")\n",
    "\n",
    "dimension_reduction = PCA(2)\n",
    "xy = dimension_reduction.fit_transform(datafit)\n",
    "ax.scatter(xy[:,0], xy[:,1], c=cluster_model_1.labels_, cmap=\"Set1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 10))\n",
    "ax.imshow(sub_image, cmap=\"Greys_r\", vmax = np.percentile(sub_image, 98), vmin = np.percentile(sub_image, 2))\n",
    "not_classified = sublattice.atom_positions[np.invert(off_edge)]\n",
    "ax.scatter(filtered_coordinates[:,0], filtered_coordinates[:,1], s=20, c=cluster_model_1.labels_, cmap=\"Set1\")\n",
    "ax.scatter(not_classified[:,0], not_classified[:,1], marker=\"s\", s=10, c=\"black\")\n",
    "ax.set_xlim(0,sub_image.shape[1])\n",
    "ax.set_ylim(0,sub_image.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close(\"all\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now color in the picture according to the nearest cluster index to segment it on a pixel by pixel basis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = cKDTree(filtered_coordinates)\n",
    "Y, X = np.mgrid[0:sub_image.shape[0], 0:sub_image.shape[1]]\n",
    "coords = np.stack([X.ravel(), Y.ravel()]).T\n",
    "nearest_point = tree.query(coords)[1]\n",
    "classification = cluster_model_1.labels_[nearest_point]\n",
    "labeled_image = classification.reshape(sub_image.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 10))\n",
    "ax.imshow(sub_image, cmap=\"Greys_r\", vmax = np.percentile(sub_image, 98), vmin = np.percentile(sub_image, 2))\n",
    "ax.imshow(labeled_image, cmap = \"Set1\", alpha=0.5)\n",
    "ax.axis(\"off\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close(\"all\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further cleaning of the image and more playing with the clustering may be able to improve the segmentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}