PCA-2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T01:48:57.409496Z",
     "iopub.status.busy": "2020-05-13T01:48:57.409309Z",
     "iopub.status.idle": "2020-05-13T01:48:59.326598Z",
     "shell.execute_reply": "2020-05-13T01:48:59.326130Z",
     "shell.execute_reply.started": "2020-05-13T01:48:57.409476Z"
=======
   "execution_count": 108,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:08:20.690397Z",
     "iopub.status.busy": "2020-05-12T23:08:20.690056Z",
     "iopub.status.idle": "2020-05-12T23:08:20.695626Z",
     "shell.execute_reply": "2020-05-12T23:08:20.694959Z",
     "shell.execute_reply.started": "2020-05-12T23:08:20.690357Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.io\n",
    "from scipy.ndimage import zoom\n",
    "from pathlib import Path\n",
    "from stl import mesh\n",
    "from mayavi import mlab\n",
    "import quaternion\n",
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T01:48:59.327379Z",
     "iopub.status.busy": "2020-05-13T01:48:59.327253Z",
     "iopub.status.idle": "2020-05-13T01:48:59.330083Z",
     "shell.execute_reply": "2020-05-13T01:48:59.329361Z",
     "shell.execute_reply.started": "2020-05-13T01:48:59.327362Z"
=======
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T22:21:34.345430Z",
     "iopub.status.busy": "2020-05-12T22:21:34.345265Z",
     "iopub.status.idle": "2020-05-12T22:21:34.348083Z",
     "shell.execute_reply": "2020-05-12T22:21:34.347547Z",
     "shell.execute_reply.started": "2020-05-12T22:21:34.345412Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "# from IPython.core.debugger import set_trace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installation: https://lukemshepherd.github.io/vox/install/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "# Classes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Class manages the loading of voxel and stl files and the methods that manipulate the bones  "
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T01:48:59.331915Z",
     "iopub.status.busy": "2020-05-13T01:48:59.331692Z",
     "iopub.status.idle": "2020-05-13T01:48:59.369360Z",
     "shell.execute_reply": "2020-05-13T01:48:59.368583Z",
     "shell.execute_reply.started": "2020-05-13T01:48:59.331884Z"
