{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How To Use HAPPy:\n",
    "\n",
    "### Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First import matplotlib (for plotting) and skan\n",
    "from matplotlib import pyplot as plt\n",
    "from skan import draw\n",
    "import numpy as np\n",
    "from skimage import exposure\n",
    "\n",
    "# Then import the radial hydride packagess\n",
    "from HAPPY import import_image\n",
    "from HAPPY import cropping_functions as crop\n",
    "from HAPPY import plot_functions as plt_f\n",
    "from HAPPY import radial_hydride_fraction as RHF\n",
    "from HAPPY import branching as branch\n",
    "from HAPPY import crack_path as cp\n",
    "from HAPPY import image_processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing Image\n",
    "\n",
    "* First, import the image using the `import_image` command. Transpose\n",
    "  the image if necessary using the `transpose` argument to make the\n",
    "  radial direction vertical.\n",
    "\n",
    "* The `cropImage` function applies a rectangular crop to the image to\n",
    "  remove scale bars, or if you have a specific rectangular region you\n",
    "  want to look at.\n",
    "\n",
    "* Input Scale Bar Value in Scale_Bar_Micron_Value and\n",
    "  Pixels_In_Scale_Bar, the scale bar will then be calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load image\n",
    "original_image = import_image.image(image_path ='data/520-5b.png', transpose = False)\n",
    "cropped_image = crop.cropImage(original_image, crop_bottom=50, crop_top=0, crop_left=0, crop_right=0)\n",
    "crop1 = cropped_image\n",
    "# Input the value of the scale bar in microns\n",
    "Scale_Bar_Micron_Value = 100\n",
    "#Input how many pixels are in your scale bar\n",
    "Pixels_In_Scale_Bar = 165.5\n",
    "Scale_Bar_Value_In_Meters = Scale_Bar_Micron_Value*(1e-6)\n",
    "scale = Scale_Bar_Value_In_Meters/Pixels_In_Scale_Bar\n",
    "scale_um = scale*1e6\n",
    "location = 'lower right'\n",
    "\n",
    "\n",
    "# Plot image\n",
    "plt_f.plot(img=cropped_image, title='Loaded image',scale=scale, location=location)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Additional Cropping\n",
    "\n",
    "The second crop function is `cropping_tube`, which should be used if\n",
    "the micrograph is curved and removes black pats of the image which are\n",
    "not the tube. A crop_param of around 0.1-0.2 is reccomended."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Crop tube\n",
    "cropped_image, crop_threshold = crop.cropping_tube(cropped_image,\n",
    "                                                   crop_param = 0.2, size_param = 1000, dilation_param = 10)\n",
    "\n",
    "# Plot comparison\n",
    "plt_f.plot_comparison(crop1, 'Original image crop', cropped_image, 'Tube crop',scale=scale,\n",
    "                   location=location)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Processing\n",
    "\n",
    "Grain contast or uneven lighting can be minimised through the\n",
    "application of a gaussian blur in the `minimize_grain_contrast`\n",
    "function. A value of 10 seems to work for most cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove grain contrast\n",
    "removed_grains = image_processing.minimize_grain_contrast(cropped_image, sigma = 10)\n",
    "\n",
    "# Plot image\n",
    "plt_f.plot(img=cropped_image, title='Minimised grain contrast', scale=scale, location=location)\n",
    "\n",
    "# Plot the histogram for removed grains so that we can see where we should threshold\n",
    "histogram = plt_f.plot_hist(removed_grains)\n",
    "\n",
    "# Print an approximate threshold value which should work well\n",
    "print('Approximate threshold: {0:.3f}'.format(\n",
    "    2*np.nanmedian(removed_grains)-np.nanpercentile(removed_grains, 90)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thresholding\n",
    "\n",
    "After this, the image is thresholded using the `simple_threshold`\n",
    "function. The threshold value should be set using the `threshold`\n",
    "argument. Small features, less than a given size in microns\n",
    "`small_obj` can optionally be removed. Note it is important not too\n",
    "over threshold the image, guidance of a value to threshold is shown\n",
    "above and can be determined by investigating the histograms plotted\n",
    "above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply threshold\n",
    "thres = image_processing.simple_threshold(removed_grains,scale_um, crop_threshold,\n",
    "                                          threshold = 0.98, small_obj = 40)\n",
    "\n",
    "# Plot the thresholded image and compare it to the original image:\n",
    "plt_f.plot_comparison(cropped_image, 'Original Image', thres,'Thresholded Image', scale=scale,location=location)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to perform the hough line transform `hough_rad`\n",
    "there are a few input parameters that should be considered: -\n",
    "`num_peaks`: should be changed dependent on the type of micrograph, if\n",
    "your hydrides are straight and not very interconnected a small value of\n",
    "around 2 is good, if in one box, there are many branches that need to be\n",
    "picked up, this value should be increased accordingly to a value of 5 or\n",
    "more. - `min_dist`, `min_angle` and `val` are pre-set and seem to\n",
    "work for most cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply Hough transform\n",
    "angle_list,len_list = RHF.hough_rad(thres, num_peaks=2, scale=scale, location=location)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Non weighted radial hydride fraction\n",
    "radial, circumferential = RHF.RHF_no_weighting_factor(angle_list, len_list)\n",
    "\n",
    "print('The non-weighted RHF  is {0:.4f}'.format(radial))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Weighted Radial Hydride Fraction\n",
    "fraction = RHF.weighted_RHF_calculation(angle_list, len_list)\n",
    "\n",
    "print('The weighted RHF is: {0:.4f}'.format(fraction))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Other Methods for Radial Hydride Fraction Calculation\n",
    "\n",
    "Here all four different RHF calculation methods are shown in the graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#chu radial hydride calculation\n",
    "deg_angle_list = np.rad2deg(angle_list)\n",
    "\n",
    "radial_list_chu=[]\n",
    "circum_list_chu = []\n",
    "\n",
    "for k in deg_angle_list:\n",
    "    if (k>0 and k<40) or (k>-40 and k<0) :\n",
    "        radial_list_chu.append(len_list)\n",
    "    elif (k>50 and k<90) or (k>-90 and k<-50):\n",
    "        circum_list_chu.append(len_list)\n",
    "\n",
    "\n",
    "rad_hyd_chu = np.sum(radial_list_chu)\n",
    "cir_hyd_chu = np.sum(circum_list_chu)\n",
    "\n",
    "\n",
    "RHFChu = rad_hyd_chu/(rad_hyd_chu+cir_hyd_chu)\n",
    "\n",
    "\n",
    "#RHF 40 deg\n",
    "radial_list_40=[]\n",
    "circum_list_40 = []\n",
    "\n",
    "for k in deg_angle_list:\n",
    "    if (k>0 and k<40) or (k>-40 and k<0) :\n",
    "        radial_list_40.append(len_list)\n",
    "    elif (k>=40 and k<90) or (k>-90 and k<=-40):\n",
    "        circum_list_40.append(len_list)\n",
    "\n",
    "\n",
    "rad_hyd_40 = np.sum(radial_list_40)\n",
    "cir_hyd_40 = np.sum(circum_list_40)\n",
    "\n",
    "\n",
    "RHF40 = rad_hyd_40/(rad_hyd_40+cir_hyd_40)\n",
    "\n",
    "import pandas as pd\n",
    "# intialise data of lists.\n",
    "data = {\"RHF\": [RHF40,radial,fraction,RHFChu]\n",
    "       }\n",
    "\n",
    "# Create DataFrame\n",
    "df = pd.DataFrame(data,index=[\"40 Degrees\", \"45 Degrees\", \"Weighted\", \"Chu\"])\n",
    "display(df)\n",
    "\n",
    "#d = {\"one\": [1.0, 2.0, 3.0, 4.0], \"two\": [4.0, 3.0, 2.0, 1.0]}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mean Hydride Length\n",
    "\n",
    "Code for determining the MHL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import ndimage\n",
    "\n",
    "hydride_len = []\n",
    "label, num_features = ndimage.label(thres > 0.1)\n",
    "slices = ndimage.find_objects(label)\n",
    "for feature in np.arange(num_features):\n",
    "    hydride_len.append(scale_um*label[slices[feature]].shape[1])\n",
    "\n",
    "#print(hydride_len)\n",
    "print(np.mean(hydride_len))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Branch Length Fraction\n",
    "\n",
    "Here we want to determine the extent of branching within the\n",
    "microstrucutre, this is done in two ways: - In image form where the\n",
    "branches are coloured differently to the main hydride - BLF the length\n",
    "fraction of branches with respect to the toatal length of all hydrides\n",
    "in the microstrucutre"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the branch length fraction\n",
    "skel,is_main,BLF = branch.branch_classification(thres);\n",
    "\n",
    "\n",
    "# Plot branching image\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax = draw.overlay_skeleton_2d_class(\n",
    "    skel,\n",
    "    skeleton_color_source=lambda s: is_main,\n",
    "    skeleton_colormap='spring',\n",
    "    axes=ax\n",
    "     )\n",
    "\n",
    "plt.axis('off')\n",
    "plt.title('Branched hydrides')\n",
    "#plt_f.addScaleBar(ax[0], scale=scale, location=location)\n",
    "plt_f.addArrows(ax[0])\n",
    "\n",
    "print('The BLF is: {0:.4f}'.format(BLF))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Crack Path\n",
    "\n",
    "Here we want to determine potential crack paths through the\n",
    "microstrucutre, we input the thresholded image `thres`. After running\n",
    "once, the area around that path (radius set with `kernel_size`) is\n",
    "discounted, then the process is repeated `num_runs` times. Here the\n",
    "`distance_weight` makes moving in the circumferential direction more\n",
    "costly, note when comparing different micrographs, ensure that this\n",
    "parameter it is kept constant. We reccomend a weighting of 1.5 and a\n",
    "kernel size of 20."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determing potential crack paths\n",
    "edist, path_list, cost_list = cp.det_crack_path(thres, crop_threshold, num_runs=5, kernel_size=20,distance_weight=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot possible crack paths\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "list_costs = []\n",
    "\n",
    "for n, (p, c) in enumerate(zip(path_list, cost_list)):\n",
    "\n",
    "    im = ax.imshow(thres, cmap='gray')\n",
    "\n",
    "    #if n==0:\n",
    "      #  plt.colorbar(im,fraction=0.03, pad=0.01)\n",
    "    ax.scatter(p[:,1], p[:,0], s=10, alpha=0.1)\n",
    "    ax.text(p[-1][1], p[-1][0], s=str(n), c='w', bbox=dict(facecolor='black', edgecolor='black'))\n",
    "    plt.axis('off')\n",
    "    print('Run #{0}\\tCost = {1:.2f}'.format(n,c))\n",
    "    list_costs.append(c)\n",
    "\n",
    "plt_f.addScaleBar(ax, scale=scale, location=location)\n",
    "plt_f.addArrows(ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Histograms for plotting the costs of each path\n",
    "plt.hist(list_costs, bins=5, cumulative = True, color = \"cornflowerblue\", ec=\"cornflowerblue\", label = \"Cumulative Distribution Function\")\n",
    "plt.hist(list_costs, bins=5, color = \"lightpink\", ec=\"lightpink\", label = \"Normal Histogram\")\n",
    "plt.legend()\n",
    "plt.xlabel('Cost', fontsize=\"12\")\n",
    "plt.ylabel('Frequency',fontsize=\"12\")\n",
    "plt.title('Paths of Lowest Cost', fontweight=\"bold\", fontsize=\"15\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can chose to skeletonize the image if you want, not reccomended\n",
    "unless there are too many hydrides to be able to distinguish between\n",
    "them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.morphology import skeletonize\n",
    "skeletonised = skeletonize(thres)\n",
    "plt.imshow(skeletonised,cmap='gray')\n",
    "plt.axis('off')\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 4
}