skeleton_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import skan\n",
    "from scipy import ndimage as ndi\n",
    "from tifffile import imsave\n",
    "\n",
    "def dendrite_pruning(skeleton_image,distance_transform):\n",
    "    n_path=-1\n",
    "    pruned_image=np.zeros([3,np.shape(skeleton_image)[0],np.shape(skeleton_image)[1],np.shape(skeleton_image)[2]])\n",
    "    skeleton_image[distance_transform==1] = 0\n",
    "    while n_path!=0:\n",
    "        skeleton_object = skan.Skeleton(skeleton_image)\n",
    "        summary=skan.summarize(skeleton_object)\n",
    "        paths=skeleton_object.paths_list()\n",
    "        \n",
    "        path_num=len(paths)\n",
    "        print(path_num,n_path)\n",
    "        if(path_num==n_path):\n",
    "            pruned_image[0,:,:,:] = skeleton_image\n",
    "            pruned_image[2,:,:,:] = distance_transform\n",
    "            for i in range(len(src)):\n",
    "                pruned_image[1,coordinates[src[i],0],coordinates[src[i],1],coordinates[src[i],2]]=euclidean_distance[i]\n",
    "                pruned_image[1,coordinates[dst[i],0],coordinates[dst[i],1],coordinates[dst[i],2]]=euclidean_distance[i]\n",
    "            imsave('images/pruned_skeleton.tif', np.einsum('ijkl->lijk', pruned_image))\n",
    "            return skeleton_object\n",
    "        else:\n",
    "            n_path=len(paths)\n",
    "        \n",
    "        t1_branch = np.where(summary['branch-type']==1)[0]\n",
    "        t2_branch = np.where(summary['branch-type']==2)[0]\n",
    "        src=[i[0] for i in paths]\n",
    "        dst=[i[-1] for i in paths]\n",
    "\n",
    "        coordinates = np.array((skeleton_object.coordinates),dtype=int)\n",
    "        euclidean_distance = np.square(coordinates[src,0]-coordinates[dst,0]) +  np.square(coordinates[src,1]-coordinates[dst,1]) +  np.square(coordinates[src,2]-coordinates[dst,2])\n",
    "        local_distance_dst = distance_transform[coordinates[dst,0],coordinates[dst,1],coordinates[dst,2]]\n",
    "        local_distance_src = distance_transform[coordinates[src,0],coordinates[src,1],coordinates[src,2]]\n",
    "\n",
    "        skeleton_image=np.zeros(np.shape(skeleton_image))\n",
    "        for i in t1_branch:\n",
    "            if(euclidean_distance[i] < 2*(local_distance_dst[i]**2) or euclidean_distance[i] < 2*(local_distance_src[i]**2) or src[i]==dst[i]):\n",
    "                continue\n",
    "            path=skeleton_object.path(i)\n",
    "            #pruned_image[1,coordinates[path[0],0],coordinates[path[0],1],coordinates[path[0],2]] = pruned_image[1,coordinates[path[0],0],coordinates[path[0],1],coordinates[path[0],2]]+1\n",
    "            #pruned_image[1,coordinates[path[-1],0],coordinates[path[-1],1],coordinates[path[-1],2]] = pruned_image[1,coordinates[path[-1],0],coordinates[path[-1],1],coordinates[path[-1],2]]+1\n",
    "            for j in path:\n",
    "                skeleton_image[coordinates[j,0],coordinates[j,1],coordinates[j,2]] = 1\n",
    "\n",
    "        for i in t2_branch:\n",
    "            if(src[i]==dst[i]):\n",
    "                continue\n",
    "            path=skeleton_object.path(i)\n",
    "            #pruned_image[1,coordinates[path[0],0],coordinates[path[0],1],coordinates[path[0],2]] = pruned_image[1,coordinates[path[0],0],coordinates[path[0],1],coordinates[path[0],2]]+1\n",
    "            #pruned_image[1,coordinates[path[-1],0],coordinates[path[-1],1],coordinates[path[-1],2]] = pruned_image[1,coordinates[path[-1],0],coordinates[path[-1],1],coordinates[path[-1],2]]+1\n",
    "            for j in path:\n",
    "                skeleton_image[coordinates[j,0],coordinates[j,1],coordinates[j,2]] = 2\n",
    "\n",
    "\n",
    "def soma_prunning(shape,skeleton_object,distance_transform):\n",
    "    n_path = 0\n",
    "    dleta_n_path = -1\n",
    "    skeleton_image = np.zeros(shape,dtype=bool)\n",
    "    while dleta_n_path!=0:\n",
    "        summary=skan.summarize(skeleton_object)\n",
    "        dleta_n_path = n_path - len(summary.index)\n",
    "        print(len(summary.index))\n",
    "\n",
    "        #if(len(summary.index)==2):\n",
    "        #    print(summary)\n",
    "        #    imsave('images/soma_skeleton.tif',np.einsum('ijk->kij', skeleton_image))\n",
    "        n_path = len(summary.index)\n",
    "        if(len(summary.index)<=1 or dleta_n_path==0):\n",
    "            #np.save('np arrays/soma_skeleton', skeleton_image)\n",
    "            return skeleton_image\n",
    "        coordinates = np.array((skeleton_object.coordinates),dtype=int)\n",
    "\n",
    "        skeleton_image = np.zeros(np.shape(skeleton_image),dtype=bool)\n",
    "        nodes = np.zeros(len(skeleton_object.coordinates))\n",
    "        paths = skeleton_object.paths_list()\n",
    "        for i in range(len(paths)):\n",
    "            nodes[paths[i][0]] = True\n",
    "            nodes[paths[i][-1]] = True\n",
    "        nodes_index = np.where(nodes==True)[0]\n",
    "        for node in nodes_index:\n",
    "            type1_branches=[]\n",
    "            n_type2_branches=0\n",
    "            n_type3_branches=0\n",
    "            euclidean_dist=[]\n",
    "            for i in summary.index:\n",
    "                if(paths[i][0]==node or paths[i][-1]==node):\n",
    "                    if(summary['branch-type'][i]==2):\n",
    "                        n_type2_branches = n_type2_branches + 1\n",
    "                    elif(summary['branch-type'][i]==1):\n",
    "                        type1_branches.append(i)\n",
    "                        euclidean_dist.append(summary['euclidean-distance'][i])\n",
    "                    elif(summary['branch-type'][i]==3):\n",
    "                        n_type3_branches=n_type3_branches+1\n",
    "            n_type1_branches=len(type1_branches)\n",
    "            if(n_type1_branches==1 and n_type2_branches==1):\n",
    "                path = paths[type1_branches[0]]\n",
    "                for j in range(len(path)):\n",
    "                    skeleton_image[coordinates[path[j],0],coordinates[path[j],1],coordinates[path[j],2]] = True \n",
    "            elif(n_type1_branches==2 and n_type2_branches==0):\n",
    "                i_max = np.argmax(euclidean_dist)\n",
    "                path = paths[type1_branches[i_max]]\n",
    "                for j in range(len(path)):\n",
    "                    skeleton_image[coordinates[path[j],0],coordinates[path[j],1],coordinates[path[j],2]] = True  \n",
    "\n",
    "                type1_branches.pop(i_max)\n",
    "                euclidean_dist.pop(i_max)\n",
    "                i_max = np.argmax(euclidean_dist)\n",
    "                path = paths[type1_branches[i_max]]\n",
    "                for j in range(len(path)):\n",
    "                    skeleton_image[coordinates[path[j],0],coordinates[path[j],1],coordinates[path[j],2]] = True\n",
    "            elif(n_type2_branches>1):\n",
    "                for i_branch in type1_branches:\n",
    "                    path = paths[i_branch]\n",
    "                    for j in range(len(path)):\n",
    "                        skeleton_image[coordinates[path[j],0],coordinates[path[j],1],coordinates[path[j],2]] = True\n",
    "            elif((n_type1_branches>1 and n_type2_branches==1) or (n_type1_branches>0 and n_type3_branches!=0)):\n",
    "                i_max = np.argmax(euclidean_dist)\n",
    "                path = paths[type1_branches[i_max]]\n",
    "                for j in range(len(path)):\n",
    "                    skeleton_image[coordinates[path[j],0],coordinates[path[j],1],coordinates[path[j],2]] = True\n",
    "\n",
    "        for i in summary.index:\n",
    "            if(summary['branch-type'][i]==2):\n",
    "                if(paths[i][0]==paths[i][-1]):\n",
    "                    print('cycle found in soma skeleton')\n",
    "                    return skeleton_image\n",
    "                for j in range(len(paths[i])):\n",
    "                    skeleton_image[coordinates[paths[i][j],0],coordinates[paths[i][j],1],coordinates[paths[i][j],2]] = True\n",
    "        skeleton_object = skan.Skeleton(skeleton_image)\n",
    "\n",
    "\n",
    "def set_biophysics(section,type):\n",
    "    section.insert('pas')\n",
    "    section.Ra = 100\n",
    "    section.cm = 1\n",
    "    section.e_pas=-90\n",
    "    section.insert('Ca_LVAst')\n",
    "    section.insert('Ca_HVA')\n",
    "    section.insert('SK_E2')\n",
    "    section.insert('SKv3_1')\n",
    "    section.insert('Im')\n",
    "    section.insert('NaTa_t')\n",
    "    section.insert('CaDynamics_E2')\n",
    "    section.insert('Ih')\n",
    "    section.insert('K_Tst')\n",
    "    section.insert('K_Pst')\n",
    "    if type=='axon':\n",
    "        section.insert('Nap_Et2')\n",
    "        for s,seg in enumerate(section):\n",
    "            seg.ek = -85\n",
    "            seg.ena = 50\n",
    "            seg.decay_CaDynamics_E2 = 277.300774\n",
    "            seg.gamma_CaDynamics_E2 = 0.000525\n",
    "            seg.gSK_E2bar_SK_E2 = 0.000047\n",
    "            seg.gSKv3_1bar_SKv3_1 = 0.473799\n",
    "            seg.gNaTa_tbar_NaTa_t = 3.89618\n",
    "            seg.gImbar_Im = 0.013322\n",
    "            seg.g_pas = 3e-5\n",
    "            seg.gCa_HVAbar_Ca_HVA = 0.000222\n",
    "            seg.gCa_LVAstbar_Ca_LVAst = 0.000813\n",
    "            seg.gIhbar_Ih = 0.0001\n",
    "            seg.gK_Tstbar_K_Tst = 0.077274\n",
    "            seg.gK_Pstbar_K_Pst = 0.188851\n",
    "            seg.gNap_Et2bar_Nap_Et2 = 0.005834\n",
    "    elif type=='soma':\n",
    "        section.insert('Nap_Et2')\n",
    "        for i,seg in enumerate(soma):\n",
    "            seg.ek = -85\n",
    "            seg.ena = 50\n",
    "            seg.gIhbar_Ih = 0.0002\n",
    "            seg.g_pas = 0.0000344\n",
    "            seg.decay_CaDynamics_E2 = 486.0\n",
    "            seg.gamma_CaDynamics_E2 = 0.000549\n",
    "            seg.gCa_LVAstbar_Ca_LVAst = 0.00432\n",
    "            seg.gCa_HVAbar_Ca_HVA = 0.000567\n",
    "            seg.gSKv3_1bar_SKv3_1 = 0.766\n",
    "            seg.gSK_E2bar_SK_E2 = 0.0556\n",
    "            seg.gK_Tstbar_K_Tst = 0.0326\n",
    "            seg.gK_Pstbar_K_Pst = 0.000547\n",
    "            seg.gNap_Et2bar_Nap_Et2 = 0.00496\n",
    "            seg.gNaTa_tbar_NaTa_t = 1.71\n",
    "            seg.gImbar_Im = 0.000008\n",
    "    elif type=='dend':\n",
    "        for s,seg in enumerate(section):\n",
    "            seg.ek = -85\n",
    "            seg.ena = 50\n",
    "            seg.decay_CaDynamics_E2 = 88.9\n",
    "            seg.gamma_CaDynamics_E2 = 0.0005\n",
    "            seg.gSK_E2bar_SK_E2 = 0.0186*10\n",
    "            seg.gSKv3_1bar_SKv3_1 = 0.00298*10\n",
    "            seg.gNaTa_tbar_NaTa_t = 0.0211\n",
    "            seg.gImbar_Im = 0.00006\n",
    "            seg.g_pas = 0.0000447\n",
    "            seg.gCa_HVAbar_Ca_HVA = 0.0000701\n",
    "            seg.gCa_LVAstbar_Ca_LVAst = 0.0001419540000\n",
    "            seg.gIhbar_Ih = 0.0001\n",
    "            seg.gK_Tstbar_K_Tst = 0.0326\n",
    "            seg.gK_Pstbar_K_Pst = 0.000547\n",
    "    else:\n",
    "        print('undefined branch type')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "21\n",
      "15\n",
      "13\n",
      "11\n",
      "9\n",
      "7\n",
      "5\n",
      "3\n",
      "3\n"
     ]
    }
   ],
   "source": [
    "distance_transform = ndi.distance_transform_edt(np.load('np arrays/foreground.npy'))\n",
    "soma_skeleton_image = np.load('np arrays/soma_skeleton.npy') #dendrites\n",
    "skeleton_object = skan.Skeleton(soma_skeleton_image)\n",
    "\n",
    "soma_image = np.zeros([3,np.shape(soma_skeleton_image)[0],np.shape(soma_skeleton_image)[1],np.shape(soma_skeleton_image)[2]])\n",
    "soma_image[0,:,:,:] = soma_skeleton_image\n",
    "soma_image[1,:,:,:] = soma_prunning(np.shape(soma_skeleton_image),skeleton_object,distance_transform)\n",
    "soma_image[2,:,:,:] = np.load('np arrays/soma_mask.npy')\n",
    "imsave('images/soma_skeleton.tif',np.einsum('ijkl->lijk', soma_image))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma_skeleton = skan.Skeleton(soma_image[1,:,:,:])\n",
    "soma_path = soma_skeleton.path(0)\n",
    "coordinates = soma_skeleton.coordinates\n",
    "soma_coordinates = np.array([coordinates[soma_path,0],coordinates[soma_path,1],coordinates[soma_path,2]],dtype=int)\n",
    "soma_diameter = np.array(distance_transform[soma_coordinates[0,:],soma_coordinates[1,:],soma_coordinates[2,:]])\n",
    "soma_coordinates = (soma_coordinates.astype(float))\n",
    "soma_diameter = (soma_diameter.astype(float)*1.5)\n",
    "\n",
    "plot_diameter = soma_diameter.copy()\n",
    "plot_coordinates = soma_coordinates.copy()\n",
    "init_potential = np.zeros(len(soma_diameter))+20\n",
    "connections = np.zeros([len(plot_diameter)-1,2])\n",
    "connections[:,0]=np.linspace(1,len(plot_diameter)-1,len(plot_diameter)-1)\n",
    "connections[:,1]=np.linspace(0,len(plot_diameter)-2,len(plot_diameter)-1)\n",
    "del soma_skeleton_image,skeleton_object,soma_image\n",
    "from neuron import h\n",
    "import numpy as np\n",
    "from neuron import gui\n",
    "\n",
    "h.nrn_load_dll('biophysical_mods/nrnmech.dll')\n",
    "soma = h.Section(name='soma')\n",
    "soma.nseg = int(len(soma_path)) + 1\n",
    "h.pt3dclear(sec=soma)\n",
    "h.pt3dadd(h.Vector(soma_coordinates[0,:]), h.Vector(soma_coordinates[1,:]), h.Vector(soma_coordinates[2,:]), h.Vector(soma_diameter), sec=soma)\n",
    "set_biophysics(soma,'soma')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6672 -1\n",
      "4237 6672\n",
      "4120 4237\n",
      "4109 4120\n",
      "4099 4109\n",
      "4099 4099\n"
     ]
    }
   ],
   "source": [
    "distance_transform = ndi.distance_transform_edt(np.load('np arrays/foreground.npy'))\n",
    "skeleton_image = np.array(np.load('np arrays/skeleton_image.npy'),dtype=int) #dendrites\n",
    "skeleton_object = dendrite_pruning(skeleton_image,distance_transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cross detection\n",
    "coordinates = np.array((skeleton_object.coordinates),dtype=int)\n",
    "nodes = np.zeros(len(coordinates))\n",
    "paths = skeleton_object.paths_list()\n",
    "for i in range(len(paths)):\n",
    "    nodes[paths[i][0]] = 1\n",
    "    nodes[paths[i][-1]] = 1\n",
    "nodes_index = np.where(nodes==1)\n",
    "\n",
    "for id in nodes_index[0]:\n",
    "    x = coordinates[id,0]\n",
    "    y = coordinates[id,1]\n",
    "    z = coordinates[id,2]\n",
    "    local_foreground_range = (distance_transform[x,y,z]*1.5)**2\n",
    "    bifercations = []\n",
    "    orientation = []\n",
    "    inverse = []\n",
    "    for i in range(len(paths)):\n",
    "        if(paths[i][0]==id):\n",
    "            the_other_end = -1\n",
    "            local_path = paths[i]\n",
    "            inverse.append(False)\n",
    "        elif(paths[i][-1]==id):\n",
    "            the_other_end = 0\n",
    "            local_path = paths[i][::-1]\n",
    "            inverse.append(True)\n",
    "        else:\n",
    "            continue\n",
    "        \n",
    "        bifercations.append(i)\n",
    "        short_branch = True\n",
    "        for coord in local_path:\n",
    "            vx = coordinates[coord,0]-x\n",
    "            vy = coordinates[coord,1]-y\n",
    "            vz = coordinates[coord,2]-z\n",
    "            if(vx**2+vy**2+vz**2>local_foreground_range):\n",
    "                orientation.append([vx,vy,vz])\n",
    "                short_branch = False\n",
    "                break\n",
    "        if(short_branch):\n",
    "            orientation.append([coordinates[paths[i][the_other_end],0] - x,\n",
    "            coordinates[paths[i][the_other_end],1] - y,\n",
    "            coordinates[paths[i][the_other_end],2] - z])\n",
    "    N_bifercation = len(bifercations)\n",
    "    if(N_bifercation<4):\n",
    "        if(N_bifercation==2):\n",
    "            if (inverse[0] and (not inverse[1])):\n",
    "                paths[bifercations[0]].extend(paths[bifercations[1]][1:])\n",
    "                paths.pop(bifercations[1])\n",
    "            elif (inverse[0] and inverse[1]):\n",
    "                paths[bifercations[0]].extend(paths[bifercations[1]][-2::-1])\n",
    "                paths.pop(bifercations[1])\n",
    "            elif ((not inverse[0]) and inverse[1]):\n",
    "                paths[bifercations[1]].extend(paths[bifercations[0]][1:])\n",
    "                paths.pop(bifercations[0])\n",
    "            elif ((not inverse[0]) and (not inverse[1])):\n",
    "                paths.append([*paths[bifercations[0]][-1:0:-1],*paths[bifercations[1]]])\n",
    "                paths.pop(bifercations[1])\n",
    "                paths.pop(bifercations[0])\n",
    "        continue\n",
    "    cross_products=np.zeros([N_bifercation,N_bifercation])-1\n",
    "    for r in range(N_bifercation-1):\n",
    "        for c in range(r+1,N_bifercation,1):\n",
    "            cross_products[r,c]=np.abs(np.dot(orientation[c],orientation[r]))\n",
    "    left_branch = 0\n",
    "    if(N_bifercation%2==1):\n",
    "        left_branch=3\n",
    "    to_be_poped=[]\n",
    "    while np.sum(cross_products>0)>left_branch:\n",
    "        r,c = np.unravel_index(np.argmax(cross_products),cross_products.shape)\n",
    "        cross_products[r,:] = -1\n",
    "        cross_products[:,c] = -1\n",
    "        if (inverse[r] and (not inverse[c]) and (not paths[bifercations[r]][0]==paths[bifercations[c]][-1])):\n",
    "            paths[bifercations[r]].extend(paths[bifercations[c]][1:])\n",
    "            to_be_poped.append(c)\n",
    "        elif (inverse[r] and inverse[c] and (not paths[bifercations[r]][0]==paths[bifercations[c]][0])):\n",
    "            paths[bifercations[r]].extend(paths[bifercations[c]][-2::-1])\n",
    "            to_be_poped.append(c)\n",
    "        elif ((not inverse[r]) and inverse[c] and (not paths[bifercations[r]][-1]==paths[bifercations[c]][0])):\n",
    "            paths[bifercations[c]].extend(paths[bifercations[r]][1:])\n",
    "            to_be_poped.append(c)\n",
    "        elif ((not inverse[r]) and (not inverse[c]) and (not paths[bifercations[r]][-1]==paths[bifercations[c]][-1])):\n",
    "            paths.append([*paths[bifercations[r]][-1:0:-1],*paths[bifercations[c]]])\n",
    "            if(c<r):\n",
    "                print('error')\n",
    "            to_be_poped.append(c)\n",
    "            to_be_poped.append(r)\n",
    "        else:\n",
    "            print(\"error: path mismatch or loop occurred\")\n",
    "    to_be_poped.sort(reverse=True)\n",
    "    for b in to_be_poped:\n",
    "        paths.pop(bifercations[b])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#cross test\n",
    "nodes = np.zeros(len(coordinates))\n",
    "for i in range(len(paths)):\n",
    "    nodes[paths[i][0]] = 1\n",
    "    nodes[paths[i][-1]] = 1\n",
    "nodes_index = np.where(nodes==1)\n",
    "\n",
    "for id in nodes_index[0]:\n",
    "    x = coordinates[id,0]\n",
    "    y = coordinates[id,1]\n",
    "    z = coordinates[id,2]\n",
    "    n=0\n",
    "    f=-1\n",
    "    s=-1\n",
    "    for i in range(len(paths)):\n",
    "        if(paths[i][0]==id or paths[i][-1]==id):\n",
    "            n=n+1\n",
    "            if(n==1):\n",
    "                f=i\n",
    "            elif(n==2):\n",
    "                s=i\n",
    "            else:\n",
    "                break\n",
    "    if(n==2):\n",
    "        print(coordinates[paths[f][0],0],coordinates[paths[f][-1],1],coordinates[paths[f][-1],2],\n",
    "        coordinates[paths[f][-1],0],coordinates[paths[f][-1],1],coordinates[paths[f][-1],2],\n",
    "        coordinates[paths[s][0],0],coordinates[paths[s][0],1],coordinates[paths[s][0],2],\n",
    "        coordinates[paths[s][-1],0],coordinates[paths[s][-1],1],coordinates[paths[s][-1],2],)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "coordinates2path_id = np.zeros([len(coordinates),3],dtype=int)-1\n",
    "direction_indicator = np.zeros([len(coordinates),3],dtype=int)\n",
    "\n",
    "def set_coordinates2path_id_map(node,i,direction):\n",
    "        if(coordinates2path_id[node,0]==-1):\n",
    "            coordinates2path_id[node,0]=i\n",
    "            direction_indicator[node,0]=direction\n",
    "        elif(coordinates2path_id[node,1]==-1):\n",
    "            coordinates2path_id[node,1]=i\n",
    "            direction_indicator[node,1]=direction\n",
    "        elif(coordinates2path_id[node,2]==-1):\n",
    "            coordinates2path_id[node,2]=i\n",
    "            direction_indicator[node,2]=direction\n",
    "        else:\n",
    "            print('error: more than 3 branches on node')\n",
    "for i in range(len(paths)):\n",
    "    set_coordinates2path_id_map(paths[i][0],i,1)\n",
    "    set_coordinates2path_id_map(paths[i][-1],i,-1)\n",
    "\n",
    "del skeleton_image,skeleton_object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#2 connection test\n",
    "if(not(len(np.where(np.sum(coordinates2path_id==-1,axis=1)==1)[0])==0)):\n",
    "    print(coordinates2path_id[np.where(np.sum(coordinates2path_id==-1,axis=1)==1)[0],:])\n",
    "    print('warning: node connects to 2 paths')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 601  637   -1]\n",
      " [3680 3694   -1]\n",
      " [3669 3695   -1]]\n",
      "warning: node connects to 2 paths\n"
     ]
    }
   ],
   "source": [
    "nodes_location = np.where(np.logical_not((direction_indicator==0).all(axis=1)))[0]\n",
    "node2path = coordinates2path_id[nodes_location,:]\n",
    "del coordinates2path_id\n",
    "direction_indicator = direction_indicator[nodes_location,:]\n",
    "nodes_corrdinates = coordinates[nodes_location,:]\n",
    "\n",
    "path2node = np.zeros([len(paths),2])-1\n",
    "for i in range(np.shape(node2path)[0]):\n",
    "    if(direction_indicator[i,0]==1):\n",
    "        path2node[node2path[i,0],0] = i\n",
    "    elif(direction_indicator[i,0]==-1):\n",
    "        path2node[node2path[i,0],1] = i\n",
    "    \n",
    "    if(direction_indicator[i,1]==1):\n",
    "        path2node[node2path[i,1],0] = i\n",
    "    elif(direction_indicator[i,1]==-1):\n",
    "        path2node[node2path[i,1],1] = i\n",
    "\n",
    "    if(direction_indicator[i,2]==1):\n",
    "        path2node[node2path[i,2],0] = i\n",
    "    elif(direction_indicator[i,2]==-1):\n",
    "        path2node[node2path[i,2],1] = i\n",
    "\n",
    "if(not(len(np.where(np.sum(node2path==-1,axis=1)==1)[0])==0)):\n",
    "    print(node2path[np.where(np.sum(node2path==-1,axis=1)==1)[0],:])\n",
    "    print('warning: node connects to 2 paths')\n",
    "\n",
    "root_ind=np.where(np.logical_and(node2path[:,0]>-1,node2path[:,1]==-1))[0]\n",
    "root = node2path[root_ind]\n",
    "# root_coordinates = nodes_corrdinates[root_ind] #use soma mask to check if root is connected to soma, trade huge storage for a bit of speed\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dendrite root reconstruction\n",
    "from scipy import ndimage as ndi\n",
    "soma_mask = ndi.binary_dilation(np.load('np arrays/soma_mask.npy'),iterations=5)\n",
    "dend = [h.Section(name='dend[%d]' % ind) for ind in range(len(paths))]\n",
    "path_front_list = []\n",
    "path_front_direction = []\n",
    "path_front_array = node2path[root_ind,0]\n",
    "reconstructed = np.zeros(len(paths),dtype=bool)\n",
    "n_path_per_node = []\n",
    "last_dend = []\n",
    "section_index=[]\n",
    "is_axon = np.zeros(len(paths),dtype=bool)\n",
    "\n",
    "src = [distance_transform[coordinates[paths[node2path[ind,0]][1],0],coordinates[paths[node2path[ind,0]][1],1],coordinates[paths[node2path[ind,0]][1],2]] for ind in root_ind] \n",
    "dst = [distance_transform[coordinates[paths[node2path[ind,0]][-1],0],coordinates[paths[node2path[ind,0]][-1],1],coordinates[paths[node2path[ind,0]][-1],2]] for ind in root_ind] \n",
    "is_axon[np.argmax([src,dst])%len(root_ind)] = True\n",
    "del src,dst\n",
    "\n",
    "\n",
    "for i,ind in enumerate(root_ind):\n",
    "    if soma_mask[coordinates[nodes_location[ind],0],coordinates[nodes_location[ind],1],coordinates[nodes_location[ind],2]]==False:\n",
    "        continue\n",
    "    i_path = node2path[ind,0]\n",
    "    if reconstructed[i_path]==True:\n",
    "        continue\n",
    "    if not len(section_index)==len(np.where(reconstructed==True)[0]):\n",
    "        print(len(section_index), len(np.where(reconstructed==True)[0]),i,ind)\n",
    "    reconstructed[i_path] = True\n",
    "    the_other_end = 0\n",
    "    if(direction_indicator[ind,0]==1):\n",
    "        path = np.array(paths[i_path],dtype=int)\n",
    "        the_other_end = path2node[i_path,1]\n",
    "    elif(direction_indicator[i,0]==-1):\n",
    "        path = np.array(paths[i_path][-1::-1],dtype=int)\n",
    "        the_other_end = path2node[i_path,0]\n",
    "    n_path = 0\n",
    "    for j,n in enumerate(node2path[int(the_other_end),:]):\n",
    "            if(n==-1):\n",
    "                continue\n",
    "            n_path = n_path + 1\n",
    "            path_front_list.append(n)\n",
    "            path_front_direction.append(direction_indicator[int(the_other_end),j])\n",
    "    n_path_per_node.append(n_path)\n",
    "    \n",
    "    path_coordinates = np.array([coordinates[path,0],coordinates[path,1],coordinates[path,2]],dtype=int)\n",
    "    diameter = distance_transform[path_coordinates[0,:],path_coordinates[1,:],path_coordinates[2,:]]\n",
    "    path_coordinates = (path_coordinates.astype(float))\n",
    "    diameter = (np.array(diameter).astype(float))\n",
    "    # mayavi plotting\n",
    "    dend_connections = np.zeros([len(diameter)-1,2])\n",
    "    dend_connections[:,0]=np.linspace(1,len(diameter)-1,len(diameter)-1) + len(plot_diameter)\n",
    "    dend_connections[:,1]=np.linspace(0,len(diameter)-2,len(diameter)-1) + len(plot_diameter)\n",
    "    connections = np.append(connections,dend_connections,axis=0)\n",
    "    plot_coordinates = np.append(plot_coordinates,path_coordinates,axis=1)\n",
    "    plot_diameter = np.append(plot_diameter,diameter)\n",
    "    init_potential = np.append(init_potential,np.zeros(len(diameter)))\n",
    "    section_index.append(i_path)\n",
    "    \n",
    "    #h.pt3dclear(sec=dend[i_path])\n",
    "    dend[i_path].nseg = int(len(diameter)) + 1\n",
    "    dend[i_path].Ra = 100\n",
    "    path_coordinates[0,:] = path_coordinates[0,:]-path_coordinates[0,0]\n",
    "    path_coordinates[1,:] = path_coordinates[1,:]-path_coordinates[1,0]\n",
    "    path_coordinates[2,:] = path_coordinates[2,:]-path_coordinates[2,0]\n",
    "    # add morphology\n",
    "    h.pt3dadd(h.Vector(path_coordinates[0,:]), h.Vector(path_coordinates[1,:]), h.Vector(path_coordinates[2,:]), h.Vector(diameter), sec=dend[i_path])\n",
    "    dist = (soma.x3d(0)-dend[i_path].x3d(0))**2 + (soma.y3d(0)-dend[i_path].y3d(0))**2 + (soma.z3d(0)-dend[i_path].z3d(0))**2\n",
    "    connect_location = 0\n",
    "    for arc in range(soma.n3d()):\n",
    "        dist0 = (soma.x3d(arc)-dend[i_path].x3d(0))**2 + (soma.y3d(arc)-dend[i_path].y3d(0))**2 + (soma.z3d(arc)-dend[i_path].z3d(0))**2\n",
    "        if(dist>dist0):\n",
    "            connect_location = arc\n",
    "            dist = dist0\n",
    "    last_dend.append(dend[i_path])\n",
    "    dend[i_path].connect(soma(float(connect_location)/soma.n3d()))\n",
    "\n",
    "    # set biophysics\n",
    "    if(is_axon[i_path]):\n",
    "        set_biophysics(dend[i_path],'axon')\n",
    "    else:\n",
    "        set_biophysics(dend[i_path],'dend')\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dendrite reconstruction\n",
    "\n",
    "while len(path_front_list)>0:\n",
    "    if not len(section_index)==len(np.where(reconstructed==True)[0]):\n",
    "        print(len(section_index), len(np.where(reconstructed==True)[0]),i,ind)\n",
    "    i_path = path_front_list[0]\n",
    "    if(reconstructed[i_path] == True):\n",
    "        path_front_list.pop(0)\n",
    "        path_front_direction.pop(0)\n",
    "        n_path_per_node[0] = n_path_per_node[0]-1\n",
    "        if(n_path_per_node[0]==0):\n",
    "            n_path_per_node.pop(0)\n",
    "            last_dend.pop(0)\n",
    "        continue\n",
    "    reconstructed[i_path] = True\n",
    "    the_other_end = 0\n",
    "\n",
    "    if(path_front_direction[0]==1):\n",
    "        path = paths[i_path]\n",
    "        the_other_end = path2node[i_path,1]\n",
    "    elif(path_front_direction[0]==-1):\n",
    "        path = paths[i_path][-1::-1]\n",
    "        the_other_end = path2node[i_path,0]\n",
    "            \n",
    "    n_path = 0\n",
    "    for j,n in enumerate(node2path[int(the_other_end),:]):\n",
    "            if(n==-1):\n",
    "                continue\n",
    "            n_path = n_path + 1\n",
    "            path_front_list.append(n)\n",
    "            path_front_direction.append(direction_indicator[int(the_other_end),j])\n",
    "    n_path_per_node.append(n_path)\n",
    "\n",
    "    path_coordinates = np.array([coordinates[path,0],coordinates[path,1],coordinates[path,2]],dtype=int)\n",
    "    diameter = distance_transform[path_coordinates[0,:],path_coordinates[1,:],path_coordinates[2,:]]\n",
    "    path_coordinates = (path_coordinates.astype(float))\n",
    "    diameter = (np.array(diameter).astype(float))\n",
    "    # mayavi plotting\n",
    "    dend_connections = np.zeros([len(diameter)-1,2])\n",
    "    dend_connections[:,0]=np.linspace(1,len(diameter)-1,len(diameter)-1) + len(plot_diameter)\n",
    "    dend_connections[:,1]=np.linspace(0,len(diameter)-2,len(diameter)-1) + len(plot_diameter)\n",
    "    connections = np.append(connections,dend_connections,axis=0)\n",
    "    plot_coordinates = np.append(plot_coordinates,path_coordinates,axis=1)\n",
    "    plot_diameter = np.append(plot_diameter,diameter)\n",
    "    init_potential = np.append(init_potential,np.zeros(len(diameter))-70)\n",
    "    section_index.append(i_path)\n",
    "    #h.pt3dclear(sec=dend[i_path])\n",
    "    dend[i_path].nseg = int(len(diameter)) + 1\n",
    "    dend[i_path].Ra = 100\n",
    "\n",
    "    # add morphology\n",
    "    path_coordinates[0,:] = path_coordinates[0,:]-path_coordinates[0,0]\n",
    "    path_coordinates[1,:] = path_coordinates[1,:]-path_coordinates[1,0]\n",
    "    path_coordinates[2,:] = path_coordinates[2,:]-path_coordinates[2,0]\n",
    "    h.pt3dadd(h.Vector(path_coordinates[0,:]), h.Vector(path_coordinates[1,:]), h.Vector(path_coordinates[2,:]), h.Vector(diameter), sec=dend[i_path])\n",
    "\n",
    "    last_dend.append(dend[i_path])\n",
    "    dend[i_path].connect(last_dend[0])\n",
    "\n",
    "    path_front_list.pop(0)\n",
    "    path_front_direction.pop(0)\n",
    "    n_path_per_node[0] = n_path_per_node[0]-1\n",
    "    if(n_path_per_node[0]==0):\n",
    "        n_path_per_node.pop(0)\n",
    "        last_dend.pop(0)\n",
    "    \n",
    "    # set biophysics\n",
    "    if(is_axon[int((last_dend[0].name())[5:-1])]):\n",
    "        is_axon[i_path] = True\n",
    "        set_biophysics(dend[i_path],'axon')\n",
    "    else:\n",
    "        set_biophysics(dend[i_path],'dend')\n",
    "unconstructed = np.where(reconstructed==False)[0]\n",
    "for i in unconstructed:\n",
    "    h.delete_section(sec=dend[i])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\ns = h.Shape()\\ns.show(False)\\n\\nv_vec = h.Vector()        # Membrane potential vector\\nt_vec = h.Vector()        # Time stamp vector\\nv_vec.record(soma(0.5)._ref_v)\\nt_vec.record(h._ref_t)\\nsimdur = 25.0\\nh.load_file('stdrun.hoc')\\nh.finitialize(-70)\\nh.continuerun(40)\\nh.topology()\\n\""
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''\n",
    "s = h.Shape()\n",
    "s.show(False)\n",
    "\n",
    "v_vec = h.Vector()        # Membrane potential vector\n",
    "t_vec = h.Vector()        # Time stamp vector\n",
    "v_vec.record(soma(0.5)._ref_v)\n",
    "t_vec.record(h._ref_t)\n",
    "simdur = 25.0\n",
    "h.load_file('stdrun.hoc')\n",
    "h.finitialize(-70)\n",
    "h.continuerun(40)\n",
    "h.topology()\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAANQElEQVR4nO3df6zdd13H8edLCpsbdgU6HftxLf7AhCmZ27UhC2RLV8ZSEYg/5jT+ICZrSDTZZoiONCMsgZh1U5iRxDRbDMuQKU6iAWa1iyaMWOCW/XJj2G0ynJ3M+WM4zIDB2z/ut3rZzu398b239/bd5yM56fecz/d8v59PTvLs6fec25uqQpLU03et9QQkSavHyEtSY0Zekhoz8pLUmJGXpMY2rPUE5tq8eXNt2bJlrachSceUAwcOPFVVp04aW1eR37JlCzMzM2s9DUk6piR5bL4xL9dIUmNGXpIaM/KS1JiRl6TGjLwkNWbkJakxIy9JjY36nnySPwF+ZLi7CfivqjonyYuBm4Bzh3PcUlW/M+ZckqSlGxX5qvr5w9tJfhd4erj7c8AJVfVjSU4CHkzykar60pjzSZKWZkV+4jVJgEuBbcNDBZycZAPw3cA3gK+uxLkkSYu3Utfk3wB8paoODvf/DPga8ATwZeCGqvqPFTqXJGmRFnwnn2QfcNqEoV1V9RfD9i8AH5kzthX4FnA68DLgU0n2VdWjE46/E9gJMDU1tbTZS5KOaMHIV9X2I40Pl2R+GjhvzsO/CPxVVX0TeDLJp4Fp4AWRr6o9wB6A6elpf+GsJK2glbhcsx14qKoen/PYl4FtmXUy8DrgoRU4lyRpCVYi8pfxnZdqAD4IvBT4B+BzwB9V1X0rcC5J0hKM/nZNVb19wmPPMPs1SknSGvInXiWpMSMvSY0ZeUlqzMhLUmNGXpIaM/KS1JiRl6TGjLwkNWbkJakxIy9JjRl5SWrMyEtSY0Zekhoz8pLUmJGXpMaMvCQ1ZuQlqTEjL0mNGXlJaszIS1JjRl6SGjPyktSYkZekxoy8JDVm5CWpMSMvSY0ZeUlqzMhLUmNGXpIaM/KS1JiRl6TGjLwkNWbkJakxIy9JjRl5SWrMyEtSY0Zekhoz8pLUmJGXpMaMvCQ1NirySc5Jsj/JPUlmkmydM/auJA8n+WKSN42fqiRpqTaMfP5u4NqquiPJjuH+hUleA1wGnA2cDuxL8uqq+tbI80mSlmDs5ZoCNg7bpwCHhu23ArdV1der6p+Ah4GtE54vSVpFY9/JXwnsTXIDs39hnD88fgawf85+jw+PvUCSncBOgKmpqZHTkSTNtWDkk+wDTpswtAu4CLiqqm5PcilwM7AdyIT9a9Lxq2oPsAdgenp64j6SpOVZMPJVtX2+sSS3AFcMdz8K3DRsPw6cNWfXM/n/SzmSpKNk7DX5Q8AFw/Y24OCw/ZfAZUlOSPIq4IeBz448lyRpicZek78cuDHJBuBZhmvrVfVAkj8FHgSeA37db9ZI0tE3KvJVdRdw3jxj7wPeN+b4kqRx/IlXSWrMyEtSY0Zekhoz8pLUmJGXpMaMvCQ1ZuQlqTEjL0mNGXlJaszIS1JjRl6SGjPyktSYkZekxoy8JDVm5CWpMSMvSY0ZeUlqzMhLUmNGXpIaM/KS1JiRl6TGjLwkNWbkJakxIy9JjRl5SWrMyEtSY0Zekhoz8pLUmJGXpMaMvCQ1ZuQlqTEjL0mNGXlJaszIS1JjRl6SGjPyktSYkZekxoy8JDVm5CWpsVGRT3JOkv1J7kkyk2Tr8PgbkxxIcv/w57aVma4kaSk2jHz+buDaqrojyY7h/oXAU8BPVdWhJD8K7AXOGHkuSdISjY18ARuH7VOAQwBVdfecfR4ATkxyQlV9feT5JElLMDbyVwJ7k9zA7KWf8yfs8zPA3fMFPslOYCfA1NTUyOlIkuZaMPJJ9gGnTRjaBVwEXFVVtye5FLgZ2D7nuWcD1wEXz3f8qtoD7AGYnp6uJc1eknREC0a+qrbPN5bkFuCK4e5HgZvmjJ0JfAz4lap6ZOQ8JUnLMPYrlIeAC4btbcBBgCSbgE8A76qqT488hyRpmcZek78cuDHJBuBZhmvrwG8APwRck+Sa4bGLq+rJkeeTJC3BqMhX1V3AeRMefy/w3jHHliSN50+8SlJjRl6SGjPyktSYkZekxoy8JDVm5CWpMSMvSY0ZeUlqzMhLUmNGXpIaM/KS1JiRl6TGjLwkNWbkJakxIy9JjRl5SWrMyEtSY0Zekhoz8pLUmJGXpMaMvCQ1ZuQlqTEjL0mNGXlJaszIS1JjRl6SGjPyktSYkZekxoy8JDVm5CWpMSMvSY0ZeUlqzMhLUmNGXpIaM/KS1JiRl6TGjLwkNWbkJakxIy9JjY2KfJJzkuxPck+SmSRbnzc+leSZJO8cN01J0nKMfSe/G7i2qs4B3j3cn+v9wB0jzyFJWqYNI59fwMZh+xTg0OGBJG8DHgW+NvIckqRlGhv5K4G9SW5g9l8F5wMkORn4beCNwBEv1STZCewEmJqaGjkdSdJcC0Y+yT7gtAlDu4CLgKuq6vYklwI3A9uBa4H3V9UzSY54/KraA+wBmJ6erqVNX5J0JKlafleTPA1sqqrKbM2frqqNST4FnDXstgn4NvDuqvqDIx1venq6ZmZmlj0fSToeJTlQVdOTxsZerjkEXAD8HbANOAhQVW+Yc/L3AM8sFHhJ0sobG/nLgRuTbACeZbi2LklaH0ZFvqruAs5bYJ/3jDmHJGn5/IlXSWrMyEtSY0Zekhoz8pLUmJGXpMaMvCQ1ZuQlqTEjL0mNGXlJaszIS1JjRl6SGjPyktSYkZekxoy8JDVm5CWpMSMvSY0ZeUlqzMhLUmNGXpIaM/KS1JiRl6TGjLwkNWbkJakxIy9JjRl5SWrMyEtSY0Zekhoz8pLUmJGXpMaMvCQ1ZuQlqTEjL0mNGXlJaszIS1JjRl6SGjPyktSYkZekxoy8JDVm5CWpsVGRT3JOkv1J7kkyk2TrnLHXJvn7JA8kuT/JieOnK0laig0jn78buLaq7kiyY7h/YZINwK3AL1fVvUleAXxz5LkkSUs0NvIFbBy2TwEODdsXA/dV1b0AVfXvI88jSVqGsZG/Etib5AZmL/2cPzz+aqCS7AVOBW6rqt0jzyVJWqIFI59kH3DahKFdwEXAVVV1e5JLgZuB7cNxXw/8BPA/wJ1JDlTVnROOvxPYCTA1NbXcdUiSJkhVLf/JydPApqqqJAGerqqNSS4DLqmqtw/7XQM8W1XXH+l409PTNTMzs+z5SNLxaHgTPT1pbOxXKA8BFwzb24CDw/Ze4LVJTho+hL0AeHDkuSRJSzT2mvzlwI1DyJ9luOxSVf+Z5PeAzzH74ewnq+oTI88lSVqiUZGvqruA8+YZu5XZr1FKktaIP/EqSY0ZeUlqzMhLUmNGXpIaM/KS1JiRl6TGjLwkNTbqvzVYaUn+DXhsreexDJuBp9Z6EkeZaz4+HG9rPlbX+/1VdeqkgXUV+WNVkpn5/t+Irlzz8eF4W3PH9Xq5RpIaM/KS1JiRXxl71noCa8A1Hx+OtzW3W6/X5CWpMd/JS1JjRl6SGjPyi5Tk5Un+JsnB4c+XzbPfJUm+mOThJFdPGH9nkkqyefVnPc7YNSe5PslDSe5L8rEkm47a5JdgEa9Zkvz+MH5fknMX+9z1arlrTnJWkr9N8oUkDyS54ujPfnnGvM7D+IuS3J3k40dv1iugqrwt4gbsBq4etq8Grpuwz4uAR4AfAF4C3Au8Zs74Wcz+asTHgM1rvabVXjNwMbBh2L5u0vPX+rbQazbsswO4AwjwOuAzi33ueryNXPMrgXOH7e8B/rH7mueM/ybwx8DH13o9S7n5Tn7x3gp8aNj+EPC2CftsBR6uqker6hvAbcPzDns/8FvM/krEY8GoNVfVX1fVc8N++4EzV3e6y7LQa8Zw/5aatR/YlOSVi3zuerTsNVfVE1X1eYCq+m/gC8AZR3PyyzTmdSbJmcBPAjcdzUmvBCO/eN9XVU8ADH9+74R9zgD+ec79x4fHSPIW4F+q6t7VnugKGrXm5/k1Zt8lrTeLmf98+yx27evNmDX/nyRbgB8HPrPyU1xxY9f8AWbfoH17lea3asb+Iu9WkuwDTpswtGuxh5jwWCU5aTjGxcud22pZrTU/7xy7gOeADy9tdkfFgvM/wj6Lee56NGbNs4PJS4HbgSur6qsrOLfVsuw1J3kz8GRVHUhy4UpPbLUZ+Tmqavt8Y0m+cvifq8M/4Z6csNvjzF53P+xM4BDwg8CrgHuTHH7880m2VtW/rtgClmEV13z4GL8KvBm4qIYLm+vMEee/wD4vWcRz16MxaybJi5kN/Ier6s9XcZ4racyafxZ4S5IdwInAxiS3VtUvreJ8V85afyhwrNyA6/nODyF3T9hnA/Aos0E//OHO2RP2+xLHxgevo9YMXAI8CJy61ms5whoXfM2YvRY79wO5zy7l9V5vt5FrDnAL8IG1XsfRWvPz9rmQY+yD1zWfwLFyA14B3AkcHP58+fD46cAn5+y3g9lvHDwC7JrnWMdK5EetGXiY2Wuc9wy3P1zrNc2zzhfMH3gH8I5hO8AHh/H7gemlvN7r8bbcNQOvZ/Yyx31zXtcda72e1X6d5xzjmIu8/62BJDXmt2skqTEjL0mNGXlJaszIS1JjRl6SGjPyktSYkZekxv4Xzy8a6ShKN6sAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "from mayavi import mlab\n",
    "mlab.options.offscreen = False\n",
    "mlab.figure(size=(1680,1050))\n",
    "\n",
    "stims=[]\n",
    "'''\n",
    "stims.append(h.epsp(soma(0.5)))\n",
    "stims[0].onset = 0.05\n",
    "stims[0].imax = 20\n",
    "'''\n",
    "\n",
    "for i in section_index[::300]:\n",
    "    print()\n",
    "    s = h.epsp(dend[i](0.5))\n",
    "    s.onset = 0.1\n",
    "    s.imax = 5\n",
    "    stims.append(s)\n",
    "\n",
    "g = h.graph()\n",
    "v_vec = h.Vector()        # Membrane potential vector\n",
    "t_vec = h.Vector()        # Time stamp vector\n",
    "v_vec.record(dend[section_index[1000]](0.5)._ref_v)\n",
    "t_vec.record(h._ref_t)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "h.load_file('stdrun.hoc')\n",
    "h.finitialize(-82)\n",
    "h.dt=0.01\n",
    "plot_coordinates[0,:] = plot_coordinates[0,:] - np.median(plot_coordinates[0,:])\n",
    "plot_coordinates[1,:] = plot_coordinates[1,:] - np.median(plot_coordinates[1,:])\n",
    "plot_coordinates[2,:] = plot_coordinates[2,:] - np.median(plot_coordinates[2,:])\n",
    "points = mlab.pipeline.scalar_scatter(plot_coordinates[0,:], plot_coordinates[1,:], plot_coordinates[2,:], plot_diameter)\n",
    "dataset = points.mlab_source.dataset\n",
    "dataset.point_data.get_array(0).name = 'diameter'\n",
    "dataset.lines = connections\n",
    "dataset.point_data.update()\n",
    "\n",
    "# The tube\n",
    "src = mlab.pipeline.set_active_attribute(points, point_scalars='diameter')\n",
    "stripper = mlab.pipeline.stripper(src)\n",
    "tube = mlab.pipeline.tube(stripper, tube_sides = 6, tube_radius = 1)\n",
    "tube.filter.capping = True\n",
    "#        tube.filter.use_default_normal = False\n",
    "tube.filter.vary_radius = 'vary_radius_by_absolute_scalar'\n",
    "tube.children = [] # Removing the old ones \n",
    "\n",
    "array_id = dataset.point_data.add_array(init_potential)\n",
    "dataset.point_data.get_array(array_id).name = 'voltage'\n",
    "dataset.point_data.update()\n",
    "src2 = mlab.pipeline.set_active_attribute(tube, \n",
    "                                            point_scalars='voltage')\n",
    "surf = mlab.pipeline.surface(src2, colormap='gist_heat')\n",
    "\n",
    "cb=mlab.colorbar(object=surf,nb_labels=10)\n",
    "lut=surf.module_manager.scalar_lut_manager.lut.table.to_array()\n",
    "lut[:,3]=np.linspace(0,255,256)\n",
    "lut[0:10,3]=np.zeros(10)+10\n",
    "\n",
    "surf.module_manager.scalar_lut_manager.lut.table = lut\n",
    "surf.module_manager.scalar_lut_manager.use_default_range=False\n",
    "membrane_potential = np.zeros(len(plot_diameter))\n",
    "#mlab.show()\n",
    "\n",
    "\n",
    "@mlab.animate(delay=100)\n",
    "def anim():\n",
    "    surf.actor.actor.rotate_y(40)\n",
    "    surf.actor.actor.rotate_x(-30)\n",
    "    surf.actor.actor.rotate_z(30)\n",
    "    surf.actor.actor.rotate_x(-90)\n",
    "    #surf.actor.actor.rotate_y(45)\n",
    "    #f = mlab.gcf()\n",
    "\n",
    "    for i in range(1500):\n",
    "        \n",
    "        surf.actor.actor.rotate_x(0.2)        \n",
    "        \n",
    "        h.fadvance()\n",
    "        \n",
    "        i_seg = 0\n",
    "        N = soma.n3d()\n",
    "        print(i,h.t)\n",
    "        \n",
    "        for seg in range(soma.n3d()):\n",
    "            membrane_potential[i_seg]=soma(soma.arc3d(seg)/soma.L).v\n",
    "            i_seg = i_seg + 1\n",
    "        \n",
    "        \n",
    "        for i_sec in section_index:\n",
    "            N=N+dend[i_sec].n3d()\n",
    "            for seg in range(dend[i_sec].n3d()):\n",
    "                membrane_potential[i_seg]=dend[i_sec](dend[i_sec].arc3d(seg)/dend[i_sec].L).v\n",
    "                i_seg = i_seg + 1\n",
    "        \n",
    "        dataset.point_data.remove_array(1)\n",
    "        dataset.point_data.update()\n",
    "        \n",
    "        # Adding the new scalar\n",
    "        dataset.point_data.add_array(membrane_potential)\n",
    "        dataset.point_data.get_array(1).name = 'voltage'\n",
    "        dataset.point_data.update()\n",
    "        #mlab.draw()\n",
    "        \n",
    "        # Updating the dataset\n",
    "        #ms = surf.mlab_source\n",
    "        #ms.m_data.update()\n",
    "\n",
    "        #surf.mlab_source.voltage = np.ones(10)*i*200\n",
    "        #f.scene.render()\n",
    "        mlab.process_ui_events()\n",
    "        \n",
    "        yield\n",
    "        mlab.savefig('video/screenshots/image'+str(i)+'.png')\n",
    "        \n",
    "anim()\n",
    "mlab.show()\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(t_vec, v_vec)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "import glob\n",
    "\n",
    "frameSize = (1680, 1007)\n",
    "\n",
    "out = cv2.VideoWriter('video/output_video.avi',cv2.VideoWriter_fourcc(*'DIVX'), 60, frameSize)\n",
    "for i in range(1500):\n",
    "    img = cv2.imread('video/screenshots/image'+str(i)+'.png')\n",
    "    out.write(img)\n",
    "\n",
    "out.release()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.8.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "21bfaabeb5f90dcb35f03e41cd0a17fbadd35308efdf98a46f30b833cf085000"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}