=======
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T22:21:36.354835Z",
     "iopub.status.busy": "2020-05-12T22:21:36.354634Z",
     "iopub.status.idle": "2020-05-12T22:21:36.376100Z",
     "shell.execute_reply": "2020-05-12T22:21:36.375506Z",
     "shell.execute_reply.started": "2020-05-12T22:21:36.354811Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "class bone:\n",
    "\n",
    "    filter_level = 0.001\n",
    "    default_color = (0.7, 1, 1)\n",
    "\n",
    "    def __init__(self, data, dtype):\n",
    "        \"\"\"\n",
    "        Performs calculations on the voxel array objects    \n",
    "        array (np.array): binary voxel object)      \n",
    "        filter_level (int/float): sets the threshold level for \n",
    "        what is considered a voxel. Everything below filter level is\n",
    "        rounded to 0, everything above rounded to 1 (ie voxel)\n",
    "        \"\"\"\n",
    "        \n",
    "        self.dtype = dtype\n",
    "        self.data = data\n",
    "        \n",
    "        self.get_xyz()\n",
    "\n",
    "                \n",
    "        \n",
    "    def get_xyz(self):\n",
    "        \"\"\"Convert 3D voxel array or STL to xyz coordinates.\n",
    "\n",
    "        filter_level (int/float): (inherited from `bone` class) sets the threshold level for \n",
    "        what is considered a voxel. Everything below filter level is\n",
    "        rounded to 0, everything above rounded to 1 (ie voxel)\n",
    "\n",
    "        returns: \n",
    "            np.array( [n x 3] )\"\"\"\n",
    "\n",
    "        \n",
    "        if self.dtype == 'voxel':\n",
    "        \n",
    "            # Everything above filter level is converted to 1\n",
    "            filtered_array = np.where(self.data < self.filter_level, 0, 1)\n",
    "\n",
    "            # records coordiates where there is a 1\n",
    "            x, y, z = np.where(filtered_array == 1)\n",
    "\n",
    "            self.xyz = np.array([x, y, z]).T\n",
    "            \n",
    "            \n",
    "        elif self.dtype == 'stl':\n",
    "            self.xyz = np.concatenate((self.data.v0,\n",
    "                                         self.data.v1,\n",
    "                                         self.data.v2), axis=0)\n",
    "\n",
    "                \n",
    "    def get_pca(self):\n",
    "        \"\"\"PCA on the xyz points array\n",
    "\n",
    "            xyz(np.array): n x 3 array of xyz coordinates\n",
    "\n",
    "            returns:    self.pc1\n",
    "                        self.pc2\n",
    "                        self.pc3\"\"\"\n",
    "\n",
    "        pca = PCA(svd_solver='full')\n",
    "        pca.fit(self.xyz)\n",
    "         \n",
    "        self.pca_list = pca.components_\n",
    "        self.pc1 = pca.components_[0]\n",
    "        self.pc2 = pca.components_[1]\n",
    "        self.pc3 = pca.components_[2]\n",
    "\n",
    "   \n",
    "    \n",
    "    @property\n",
    "    def mean(self):\n",
    "        \"\"\"The mean of the xyz atriube \n",
    "            returns:\n",
    "            tupple (mean_of_x, mean_of_y ,mean_of_z)\"\"\"\n",
    "    \n",
    "        return (np.mean(self.xyz[:, 0]), np.mean(self.xyz[:, 1]), np.mean(self.xyz[:, 2]))\n",
    "\n",
    "\n",
    "    def center_to_origin(self):\n",
    "        \"\"\" sets the mean of the bone to 0,0,0\"\"\"\n",
    "\n",
    "        # set transformation (tfm) value\n",
    "        self.tfm = self.mean\n",
    "\n",
    "        self.xyz = self.xyz - self.mean\n",
    "\n",
    "    def reset_position(self):\n",
    "        \"\"\" resets the position of the bone to its orginal one\"\"\"\n",
    "        self.xyz = self.xyz + self.tfm\n",
    "\n",
    "\n",
    "    def plot(self, user_color=None, mesh=False, PCA_inv=False, PCA=True):\n",
    "        \"\"\" Plot voxels with optional PCA, and colours\n",
    "        \n",
    "            user_color (tupple): RGB color of the bone where 1 is maxium\n",
    "                                    eg: red = (1,0,0)\n",
    "                                    \n",
    "            PCA (boolean): plots the PCAs of the voxel\n",
    "            \n",
    "            PCA_inv (boolean): plots the inverse of each PCA so the axes go in both directions\n",
    "        \"\"\"\n",
    "                \n",
    "        if hasattr(self, 'pc1') is False:\n",
    "            self.get_pca()\n",
    "        \n",
    "    \n",
    "        if user_color is None:\n",
    "            user_color = self.default_color\n",
    "            \n",
    "            \n",
    "            \n",
    "        if mesh is False:\n",
    "            #plots points\n",
    "            mlab.points3d(self.xyz[:, 0],\n",
    "                          self.xyz[:, 1],\n",
    "                          self.xyz[:, 2],\n",
    "                          mode = \"cube\",\n",
    "                          color= user_color,\n",
    "                          scale_factor = 1)  \n",
    "            \n",
    "        else:\n",
    "             mlab.mesh(self.data.x, self.data.y, self.data.z)\n",
    "        \n",
    "\n",
    "        def quiver_pca(n,i):\n",
    "            mlab.quiver3d(*self.mean, *(getattr(self,f'pc{n}')*i),\n",
    "                                  line_width=6,\n",
    "                                  scale_factor=100/n,\n",
    "                                  color=c)\n",
    "\n",
    "        for n in range(1,4):\n",
    "            #sets color: red = pc1, blue = pc2, green = pc3\n",
    "            c = [0,0,0]\n",
    "            c[n-1] = 1\n",
    "            c = tuple(c)\n",
    "            \n",
    "            # plots pca arrows\n",
    "            if PCA is True:\n",
    "                quiver_pca(n,1)\n",
    "                \n",
    "            #plots the pca *-1      \n",
    "            if PCA_inv is True:\n",
    "                quiver_pca(n,-1)\n",
    "            \n",
    "        \n",
    "    def scale(self, n, algo = 'constant'):\n",
    "        \"\"\" up-scales the bone by n\n",
    "        \n",
    "            n: scale factor \n",
    "            \n",
    "            algo: method of upscaling array \n",
    "            scipy.ndimagezoom(mode=...)\"\"\"\n",
    "        \n",
    "        self.data = zoom(self.data, (n, n, n), mode=algo)\n",
    "        \n",
    "        #update xyz\n",
    "        self.get_xyz()\n",
    "            \n",
    "            \n",
    "    def xyz_to_array(self,array_dim=(256,256,256)):\n",
    "        \"\"\" Converts xyz coordinates to numpy voxel array\"\"\"\n",
    "\n",
    "        #empty array\n",
    "        vx_array = np.zeros(array_dim,dtype=bool)\n",
    "    \n",
    "        #for every xyz coord: if int(+- tolarance) write to array\n",
    "        for i in self.xyz:\n",
    "            if np.allclose(i, np.around(i), rtol= 0.5, equal_nan=True):\n",
    "                vx_array[tuple(np.around(i).astype(int))] = True\n",
    "                \n",
    "        x = np.count_nonzero(vx_array)/self.xyz.shape[0]\n",
    "        \n",
    "        print(f'{x*100}% reconstructed')\n",
    "        \n",
    "        return vx_array\n",
    "\n",
    "\n",
    "    @classmethod\n",
    "    def from_matlab_path(cls, matlab_file):\n",
    "        \"\"\"Imports matlab file drectly\n",
    "\n",
    "           path: path object/string \n",
    "\n",
    "           retruns np.array (n x n x n )\"\"\"\n",
    "\n",
    "        matlab_file = Path(matlab_file)\n",
    "\n",
    "        matlab_object = scipy.io.loadmat(matlab_file)\n",
    "        obj = matlab_object.keys()\n",
    "        obj = list(obj)\n",
    "        data = matlab_object[obj[-1]]\n",
    "\n",
    "        return cls(data, dtype=\"voxel\")\n",
    "\n",
    "    @classmethod\n",
    "    def from_stl_path(cls, stl_file):\n",
    "        \"\"\"Imports stl file drectly\n",
    "\n",
    "       path: path object/string \n",
    "\n",
    "       retruns np.array (n x n x n )\"\"\"\n",
    "\n",
    "        stl_file = Path(stl_file)\n",
    "\n",
    "        data = mesh.Mesh.from_file(stl_file)\n",
    "\n",
    "        return cls(data, dtype=\"stl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maths functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## mag"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T01:48:59.370555Z",
     "iopub.status.busy": "2020-05-13T01:48:59.370384Z",
     "iopub.status.idle": "2020-05-13T01:48:59.377848Z",
     "shell.execute_reply": "2020-05-13T01:48:59.377228Z",
     "shell.execute_reply.started": "2020-05-13T01:48:59.370533Z"
