diff --git a/.gitignore b/.gitignore index d948def..17c48f4 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ access-token.txt local-build.sh +OSS_platform/**/__pycache__ + diff --git a/OSS-DBS - tutorial.pdf b/OSS-DBS - tutorial.pdf index f3c81c7..3ba5b4b 100644 Binary files a/OSS-DBS - tutorial.pdf and b/OSS-DBS - tutorial.pdf differ diff --git a/OSS_platform/Axon_files/axon.py b/OSS_platform/Axon_files/axon.py index 1c59c33..4f12931 100644 --- a/OSS_platform/Axon_files/axon.py +++ b/OSS_platform/Axon_files/axon.py @@ -138,21 +138,21 @@ def get_axonparams(self): 'total_nodes' int Total number of compartments 'ranvier_nodes' int - Number of node of Ranvier compartments + Number of node of Ranvier compartments 'para1_nodes' int - Number of first paranodal compartments + Number of first paranodal compartments (MYSA) 'para2_nodes' int - Number of second paranodal compartments + Number of second paranodal compartments (FLUT) 'inter_nodes' int - Number of internodal compartments + Number of internodal compartments (STIN) 'ranvier_length' float Length of the Node of Ranvier 'para1_length' float - Length of the first paranodal compartments + Length of the first paranodal compartments (MYSA) 'para2_length' float - Length of the second paranodal compartments + Length of the second paranodal compartments (FLUT) 'inter_length' float - Length of the internodal compartments + Length of the internodal compartments (STIN) 'deltax' float Length from a node of Ranvier to the next 'fiberD' float @@ -162,9 +162,9 @@ def get_axonparams(self): 'node_diameter': float Diameter of the nodes of Ranvier 'para1_diameter': float - Diameter of the first paranodal compartments + Diameter of the first paranodal compartments (MYSA) 'para2_diameter': float, - Diameter of the second paranodal compartments + Diameter of the second paranodal compartments (FLUT) 'condg': float Axial conductivity 'lamellas': int @@ -191,23 +191,40 @@ def __create_ref_nodes(self): ncomp = nstin+nmysa+nflut+1 ref_nodes = np.zeros(n,) - for i in range(0, n): - j = i % ncomp - if j == 0 or j == 1: # mysa <-> ranvier node - if i == 0: # first of all nodes - ref_nodes[i] = l_ranvier/2. - else: - ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_ranvier/2. - elif (j > 1 and j < (int(nmysa/2)+1)) or (j > (int(nmysa/2)+nflut+nstin) and j <(nmysa+nflut+nstin)): # mysa <-> mysa node - ref_nodes[i] = ref_nodes[i-1]+l_para1 - elif j == (int(nmysa/2)+1) or j == (int(nmysa/2)+nflut+nstin): # mysa <-> flut node - ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_para2/2. - elif (j > (int(nmysa/2)+1) and j < (int(nmysa/2)+int(nflut/2)+1)) or (j > (int(nmysa/2)+int(nflut/2)+nstin) and j < (int(nmysa/2)+nflut+nstin)): # flut <-> flut node - ref_nodes[i] = ref_nodes[i-1]+l_para2 - elif j == (int(nmysa/2)+int(nflut/2)+1) or j == (int(nmysa/2)+int(nflut/2)+nstin): # flut <-> stin node - ref_nodes[i] = ref_nodes[i-1]+l_para2/2.+l_inter/2. - else: # stin <-> stin node - ref_nodes[i] = ref_nodes[i-1]+l_inter + + # An easier implementation based on the structure of the internodals segments + structure = 'RMF'+6*'S'+'FM' # Node of Rnavier + internodal structure + + # Let's rename the lengths, just for convenience + L_R = l_ranvier + L_F = l_para2 + L_M = l_para1 + L_S = l_inter + + compartment = [0] # intial node of Ranvier + for i in range(n-1): + l1 = eval('L_'+structure[i%len(structure)]) # current segment length + l2 = eval('L_'+structure[(i+1)%len(structure)]) # next segment length + compartment.append(compartment[-1]+(l1+l2)/2) + + ref_nodes = np.array(compartment) + # for i in range(0, n): + # j = i % ncomp + # if j == 0 or j == 1: # mysa <-> ranvier node + # if i == 0: # first of all nodes + # ref_nodes[i] = l_ranvier/2. + # else: + # ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_ranvier/2. + # elif (j > 1 and j < (int(nmysa/2)+1)) or (j > (int(nmysa/2)+nflut+nstin) and j <(nmysa+nflut+nstin)): # mysa <-> mysa node + # ref_nodes[i] = ref_nodes[i-1]+l_para1 + # elif j == (int(nmysa/2)+1) or j == (int(nmysa/2)+nflut+nstin): # mysa <-> flut node + # ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_para2/2. + # elif (j > (int(nmysa/2)+1) and j < (int(nmysa/2)+int(nflut/2)+1)) or (j > (int(nmysa/2)+int(nflut/2)+nstin) and j < (int(nmysa/2)+nflut+nstin)): # flut <-> flut node + # ref_nodes[i] = ref_nodes[i-1]+l_para2 + # elif j == (int(nmysa/2)+int(nflut/2)+1) or j == (int(nmysa/2)+int(nflut/2)+nstin): # flut <-> stin node + # ref_nodes[i] = ref_nodes[i-1]+l_para2/2.+l_inter/2. + # else: # stin <-> stin node + # ref_nodes[i] = ref_nodes[i-1]+l_inter return ref_nodes * 1e-6 # um -> m @@ -217,6 +234,9 @@ def __create_nodes(self): self.__ref_nodes[-1]/2 else: ref_nodes = self.__ref_nodes + + # The reference nodes are not 100% symmetric, i.e. sometimes ref_nodes[i]!=ref_nodes[-i] + # due to the fitting option # standard position(along x axis in [x,y,z]) nodes = np.zeros((len(ref_nodes),3)) nodes[:,0] = ref_nodes @@ -256,17 +276,17 @@ def __create_axon_parameters(self, diameter): axon_parameters = copy.deepcopy(self.AXONPARAMETERS[diameter]) - nranvier = axon_parameters['ranvier_nodes'] - nstin = int(float(axon_parameters['inter_nodes'])/(nranvier-1)) - nmysa = int(float(axon_parameters['para1_nodes'])/(nranvier-1)) - nflut = int(float(axon_parameters['para2_nodes'])/(nranvier-1)) + nranvier = axon_parameters['ranvier_nodes'] #21 + nstin = int(float(axon_parameters['inter_nodes'])/(nranvier-1)) # 6 + nmysa = int(float(axon_parameters['para1_nodes'])/(nranvier-1)) # 2 + nflut = int(float(axon_parameters['para2_nodes'])/(nranvier-1)) # 2 axon_parameters['fiberD']=diameter axon_parameters['total_nodes'] = axon_parameters['ranvier_nodes'] + \ axon_parameters['para1_nodes']+axon_parameters['para2_nodes'] + \ - axon_parameters['inter_nodes'] + axon_parameters['inter_nodes'] #221 axon_parameters['inter_length'] = (axon_parameters['deltax'] - axon_parameters['ranvier_length'] - nmysa*axon_parameters['para1_length'] - - nflut*axon_parameters['para2_length'])/float(nstin) - return axon_parameters \ No newline at end of file + nflut*axon_parameters['para2_length'])/float(nstin) #70.5 + return axon_parameters diff --git a/OSS_platform/CSF_refinement_new.py b/OSS_platform/CSF_refinement_new.py index 51e9177..d8e0363 100644 --- a/OSS_platform/CSF_refinement_new.py +++ b/OSS_platform/CSF_refinement_new.py @@ -304,6 +304,7 @@ def Refine_CSF(MRI_param,DTI_param,Scaling,Domains,Field_calc_param,rel_div,CSF_ z_neuron_min=min(Vertices_neur[:,2]) #go over all voxels and check whether it contains CSF and intersect with the mesh + eps = 0.000000001 for x_coord in x_vect: for y_coord in y_vect: for z_coord in z_vect: @@ -314,9 +315,9 @@ def Refine_CSF(MRI_param,DTI_param,Scaling,Domains,Field_calc_param,rel_div,CSF_ if (x_pos<=x_neuron_max+CSF_ref_add and x_pos>=x_neuron_min-CSF_ref_add and y_pos<=y_neuron_max+CSF_ref_add and y_pos>=y_neuron_min-CSF_ref_add and z_pos<=z_neuron_max+CSF_ref_add and z_pos>=z_neuron_min-CSF_ref_add): - xv_mri=int((x_coord)/MRI_param.x_vox_size-0.000000001) #defines number of steps to get to the voxels containing x[0] coordinate - yv_mri=(int((y_coord)/MRI_param.y_vox_size-0.000000001))*MRI_param.M_x #defines number of steps to get to the voxels containing x[0] and x[1] coordinates - zv_mri=(int((z_coord)/MRI_param.z_vox_size-0.000000001))*MRI_param.M_x*MRI_param.M_y #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates + xv_mri= x_coord//(MRI_param.x_vox_size-eps) #defines number of steps to get to the voxels containing x[0] coordinate + yv_mri=(y_coord//(MRI_param.y_vox_size-eps))*MRI_param.M_x #defines number of steps to get to the voxels containing x[0] and x[1] coordinates + zv_mri=(z_coord//(MRI_param.z_vox_size-eps))*MRI_param.M_x*MRI_param.M_y #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates glob_index=int(xv_mri+yv_mri+zv_mri) pnt=Point(x_pos,y_pos,z_pos) diff --git a/OSS_platform/Dict_corrector.py b/OSS_platform/Dict_corrector.py index dfbb19b..b218b6f 100644 --- a/OSS_platform/Dict_corrector.py +++ b/OSS_platform/Dict_corrector.py @@ -4,9 +4,10 @@ @author: konstantin """ +import numpy as np #converts from user friendly units to the platform's format - + def rearrange_Inp_dict(dict_from_GUI): #fitst, change empty lines to 0 @@ -14,7 +15,6 @@ def rearrange_Inp_dict(dict_from_GUI): if dict_from_GUI[key]=='' or dict_from_GUI[key]=='0': dict_from_GUI[key]=0 - import numpy as np #second, some angles will be changed to rads if dict_from_GUI['Global_rot'] == 1: dict_from_GUI['alpha_array_glob'] =[i*np.pi/(180.0) for i in dict_from_GUI['alpha_array_glob']] diff --git a/OSS_platform/FEM_in_spectrum.py b/OSS_platform/FEM_in_spectrum.py index d528a60..0d8867d 100644 --- a/OSS_platform/FEM_in_spectrum.py +++ b/OSS_platform/FEM_in_spectrum.py @@ -198,8 +198,12 @@ class Conductivity : public dolfin::Expression if int(sine_freq)==int(signal_freq): c_unscaled = CompiledExpression(compile_cpp_code(conductivity_code).Conductivity(), - c00=c_unscaled00, c01=c_unscaled01, c02=c_unscaled02, c11=c_unscaled11, c12=c_unscaled12, c22=c_unscaled22, degree=0) - tensor= as_tensor([[c_unscaled[0], c_unscaled[1], c_unscaled[2]], [c_unscaled[1], c_unscaled[3], c_unscaled[4]],[c_unscaled[2],c_unscaled[4],c_unscaled[5]]]) + c00=c_unscaled00, c01=c_unscaled01, c02=c_unscaled02, + c11=c_unscaled11, c12=c_unscaled12, c22=c_unscaled22, + degree=0) + tensor = as_tensor([[c_unscaled[0], c_unscaled[1], c_unscaled[2]], + [c_unscaled[1], c_unscaled[3], c_unscaled[4]], + [c_unscaled[2],c_unscaled[4],c_unscaled[5]]]) f_vector_repr=project(tensor,TensorFunctionSpace(mesh, "Lagrange", 1),solver_type="cg", preconditioner_type="amg") file=File('Tensors/Ellipsoids_unscaled_at_'+str(signal_freq)+'_Hz.pvd') file< 3.0: # if larger than 3.0, 6 central compartments are required, otherwise 3 - if loc_inx == 0: - l_step = (nr["para1_length"] + nr["ranvier_length"]) / 2000 #switch to mm from µm - if loc_inx == 1 or loc_inx == 11: - l_step = (nr["ranvier_length"] + nr["para1_length"]) / 2000 #switch to mm from µm - if loc_inx == 2 or loc_inx == 10: - l_step = (nr["para1_length"] + nr["para2_length"]) / 2000 #switch to mm from µm - if loc_inx == 3 or loc_inx == 9: - l_step = (nr["para2_length"] + nr["inter_length"]) / 2000 #switch to mm from µm - if loc_inx == 4 or loc_inx == 5 or loc_inx == 6 or loc_inx == 7 or loc_inx == 8: - l_step = nr["inter_length"] / 1000 #switch to mm from µm - else: - if loc_inx == 0: - l_step = (nr["para1_length"]+nr["ranvier_length"]) / 2000 #switch to mm from µm - if loc_inx == 1 or loc_inx == 8: - l_step = (nr["ranvier_length"] + nr["para1_length"]) / 2000 #switch to mm from µm - if loc_inx == 2 or loc_inx == 7: - l_step = (nr["para1_length"] + nr["para2_nodes"]) / 2000 #switch to mm from µm - if loc_inx == 3 or loc_inx == 6: - l_step = (nr["para2_nodes"] + nr["inter_length"]) / 2000 #switch to mm from µm - if loc_inx == 4 or loc_inx == 5: - l_step = nr["inter_length"] / 1000 #switch to mm from µm + # this implementation doesnt + # for inx in range(1,self.pattern['num_segments']): + # if self.pattern['fiber_diameter'] > 3.0: # if larger than 3.0, 6 central compartments are required, otherwise 3 + # if loc_inx == 0: + # l_step = (nr["para1_length"] + nr["ranvier_length"]) / 2000 #switch to mm from µm + # if loc_inx == 1 or loc_inx == 11: + # l_step = (nr["ranvier_length"] + nr["para1_length"]) / 2000 #switch to mm from µm + # if loc_inx == 2 or loc_inx == 10: + # l_step = (nr["para1_length"] + nr["para2_length"]) / 2000 #switch to mm from µm + # if loc_inx == 3 or loc_inx == 9: + # l_step = (nr["para2_length"] + nr["inter_length"]) / 2000 #switch to mm from µm + # if loc_inx == 4 or loc_inx == 5 or loc_inx == 6 or loc_inx == 7 or loc_inx == 8: + # l_step = nr["inter_length"] / 1000 #switch to mm from µm + # else: + # if loc_inx == 0: + # l_step = (nr["para1_length"]+nr["ranvier_length"]) / 2000 #switch to mm from µm + # if loc_inx == 1 or loc_inx == 8: + # l_step = (nr["ranvier_length"] + nr["para1_length"]) / 2000 #switch to mm from µm + # if loc_inx == 2 or loc_inx == 7: + # l_step = (nr["para1_length"] + nr["para2_nodes"]) / 2000 #switch to mm from µm + # if loc_inx == 3 or loc_inx == 6: + # l_step = (nr["para2_nodes"] + nr["inter_length"]) / 2000 #switch to mm from µm + # if loc_inx == 4 or locranvier_nodes_inx == 5: + # l_step = nr["inter_length"] / 1000 #switch to mm from µm - loc_inx += 1 - self.pattern['Array_coord_pattern_O'][inx,:] = [0.0, self.pattern['Array_coord_pattern_O'][inx-1,1] + l_step , 0.0] # pattern along Y-axis + # loc_inx += 1 + # self.pattern['Array_coord_pattern_O'][inx,:] = [0.0, self.pattern['Array_coord_pattern_O'][inx-1,1] + l_step , 0.0] # pattern along Y-axis + + # if inx % n_comp == 0: + # loc_inx = 0 + + # An easier implementation based on the structure of the internodals segments: + # Ranvier - MYSA - FLUT - (6x or 3x) STIN - FLUT - MYSA - Ranvier - if inx % n_comp == 0: - loc_inx = 0 + structure = 'RMF'+3*'S'+'FM' # Node of Rnavier + internodal structure + if self.pattern['fiber_diameter'] > 3.0: # if larger than 3.0, 6 central compartments are required, otherwise 3 + structure = 'RMF'+6*'S'+'FM' # Node of Rnavier + internodal structure + + # Let's rename the lengths, just for convenience + L_R = nr['ranvier_length'] + L_M = nr['para1_length'] + L_F = nr['para2_length'] + L_S = nr['inter_length'] + + y_pattern = [0] # intial node of Ranvier + for i in range(self.pattern['num_segments']-1): + l1 = eval('L_'+structure[i%len(structure)]) # current segment length + l2 = eval('L_'+structure[(i+1)%len(structure)]) # next segment length + y_pattern.append(y_pattern[-1]+(l1+l2)/2) + y_pattern = np.array(y_pattern) + y_pattern -= y_pattern.mean() # center the pattern elif self.neuron_model == 'Reilly2016': - l_step = nr["deltax"] / 2000.0 # linear change of the internodal distance from 1 to 2 mm - for inx in range(1,self.pattern['num_segments']): - self.pattern['Array_coord_pattern_O'][inx,:] = [0.0, self.pattern['Array_coord_pattern_O'][inx-1,1] + l_step , 0.0] # pattern along Y-axis + l_step = nr["deltax"] / 2 # linear change of the internodal distance from 1 to 2 mm + y_pattern = np.arange(0,(1+self.pattern['num_segments'])*l_step, l_step) + y_pattern -= y_pattern.mean() # center the pattern + + # for inx in range(1,self.pattern['num_segments']): + # self.pattern['Array_coord_pattern_O'][inx,:] = [0.0, self.pattern['Array_coord_pattern_O'][inx-1,1] + l_step , 0.0] # pattern along Y-axis else: print("The neuron model is not implemented, exiting") raise SystemExit - + + self.pattern['Array_coord_pattern_O'][:,1] = np.array(y_pattern)*1e-3 # µm -> mm self.pattern['Array_coord_pattern_O'] = np.round(self.pattern['Array_coord_pattern_O'],8) np.savetxt('Neuron_model_arrays/'+str(self.pattern['name']), self.pattern['Array_coord_pattern_O'], delimiter=" ") @@ -327,23 +364,41 @@ def allocate_neurons_PO(self): self.ROI_radius=max(distance.cdist(self.custom_coord_PO, self.impl_coords_PO, 'euclidean'))[0] elif self.Type == 'VTA': # if placement in the ordered array (as for VTA) - - total_number_of_compartments = self.pattern['num_segments'] * self.VTA_seeds_PO.shape[0] - self.VTA_coord_PO = np.zeros((total_number_of_compartments,3),float) + if self.VTA_type=='cylinder': + generate_nodes_around_lead(MRI_shift= self.MRI_first_voxel, + Nx= self.VTA_structure['Nx'], + Nz_p= self.VTA_structure['Nz_p'], Nz_n= self.VTA_structure['Nz_n'], + dx= self.VTA_structure['dx'], dz=self.VTA_structure['dz'], + xmin= self.VTA_structure['xmin'], z0=self.VTA_structure['z0'], + n_dirs= self.VTA_structure['n_dirs']) + + # ROI is the extent of the array. Thus, it is defined as sqrt[x_max^2 + z^max^2] + xmax = self.VTA_structure['Nx']*self.VTA_structure['dx'] + self.VTA_structure['xmin'] + zmax = max(self.VTA_structure['Nz_p'],self.VTA_structure['Nz_n'])*self.VTA_structure['dz'] +\ + self.VTA_structure['z0'] + self.ROI_radius = np.sqrt(xmax**2+zmax**2) + + self.VTA_coord_PO = np.loadtxt('Neuron_model_arrays/All_neuron_models.csv', delimiter=' ') + + + else: + total_number_of_compartments = self.pattern['num_segments'] * self.VTA_seeds_PO.shape[0] + self.VTA_coord_PO = np.zeros((total_number_of_compartments,3),float) + + loc_ind=0 - loc_ind=0 - - for inx in range(self.VTA_seeds_PO.shape[0]): #goes through all seeding (central) nodes of neurons (axons) - angle_X = self.VTA_structure["X_angles_glob"][int(inx/self.N_models_in_plane)] - angle_Y = self.VTA_structure["Y_angles_glob"][int(inx/self.N_models_in_plane)] - angle_Z = self.VTA_structure["Z_angles_glob"][int(inx/self.N_models_in_plane)] - self.VTA_coord_PO[loc_ind:loc_ind + self.pattern['num_segments'],:] = self.place_neuron(angle_X,angle_Y,angle_Z,self.VTA_seeds_PO[inx,0],self.VTA_seeds_PO[inx,1],self.VTA_seeds_PO[inx,2]) - - loc_ind = loc_ind + self.pattern['num_segments'] + for inx in range(self.VTA_seeds_PO.shape[0]): #goes through all seeding (central) nodes of neurons (axons) + angle_X = self.VTA_structure["X_angles_glob"][int(inx/self.N_models_in_plane)] + angle_Y = self.VTA_structure["Y_angles_glob"][int(inx/self.N_models_in_plane)] + angle_Z = self.VTA_structure["Z_angles_glob"][int(inx/self.N_models_in_plane)] + self.VTA_coord_PO[loc_ind:loc_ind + self.pattern['num_segments'],:] = self.place_neuron(angle_X,angle_Y,angle_Z,self.VTA_seeds_PO[inx,0],self.VTA_seeds_PO[inx,1],self.VTA_seeds_PO[inx,2]) - np.savetxt('Neuron_model_arrays/All_neuron_models.csv', self.VTA_coord_PO , delimiter=" ") - self.ROI_radius = max(distance.cdist(self.VTA_coord_PO, self.impl_coords_PO, 'euclidean'))[0] - + loc_ind = loc_ind + self.pattern['num_segments'] + + self.ROI_radius = max(distance.cdist(self.VTA_coord_PO, self.impl_coords_PO, 'euclidean'))[0] + np.savetxt('Neuron_model_arrays/All_neuron_models.csv', self.VTA_coord_PO , delimiter=" ") + + elif self.Type == 'Imported': # we already have everything in self.imp_proc_coord_PO #self.imp_proc_coord_PO # here already merged from all projections !!! (extract each and @@ -554,7 +609,7 @@ def adjust_neuron_models(self,Domains,MRI_param): points_csf, points_encap, points_outside = (0,0,0) - # now we will folter out unphysiological neurins + # now we will filter out unphysiological neurons mesh = Mesh("Meshes/Mesh_unref.xml") subdomains_assigned = MeshFunction('size_t',mesh,"Meshes/Mesh_unref_physical_region.xml") @@ -574,6 +629,7 @@ def adjust_neuron_models(self,Domains,MRI_param): submesh_encup = SubMesh(mesh, subdomains_enc, 1) inx = 0 + #import pdb if self.Type != 'Imported': while inx < Array_coord.shape[0]: @@ -595,7 +651,7 @@ def adjust_neuron_models(self,Domains,MRI_param): check2_1=(voxel_array_CSF_shifted[:,0]>=Array_coord[inx,0]) # we could just check the sign check2_2=(voxel_array_CSF_shifted[:,1]>=Array_coord[inx,1]) check2_3=(voxel_array_CSF_shifted[:,2]>=Array_coord[inx,2]) - + check3=np.logical_and(np.logical_and(check1_1,check2_1),np.logical_and(np.logical_and(check1_2,check2_2),np.logical_and(check1_3,check2_3))) a=np.where((check3 == (True))) if str(a)!='(array([], dtype=int64),)': diff --git a/OSS_platform/Neuron_models_arangement_new.py b/OSS_platform/Neuron_models_arangement_new.py index de87ce1..2eb9efb 100644 --- a/OSS_platform/Neuron_models_arangement_new.py +++ b/OSS_platform/Neuron_models_arangement_new.py @@ -65,6 +65,13 @@ def rotate_globally(x_pos,y_pos,z_pos,alpha,beta,gamma): return Point_coords def create_structured_array(d): + """ + makes a pattern of coordinates, centered on the (0, 0, 0), + and rotates them based on the euler angles specified in the + input dictionary. + + Needs a rewrite! + """ x_coord=[] y_coord=[] @@ -84,7 +91,10 @@ def create_structured_array(d): for inx_angle in range(len(d["alpha_array_glob"])): for inx in range(len(x_coord)): - Rotated_point=rotate_globally(x_coord[inx],y_coord[inx],z_coord[inx],d["alpha_array_glob"][inx_angle],d["beta_array_glob"][inx_angle],d["gamma_array_glob"][inx_angle]) + Rotated_point=rotate_globally(x_coord[inx], y_coord[inx], z_coord[inx], + d["alpha_array_glob"][inx_angle], + d["beta_array_glob"][inx_angle], + d["gamma_array_glob"][inx_angle]) if not('Array_coord_str' in locals()): Array_coord_str=Rotated_point @@ -97,26 +107,29 @@ def create_structured_array(d): return (x_coord,y_coord,z_coord) -# takes a pattern model, rotates it around X,Y,Z and shifts its center to (x_loc,y_loc,z_loc) + def place_neuron(Arr_coord,alpha,beta,gamma,x_loc,y_loc,z_loc): - + """ + Arr_coord is the neural model pattern composed of several nodes. It's rotated by the euler + angles and then shifted by (x_loc, y_loc, z_loc). + """ Neuron_coord=np.zeros((Arr_coord.shape[0],3),float) alpha_matrix=np.array([[1,0,0], - [0,cos(alpha),-1*sin(alpha)], + [0,cos(alpha),-sin(alpha)], [0,sin(alpha),cos(alpha)]]) beta_matrix=np.array([[cos(beta),0,sin(beta)], [0,1,0], - [-1*sin(beta),0,cos(beta)]]) + [-sin(beta),0,cos(beta)]]) - gamma_matrix=np.array([[cos(gamma),-1*sin(gamma),0], + gamma_matrix=np.array([[cos(gamma),-sin(gamma),0], [sin(gamma),cos(gamma),0], [0,0,1]]) for inx in range(Neuron_coord.shape[0]): - xyz_alpha=np.array([Arr_coord[inx,0],Arr_coord[inx,1],Arr_coord[inx,2]]) + xyz_alpha=np.array([Arr_coord[inx,0], Arr_coord[inx,1], Arr_coord[inx,2]]) [Neuron_coord[inx,0],Neuron_coord[inx,1],Neuron_coord[inx,2]]=alpha_matrix.dot(xyz_alpha) xyz_beta=np.array([Neuron_coord[inx,0],Neuron_coord[inx,1],Neuron_coord[inx,2]]) @@ -214,9 +227,14 @@ def generate_pattern_model(name,N_Ranv,axon_param,Axon_Model_Type): # builds or resaves the full neuron array in the positive octant, computes its extent (distance between the implantation point and the furtherst compartment). #def get_neuron_models_dims(Full_model_ready,X_imp,Y_imp,Z_imp,Neuron_param,axon_param,plane_rot): -def get_neuron_models_dims(Full_model_ready,X_imp,Y_imp,Z_imp,Neuron_param,plane_rot): +def get_neuron_models_dims(Full_model_ready,X_imp,Y_imp,Z_imp,Neuron_param,plane_rot): + """ + If the full_model is not ready loads the model pattern, rotates and translates + it to the center points it should be on. + """ + start_neuron_models=time.clock() - + if Full_model_ready==0: #[ranvier_nodes, para1_nodes, para2_nodes, inter_nodes]=(axon_param[:4]) if Neuron_param.name != 'pattern.csv': #get a pattern from the externally provided file @@ -240,7 +258,11 @@ def get_neuron_models_dims(Full_model_ready,X_imp,Y_imp,Z_imp,Neuron_param,plane for alpha_angle in Neuron_param.alpha_array: for beta_angle in Neuron_param.beta_array: for gamma_angle in Neuron_param.gamma_array: - New_model_coord=place_neuron(Array_coord_temp,alpha_angle,beta_angle,gamma_angle,Neuron_param.x_loc[inx],Neuron_param.y_loc[inx],Neuron_param.z_loc[inx]) + New_model_coord=place_neuron(Array_coord_temp, # the pattern + alpha_angle, beta_angle, gamma_angle, + Neuron_param.x_loc[inx], + Neuron_param.y_loc[inx], + Neuron_param.z_loc[inx]) if not('Array_coord' in locals()): Array_coord=New_model_coord else: @@ -248,7 +270,11 @@ def get_neuron_models_dims(Full_model_ready,X_imp,Y_imp,Z_imp,Neuron_param,plane if plane_rot==1: # if placement in the ordered array (as for VTA) for inx in range(len(Neuron_param.x_loc)): #goes through all seeding (central) nodes of neurons (axons) - New_model_coord=place_neuron(Array_coord_temp,Neuron_param.alpha_array_global[int(inx/Neuron_param.N_models_in_plane)],Neuron_param.beta_array_global[int(inx/Neuron_param.N_models_in_plane)],Neuron_param.gamma_array_global[int(inx/Neuron_param.N_models_in_plane)],Neuron_param.x_loc[inx],Neuron_param.y_loc[inx],Neuron_param.z_loc[inx]) + New_model_coord=place_neuron(Array_coord_temp, + Neuron_param.alpha_array_global[int(inx/Neuron_param.N_models_in_plane)], + Neuron_param.beta_array_global[int(inx/Neuron_param.N_models_in_plane)], + Neuron_param.gamma_array_global[int(inx/Neuron_param.N_models_in_plane)], + Neuron_param.x_loc[inx], Neuron_param.y_loc[inx], Neuron_param.z_loc[inx]) if not('Array_coord' in locals()): Array_coord=New_model_coord else: @@ -632,15 +658,19 @@ def build_neuron_models(d,MRI_param): if d["Global_rot"]==1: # for VTA - #center node of the center neuron is located at the seed (the seeding coordinate is in the positive octant space) + # center node of the center neuron is located at the seed (the seeding coordinate is in the positive octant space) x_seeding=MRI_param.x_shift+d["x_seed"] y_seeding=MRI_param.y_shift+d["y_seed"] z_seeding=MRI_param.z_shift+d["z_seed"] #number of models, seeded in one plane. This parameter is 0 by default in Neuron_info class N_models_in_plane=(d["x_steps"]+1)*(d["y_steps"]+1)*(d["z_steps"]+1) - - [X_coord_old,Y_coord_old,Z_coord_old]=create_structured_array(d) + + # coordinates of the axon centers. Centered around (0,0,0) + [X_coord_old,Y_coord_old,Z_coord_old]=create_structured_array(d) + + # translates the coordinates of the centers to the positive octant + # but why no broadcasting? why changing to list? X_coord=[xs+x_seeding for xs in X_coord_old] Y_coord=[ys+y_seeding for ys in Y_coord_old] Z_coord=[zs+z_seeding for zs in Z_coord_old] @@ -685,8 +715,13 @@ def build_neuron_models(d,MRI_param): Array_coord_load=Array_coord_load_get.values n_segments=Array_coord_load.shape[0] #all segments should be in the pattern model del Array_coord_load - - Neuron_param=Neuron_info(pattern_model_name,X_coord,Y_coord,Z_coord,d["YZ_angles"],d["ZX_angles"],d["XY_angles"],d["alpha_array_glob"],d["beta_array_glob"],d["gamma_array_glob"],N_models_in_plane) + + # a class that is used as a dictionary only!!! + Neuron_param=Neuron_info(pattern_model_name, X_coord, Y_coord, Z_coord, + d["YZ_angles"], d["ZX_angles"], d["XY_angles"], + d["alpha_array_glob"], d["beta_array_glob"], d["gamma_array_glob"], + N_models_in_plane) + with open('Neuron_model_arrays/Neuron_param_class.file', "wb") as f: pickle.dump(Neuron_param, f, pickle.HIGHEST_PROTOCOL) @@ -697,13 +732,18 @@ def build_neuron_models(d,MRI_param): "definition of ROI is required for adaptive mesh refinement. ROI is a sphere encompassing all neuron models" #ROI_radius=get_neuron_models_dims(d["Neuron_model_array_prepared"],Xt_new,Yt_new,Zt_new,Neuron_param,param_axon,d["Global_rot"]) - ROI_radius=get_neuron_models_dims(d["Neuron_model_array_prepared"],Xt_new,Yt_new,Zt_new,Neuron_param,d["Global_rot"]) + ROI_radius=get_neuron_models_dims(d["Neuron_model_array_prepared"], + Xt_new, Yt_new,Zt_new, + Neuron_param, d["Global_rot"]) if d["Axon_Model_Type"] == 'McIntyre2002': - Neuron_model_misc=np.array([nr["ranvier_nodes"], nr["para1_nodes"], nr["para2_nodes"], nr["inter_nodes"], nr["ranvier_length"], nr["para1_length"], nr["para2_length"], nr["inter_length"], nr["deltax"],d["diam_fib"],d["n_Ranvier"],ROI_radius,n_segments]) + Neuron_model_misc=np.array([nr["ranvier_nodes"], nr["para1_nodes"], nr["para2_nodes"], nr["inter_nodes"], + nr["ranvier_length"], nr["para1_length"], nr["para2_length"], nr["inter_length"], + nr["deltax"],d["diam_fib"],d["n_Ranvier"],ROI_radius,n_segments]) elif d["Axon_Model_Type"] == 'Reilly2016': - Neuron_model_misc=np.array([d["n_Ranvier"], 0, 0, d["n_Ranvier"]-1, 0.0, 0.0, 0.0, 0.0, internodel_length,d["diam_fib"],d["n_Ranvier"],ROI_radius,n_segments]) + Neuron_model_misc=np.array([d["n_Ranvier"], 0, 0, d["n_Ranvier"]-1, 0.0, 0.0, 0.0, 0.0, + internodel_length,d["diam_fib"],d["n_Ranvier"],ROI_radius,n_segments]) elif d["Axon_Model_Type"] == 'Gillies2006': Neuron_model_misc=np.array([0, 0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0,0,ROI_radius,n_segments]) diff --git a/OSS_platform/Tissue_marking_new.py b/OSS_platform/Tissue_marking_new.py index ec9a24b..9e72d15 100644 --- a/OSS_platform/Tissue_marking_new.py +++ b/OSS_platform/Tissue_marking_new.py @@ -50,6 +50,7 @@ def get_cellmap(mesh,subdomains_assigned,Domains,MRI_param,default_material): subdomains.set_all(0) #subdomains = MeshFunctionSizet(mesh, 3, 0) + eps = 0.000000001 for cell in cells(mesh): x_coord=cell.midpoint().x() @@ -58,16 +59,16 @@ def get_cellmap(mesh,subdomains_assigned,Domains,MRI_param,default_material): z_coord=cell.midpoint().z() - xv_mri=int((x_coord)/voxel_size_x-0.000000001) #defines number of steps to get to the voxels containing x[0] coordinate - yv_mri=(int((y_coord)/voxel_size_y-0.000000001))*Mx_mri #defines number of steps to get to the voxels containing x[0] and x[1] coordinates - zv_mri=(int((z_coord)/voxel_size_z-0.000000001))*Mx_mri*My_mri #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates + xv_mri=int((x_coord)/voxel_size_x-eps) #defines number of steps to get to the voxels containing x[0] coordinate + yv_mri=(int((y_coord)/voxel_size_y-eps))*Mx_mri #defines number of steps to get to the voxels containing x[0] and x[1] coordinates + zv_mri=(int((z_coord)/voxel_size_z-eps))*Mx_mri*My_mri #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates k_mri=xv_mri+yv_mri+zv_mri #print k_mri k_mri=int(k_mri) - x_vect_index=(int((x_coord)/voxel_size_x-0.000000001)) - y_vect_index=(int((y_coord)/voxel_size_y-0.000000001)) - z_vect_index=(int((z_coord)/voxel_size_z-0.000000001)) + x_vect_index=(int((x_coord)/voxel_size_x-eps)) + y_vect_index=(int((y_coord)/voxel_size_y-eps)) + z_vect_index=(int((z_coord)/voxel_size_z-eps)) if x_vect_index>=Mx_mri or y_vect_index>=My_mri or z_vect_index>=Mz_mri or x_coord<0.0 or y_coord<0.0 or z_coord<0.0: x_vect_index,y_vect_index,z_vect_index=0,0,0 #cell is outside MRI data, assign first voxel of MRI, which will not fullfil the following condition, and the cell will be remain marked with 0 (default) @@ -79,18 +80,10 @@ def get_cellmap(mesh,subdomains_assigned,Domains,MRI_param,default_material): '''First, find the corresponding voxel in MRI, then check the value of the voxel and assign a corresponding conductivity''' if (round(y_vect[y_vect_index]-y_coord,8)<=voxel_size_y and round(y_vect[y_vect_index]-y_coord,8)>=0 and (round(x_vect[x_vect_index]-x_coord,8)<=voxel_size_x and round(x_vect[x_vect_index]-x_coord,8)>=0) and round(z_vect[z_vect_index]-z_coord)<=voxel_size_z and round(z_vect[z_vect_index]-z_coord,8)>=0): - if int(Tissue_array[k_mri])==3: - #cell.mark(subdomains, 1) - subdomains[cell]=3 - - if int(Tissue_array[k_mri])==2: - subdomains[cell]=2 - - if int(Tissue_array[k_mri])==1: - subdomains[cell]=1 - if int(Tissue_array[k_mri])==0: subdomains[cell]=default_material + else: + subdomains[cell]=int(Tissue_array[k_mri]) #old way # if Domains.Float_contacts!=-1: # subdomains.array()[subdomains_assigned.array()==Domains.Float_contacts]=5 #5 is index for float contact (very high cond and perm) @@ -199,18 +192,17 @@ def get_cellmap_tensors(mesh,subdomains_assigned,Domains,MRI_param,DTI_param,def c12 = MeshFunction("double", mesh, 3, 0.0) c22 = MeshFunction("double", mesh, 3, 1.0) - + eps = 0.000000001 for cell in cells(mesh): x_coord=cell.midpoint().x() - y_coord=cell.midpoint().y() z_coord=cell.midpoint().z() - xv_mri=int(x_coord/voxel_size_x-0.000000001) #defines number of steps to get to the voxels containing x[0] coordinate - yv_mri=xv_mri+(int(y_coord/voxel_size_y-0.000000001))*Mx_mri #defines number of steps to get to the voxels containing x[0] and x[1] coordinates - zv_mri=yv_mri+(int(z_coord/voxel_size_z-0.000000001))*Mx_mri*My_mri #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates + xv_mri=x_coord//(voxel_size_x-eps) #defines number of steps to get to the voxels containing x[0] coordinate + yv_mri=xv_mri+ (y_coord//(voxel_size_y-eps))*Mx_mri #defines number of steps to get to the voxels containing x[0] and x[1] coordinates + zv_mri=yv_mri+ (z_coord//(voxel_size_z-eps))*Mx_mri*My_mri #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates k_mri=zv_mri k_mri=int(k_mri) @@ -218,18 +210,18 @@ def get_cellmap_tensors(mesh,subdomains_assigned,Domains,MRI_param,DTI_param,def if k_mri>=Mx_mri*My_mri*Mz_mri or k_mri<0: k_mri=0 #cell is outside MRI data, assign first voxel of MRI, which will not fullfil the following condition, and the cell will be remain marked with 0 (default) - x_vect_index=(int((x_coord)/voxel_size_x-0.000000001)) - y_vect_index=(int((y_coord)/voxel_size_y-0.000000001)) - z_vect_index=(int((z_coord)/voxel_size_z-0.000000001)) + x_vect_index=(int((x_coord)/voxel_size_x-eps)) + y_vect_index=(int((y_coord)/voxel_size_y-eps)) + z_vect_index=(int((z_coord)/voxel_size_z-eps)) if x_vect_index>=Mx_mri or y_vect_index>=My_mri or z_vect_index>=Mz_mri or x_coord<0.0 or y_coord<0.0 or z_coord<0.0: x_vect_index,y_vect_index,z_vect_index=0,0,0 #cell is outside MRI data, assign first voxel of MRI, which will not fullfil the following condition, and the cell will be remain marked with 0 (default) - xv_dti=int((x_coord-DTI_param.x_start)/voxel_size_x_DTI-0.00000001) - yv_dti=xv_dti+(int((y_coord-DTI_param.y_start)/voxel_size_y_DTI-0.00000001))*Mx_dti - zv_dti=yv_dti+(int((z_coord-DTI_param.z_start)/voxel_size_z_DTI-0.00000001))*Mx_dti*My_dti + xv_dti=(x_coord-DTI_param.x_start)//(voxel_size_x_DTI-eps) + yv_dti=xv_dti+((y_coord-DTI_param.y_start)//(voxel_size_y_DTI-eps))*Mx_dti + zv_dti=yv_dti+((z_coord-DTI_param.z_start)//(voxel_size_z_DTI-eps))*Mx_dti*My_dti k_dti=zv_dti k_dti=int(k_dti) @@ -238,9 +230,9 @@ def get_cellmap_tensors(mesh,subdomains_assigned,Domains,MRI_param,DTI_param,def k_dti=0 #cell is outside DTI data, assign first voxel of DTI, which will not fullfil the following condition, and the cell will be considered as isotropic ''' otherwise comment it out''' - x_vect_DTI_index=(int((x_coord-DTI_param.x_start)/voxel_size_x_DTI-0.000000001)) - y_vect_DTI_index=(int((y_coord-DTI_param.y_start)/voxel_size_y_DTI-0.000000001)) - z_vect_DTI_index=(int((z_coord-DTI_param.z_start)/voxel_size_z_DTI-0.000000001)) + x_vect_DTI_index=(int((x_coord-DTI_param.x_start)/voxel_size_x_DTI-eps)) + y_vect_DTI_index=(int((y_coord-DTI_param.y_start)/voxel_size_y_DTI-eps)) + z_vect_DTI_index=(int((z_coord-DTI_param.z_start)/voxel_size_z_DTI-eps)) if x_vect_DTI_index>=Mx_dti or y_vect_DTI_index>=My_dti or z_vect_DTI_index>=Mz_dti or x_coord<0.0 or y_coord<0.0 or z_coord<0.0 or z_vect_DTI_index<0 or y_vect_DTI_index<0 or x_vect_DTI_index<0: x_vect_DTI_index,y_vect_DTI_index,z_vect_DTI_index=0,0,0 #cell is outside MRI data, assign first voxel of MRI, which will not fullfil the following condition, and the cell will be remain marked with 0 (default) diff --git a/OSS_platform/VTA/cylindrical_array.py b/OSS_platform/VTA/cylindrical_array.py new file mode 100644 index 0000000..507f092 --- /dev/null +++ b/OSS_platform/VTA/cylindrical_array.py @@ -0,0 +1,203 @@ +import numpy as np +import pandas as pd +from scipy.spatial.transform import Rotation as R + +def generate_nodes_around_lead(MRI_shift, + Nx, Nz_p, Nz_n, + dx, dz, xmin, z0, + n_dirs, + produce_centers=True, + alpha=None, gamma=None, order='x_first', + test=True, + input_dict = 'GUI_inp_dict.py', + target_path = 'Neuron_model_arrays/', + pattern_path= 'Neuron_model_arrays/default_pattern.csv'): + + """ + Generated a grid of nodes surrounding the lead cylindrically, in a system + of coordinates (S) attached to the lead (center at the implantation + location, z+ at the direction of the lead tip toward the skull). A 2D grid + will be generated in x-z surface. Extra surfaces will be generated by + revolving the seeding surface along the z-direction with uniform angular + seperation. + + elementary arguments: + ===================== + input_dict: path to the input dictionary + target_path: save directort + pattern_path: path to the pattern that encapsulates the axon model nodes + produce_centers: if true, an extra csv file containing the location of the + centeral points in axons will be generated. These central + points may or may not correspond to nodes of Ranvier + (depending on the model and the length of the axon.) + + Geometrical inputs: + =================== + Nx: number of seeds in x+ and x- direction + Nz_p: number of seeds in z+ direction (starts at the implantation level) + Nz_n: number of seeds in z- direction (starts at the implantation level) + dx: distance between seeds in x direction + dy: distance between seeds in y direction + xmin: minimum seeding distance lateral to the lead (in S coordinates) + z0: seeding starting level from the tip of the lead upward in S coordinates (mm) + n_dirs:number of directions for VTA calculations. `n_dir` surfaces with a + uniform angular seperation and the same seeding pattern will be placed + around the lead. E.g., 1 corresponds to one surface, faced toward the + y-axis in (S coordinate system); 3 corresponds to three surfaces, one + faced toward the y-axis and the other two each with 60 degree seperation + from each other. In general the angle between the planes is `180/n_dirs`. + + advanced inputs: + ================ + In S coordinates, seeding points are in x-z plane and the axon nodes are + place perpendicular to this plane, i.e. y direction. It's however possible + to rotate the axons such that they are not necessarily in y-direction. To + do so, one can specify rotation angles `alpha` and `gamma`, around the x + and z axis, respectively. Note that since the axons are along the y-axis + at first, a rotation along y doesn't make any difference. This will not + impose any limitation though. All possible directions can be specified by + two directions. + + In a spherical coordiante system that shares the same x,y, z axes with S + coordinates system we have: + + alpha = theta - 90 + gamma = phi - 90 + + Look at here for the definition of the (theta, phi) in spherical + coordinates (physics convention): + + https://en.wikipedia.org/wiki/File:3D_Spherical.svg + + alpha: rotation along the x direction (in degree) + gamma: rotation along the y direction (in degree) + order: `x_first` will first rotate along the x and the z (default). The + alternative is `z_first` + + + output: + ======= + nodes: an array which contains the (x,y,z) of all the nodes on the grid. + Since each node can be fully specified by its direction_id, seeding + index in x and z, and its location along the axon, these (x,y,z) + coordinates are augmented by a 4-tuple to address the index of the + node. Each entry thus has a form of (dir_id, i, j, k, x,y,z) with + `i, j`, and `k` standing for the index in x,z, and y (on the axon) + and`dir_id` the direction id. `i,j` start from the least value of + corresponding x,z seeding point. + + axons: an array, almost with the same structure as the nodes, except that + the coordinates are averaged alongside the axon (`k` direction) for + each index `(dir_id, i,j)`. + """ + + import importlib.util + spec = importlib.util.spec_from_file_location("d",input_dict) + foo = importlib.util.module_from_spec(spec) + spec.loader.exec_module(foo) + d = foo.d + + ## Constructing the S coordinate system + origin = np.array([d['Implantation_coordinate_X'], + d['Implantation_coordinate_Y'], + d['Implantation_coordinate_Z']]) + + # z-direction + zhat = np.array([d['Second_coordinate_X']-d['Implantation_coordinate_X'], + d['Second_coordinate_Y']-d['Implantation_coordinate_Y'], + d['Second_coordinate_Z']-d['Implantation_coordinate_Z']]) + zhat /= np.linalg.norm(zhat) + + # x-direction + xhat = np.array([1,0,0]) # an arbitrary direction + xhat = xhat -zhat*np.dot(xhat,zhat) # Gram-schmit + xhat /= np.linalg.norm(xhat) + + # y-direction + yhat = np.cross(zhat, xhat) + + + ## Constructing the grid structure + xc = np.arange(xmin, xmin+Nx*dx, dx) + xc = np.sort(np.append(-xc, xc)) + zc = np.arange(z0-Nz_n*dz, z0+(1+Nz_p)*dz, dz) + + # y direction is based on the neuron model + pattern = np.loadtxt(pattern_path, delimiter=' ') + pattern = (pattern-pattern[::-1])/2 # enforce symmetry (a bug in axon model) + + + if ((alpha==None) & (gamma==None)): + y_nodes = pattern[:,1] + x,y,z = np.meshgrid(xc, y_nodes, zc, indexing='xy') + ref = yhat.copy() # reference direction + + else: + # convert None to zeros + if alpha==None: + alpha=0. + if gamma==None: + gamma=0. + + rx = R.from_euler('x', alpha, degrees=True).as_dcm()#.as_matrix() + rz = R.from_euler('z', gamma, degrees=True).as_dcm()#.as_matrix() + if order=='x_first': + Rtot = rz @ rx + else: + Rtot = rx @ rz + + # Only the pattern are rotated here. The nodes, originally on y-axis, + # will have a general 3D coordinates after rotation. To braodcast it + # to all seed centers, we add the x, and z components after making a + # meshgrid with the xc, rotated_y, and zc. + pattern_= (Rtot @ pattern.T).T + x_nodes = pattern_[:,0] + y_nodes = pattern_[:,1] + z_nodes = pattern_[:,2] + + x,y,z = np.meshgrid(xc, y_nodes, zc, indexing='xy') + x = (x.transpose(1,2,0) + x_nodes).transpose(2,0,1) + z = (z.transpose(1,2,0) + z_nodes).transpose(2,0,1) + ref = (Rtot @ yhat.T).T + + + ## Constructing the nodes array + dir_id = [0]*len(x.flatten()) # first direction_id = 0 + k,i,j=np.indices(x.shape) + nodes = np.stack([x.flatten(), y.flatten(), z.flatten(), dir_id,\ + i.flatten(), j.flatten(), k.flatten()], axis=-1) + + ## Rovolving + for n in range(1,n_dirs): + theta = n/(n_dirs)*np.pi + Rrev = R.from_euler('z', theta).as_dcm()#.as_matrix() + revolved = np.dot(np.array([x, y, z]).transpose(1,2,3,0), Rrev.T) + x_,y_,z_ = revolved.transpose(3,0,1,2) + + dir_id = [n]*len(x_.flatten()) + k_,i_,j_=np.indices(x_.shape) + nodes_ = np.stack([x_.flatten(), y_.flatten(), z_.flatten(), dir_id, + i_.flatten(), j_.flatten(), k_.flatten()],axis=-1) + nodes = np.vstack((nodes,nodes_)) + + ## transferring to the lab coordinates + A = np.array([xhat,yhat,zhat]).T + nodes[:,:3] = nodes[:,:3] @ A.T # rotate + nodes[:,:3] += origin # translate + nodes[:,:3] -= MRI_shift # shift the positive octant + + df_nodes = pd.DataFrame(data=nodes, + columns=['x','y','z','dir_id','i','j','k']) + df_nodes = df_nodes.set_index(['dir_id','i','j','k']).sort_index() + df_nodes.to_csv(target_path+'nodes.csv',float_format='%.5f') + + if produce_centers: + df_axons = df_nodes.groupby(['dir_id','i','j']).mean() + df_axons.to_csv(target_path+'axons.csv',float_format='%.5f') + + # Saving the arrays as csv with correct order + np.savetxt(target_path+'All_neuron_models.csv', + df_nodes[['x','y','z']].values, delimiter=' ') + + if test: + return nodes, df_nodes, df_axons, xhat, yhat, zhat diff --git a/OSS_platform/VTA/visualize.py b/OSS_platform/VTA/visualize.py new file mode 100644 index 0000000..dc81e97 --- /dev/null +++ b/OSS_platform/VTA/visualize.py @@ -0,0 +1,254 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + +import os +path = os.getcwd() +os.chdir(path) + +import importlib.util +spec = importlib.util.spec_from_file_location("d",path+'/GUI_inp_dict.py') +foo = importlib.util.module_from_spec(spec) +spec.loader.exec_module(foo) +d = foo.d + +def plot_vta_figs(path_target = path+'/Images/', + path_nodes = path+'/Neuron_model_arrays/All_neuron_models.csv', + path_pattern = path+'/Neuron_model_arrays/default_pattern.csv', + path_activation = path+'/Field_solutions/Activation/Network_status.npy', + activate_close_axons=False): + + kwargs = {'path_target': path_target, 'path_nodes':path_nodes, 'path_pattern': path_pattern, + 'path_activation': path_activation, 'activate_close_axons': activate_close_axons} + + if d['VTA_type']=='cylinder': + plot_cylindircal_vta(**kwargs) + else: + plot_grid_vta(**kwargs) + + + +def plot_cylindircal_vta(path_activation, path_target, activate_close_axons, **kwargs): + """ + Visualizes the VTA cross-section for each surface. + """ + + # reading the input dict + #import importlib.util + #spec = importlib.util.spec_from_file_location("d",input_dict) + #foo = importlib.util.module_from_spec(spec) + #spec.loader.exec_module(foo) + #d = foo.d + + Nx = d['Nx'] + Nz_p = d['Nz_p'] + Nz_n = d['Nz_n'] + dx = d['dx'] + dz = d['dz'] + xmin = d['xmin'] + z0 = d['z0'] + n_dirs = d['n_dirs'] + + # constructing the grid + xc = np.arange(xmin, xmin+Nx*dx, dx) + xc = np.sort(np.append(-xc, xc)) + zc = np.arange(z0-Nz_n*dz, z0+(1+Nz_p)*dz, dz) + x,z = np.meshgrid(xc, zc) + + # reading the activation status + activations = np.load(path_activation) + if activate_close_axons: + activations= np.abs(activations) + + axon_per_plane = len(xc)*len(zc) + for n in range(n_dirs): + vta_plane = activations[n*axon_per_plane: (n+1)*axon_per_plane] + vta_plane = vta_plane.reshape(len(xc), len(zc)) + + plt.figure() + vmin = vmin = (activate_close_axons-1) + ax = plt.pcolormesh(xc, zc, vta_plane.T, shading='nearest', cmap='coolwarm', + edgecolors='gray', linewidths=0.25, vmin=vmin, vmax=1) + plt.colorbar(ax,fraction=0.03,label='activation state') + ax = plt.gca() + ax.set_aspect('equal') + plt.title('VTA profile prependicular to the direction '+ str(n)) + plt.tight_layout() + plt.savefig(path_target+'Cyl_VTA_dir'+str(n)+'.png', dpi=300, bbox_inches='tight') + plt.close() + + +def plot_grid_vta(path_nodes, path_pattern, path_activation, path_target, + activate_close_axons, **kwargs): + """ + extracts axons directions and the centers of the axon arrays from the + `All_neuron_models.csv`. Makes png (and possibly) nifti files to show VTA + orthogonal to each direction. + + If axon centers are seeded in a plane (which will be recognized from the + input dictionary) then just png files will be generated. Otherwise, for + each direction, a nifti file is also produced and saved in the same directory. + + arguments: + ---------- + path_target: loaction to save VTA cross-section images and niftis + path_input_dict: path of the input dictionary (.py) + path_nodes: path of the full neuron array model (.csv) file + path_pattern: path of the neural model pattern (.csv) file + path_activation: path of activation status (.npy) files + activate_close_axons: True if removed axons are assumed activated (bool) + + """ + import nibabel as nib + from scipy.spatial.transform import Rotation as R + + #import importlib.util + #spec = importlib.util.spec_from_file_location("d",path_input_dict) + #foo = importlib.util.module_from_spec(spec) + #spec.loader.exec_module(foo) + #d = foo.d + + eps = 1e-4 + nodes = np.loadtxt(path_nodes, delimiter=' ') + pattern = np.loadtxt(path_pattern, delimiter=' ') + + ## Preprocessing ## + # N is number of nodes in a direction (not num. segments!) + Nx = d['x_steps']+1 + Ny = d['y_steps']+1 + Nz = d['z_steps']+1 + + # dx is the length of segments + dx = d['x_step'] + dy = d['y_step'] + dz = d['z_step'] + + # electrode location + elec_loc = np.array([d['Implantation_coordinate_X'], + d['Implantation_coordinate_Y'], + d['Implantation_coordinate_Z']]) + + # eurler angles + alphas = d['alpha_array_glob'] + betas = d['beta_array_glob'] + gammas = d['gamma_array_glob'] + + # some good numbers about neuron array + n_dir = len(d["alpha_array_glob"]) # number of unique directions + n_axons_per_dir = Nx*Ny*Nz + n_axons = n_axons_per_dir*n_dir + n_nodes_per_axon = len(pattern)#len(nodes)//n_axons # total node per axon (ranvier and non-ranvier) + + # center axons (from positive octant) + CEN = nodes[:,:3].mean(axis=0) + nodes[:,:3] -= CEN + + ## Extracting centers and directions ## + + # information about each direction will be stored in a dict + directions_info = {} + + # information about each axons in a dataframe + centers = [] + direction_ids = [] + for axon_id in range(n_axons): + vec = np.array(nodes[(axon_id+1)*n_nodes_per_axon-1, :3]-\ + nodes[axon_id*n_nodes_per_axon, :3]) # axon vector + cen = vec/2. + nodes[axon_id*n_nodes_per_axon, :3] + centers.append(cen) # the center location + + # save direction if new + if axon_id%n_axons_per_dir==0: + id_ = axon_id//n_axons_per_dir + info = {'direction': np.round(np.array(vec/np.linalg.norm(vec)),5), + 'alpha': alphas[id_], + 'beta': betas[id_], + 'gamma': gammas[id_]} + directions_info[id_] = info + direction_ids += [id_]*n_axons_per_dir + + centers = np.round(np.array(centers), 5) + activations = np.load(path_activation) + data = {'center': list(centers), 'active':activations, 'dir_id': direction_ids} + df = pd.DataFrame(data = data) + + ## Making nii/png files ## + # centeral node + xc = np.arange(-(Nx-1)*dx/2, Nx*dx/2, dx) + yc = np.arange(-(Ny-1)*dy/2, Ny*dy/2, dy) + zc = np.arange(-(Nz-1)*dz/2, Nz*dz/2, dz) + centers = np.meshgrid(xc,yc,zc, indexing='ij') + +# Max = np.stack(df.center.values).max() +# Min = np.stack(df.center.values).min() +# M = max(Max, abs(Min)) + + for dir_id in sorted(directions_info.keys()): + + VTA = np.zeros_like(centers[0]) #shape = (Nx,Ny,Nz) + + alpha = directions_info[dir_id]['alpha'] + beta = directions_info[dir_id]['beta'] + gamma = directions_info[dir_id]['gamma'] + + rx = R.from_euler('x', alpha, degrees=True).as_dcm()#.as_matrix() + ry = R.from_euler('y', beta, degrees=True).as_dcm()#.as_matrix() + rz = R.from_euler('z', gamma, degrees=True).as_dcm()#.as_matrix() + R_tot = np.dot(np.dot(rz,ry), rx) + + # rotate the array + x,y,z = np.dot(R_tot, np.array(centers).transpose(1,2,0,3)) + + # make a ref direction (+y) +# ref = np.array([[0,0],[0,M*0.8],[0,0]]) +# ref = np.dot(R_tot,ref) + + # finding the activated cells + + active_cond = (df.dir_id ==dir_id) & (df.active!=0) + actives = np.stack(df[active_cond].center) + status = df[active_cond].active.reset_index(drop=True) + if activate_close_axons: + status = np.abs(status) + + for idx, row in enumerate(actives): + x0,y0,z0 = row + + cond_x = abs(x-x0)1: + aff_ = np.zeros((4,4)) + scl = np.diag([dx, dy, dz]) + aff_[:3,:3] = scl + aff_[:3,3] = np.array([center.min() for center in centers])+ elec_loc + nii = nib.nifti1.Nifti1Image(VTA.astype('int8'), aff_,) + nib.save(nii, path_target+'VTA4nii_dir'+str(dir_id)+'.nii') + + + VTA = VTA[:,0,:] # the section perpendicular to the direction + + vmin = (activate_close_axons-1) + pcm = plt.pcolormesh(xc, zc, VTA, shading='nearest', cmap='coolwarm', + edgecolors='gray', linewidths=1, vmin=vmin, vmax=1, + in_layout=False) + + ax = plt.gca() + ax.set_title('VTA profile prependicular to the direction '+ str(dir_id)) + ax.set_aspect('equal') + + divider = make_axes_locatable(plt.gca()) + cax = divider.append_axes("right", "5%", pad="3%") + plt.colorbar(pcm, cax=cax) + + plt.tight_layout() + plt.savefig(path_target+'Grid_VTA_dir'+str(dir_id)+'.png', dpi=300, + bbox_inches='tight') + plt.close() + +