{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Radial Metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "# from project_heart.lv import LV\n",
    "from project_heart.utils import set_jupyter_backend\n",
    "from project_heart.enums import *\n",
    "set_jupyter_backend(\"pythreejs\")\n",
    "import logging\n",
    "import numpy as np\n",
    "float_formatter = \"{:.5f}\".format\n",
    "np.set_printoptions(formatter={'float_kind':float_formatter})\n",
    "\n",
    "from project_heart.examples import get_lv_ideal\n",
    "from project_heart.examples import get_lv_typeA\n",
    "\n",
    "lv_ideal = get_lv_ideal(Path(\"../../_static/sample_files/ideal_linear_pressure_increase copy.xplt\"))\n",
    "lv_typeA = get_lv_typeA(\n",
    "    Path(\"../../_static/sample_files/lv_typeA_hex.vtk\"),\n",
    "    Path(\"../../_static/sample_files/sample_displacement_lv_typeA_hex_states_linear_press_incr.pbz2\")\n",
    "    )\n",
    "\n",
    "sample_spk_typeA = lv_typeA.get_speckles(spk_name=\"SAMPLE\", spk_group=\"epi\", spk_collection=\"SAMPLE\")[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radial metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our implementation, 'radial metrics' are defined as any metric that relates a point with the left ventricle's longitudinal axis. In a clinical setting, these metrics are used to compute radial shortening and radial strain. However, it is often unclear how to proper establish its calculation as different reference points related to the longitudinal axis can be selected. \n",
    "For instance, we can use a 'center' point located at a specified 'height' of the longitudinal axis and compute the vector length from this reference to a point on the endocardium surface. On the other hand, we can compute the perpendicular distance between a point on the endocardium and the longitudinal axis. \n",
    "With this in mind, we provided two distinct metrics for radial measurements: 'Radial Distance' and 'Radial Length'. \n",
    "\n",
    "In addition, our implementation consider the usage of [Speckles](../basic_definitions/speckles.ipynb), which helps to approximate a clinical setting, facilitates the computation for certain metrics and helps to minimizing errors due to noise in geometry by averaging values based on local regions. See [docs](../basic_definitions/speckles.ipynb) for details.\n",
    "\n",
    "Moreover, the final metric is computed by applying [reduction](./radial_metrics.ipynb) accross values computed from individual speckles. The [default reduction method](./radial_metrics.ipynb) is the 'mean' value for all selected speckles. For instance, if we apply this computation accross speckles at endocardium, we will have the mean values for the endocardium region. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For documentation purposes, we will be using 'Sample' speckles at the Epicardium:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_spk_ideal = lv_ideal.get_speckles(spk_name=\"SAMPLE\", spk_group=\"epi\", spk_collection=\"SAMPLE\")[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Radial Distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Definition**: Perpendicular distance from a given point to normal vector (defined from longitudinal line)\n",
    "\n",
    "The 'Radial distance' is computed as the perpendicular distance (shortest euclidean distance) between a point (based on speckles) at a given region of the geometry and the longitudinal axis. Considering tha such axis is aligned with the Z axis, this metric can be thought as the '2D radius' from a top view (xy plane)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rd = lv_ideal.plot_speckles_radial_distances(sample_spk_ideal, t=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*rd represents the average radial distance among the selected speckles. We will use it to compare with our metric computation.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rds = lv_ideal.radial_distance(\n",
    "            sample_spk_ideal,         # selected speckles\n",
    "            approach=\"moving_vector\", # choice of reference LA axis\n",
    "            method=\"mean\",            # reduction method (default)\n",
    "            t=0.0,                    # return value at given timestep\n",
    "            recompute=True,           # forces recomputation (ignores previous computations)\n",
    "            log_level=logging.ERROR, \n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check if our algorithm implementation matches the expected result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Sample radial distance:\", rd)\n",
    "print(\"Reduced radial distance by internal algorithm:\", rds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fixed vs Moving vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During a cardiac cycle, contractions may cause the longitudinal axis of the left ventricle to change. This occurs spontaneously owing to the fiber composition of the heart tissue and affects the radial distance value. \n",
    "\n",
    "Taking this into account, we have two approaches regarding radial distance computation: ***fixed_vector*** and ***moving_vetor***. The first approach uses the reference configuration as the basis vector for computation, while the second accounts for possible oscillations that might occur thoughout timesteps by using variable positions for apex and base reference nodes (thus creating a 'new' longitudinal axis for each timestep).\n",
    "\n",
    "Here is an example of how these two approaches can impact the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rd_mv = lv_typeA.plot_speckles_radial_distances(sample_spk_typeA, t=0.15, \n",
    "                                             approach=\"moving_vector\", \n",
    "                                             plot_kwargs={\"style\": None})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rd_fv = lv_typeA.plot_speckles_radial_distances(sample_spk_typeA, t=0.15, \n",
    "                                             approach=\"fixed_vector\", \n",
    "                                             plot_kwargs={\"style\": None})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the numerical comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"with moving vector:\", rd_mv)\n",
    "print(\"with fixed vector:\", rd_fv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in [1](https://core.ac.uk/reader/53744593?utm_source=linkout), the base reference node should be estimated at mid-point the 'right' and 'left' base; transposing this to 3 dimensions, the base should be estimated as the centroid of top-most nodes of endocardium. Our implementation takes this into consideration, see [geo references](../basic_definitions/geo_references.ipynb) for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Radial Length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Definition**: Vector magnitude from a given point to a point located at the longitudinal axis (not necessarily the perpendicular distance)\n",
    "\n",
    "The 'Radial Length' is computed using the vector maginitude from a point at given region of the geometry and a corresponding 'center' located along the longitudinal axis. This reference 'center' is defined based on the longitudinal axis at an average height of the query points. If only one point is in the query, then this metric should reflect the same value as 'radial distance'.\n",
    "The purpose of this metric is to capture not only radial characteristics in a single plane (perpendicular to geometry's normal), but also curvature information in the given region. \n",
    "\n",
    "Note that this metric is intended to be used with speckles, therefore, the reference 'center' will be the corresponding point along the longitudinal axis at the mean heigh of the given speckles. \n",
    "\n",
    "In addition, as we base the 'height' on Z axis *(for now)*, this algorithm is constrained on the geometry's longitudinal axis being aligned with the Z axis.\n",
    "\n",
    "Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rl = lv_ideal.plot_speckles_radial_lengths(sample_spk_ideal, t=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rls = lv_ideal.radial_length(\n",
    "            sample_spk_ideal,          # selected speckles\n",
    "            approach=\"moving_centers\", # choice of reference LA axis\n",
    "            method=\"mean\",             # reduction method (default)\n",
    "            t=0.0,                     # return value at given timestep\n",
    "            recompute=True,            # forces recomputation (ignores previous computations)\n",
    "            log_level=logging.ERROR, \n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Sample radial length:\", rl)\n",
    "print(\"Reduced radial length by internal algorithm:\", rls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fixed vs moving centers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly from radial distance, we implemented two approaches for radial length: ***fixed_centers*** and ***moving_centers***. The first uses the reference center when speckles were created, while the second recomputes the centers based on the 'moving' longitudinal axis. Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rl_mv = lv_typeA.plot_speckles_radial_lengths(sample_spk_typeA, t=0.15, \n",
    "                                             approach=\"moving_centers\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rl_fv = lv_typeA.plot_speckles_radial_lengths(sample_spk_typeA, t=0.15, \n",
    "                                             approach=\"fixed_centers\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"with moving centers:\", rl_mv)\n",
    "print(\"with fixed centers:\", rl_fv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing Radial Metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing thickness between speckles and their effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "def plot_all_radial_metrics(lv, spks, title=\"\"):\n",
    "        timestep = lv.timesteps()\n",
    "        mask = timestep >=0.1\n",
    "        \n",
    "        radial_dist_1 = lv.radial_distance(spks, \n",
    "                                        reduce_by={\"group\", \"name\"},\n",
    "                                        approach=\"fixed_vector\",\n",
    "                                        recompute=True, \n",
    "                                        log_level=logging.ERROR)[mask]\n",
    "        radial_dist_2 = lv.radial_distance(spks, \n",
    "                                        reduce_by={\"group\", \"name\"},\n",
    "                                        approach=\"moving_vector\",\n",
    "                                        #  ql=\"\",\n",
    "                                        recompute=True, \n",
    "                                        log_level=logging.ERROR)[mask]\n",
    "        radial_length_1 = lv.radial_length(spks, \n",
    "                                        recompute=True, \n",
    "                                        reduce_by={\"group\", \"name\"},\n",
    "                                        approach=\"fixed_centers\",\n",
    "                                        log_level=logging.ERROR)[mask]\n",
    "        radial_length_2 = lv.radial_length(spks, \n",
    "                                        recompute=True, \n",
    "                                        reduce_by={\"group\", \"name\"},\n",
    "                                        approach=\"moving_centers\",\n",
    "                                        log_level=logging.ERROR)[mask]\n",
    "        timestep = timestep[mask] - 0.1\n",
    "        data = np.hstack([\n",
    "                        timestep.reshape((-1,1)), \n",
    "                        radial_dist_1.reshape((-1,1)), \n",
    "                        radial_dist_2.reshape((-1,1)),\n",
    "                        radial_length_1.reshape((-1,1)), \n",
    "                        radial_length_2.reshape((-1,1))\n",
    "                        ])\n",
    "        \n",
    "        import pandas as pd\n",
    "        \n",
    "\n",
    "        df = pd.DataFrame(data,columns=[\n",
    "                \"timestep\", \n",
    "                \"radial_distance_fixed_vector\", \n",
    "                \"radial_distance_moving_vector\", \n",
    "                \"radial_length_fixed_centers\",\n",
    "                \"radial_length_moving_centers\",\n",
    "                ])\n",
    "        df.plot(x=\"timestep\", figsize=(10,5), grid=True, title=title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we are comparing the effect of using thick, thin or mid-size speckes in terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thick_spks = lv_ideal.get_speckles(spk_name=\"THICK\", spk_group=\"endo\", spk_collection=\"SAMPLE\")\n",
    "mid_spks = lv_ideal.get_speckles(spk_name=\"MID\", spk_group=\"endo\", spk_collection=\"SAMPLE\")\n",
    "thin_spks = lv_ideal.get_speckles(spk_name=\"THIN\", spk_group=\"endo\", spk_collection=\"SAMPLE\")\n",
    "\n",
    "plot_all_radial_metrics(lv_ideal, thick_spks, \"Ideal, using THICK speckles.\")\n",
    "plot_all_radial_metrics(lv_ideal, mid_spks, \"Ideal, using MID speckles.\")\n",
    "plot_all_radial_metrics(lv_ideal, thin_spks, \"Ideal, using THIN speckles.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same comparison, but on a typeA geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thick_spks = lv_typeA.get_speckles(spk_name=\"THICK\", spk_group=\"endo\", spk_collection=\"SAMPLE\")\n",
    "mid_spks = lv_typeA.get_speckles(spk_name=\"MID\", spk_group=\"endo\", spk_collection=\"SAMPLE\")\n",
    "thin_spks = lv_typeA.get_speckles(spk_name=\"THIN\", spk_group=\"endo\", spk_collection=\"SAMPLE\")\n",
    "\n",
    "plot_all_radial_metrics(lv_typeA, thick_spks, \"Ideal, using THICK speckles.\")\n",
    "plot_all_radial_metrics(lv_typeA, mid_spks, \"Ideal, using MID speckles.\")\n",
    "plot_all_radial_metrics(lv_typeA, thin_spks, \"Ideal, using THIN speckles.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing Radial metrics at different locations of the geometry:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compare how the metric behaves when using different speckles at different regions of the geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_all_radial_metrics(lv_ideal, sample_spk_ideal, \"Ideal, using SAMPLE speckles.\")\n",
    "plot_all_radial_metrics(lv_typeA, sample_spk_typeA, \"TypeA, using SAMPLE speckles.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circ_spks = lv_ideal.get_speckles(spk_name=\"base\", spk_group=\"endo\", spk_collection=\"circ-6\")\n",
    "plot_all_radial_metrics(lv_ideal, circ_spks, \"Ideal, using all 'endo' at base speckles.\")\n",
    "circ_spks = lv_typeA.get_speckles(spk_name=\"base\", spk_group=\"endo\", spk_collection=\"circ-6\")\n",
    "plot_all_radial_metrics(lv_typeA, circ_spks, \"TypeA, using all 'endo' at base speckles.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circ_spks = lv_ideal.get_speckles(spk_group=\"endo\", spk_collection=\"circ-6\")\n",
    "plot_all_radial_metrics(lv_ideal, circ_spks, \"Ideal, using all 'endo' speckles.\")\n",
    "circ_spks = lv_typeA.get_speckles(spk_group=\"endo\", spk_collection=\"circ-6\")\n",
    "plot_all_radial_metrics(lv_typeA, circ_spks, \"TypeA, using all 'endo' speckles.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circ_spks = lv_ideal.get_speckles(spk_collection=\"circ-6\")\n",
    "plot_all_radial_metrics(lv_ideal, circ_spks, \"Ideal, using all 'circ-6' speckles.\")\n",
    "circ_spks = lv_typeA.get_speckles(spk_collection=\"circ-6\")\n",
    "plot_all_radial_metrics(lv_typeA, circ_spks, \"TypeA, using all 'circ-6' speckles.\")"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "a4b270fe5129ef310552b5197d637a230d47119f1e46f0b92ff32b1c5f5c5425"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 ('project_heart')",
   "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.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}