=======
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T22:21:41.023262Z",
     "iopub.status.busy": "2020-05-12T22:21:41.023095Z",
     "iopub.status.idle": "2020-05-12T22:21:41.026196Z",
     "shell.execute_reply": "2020-05-12T22:21:41.025624Z",
     "shell.execute_reply.started": "2020-05-12T22:21:41.023244Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "def mag(v):\n",
    "    \"\"\" Finds magnitude of vector\n",
    "\n",
    "        v (np.array): vector\"\"\"\n",
    "    return math.sqrt(np.dot(v, v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## angle\n",
    "\n",
    "https://uk.mathworks.com/matlabcentral/answers/101590-how-can-i-determine-the-angle-between-two-vectors-in-matlab"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T01:48:59.378794Z",
     "iopub.status.busy": "2020-05-13T01:48:59.378645Z",
     "iopub.status.idle": "2020-05-13T01:48:59.384447Z",
     "shell.execute_reply": "2020-05-13T01:48:59.383881Z",
     "shell.execute_reply.started": "2020-05-13T01:48:59.378775Z"
=======
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T22:21:42.625733Z",
     "iopub.status.busy": "2020-05-12T22:21:42.625139Z",
     "iopub.status.idle": "2020-05-12T22:21:42.642678Z",
     "shell.execute_reply": "2020-05-12T22:21:42.640279Z",
     "shell.execute_reply.started": "2020-05-12T22:21:42.625659Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "def angle(v1, v2):\n",
    "    \"\"\" Finds angel between 2 vectors\n",
    "\n",
    "    returns: ang , v1\"\"\"\n",
    "    \n",
    "    try:\n",
    "       \n",
    "        ang = math.atan2(np.linalg.norm(np.cross(v1,v2)),np.dot(v1,v2))\n",
    "                         \n",
    "        if ang > math.pi/2:\n",
    "            v1 = -v1\n",
    "            ang = math.atan2(np.linalg.norm(np.cross(v1,v2)),np.dot(v1,v2))\n",
    "            \n",
    "            print(f'{ang} PC inverted')\n",
    "\n",
    "        else:\n",
    "            print(f'{ang} no invert')\n",
    "\n",
    "    except:\n",
    "        #vang = 0\n",
    "        print(f'ERROR: vectors v1= {v1}, v2= {v2}')\n",
    "        ang = 'ERROR'\n",
    "\n",
    "    return ang, v1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bones are rotated with quaternions.\n",
    "\n",
    "The angle between the two PC1 vectors is taken. The object is then rotated (by a quaternion) around the cross product between the PC1 vectors.\n",
    "\n",
    "The new angles between the next PCs are calculates and the process is repeated for the other PCs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## quaternion_rotation"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T01:48:59.385457Z",
     "iopub.status.busy": "2020-05-13T01:48:59.385291Z",
     "iopub.status.idle": "2020-05-13T01:48:59.392954Z",
     "shell.execute_reply": "2020-05-13T01:48:59.392439Z",
     "shell.execute_reply.started": "2020-05-13T01:48:59.385436Z"
=======
   "execution_count": 97,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:23.141275Z",
     "iopub.status.busy": "2020-05-12T23:04:23.141026Z",
     "iopub.status.idle": "2020-05-12T23:04:23.145718Z",
     "shell.execute_reply": "2020-05-12T23:04:23.145049Z",
     "shell.execute_reply.started": "2020-05-12T23:04:23.141245Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "def quaternion_rotation_from_angle(v, c_axis, theta):\n",
    "    \n",
    "    \"\"\"\n",
    "    rotates vector around axis by theta\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    rotation_axis = np.array([0.] + c_axis)\n",
    "    axis_angle = (theta*0.5) * rotation_axis/np.linalg.norm(rotation_axis)\n",
    "\n",
    "    vec = np.quaternion(*v)\n",
    "\n",
    "    # quaternion from exp of axis angle\n",
    "    qlog = np.quaternion(*axis_angle)\n",
    "    q = np.exp(qlog)\n",
    "\n",
    "    # double cover quaternion rotation\n",
    "    v_prime = q * vec * np.conjugate(q)\n",
    "\n",
    "    return v_prime.imag , q"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T01:48:59.393835Z",
     "iopub.status.busy": "2020-05-13T01:48:59.393689Z",
     "iopub.status.idle": "2020-05-13T01:48:59.399968Z",
     "shell.execute_reply": "2020-05-13T01:48:59.399583Z",
     "shell.execute_reply.started": "2020-05-13T01:48:59.393816Z"
=======
   "execution_count": 98,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:23.538472Z",
     "iopub.status.busy": "2020-05-12T23:04:23.538308Z",
     "iopub.status.idle": "2020-05-12T23:04:23.541290Z",
     "shell.execute_reply": "2020-05-12T23:04:23.540729Z",
     "shell.execute_reply.started": "2020-05-12T23:04:23.538454Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "def quaternion_rotation_from_quaternion(v, q):\n",
<<<<<<< HEAD
    "    \"\"\"\n",
    "    rotates vector by quaternion\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    # double cover quaternion rotation\n",
    "    vec = np.quaternion(*v)\n",
=======
    "    # double cover quaternion rotation\n",
    "    vec = quat.quaternion(*v)\n",
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    "    \n",
    "    v_prime = q * vec * np.conjugate(q)\n",
    "    \n",
    "    return v_prime.imag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## rotate"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 39,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T02:01:27.427239Z",
     "iopub.status.busy": "2020-05-13T02:01:27.427057Z",
     "iopub.status.idle": "2020-05-13T02:01:27.435706Z",
     "shell.execute_reply": "2020-05-13T02:01:27.435165Z",
     "shell.execute_reply.started": "2020-05-13T02:01:27.427218Z"
=======
   "execution_count": 99,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:24.157541Z",
     "iopub.status.busy": "2020-05-12T23:04:24.157130Z",
     "iopub.status.idle": "2020-05-12T23:04:24.171313Z",
     "shell.execute_reply": "2020-05-12T23:04:24.170596Z",
     "shell.execute_reply.started": "2020-05-12T23:04:24.157497Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "def rotate(bone_f1, bone_f2, interpolate = False, scale_factor= 2):\n",
    "\n",
    "    if interpolate is True: \n",
    "        print(f'scalling bone by {scale_factor}')\n",
    "        bone_f1.scale(scale_factor)\n",
    "    \n",
    "    # center bones too 0,0,0,\n",
    "    bone_f1.center_to_origin()\n",
    "    bone_f2.center_to_origin()\n",
    "\n",
    "    # PCA on bones\n",
    "    bone_f1.get_pca()\n",
    "    bone_f2.get_pca()\n",
    "\n",
    "    # for 1 to 3 principle conponents of the object\n",
    "    for i in range(1, 4):\n",
    "\n",
    "        # takes cross product axis\n",
    "        cross_product_axis = np.cross(\n",
    "            getattr(bone_f1, f'pc{i}'),\n",
    "            getattr(bone_f2, f'pc{i}'))\n",
    "\n",
    "        # finds angle between PCs for f1 vs f2\n",
    "        theta, vector = angle(\n",
    "            getattr(bone_f1, f'pc{i}'),\n",
    "            getattr(bone_f2, f'pc{i}'))\n",
    "\n",
    "        # sets any new values needed\n",
    "        setattr(bone_f1, f'pc{i}', vector)\n",
    "\n",
    "        # rotates each PC of the bone\n",
    "        for n in range(1, 4):\n",
    "            transformed_pc, q = quaternion_rotation_from_angle(\n",
    "                v=getattr(bone_f1, f'pc{n}'),\n",
    "                c_axis=cross_product_axis,\n",
    "                theta=theta)\n",
    "            \n",
<<<<<<< HEAD
    "            # sets new PCA\n",
    "            setattr(bone_f1, f'pc{n}', transformed_pc)\n",
    "            \n",
    "        #logs quaternion\n",
    "        setattr(bone_f1,f'q_prod_{i}', q)\n",
    "        \n",
    "    q_total = getattr(bone_f1, 'q_prod_3') * getattr(bone_f1, 'q_prod_2') * getattr(bone_f1, 'q_prod_1')\n",
    "            \n",
    "            \n",
=======
    "            if hasattr(bone_f1,'q_prod') == True:\n",
    "                bone_f1.q_prod = bone_f1.q_prod * q\n",
    "                \n",
    "            else:\n",
    "                setattr(bone_f1,'q_prod', q)\n",
    "                \n",
    "            # sets new PCA\n",
    "            setattr(bone_f1, f'pc{n}', transformed_pc)\n",
    "\n",
    "\n",
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    "    # rotates xyz array with the quaterion product\n",
    "    rotated_xyz = np.apply_along_axis(\n",
    "        quaternion_rotation_from_quaternion,1, \n",
    "        getattr(bone_f1, 'xyz'),\n",
<<<<<<< HEAD
    "        q_total)\n",
=======
    "        bone_f1.q_prod)\n",
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    "        \n",
    "    setattr(bone_f1, 'xyz', rotated_xyz)\n",
    "\n",
    "    bone_f1.reset_position()\n",
    "    bone_f2.reset_position()\n",
    "\n",
    "    # reduce bone to orginal size \n",
    "    if interpolate is True: \n",
    "        print(f'scalling bone by {1/scale_factor}')\n",
    "        bone_f1.scale(1/scale_factor)   \n",
    "        \n",
    "\n",
    "    if bone_f1.dtype is 'stl':\n",
    "\n",
    "        #update internal data\n",
    "        bone_f1.data.v0 , bone_f1.data.v1, bone_f1.data.v2 = np.array_split(bone_f1.xyz, 3)\n",
    "        bone_f1.data.update_normals()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of Angles\n",
    "\n",
    "plots at pandas data frame of the angle differences between the two bone's PCAs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## df_angles"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T01:48:59.411922Z",
     "iopub.status.busy": "2020-05-13T01:48:59.411775Z",
     "iopub.status.idle": "2020-05-13T01:48:59.420295Z",
     "shell.execute_reply": "2020-05-13T01:48:59.419600Z",
     "shell.execute_reply.started": "2020-05-13T01:48:59.411903Z"
=======
   "execution_count": 100,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:25.285283Z",
     "iopub.status.busy": "2020-05-12T23:04:25.284558Z",
     "iopub.status.idle": "2020-05-12T23:04:25.300122Z",
     "shell.execute_reply": "2020-05-12T23:04:25.298855Z",
     "shell.execute_reply.started": "2020-05-12T23:04:25.285203Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "def df_angles(bone_f1, bone_f2, degrees=False, name='UN-NAMED BONE'):\n",
    "    \"\"\"\n",
    "    Compares the PCA angles between to bones.\n",
    "\n",
    "    Input:  bone_f1 = bone in 1st position\n",
    "            bone_f2 = bone in 2nd position\n",
    "\n",
    "    Returns: pandas dataframe\n",
    "    \"\"\"\n",
    "\n",
    "    df = pd.DataFrame()\n",
    "\n",
    "    \n",
    "    # Check for PCAs\n",
    "    if hasattr(bone_f1, 'pc1') is False:\n",
    "            bone_f1.get_pca()\n",
    "    \n",
    "    if hasattr(bone_f2, 'pc1') is False:\n",
    "            bone_f2.get_pca()\n",
    "            \n",
    "    # loops over each PCA\n",
    "    for n in range(1, 4):\n",
    "        theta, _ = angle(\n",
    "            getattr(bone_f1, f'pc{n}'), getattr(bone_f2, f'pc{n}'))\n",
    "        \n",
    "        if degrees is True: theta = math.degrees(theta)\n",
    "        \n",
    "        # Sets the column names\n",
    "        df.loc[f'{name} f1: pc{n}', f'{name} f2: pc{n}'] = theta\n",
    "\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to use:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load the data that you want to use\n",
    "*you have to use the absoslute file paths*"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 40,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T02:01:31.027130Z",
     "iopub.status.busy": "2020-05-13T02:01:31.026957Z",
     "iopub.status.idle": "2020-05-13T02:01:31.646247Z",
     "shell.execute_reply": "2020-05-13T02:01:31.645731Z",
     "shell.execute_reply.started": "2020-05-13T02:01:31.027110Z"
=======
   "execution_count": 101,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:26.147150Z",
     "iopub.status.busy": "2020-05-12T23:04:26.146981Z",
     "iopub.status.idle": "2020-05-12T23:04:26.787066Z",
     "shell.execute_reply": "2020-05-12T23:04:26.786561Z",
     "shell.execute_reply.started": "2020-05-12T23:04:26.147131Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "# load data\n",
    "tibia_f2 = bone.from_matlab_path(\n",
    "    matlab_file='/home/luke/ankle_data/phantom/phantom_tibia_f2.mat')\n",
    "\n",
    "tibia_f1 = bone.from_matlab_path(\n",
    "    matlab_file='/home/luke/ankle_data/phantom/phantom_tibia_f1.mat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Or STL files*"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 41,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T02:01:34.489996Z",
     "iopub.status.busy": "2020-05-13T02:01:34.489834Z",
     "iopub.status.idle": "2020-05-13T02:01:34.492368Z",
     "shell.execute_reply": "2020-05-13T02:01:34.491854Z",
     "shell.execute_reply.started": "2020-05-13T02:01:34.489978Z"
=======
   "execution_count": 102,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:26.788129Z",
     "iopub.status.busy": "2020-05-12T23:04:26.787953Z",
     "iopub.status.idle": "2020-05-12T23:04:26.790225Z",
     "shell.execute_reply": "2020-05-12T23:04:26.789798Z",
     "shell.execute_reply.started": "2020-05-12T23:04:26.788109Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "# tibia_stl = bone.from_stl_path('/home/luke/ankle_data/stl_foot/tibia.stl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Set custom filter level (optional)*\n",
    "only for voxel arrays"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 42,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T02:01:34.881182Z",
     "iopub.status.busy": "2020-05-13T02:01:34.880947Z",
     "iopub.status.idle": "2020-05-13T02:01:34.884047Z",
     "shell.execute_reply": "2020-05-13T02:01:34.883413Z",
     "shell.execute_reply.started": "2020-05-13T02:01:34.881158Z"
=======
   "execution_count": 103,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:27.033838Z",
     "iopub.status.busy": "2020-05-12T23:04:27.033490Z",
     "iopub.status.idle": "2020-05-12T23:04:27.038022Z",
     "shell.execute_reply": "2020-05-12T23:04:27.036994Z",
     "shell.execute_reply.started": "2020-05-12T23:04:27.033797Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "#bone.filter_level = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Set custom colour for bone (optional)*"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 43,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T02:01:35.669838Z",
     "iopub.status.busy": "2020-05-13T02:01:35.669630Z",
     "iopub.status.idle": "2020-05-13T02:01:35.672919Z",
     "shell.execute_reply": "2020-05-13T02:01:35.672241Z",
     "shell.execute_reply.started": "2020-05-13T02:01:35.669813Z"
=======
   "execution_count": 104,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:28.060377Z",
     "iopub.status.busy": "2020-05-12T23:04:28.059329Z",
     "iopub.status.idle": "2020-05-12T23:04:28.072691Z",
     "shell.execute_reply": "2020-05-12T23:04:28.070751Z",
     "shell.execute_reply.started": "2020-05-12T23:04:28.060239Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "tibia_f1.default_color = (0.8, 0.3, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of angles"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 44,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T02:01:36.547454Z",
     "iopub.status.busy": "2020-05-13T02:01:36.547271Z",
     "iopub.status.idle": "2020-05-13T02:01:36.549754Z",
     "shell.execute_reply": "2020-05-13T02:01:36.549276Z",
     "shell.execute_reply.started": "2020-05-13T02:01:36.547433Z"
=======
   "execution_count": 105,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:28.801701Z",
     "iopub.status.busy": "2020-05-12T23:04:28.801531Z",
     "iopub.status.idle": "2020-05-12T23:04:28.804478Z",
     "shell.execute_reply": "2020-05-12T23:04:28.803976Z",
     "shell.execute_reply.started": "2020-05-12T23:04:28.801675Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "# df_angles(tibia_f1, tibia_f2, degrees=True, name='tibia')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rotate the bone"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 45,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T02:01:37.105349Z",
     "iopub.status.busy": "2020-05-13T02:01:37.105176Z",
     "iopub.status.idle": "2020-05-13T02:01:37.223309Z",
     "shell.execute_reply": "2020-05-13T02:01:37.222688Z",
     "shell.execute_reply.started": "2020-05-13T02:01:37.105329Z"
=======
   "execution_count": 106,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:34.011413Z",
     "iopub.status.busy": "2020-05-12T23:04:34.011077Z",
     "iopub.status.idle": "2020-05-12T23:04:34.130655Z",
     "shell.execute_reply": "2020-05-12T23:04:34.130226Z",
     "shell.execute_reply.started": "2020-05-12T23:04:34.011373Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.17471680699860007 no invert\n",
      "0.5819453559235347 no invert\n",
      "3.906788298496456e-16 no invert\n"
     ]
    }
   ],
   "source": [
    "rotate(tibia_f1, tibia_f2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 46,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T02:01:37.852373Z",
     "iopub.status.busy": "2020-05-13T02:01:37.852202Z",
     "iopub.status.idle": "2020-05-13T02:01:48.943780Z",
     "shell.execute_reply": "2020-05-13T02:01:48.942817Z",
     "shell.execute_reply.started": "2020-05-13T02:01:37.852354Z"
=======
   "execution_count": 107,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T23:04:46.639076Z",
     "iopub.status.busy": "2020-05-12T23:04:46.638907Z",
     "iopub.status.idle": "2020-05-12T23:04:49.218504Z",
     "shell.execute_reply": "2020-05-12T23:04:49.217403Z",
     "shell.execute_reply.started": "2020-05-12T23:04:46.639058Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "tibia_f2.plot()\n",
    "tibia_f1.plot()\n",
    "mlab.show()"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 47,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T02:01:48.946064Z",
     "iopub.status.busy": "2020-05-13T02:01:48.945663Z",
     "iopub.status.idle": "2020-05-13T02:01:48.951554Z",
     "shell.execute_reply": "2020-05-13T02:01:48.950565Z",
     "shell.execute_reply.started": "2020-05-13T02:01:48.946013Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "quaternion(1, 0, 0, 0)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tibia_f1.q_prod_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-13T01:44:52.627487Z",
     "iopub.status.busy": "2020-05-13T01:44:52.627230Z",
     "iopub.status.idle": "2020-05-13T01:44:52.630266Z",
     "shell.execute_reply": "2020-05-13T01:44:52.629768Z",
     "shell.execute_reply.started": "2020-05-13T01:44:52.627459Z"
=======
   "execution_count": 95,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-05-12T22:10:40.616182Z",
     "iopub.status.busy": "2020-05-12T22:10:40.615962Z",
     "iopub.status.idle": "2020-05-12T22:10:40.619046Z",
     "shell.execute_reply": "2020-05-12T22:10:40.618476Z",
     "shell.execute_reply.started": "2020-05-12T22:10:40.616148Z"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
    }
   },
   "outputs": [],
   "source": [
    "# optional mesh ploting for stl objects\n",
    "# tibia_stl.plot(mesh=True)\n",
    "# mlab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Voxel reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "85.37965945697192% reconstructed\n"
     ]
    }
   ],
   "source": [
    "vx_array = tibia_f1.xyz_to_array()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Validate the voxel array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, z = np.where(vx_array == True)\n",
    "\n",
    "recon_xyz = np.array([x, y, z]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlab.points3d(recon_xyz[:, 0],\n",
    "                      recon_xyz[:, 1],\n",
    "                      recon_xyz[:, 2],\n",
    "                      mode = \"cube\",\n",
    "                      scale_factor = 1) \n",
    "mlab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STL file out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://numpy-stl.readthedocs.io/en/latest/stl.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tibia_stl.data.save()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
<<<<<<< HEAD
   "display_name": "Python [conda env:vox]",
   "language": "python",
   "name": "conda-env-vox-py"
=======
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
>>>>>>> 35409e49361ce24494ea4506a39098928bd7d02d
  },
  "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.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "308px",
    "width": "236.667px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "65.8167px",
    "left": "1318.05px",
    "top": "110.633px",
    "width": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  },
  "toc-autonumbering": false,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}