From 3a99c0bf466e282c03f48421d12e61f310adb5c2 Mon Sep 17 00:00:00 2001 From: arashgmn Date: Sat, 15 May 2021 16:09:58 +0200 Subject: [PATCH 01/14] unstage my additional directories --- .gitignore | 2 ++ OSS-DBS - tutorial.pdf | Bin 1781190 -> 1784149 bytes 2 files changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index d948def..350dc19 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ access-token.txt local-build.sh +OSS_paltform/my_codes +OSS_paltform/extra diff --git a/OSS-DBS - tutorial.pdf b/OSS-DBS - tutorial.pdf index f3c81c7aaa7b172027e4792a03b218b3cb96da0b..3ba5b4b7405e3e7f163787bd4458d0170f079b9e 100644 GIT binary patch delta 2829 zcmb_dO>7%Q6rLDp!NmnyehL9;4n%b%lHJ+$uI;E)DoObffu>F>Ku)L}XA`@1yq3K# zX%B1yBh&+jA|!m`0s^V`N=TeI<-m=ZfAIo8aiX*FOJcaQaqG~r52p_tO&O=xf){3kQsB|Mg>Mp1ch-MQ zOm|9uCbFINUlWHz<5a!w2Q4mK0x9s5C%IQ&0$3_3T)U5J%w2AHTsw24;FSYjv~%N@ z!?UJ4Zrkv+T*nsNwbY_*Mi&>jG231Y6vUjj7I5Pgo}QYpOe++oBZSBZ%dtnfG0BU9 z3uHyajNfQCu&<2TmodlFe{Zo8YvO6K3Y=%S-m)?2*8{JPLy$*{wb>H$8CTzBc9D6A z;dt7t4m*ZV9T=-v=PW$KX7M?{!EoF7$jB>i4A$%9ZxM%$}pC7Xcc6)Fd+4?xE08NkV1 zx}u}S-#p!(W7}v4nt|xeN3|282lTJe%tghcMbNME0r8L^;*)1B?q&tgdQVCg4RR4(4+13 zg11bsmb^7yG%rVO{}_q11E2`U(T)M#5l(<`l}k_>-CTU$eHcx)5D z*zox8%9uBl5ie_VO&R;3Ht({KTr)qZf1NU8sU)Uv53@IjTS2KAtTnx5%;ap&7Ir?B z92q%%=K23Vqr>gVjg$518ywP9ppM9xdDDsMbN_lqhcv498`5_&qo-7zTtt1)Y~FWf z)BM%r{0@fp|2VV5?$F-ZUZM4V<#p9>#jF47x93K#xBs~L{dbkld^7*ddok+oswYwb oy8dD+fYgr~z+O}r+@l6i{d?AcdsQhqSKLq0N3^7K`E5!61L!iMOaK4? delta 63 zcmV~$NfCko006+@1_3$c& Date: Thu, 27 May 2021 17:29:40 +0200 Subject: [PATCH 02/14] some refactoring --- .gitignore | 4 +- OSS_platform/CSF_refinement_new.py | 7 +- OSS_platform/FEM_in_spectrum.py | 71 ++++++++++++---- OSS_platform/GUI_inp_dict.py | 36 ++++----- OSS_platform/Launcher_OSS_lite.py | 9 ++- OSS_platform/MRI_DTI_prep_new.py | 125 ++++++++++++++++------------- OSS_platform/Tissue_marking_new.py | 40 ++++----- 7 files changed, 176 insertions(+), 116 deletions(-) diff --git a/.gitignore b/.gitignore index 350dc19..ef5a27b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ access-token.txt local-build.sh -OSS_paltform/my_codes -OSS_paltform/extra +/OSS_platform/extra/ +OSS_platform/my_codes/ diff --git a/OSS_platform/CSF_refinement_new.py b/OSS_platform/CSF_refinement_new.py index aa1e79a..6d45a77 100644 --- a/OSS_platform/CSF_refinement_new.py +++ b/OSS_platform/CSF_refinement_new.py @@ -228,6 +228,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: @@ -238,9 +239,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/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<=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) @@ -199,18 +200,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 +218,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 +238,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) From a09bd83bc930d42c630fd6143e768f3f48379aae Mon Sep 17 00:00:00 2001 From: Arash Gol-Mohammadi <30870706+arashgmn@users.noreply.github.com> Date: Thu, 17 Jun 2021 19:30:54 +0200 Subject: [PATCH 03/14] Update Launcher_OSS_lite.py to resolve the conflict --- OSS_platform/Launcher_OSS_lite.py | 142 ++++++++++++++++-------------- 1 file changed, 78 insertions(+), 64 deletions(-) diff --git a/OSS_platform/Launcher_OSS_lite.py b/OSS_platform/Launcher_OSS_lite.py index 60ff463..ffab69f 100644 --- a/OSS_platform/Launcher_OSS_lite.py +++ b/OSS_platform/Launcher_OSS_lite.py @@ -79,27 +79,22 @@ def run_full_model(master_dict): if d["voxel_arr_MRI"]==0 and d["voxel_arr_DTI"]==1: print("MRI data is new, the DTI data will be reprocessed") - d["voxel_arr_DTI"]==0 + d["voxel_arr_DTI"]=0 #loading of meta data depending on the simulatation setup and state - if d["Init_neuron_model_ready"]==1: - _out = np.genfromtxt('Neuron_model_arrays/Neuron_model_misc.csv', delimiter=' ') - - [ranvier_nodes, para1_nodes, para2_nodes, inter_nodes, ranvier_length, para1_length, - para2_length, inter_length, deltax, diam_fib, n_Ranvier,ROI_radius,N_segm] = _out - - param_axon=[ranvier_nodes, para1_nodes, para2_nodes, inter_nodes, ranvier_length, - para1_length, para2_length, inter_length, deltax, diam_fib] - - if d["Neuron_model_array_prepared"]==0: - with open('Neuron_model_arrays/Neuron_param_class.file', "rb") as f: - Neuron_param = pickle.load(f) - if d["Neuron_model_array_prepared"]==1: - Neuron_param=0 #not needed for a pre-defined neuron array - if d["Name_prepared_neuron_array"][-3:]=='.h5': - n_segments_fib_diam_array=np.load('Neuron_model_arrays/Neuron_populations_misc.npy') - N_segm=n_segments_fib_diam_array[:,0] - N_segm=N_segm.astype(int) + #if d["Init_neuron_model_ready"]==1: + # [ranvier_nodes, para1_nodes, para2_nodes, inter_nodes, ranvier_length, para1_length, para2_length, inter_length, deltax, diam_fib,n_Ranvier,ROI_radius,N_segm]=np.genfromtxt('Neuron_model_arrays/Neuron_model_misc.csv', delimiter=' ') + # param_axon=[ranvier_nodes, para1_nodes, para2_nodes, inter_nodes, ranvier_length, para1_length, para2_length, inter_length, deltax, diam_fib] + + #if d["Neuron_model_array_prepared"]==0: + # with open('Neuron_model_arrays/Neuron_param_class.file', "rb") as f: + # Neuron_param = pickle.load(f) + #if d["Neuron_model_array_prepared"]==1: + # Neuron_param=0 #not needed for a pre-defined neuron array + # if d["Name_prepared_neuron_array"][-3:]=='.h5': + # n_segments_fib_diam_array=np.load('Neuron_model_arrays/Neuron_populations_misc.npy') + # N_segm=n_segments_fib_diam_array[:,0] + # N_segm=N_segm.astype(int) if d["Init_mesh_ready"]==1: with open('Meshes/Mesh_ind.file', "rb") as f: Domains = pickle.load(f) @@ -122,75 +117,94 @@ def run_full_model(master_dict): from CAD_Salome import build_brain_approx x_length,y_length,z_length=build_brain_approx(d,MRI_param) #also creates 'brain_subsitute.brep' Brain_shape_name='Brain_substitute.brep' + d["Brain_shape_name"] = Brain_shape_name if d["Brain_shape_name"]!=0: Brain_shape_name=d["Brain_shape_name"] #==================Initial neuron array generation====================# + + from Neural_array_processing import Neuron_array + N_array = Neuron_array(d, MRI_param) if d["Init_neuron_model_ready"]==0 and d["Neuron_model_array_prepared"]==0: print("----- Creating initial neuron array -----") + N_array.build_neuron_models() #builds a pattern model, if not provided, then builds a neuron array and stores in 'Neuron_model_arrays/All_neuron_models.csv' + with open('Neuron_model_arrays/Neuron_array_class.file', "wb") as f: + pickle.dump(N_array, f, pickle.HIGHEST_PROTOCOL) + elif d["Neuron_model_array_prepared"]==1 and d["Init_neuron_model_ready"]==0: + print("----- Creating initial neuron array from a provided neuron array -----") + N_array.process_external_array() # adjusts the prepared neuron array to the computational domain (only for brain substitutes!), then stores the nueron array in 'Neuron_model_arrays/All_neuron_models.csv' + with open('Neuron_model_arrays/Neuron_array_class.file', "wb") as f: + pickle.dump(N_array, f, pickle.HIGHEST_PROTOCOL) + else: + print("--- Initial neuron array was loaded\n") + with open('Neuron_model_arrays/Neuron_array_class.file', "rb") as f: + N_array = pickle.load(f) + + + #if d["Init_neuron_model_ready"]==0 and d["Neuron_model_array_prepared"]==0: + # print("----- Creating initial neuron array -----") - from Neuron_models_arangement_new import build_neuron_models - Neuron_param,param_axon,ROI_radius,N_segm=build_neuron_models(d,MRI_param) #builds a pattern model, if not provided, then builds a neuron array and stores in 'Neuron_model_arrays/All_neuron_models.csv'. Also, creates corresp. meta data. + # from Neuron_models_arangement_new import build_neuron_models + # Neuron_param,param_axon,ROI_radius,N_segm=build_neuron_models(d,MRI_param) #builds a pattern model, if not provided, then builds a neuron array and stores in 'Neuron_model_arrays/All_neuron_models.csv'. Also, creates corresp. meta data. - if d["Neuron_model_array_prepared"]==1 and d["Init_neuron_model_ready"]==0: - print("----- Creating initial neuron array from a provided neuron array -----") - from Neuron_models_arangement_new import create_meta_data_for_predefined_models, cut_models_by_domain + #if d["Neuron_model_array_prepared"]==1 and d["Init_neuron_model_ready"]==0: + # print("----- Creating initial neuron array from a provided neuron array -----") + # from Neuron_models_arangement_new import create_meta_data_for_predefined_models, cut_models_by_domain - cut_models_by_domain(d,Brain_shape_name,d["Name_prepared_neuron_array"]) #to adjust the prepared neuron array to the computational domain (only for brain substitutes!) - ROI_radius,param_axon,N_segm=create_meta_data_for_predefined_models(d,MRI_param) #shifts coordinates of the provided neuron array to the positive octant coord. and stores in 'Neuron_model_arrays/All_neuron_models.csv'. Also creates corresp. meta data. - Neuron_param=0 #This class is irrelevant if we have a full neuron array predefined + # cut_models_by_domain(d,Brain_shape_name,d["Name_prepared_neuron_array"]) #to adjust the prepared neuron array to the computational domain (only for brain substitutes!) + # ROI_radius,param_axon,N_segm=create_meta_data_for_predefined_models(d,MRI_param) #shifts coordinates of the provided neuron array to the positive octant coord. and stores in 'Neuron_model_arrays/All_neuron_models.csv'. Also creates corresp. meta data. + # Neuron_param=0 #This class is irrelevant if we have a full neuron array predefined - if d["Neuron_model_array_prepared"]==1 and d["Init_neuron_model_ready"]==1: - Neuron_param=0 #This class is irrelevant if we have a full neuron array predefined - print("--- Initial neuron array was loaded\n") + #if d["Neuron_model_array_prepared"]==1 and d["Init_neuron_model_ready"]==1: + # Neuron_param=0 #This class is irrelevant if we have a full neuron array predefined + # print("--- Initial neuron array was loaded\n") #if brain substitute is used, it will be enlarged to encompass previosly defined neuron array (if necessary) if Brain_shape_name=='Brain_substitute.brep': needs_a_rebuid=0 - if ROI_radius > (x_length/2): - d['Approximating_Dimensions'][0]=ROI_radius*2+0.1 + if N_array.ROI_radius > (x_length/2): + d['Approximating_Dimensions'][0]=N_array.ROI_radius*2+0.1 print("increasing length along x to encompass the neuron array\n") needs_a_rebuid=1 - if ROI_radius > (y_length/2): - d['Approximating_Dimensions'][1]=ROI_radius*2+0.1 + if N_array.ROI_radius > (y_length/2): + d['Approximating_Dimensions'][1]=N_array.ROI_radius*2+0.1 print("increasing length along y to encompass the neuron array\n") needs_a_rebuid=1 - if ROI_radius > (z_length/2): - d['Approximating_Dimensions'][2]=ROI_radius*2+0.1 + if N_array.ROI_radius > (z_length/2): + d['Approximating_Dimensions'][2]=N_array.ROI_radius*2+0.1 print("increasing length along z to encompass the neuron array\n") needs_a_rebuid=1 if needs_a_rebuid==1: x_length,y_length,z_length=build_brain_approx(d,MRI_param) #also creates 'brain_subsitute.brep' print("\n") - if ROI_radius > min((x_length/2),(y_length/2),(z_length/2)): - print("ROI_radius: ",ROI_radius) + if N_array.ROI_radius > min((x_length/2),(y_length/2),(z_length/2)): + print("ROI_radius: ",N_array.ROI_radius) print("ROI is still bigger than the computational domain.") raise SystemExit #===================Final geometry generation=========================# from CAD_Salome import build_final_geometry - Domains=build_final_geometry(d,MRI_param,Brain_shape_name,ROI_radius,cc_multicontact) #creates and discretizes the geometry with the implanted electrode, encapsulation layer and ROI, converts to the approp. format. The files are stored in Meshes/ + Domains=build_final_geometry(d,MRI_param,Brain_shape_name,N_array.ROI_radius,cc_multicontact) #creates and discretizes the geometry with the implanted electrode, encapsulation layer and ROI, converts to the approp. format. The files are stored in Meshes/ #===============Adjusting neuron array====================================# if d["Adjusted_neuron_model_ready"]==0: from Neuron_models_arangement_new import adjust_neuron_models - N_models=adjust_neuron_models(d,MRI_param,Domains,Neuron_param,param_axon) #subtracts neurons from the previously define All_neuron_models.csv, if the are non-physical (inside encap. layer or CSF, outside of the domain, intersect with the electrode geometry.) + if d["Init_mesh_ready"]==1: + with open('Neuron_model_arrays/Neuron_array_class.file', "rb") as f: + N_array = pickle.load(f) + + N_array.adjust_neuron_models(Domains, MRI_param) + with open('Neuron_model_arrays/Neuron_array_class.file', "wb") as f: + pickle.dump(N_array, f, pickle.HIGHEST_PROTOCOL) else: - if d["Neuron_model_array_prepared"]==1: - if d["Name_prepared_neuron_array"][-3:]=='.h5': #if imported with h5 - N_models = np.genfromtxt('Neuron_model_arrays/Adjusted_neuron_array_info.csv', delimiter=' ') - N_models=N_models.astype(int) - else: - N_models,points_csf,points_encap_and_float_contacts,points_outside=np.genfromtxt('Neuron_model_arrays/Adjusted_neuron_array_info.csv', delimiter=' ') - N_models=int(N_models) - else: - N_models,points_csf,points_encap_and_float_contacts,points_outside=np.genfromtxt('Neuron_model_arrays/Adjusted_neuron_array_info.csv', delimiter=' ') - N_models=int(N_models) - print("--- Neuron array meta data were loaded\n") - - number_of_points=int(np.sum(N_segm*N_models)) + print("--- Adjusted neuron array was loaded\n") + with open('Neuron_model_arrays/Neuron_array_class.file', "rb") as f: + N_array = pickle.load(f) + + number_of_points = int(np.sum(N_array.pattern['num_segments']*N_array.N_models)) + #number_of_points=int(np.sum(N_segm*N_models)) #subprocess.call('python Paraview_InitMesh_and_Neurons.py', shell=True) if d['Show_paraview_screenshots']==1: @@ -354,7 +368,7 @@ def run_full_model(master_dict): if d["IFFT_ready"] == 0 and d["Full_Field_IFFT"] == 1: from Full_IFFT_field_function import get_field_in_time - VTA_size = get_field_in_time(d,FR_vector_signal,Xs_signal_norm,t_vector,N_models,N_segm) # also uses data from Field_solutions_functions/ and if spectrum truncation is applied, than data from Stim_Signal/ + VTA_size = get_field_in_time(d,FR_vector_signal,Xs_signal_norm,t_vector,N_array.N_models,N_array.pattern['num_segments']) # also uses data from Field_solutions_functions/ and if spectrum truncation is applied, than data from Stim_Signal/ # if VTA_from_NEURON is enabled, will save pointwise solutions in Points_in_time/ d["IFFT_ready"] = 1 #modification of dictionary @@ -379,9 +393,9 @@ def run_full_model(master_dict): hf.close() for i in range(len(d["n_Ranvier"])): print("in ",lst_population_names[i]," population") - last_point=convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_models[i],N_segm[i],FR_vector_signal,t_vector,A,'Field_solutions/sorted_solution.csv',dif_axons=True,last_point=last_point) + last_point=convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_array.N_models[i],N_array.pattern['num_segments'][i],FR_vector_signal,t_vector,A,'Field_solutions/sorted_solution.csv',dif_axons=True,last_point=last_point) else: - convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_models,N_segm,FR_vector_signal,t_vector,A,'Field_solutions/sorted_solution.csv') + convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_array.N_models,N_array.pattern['num_segments'],FR_vector_signal,t_vector,A,'Field_solutions/sorted_solution.csv') else: print("Truncation of the obtained full solution is only for high. ampl and cutoff methods") @@ -399,9 +413,9 @@ def run_full_model(master_dict): hf.close() for i in range(len(d["n_Ranvier"])): print("in ",lst_population_names[i]," population") - last_point=convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_models[i],N_segm[i],FR_vector_signal,t_vector,A,name_sorted_solution,dif_axons=True,last_point=last_point) + last_point=convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_array.N_models[i],N_array.pattern['num_segments'][i],FR_vector_signal,t_vector,A,name_sorted_solution,dif_axons=True,last_point=last_point) else: - convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_models,N_segm,FR_vector_signal,t_vector,A,name_sorted_solution) + convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_array.N_models,N_array.pattern['num_segments'],FR_vector_signal,t_vector,A,name_sorted_solution) if (d["spectrum_trunc_method"]=='High Amplitude Method' or d["spectrum_trunc_method"]=='Cutoff Method') and d["Truncate_the_obtained_full_solution"]==0: if isinstance(d["n_Ranvier"],list): #if different populations @@ -411,9 +425,9 @@ def run_full_model(master_dict): hf.close() for i in range(len(d["n_Ranvier"])): print("in ",lst_population_names[i]," population") - last_point=convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm_new,N_models[i],N_segm[i],FR_vector_signal_new,t_vector,A,name_sorted_solution,dif_axons=True,last_point=last_point) + last_point=convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm_new,N_array.N_models[i],N_array.pattern['num_segments'][i],FR_vector_signal_new,t_vector,A,name_sorted_solution,dif_axons=True,last_point=last_point) else: - convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm_new,N_models,N_segm,FR_vector_signal_new,t_vector,A,name_sorted_solution) + convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm_new,N_array.N_models,N_array.pattern['num_segments'],FR_vector_signal_new,t_vector,A,name_sorted_solution) if d["spectrum_trunc_method"]=='Octave Band Method': if isinstance(d["n_Ranvier"],list): #if different populations @@ -423,9 +437,9 @@ def run_full_model(master_dict): hf.close() for i in range(len(d["n_Ranvier"])): print("in ",lst_population_names[i]," population") - last_point=convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_models[i],N_segm[i],FR_vector_signal,t_vector,A,name_sorted_solution,inx_st_oct=inx_start_octv,dif_axons=True,last_point=last_point) + last_point=convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_array.N_models[i],N_array.pattern['num_segments'][i],FR_vector_signal,t_vector,A,name_sorted_solution,inx_st_oct=inx_start_octv,dif_axons=True,last_point=last_point) else: - convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_models,N_segm,FR_vector_signal,t_vector,A,name_sorted_solution,inx_st_oct=inx_start_octv,dif_axons=False,last_point=0) + convolute_signal_with_field_and_compute_ifft(d,Xs_signal_norm,N_array.N_models,N_array.pattern['num_segments'],FR_vector_signal,t_vector,A,name_sorted_solution,inx_st_oct=inx_start_octv,dif_axons=False,last_point=0) @@ -454,12 +468,12 @@ def run_full_model(master_dict): os.chdir("Axon_files/") if d["Axon_Model_Type"] == 'Reilly2016': os.chdir("Reilly2016/") - last_point=N_segm[i]*N_models[i]+last_point + last_point=N_array.pattern['num_segments'][i]*N_array.N_models[i]+last_point os.chdir("..") if d["Axon_Model_Type"] == 'Reilly2016': os.chdir("..") else: - Number_of_activated=run_simulation_with_NEURON(0,-1,d["diam_fib"],1000*d["t_step"],1000.0/d["freq"],d["n_Ranvier"],N_models,d["v_init"],t_vector.shape[0],d["Ampl_scale"],d["number_of_processors"]) + Number_of_activated=run_simulation_with_NEURON(0,-1,d["diam_fib"],1000*d["t_step"],1000.0/d["freq"],d["n_Ranvier"],N_array.N_models,d["v_init"],t_vector.shape[0],d["Ampl_scale"],d["number_of_processors"]) if isinstance(d["n_Ranvier"],list) and len(d["n_Ranvier"])>1: with open(os.devnull, 'w') as FNULL: subprocess.call('python Visualization_files/Paraview_connections_activation.py', shell=True, stdout=FNULL, stderr=subprocess.STDOUT) From 0af3482008309e75ef110d5213dd6218d3de9e8a Mon Sep 17 00:00:00 2001 From: Arash Gol-Mohammadi <30870706+arashgmn@users.noreply.github.com> Date: Thu, 17 Jun 2021 19:32:15 +0200 Subject: [PATCH 04/14] Update MRI_DTI_prep_new.py to resolve the conflict --- OSS_platform/MRI_DTI_prep_new.py | 125 +++++++++++++------------------ 1 file changed, 54 insertions(+), 71 deletions(-) diff --git a/OSS_platform/MRI_DTI_prep_new.py b/OSS_platform/MRI_DTI_prep_new.py index 5ad2251..049e02e 100644 --- a/OSS_platform/MRI_DTI_prep_new.py +++ b/OSS_platform/MRI_DTI_prep_new.py @@ -49,8 +49,6 @@ def __init__(self,MRI_n,Mx,My,Mz,voxel_size_x,voxel_size_y,voxel_size_z,x_st,y_s def map_MRI(MRI_name,MRI_data_in_m,default_material,CSF_inx,WM_inx,GM_inx,from_grid_txt): # exctracts MRI data from the grid txt file (COMSOL format) start_voxel=time.clock() - # meter to mm ratio (used only if MRI is in meters) - r = 1000. if MRI_data_in_m else 1. if from_grid_txt==True: #kicks out strings with comments @@ -106,14 +104,18 @@ def map_MRI(MRI_name,MRI_data_in_m,default_material,CSF_inx,WM_inx,GM_inx,from_g tissue_array = img.get_fdata() voxel_array_temp=tissue_array.flatten('F') - #number of voxels along axes - Mx,My,Mz=(tissue_array.shape[0],tissue_array.shape[1],tissue_array.shape[2]) - - voxel_size_x=img.header.get_zooms()[0]*r - voxel_size_y=img.header.get_zooms()[1]*r - voxel_size_z=img.header.get_zooms()[2]*r - img_start_x,img_start_y,img_start_z=(img.affine[0,3]*r, img.affine[1,3]*r, img.affine[2,3]*r) + Mx,My,Mz=(tissue_array.shape[0],tissue_array.shape[1],tissue_array.shape[2]) #number of voxels along axes + if MRI_data_in_m==1: #switch to mm if the MRI data is in m + voxel_size_x=img.header.get_zooms()[0]*1000 + voxel_size_y=img.header.get_zooms()[1]*1000 + voxel_size_z=img.header.get_zooms()[2]*1000 + img_start_x,img_start_y,img_start_z=(img.affine[0,3]*1000,img.affine[1,3]*1000,img.affine[2,3]*1000) + else: + voxel_size_x=img.header.get_zooms()[0] + voxel_size_y=img.header.get_zooms()[1] + voxel_size_z=img.header.get_zooms()[2] + img_start_x,img_start_y,img_start_z=(img.affine[0,3],img.affine[1,3],img.affine[2,3]) x_arr=np.arange(img_start_x,img_start_x+voxel_size_x*Mx,voxel_size_x) y_arr=np.arange(img_start_y,img_start_y+voxel_size_y*My,voxel_size_y) @@ -157,11 +159,11 @@ def map_MRI(MRI_name,MRI_data_in_m,default_material,CSF_inx,WM_inx,GM_inx,from_g np.save('MRI_DTI_derived_data/Tissue_array_MRI', voxel_arr.astype('b'), allow_pickle=False, fix_imports=False) del voxel_arr,voxel_array_temp - #switch to mm if the MRI data is in m - x_arr[:] = [x_i*r for x_i in x_arr] - y_arr[:] = [y_i*r for y_i in y_arr] - z_arr[:] = [z_i*r for z_i in z_arr] - + if MRI_data_in_m==1: #switch to mm if the MRI data is in m + x_arr[:] = [x_i * 1000.0 for x_i in x_arr] + y_arr[:] = [y_i * 1000.0 for y_i in y_arr] + z_arr[:] = [z_i * 1000.0 for z_i in z_arr] + voxel_size_x=abs(round(x_arr[1]-x_arr[0],6)) #size of voxels along x-axis voxel_size_y=abs(round(y_arr[1]-y_arr[0],6)) voxel_size_z=abs(round(z_arr[1]-z_arr[0],6)) @@ -184,24 +186,16 @@ def map_MRI(MRI_name,MRI_data_in_m,default_material,CSF_inx,WM_inx,GM_inx,from_g return (Mx,My,Mz,round(min(x_arr),6),round(min(y_arr),6),round(min(z_arr),6),round(max(x_arr),6),round(max(y_arr),6),round(max(z_arr),6),voxel_size_x,voxel_size_y,voxel_size_z) -def map_DTI(d, DTI_name, DTI_data_in_m, from_grid_txt): - """ - exctracts Tensor data from the grid txt file (COMSOL format) - """ - +def map_DTI(d,DTI_name,DTI_data_in_m,from_grid_txt): # exctracts Tensor data from the grid txt file (COMSOL format) + start_voxel=time.clock() - # meter to mm ratio (used only if DTI is in meters) - r = 1000. if DTI_data_in_m else 1. - if from_grid_txt==True: infile = open(DTI_name,'r').readlines() with open('MRI_DTI_derived_data/Filtered_'+DTI_name,'w') as outfile: - # first, we check the length of z and y vectors (one line of - # DTI data corresponds to the same y and z coordinate, - # directions are separated by a comment) + #first, we check the length of z and y vectors (one line of DTI data corresponds to the same y and z coordinate, directions are separated by a comment) for index,line in enumerate(infile): if index==2: line=line.rstrip() @@ -305,13 +299,12 @@ def map_DTI(d, DTI_name, DTI_data_in_m, from_grid_txt): print("Error while processing DTI data, maybe not all slices were extracted") raise SystemExit - # converts DTI data to mm if it is not already - x_arr[:] = [x_i *r for x_i in x_arr] - y_arr[:] = [y_i *r for y_i in y_arr] - z_arr[:] = [z_i *r for z_i in z_arr] - - #size of voxel along xyz-axis - voxel_size_x=abs(round(x_arr[1]-x_arr[0],6)) + if DTI_data_in_m==1: #if DTI data is in m + x_arr[:] = [x_i * 1000.0 for x_i in x_arr] + y_arr[:] = [y_i * 1000.0 for y_i in y_arr] + z_arr[:] = [z_i * 1000.0 for z_i in z_arr] + + voxel_size_x=abs(round(x_arr[1]-x_arr[0],6)) #size of voxel along x-axis voxel_size_y=abs(round(y_arr[1]-y_arr[0],6)) voxel_size_z=abs(round(z_arr[1]-z_arr[0],6)) else: @@ -324,13 +317,13 @@ def map_DTI(d, DTI_name, DTI_data_in_m, from_grid_txt): if d["Brain_shape_name"]==0: print("Extracting a subset of the tensor data for the appox. volume") - res_x, res_y, res_z=(img.header.get_zooms()[0]*r, - img.header.get_zooms()[1]*r, - img.header.get_zooms()[2]*r) - img_start_x, img_start_y, img_start_z =(img.affine[0,3]*r, - img.affine[1,3]*r, - img.affine[2,3]*r) - + if DTI_data_in_m==1: + res_x,res_y,res_z=(img.header.get_zooms()[0]*1000.0,img.header.get_zooms()[1]*1000.0,img.header.get_zooms()[2]*1000.0) + img_start_x,img_start_y,img_start_z=(img.affine[0,3]*1000,img.affine[1,3]*1000,img.affine[2,3]*1000) + else: + res_x,res_y,res_z=(img.header.get_zooms()[0],img.header.get_zooms()[1],img.header.get_zooms()[2]) + img_start_x,img_start_y,img_start_z=(img.affine[0,3],img.affine[1,3],img.affine[2,3]) + start_vox_x=int((d['Implantation_coordinate_X']-img_start_x)/res_x)-int(d['Approximating_Dimensions'][0]/(2.0*res_x)) vox_window_x=int(d['Approximating_Dimensions'][0]/(res_x)) @@ -368,11 +361,17 @@ def map_DTI(d, DTI_name, DTI_data_in_m, from_grid_txt): raise SystemExit Tensor_array=np.zeros((voxel_arr_c11.shape[0],6),float) - - voxel_size_x=img.header.get_zooms()[0]*r - voxel_size_y=img.header.get_zooms()[1]*r - voxel_size_z=img.header.get_zooms()[2]*r - img_start_x,img_start_y,img_start_z=(img.affine[0,3]*r, img.affine[1,3]*r, img.affine[2,3]*r) + + if DTI_data_in_m==1: #switch to mm if the MRI data is in m + voxel_size_x=img.header.get_zooms()[0]*1000 + voxel_size_y=img.header.get_zooms()[1]*1000 + voxel_size_z=img.header.get_zooms()[2]*1000 + img_start_x,img_start_y,img_start_z=(img.affine[0,3]*1000,img.affine[1,3]*1000,img.affine[2,3]*1000) + else: + voxel_size_x=img.header.get_zooms()[0] + voxel_size_y=img.header.get_zooms()[1] + voxel_size_z=img.header.get_zooms()[2] + img_start_x,img_start_y,img_start_z=(img.affine[0,3],img.affine[1,3],img.affine[2,3]) x_arr=np.arange(img_start_x,img_start_x+voxel_size_x*Mx,voxel_size_x) y_arr=np.arange(img_start_y,img_start_y+voxel_size_y*My,voxel_size_y) @@ -392,7 +391,7 @@ def map_DTI(d, DTI_name, DTI_data_in_m, from_grid_txt): Tensor_array[:,:]=np.vstack((voxel_arr_c11[:],voxel_arr_c21[:],voxel_arr_c31[:],voxel_arr_c22[:],voxel_arr_c32[:],voxel_arr_c33[:])).T - + #'''Initialyly the data should be in ascending order''' #'''To ensure that coordinate vectors go in ascending order''' #voxel_array = voxel_array[voxel_array[:,0].argsort()] @@ -419,33 +418,21 @@ def obtain_MRI_class(inp_dict): if inp_dict["voxel_arr_MRI"]==0: # 1 if MRI data were already processed by the platform and corresp. meta data were created print("--- processing provided MRI data") if inp_dict["MRI_data_name"][-3:]=='nii' or inp_dict["MRI_data_name"][-6:]=='nii.gz': - from_grid_txt= False + [Mx_mri,My_mri,Mz_mri,x_min,y_min,z_min,x_max,y_max,z_max,MRI_voxel_size_x,MRI_voxel_size_y,MRI_voxel_size_z]=map_MRI(inp_dict["MRI_data_name"],inp_dict["MRI_in_m"],inp_dict["default_material"],inp_dict["CSF_index"],inp_dict["WM_index"],inp_dict["GM_index"],False) #will also prepare voxel_array! else: - from_grid_txt= True + [Mx_mri,My_mri,Mz_mri,x_min,y_min,z_min,x_max,y_max,z_max,MRI_voxel_size_x,MRI_voxel_size_y,MRI_voxel_size_z]=map_MRI(inp_dict["MRI_data_name"],inp_dict["MRI_in_m"],inp_dict["default_material"],inp_dict["CSF_index"],inp_dict["WM_index"],inp_dict["GM_index"],True) #will also prepare voxel_array! - [Mx_mri, My_mri, Mz_mri, x_min, y_min, z_min, x_max, y_max, z_max, - MRI_voxel_size_x, MRI_voxel_size_y, MRI_voxel_size_z]=map_MRI(inp_dict["MRI_data_name"],inp_dict["MRI_in_m"], - inp_dict["default_material"],inp_dict["CSF_index"], - inp_dict["WM_index"],inp_dict["GM_index"], from_grid_txt) #will also prepare voxel_array! - '''Save meta data for the future simulations with the current MRI data set''' - MRI_misc = np.array([Mx_mri,My_mri,Mz_mri,x_min,y_min,z_min,x_max,y_max,z_max, - MRI_voxel_size_x,MRI_voxel_size_y,MRI_voxel_size_z]) + MRI_misc=np.array([Mx_mri,My_mri,Mz_mri,x_min,y_min,z_min,x_max,y_max,z_max,MRI_voxel_size_x,MRI_voxel_size_y,MRI_voxel_size_z]) np.savetxt('MRI_DTI_derived_data/MRI_misc.csv', MRI_misc, delimiter=" ") print("--- MRI meta data were created\n") else: - [Mx_mri,My_mri, Mz_mri, - x_min, y_min, z_min, x_max, y_max, z_max, - MRI_voxel_size_x, MRI_voxel_size_y, MRI_voxel_size_z] =np.genfromtxt('MRI_DTI_derived_data/MRI_misc.csv', delimiter=' ') + [Mx_mri,My_mri,Mz_mri,x_min,y_min,z_min,x_max,y_max,z_max,MRI_voxel_size_x,MRI_voxel_size_y,MRI_voxel_size_z]=np.genfromtxt('MRI_DTI_derived_data/MRI_misc.csv', delimiter=' ') print("--- MRI meta data were loaded\n") x_shift,y_shift,z_shift=(-1*(x_min),-1*(y_min),-1*(z_min)) #shift of MRI to have it in the positive octant and start in (0,0,0) - MRI_param = MRI_info(inp_dict["MRI_data_name"], - Mx_mri, My_mri, Mz_mri, - MRI_voxel_size_x, MRI_voxel_size_y, MRI_voxel_size_z, - x_min, y_min, z_min, x_max, y_max, z_max, - x_shift, y_shift, z_shift) + MRI_param=MRI_info(inp_dict["MRI_data_name"],Mx_mri,My_mri,Mz_mri,MRI_voxel_size_x,MRI_voxel_size_y,MRI_voxel_size_z,x_min,y_min,z_min,x_max,y_max,z_max,x_shift,y_shift,z_shift) with open('MRI_DTI_derived_data/MRI_class.file', "wb") as f: pickle.dump(MRI_param, f, pickle.HIGHEST_PROTOCOL) @@ -455,17 +442,13 @@ def obtain_DTI_class(inp_dict,MRI_param): if inp_dict["voxel_arr_DTI"]==0: #1 if DTI data were already processed by the platform and corresp. meta data were created if inp_dict["DTI_data_name"][-3:]=='nii' or inp_dict["DTI_data_name"][-6:]=='nii.gz': - from_grid_txt= False + [Mx_dti,My_dti,Mz_dti,x_min_dti,y_min_dti,z_min_dti,DTI_voxel_size_x,DTI_voxel_size_y,DTI_voxel_size_z]=map_DTI(inp_dict,inp_dict["DTI_data_name"],inp_dict["MRI_in_m"],False) else: - from_grid_txt=True - - [Mx_dti, My_dti, Mz_dti, x_min_dti, y_min_dti, z_min_dti, - DTI_voxel_size_x, DTI_voxel_size_y, DTI_voxel_size_z]= map_DTI(inp_dict,inp_dict["DTI_data_name"], - inp_dict["MRI_in_m"], from_grid_txt) + [Mx_dti,My_dti,Mz_dti,x_min_dti,y_min_dti,z_min_dti,DTI_voxel_size_x,DTI_voxel_size_y,DTI_voxel_size_z]=map_DTI(inp_dict,inp_dict["DTI_data_name"],inp_dict["MRI_in_m"],True) - x_start_dti= x_min_dti-MRI_param.x_min #DTI can be shifted from the MRI origin (0,0,0) (but only to the positive direction). - y_start_dti= y_min_dti-MRI_param.y_min - z_start_dti= z_min_dti-MRI_param.z_min + x_start_dti=x_min_dti-MRI_param.x_min #DTI can be shifted from the MRI origin (0,0,0) (but only to the positive direction). + y_start_dti=y_min_dti-MRI_param.y_min + z_start_dti=z_min_dti-MRI_param.z_min if x_start_dti<0 or y_start_dti<0 or z_start_dti<0: print("DTI data is not in positive octant, exiting.") From b8d2c66c69482ba4d8a21710674efe82a127b9c2 Mon Sep 17 00:00:00 2001 From: arashgmn Date: Wed, 23 Jun 2021 18:21:15 +0200 Subject: [PATCH 05/14] corrected the Axon class for McIntyre's model --- OSS_platform/Axon_files/axon.py | 86 ++++++++++++++++++++------------- 1 file changed, 53 insertions(+), 33 deletions(-) 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 From 0fd0591ad2e7a11e9d1d26322be0290f62b351ff Mon Sep 17 00:00:00 2001 From: arashgmn Date: Wed, 23 Jun 2021 18:48:47 +0200 Subject: [PATCH 06/14] more refactoring --- OSS_platform/Neural_array_processing.py | 14 ++-- OSS_platform/Neuron_models_arangement_new.py | 78 +++++++++++++++----- OSS_platform/Tissue_marking_new.py | 14 +--- 3 files changed, 68 insertions(+), 38 deletions(-) diff --git a/OSS_platform/Neural_array_processing.py b/OSS_platform/Neural_array_processing.py index e274497..5b1cede 100644 --- a/OSS_platform/Neural_array_processing.py +++ b/OSS_platform/Neural_array_processing.py @@ -52,14 +52,13 @@ def __init__(self, input_dict, MRI_param): #for better readability, resaving t 'Dist_seeding_steps': [input_dict['x_step'], input_dict['y_step'], input_dict['z_step']], # distance between axons (center nodes) along the axes 'X_angles_glob' : input_dict['alpha_array_glob'], # list containing rotational angles around X axis, already converted to radians in Dict_corrector.py 'Y_angles_glob' : input_dict['beta_array_glob'], - 'Z_angles_glob' : input_dict['gamma_array_glob'], } + 'Z_angles_glob' : input_dict['gamma_array_glob']} else: self.Type = 'Custom' # custom allocation of neurons according to entries 'X_coord_old', ... and 'YZ_angles', ... in GUI_inp_dict.py self.custom_structure = {'Center_coordinates' : [input_dict['X_coord_old'], input_dict['Y_coord_old'], input_dict['Z_coord_old']] , # list of lists [[x1,x2,x2],[y1,y2,y3],...], usually at the electrode's tip or the last contact 'X_angles_loc' : input_dict['YZ_angles'], # list containing rotational angles around X axis. IMPORTANT: the neurons are first centered at O(0,0,0), rotated, and then translated to their center coordinates 'Y_angles_loc' : input_dict['ZX_angles'], # already converted to radians in Dict_corrector.py - 'Z_angles_loc' : input_dict['XY_angles'], - } + 'Z_angles_loc' : input_dict['XY_angles']} def rotate_globally(self,point_coords,inx_angle): @@ -164,11 +163,13 @@ def create_structured_array(self): # contain coordinates of the VTA axons centered at (0,0,0) self.N_models_in_plane=(self.VTA_structure['N_seeding_steps'][0] + 1)*(self.VTA_structure['N_seeding_steps'][1] + 1)*(self.VTA_structure['N_seeding_steps'][2] + 1) - self.N_total_of_axons = (self.VTA_structure['N_seeding_steps'][0] + 1)*(self.VTA_structure['N_seeding_steps'][1] + 1)*(self.VTA_structure['N_seeding_steps'][2] + 1)*len(self.VTA_structure["X_angles_glob"]) #+1 because we have the initial axon + self.N_total_of_axons =(self.VTA_structure['N_seeding_steps'][0] + 1)*(self.VTA_structure['N_seeding_steps'][1] + 1)*(self.VTA_structure['N_seeding_steps'][2] + 1)*len(self.VTA_structure["X_angles_glob"]) #+1 because we have the initial axon self.VTA_seeds_O = np.zeros((self.N_total_of_axons,3),float) # computed as x_start_point=0-(d["x_step"]*(d["x_steps"])/2) - start_point = [-1*self.VTA_structure['Dist_seeding_steps'][0] * self.VTA_structure['N_seeding_steps'][0] / 2 , -1*self.VTA_structure['Dist_seeding_steps'][1] * self.VTA_structure['N_seeding_steps'][1] / 2 , -1*self.VTA_structure['Dist_seeding_steps'][2] * self.VTA_structure['N_seeding_steps'][2] / 2] + start_point = [-1*self.VTA_structure['Dist_seeding_steps'][0] * self.VTA_structure['N_seeding_steps'][0] / 2 , + -1*self.VTA_structure['Dist_seeding_steps'][1] * self.VTA_structure['N_seeding_steps'][1] / 2 , + -1*self.VTA_structure['Dist_seeding_steps'][2] * self.VTA_structure['N_seeding_steps'][2] / 2] # seed axons for one direction x_one_dir=[] @@ -193,9 +194,6 @@ def create_structured_array(self): gl_ind+=1 - - - def get_neuron_morhology(self,fib_diam,N_Ranvier): # pass fib_diam and N_Ranvier explicitly, because they might be from a list """ 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 c9ef3be..9e72d15 100644 --- a/OSS_platform/Tissue_marking_new.py +++ b/OSS_platform/Tissue_marking_new.py @@ -80,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) @@ -208,7 +200,7 @@ def get_cellmap_tensors(mesh,subdomains_assigned,Domains,MRI_param,DTI_param,def z_coord=cell.midpoint().z() - xv_mri=x_coord//(voxel_size_x-eps) #defines number of steps to get to the voxels containing x[0] coordinate + 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 From 2a8ea43746a9130583cd4afa591f6a85c7cce4ba Mon Sep 17 00:00:00 2001 From: arashgmn Date: Wed, 23 Jun 2021 18:51:32 +0200 Subject: [PATCH 07/14] update gitignore to skip pycache folders --- .gitignore | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index ef5a27b..a3fcdbd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ access-token.txt local-build.sh -/OSS_platform/extra/ -OSS_platform/my_codes/ +__pycache__/ +*/*/__pychach__/ From 39be161374f1655429cdd59a4c0810aa1c664581 Mon Sep 17 00:00:00 2001 From: arashgmn Date: Wed, 23 Jun 2021 18:59:47 +0200 Subject: [PATCH 08/14] restoring to the default input_dict --- OSS_platform/GUI_inp_dict.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/OSS_platform/GUI_inp_dict.py b/OSS_platform/GUI_inp_dict.py index c9eca10..cf6aeb8 100644 --- a/OSS_platform/GUI_inp_dict.py +++ b/OSS_platform/GUI_inp_dict.py @@ -45,23 +45,23 @@ 'encap_scaling_perm': 0.8, 'pattern_model_name': '0', 'diam_fib': [5.7], - 'n_Ranvier': [2], + 'n_Ranvier': [21], 'v_init': -80.0, 'Neuron_model_array_prepared': 0, 'Name_prepared_neuron_array': '0', 'Global_rot': 1, - 'x_seed': 100, #10.929, - 'y_seed': 100, #-12.117, - 'z_seed': 100, #-7.697, + 'x_seed': 10.929, + 'y_seed': -12.117, + 'z_seed': -7.697, 'x_steps': 10, # number of nodes around x 'y_steps': 10, # number of nodes around y 'z_steps': 10, # number of nodes around z 'x_step': 1, # dx [mm] 'y_step': 1, # dy [mm] 'z_step': 1, # dz [mm] - 'alpha_array_glob': [0, 30, 0, 30, 30, 0, 30], - 'beta_array_glob': [0, 0, 0, 30, 0, 30, 30],#[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], - 'gamma_array_glob': [0, 0, 30, 0, 30, 30, 30],#[0.0, 45.0, 90.0, 135.0, 0.0, 0.0, 0.0], + 'alpha_array_glob': [0.0, 0.0, 0.0, 0.0, 45, 90, 135], + 'beta_array_glob': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + 'gamma_array_glob': [0.0, 45.0, 90.0, 135.0, 0.0, 0.0, 0.0], 'X_coord_old': 0, 'Y_coord_old': 0, 'Z_coord_old': 0, @@ -69,7 +69,7 @@ 'ZX_angles': [0], 'XY_angles': [0], 'EQS_core': 'QS', - 'Skip_mesh_refinement': 1, + 'Skip_mesh_refinement': 0, 'refinement_frequency': [130.0], 'num_ref_freqs': -1, 'rel_div_CSF': 5.0, From 9336d9f0d00a15dfdd27d866be4bd81fb62caf7b Mon Sep 17 00:00:00 2001 From: arashgmn Date: Thu, 24 Jun 2021 11:11:56 +0200 Subject: [PATCH 09/14] fixed assymetric pattern in Neural_array_processing + updated gitignore --- .gitignore | 4 +- OSS_platform/Neural_array_processing.py | 96 ++++++++++++++++--------- 2 files changed, 63 insertions(+), 37 deletions(-) diff --git a/.gitignore b/.gitignore index a3fcdbd..17c48f4 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ access-token.txt local-build.sh -__pycache__/ -*/*/__pychach__/ +OSS_platform/**/__pycache__ + diff --git a/OSS_platform/Neural_array_processing.py b/OSS_platform/Neural_array_processing.py index 5b1cede..c261d24 100644 --- a/OSS_platform/Neural_array_processing.py +++ b/OSS_platform/Neural_array_processing.py @@ -241,51 +241,77 @@ def generate_pattern(self): self.pattern['Array_coord_pattern_O'] = np.zeros((self.pattern['num_segments'],3),float) # _O refers to that fact that the pattern model is centered at O(0,0,0) # first compartment (not the middle!) - self.pattern['Array_coord_pattern_O'][0,0] = 0.0 - self.pattern['Array_coord_pattern_O'][0,1] = -0.001 * nr["deltax"] *(self.pattern['num_Ranvier'] - 1) / 2.0 # deltax is in µm, therefore scaling (OSS-DBS operates with mm) - self.pattern['Array_coord_pattern_O'][0,2] = 0.0 + # self.pattern['Array_coord_pattern_O'][0,0] = 0.0 + # self.pattern['Array_coord_pattern_O'][0,1] = -0.001 * nr["deltax"] *(self.pattern['num_Ranvier'] - 1) / 2.0 # deltax is in µm, therefore scaling (OSS-DBS operates with mm) + # self.pattern['Array_coord_pattern_O'][0,2] = 0.0 - loc_inx = 1 #because first node (with index 0) was already seeded + # loc_inx = 1 #because first node (with index 0) was already seeded if self.neuron_model == 'McIntyre2002': - 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 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 + + 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'] - if inx % n_comp == 0: - loc_inx = 0 + 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=" ") From b10263a4b2ca316c9d7e595c00c7f9b7a8ca9e42 Mon Sep 17 00:00:00 2001 From: arashgmn Date: Mon, 28 Jun 2021 19:58:49 +0200 Subject: [PATCH 10/14] cylindrical array generation + visualization of grid/cylindrical vta --- OSS_platform/GUI_inp_dict.py | 142 +++++++------- OSS_platform/VTA/cylindrical_array.py | 203 ++++++++++++++++++++ OSS_platform/VTA/visualize.py | 260 ++++++++++++++++++++++++++ 3 files changed, 537 insertions(+), 68 deletions(-) create mode 100644 OSS_platform/VTA/cylindrical_array.py create mode 100644 OSS_platform/VTA/visualize.py diff --git a/OSS_platform/GUI_inp_dict.py b/OSS_platform/GUI_inp_dict.py index cf6aeb8..6a9191e 100644 --- a/OSS_platform/GUI_inp_dict.py +++ b/OSS_platform/GUI_inp_dict.py @@ -16,73 +16,46 @@ 'Adapted_mesh_ready': 0, 'signal_generation_ready': 0, 'Parallel_comp_ready': 0, - 'Parallel_comp_interrupted': 0, + 'Parallel_comp_interrupted':0, 'IFFT_ready': 0, - 'MRI_data_name': 'icbm_avg_152_segmented.nii.gz', + 'MRI_data_name': 'segmask.nii', 'MRI_in_m': 0, 'DTI_data_name': '', - 'DTI_in_m': 0, # Boolean, True if DTI data is in meters - 'CSF_index': 1.0, - 'WM_index': 3.0, - 'GM_index': 2.0, - 'default_material': 3, + 'DTI_in_m': 0, + 'CSF_index': 3.0, + 'WM_index': 2.0, + 'GM_index': 1.0, + 'default_material': 1, 'Electrode_type': 'Medtronic3389', 'Brain_shape_name': '0', - 'x_length': 40.0, - 'y_length': 40.0, - 'z_length': 40.0, - 'Aprox_geometry_center': [10.92957028, -12.11697637, -7.69744601], - 'Implantation_coordinate_X': 10.929, - 'Implantation_coordinate_Y': -12.117, - 'Implantation_coordinate_Z': -7.697, - 'Second_coordinate_X': 10.929, - 'Second_coordinate_Y': -9.437, - 'Second_coordinate_Z': 3.697, + 'Aprox_geometry_center': [11.9418, 24.0083 ,13.7399], + 'Approximating_Dimensions': [80.0, 80.0, 80.0], + 'Implantation_coordinate_X': 7.9418, + 'Implantation_coordinate_Y': 24.0083, + 'Implantation_coordinate_Z': 10.7399, + 'Second_coordinate_X': 14.2483, + 'Second_coordinate_Y': 26.1642, + 'Second_coordinate_Z': 18.9309, 'Rotation_Z': 0.0, - 'encap_thickness': 0.20, + 'encap_thickness': 0.1, 'encap_tissue_type': 2, - 'encap_scaling_cond': 0.8, - 'encap_scaling_perm': 0.8, + 'encap_scaling_cond': 0.9, + 'encap_scaling_perm': 1.0, 'pattern_model_name': '0', + 'Axon_Model_Type': 'McIntyre2002', 'diam_fib': [5.7], - 'n_Ranvier': [21], + 'n_Ranvier': [3], 'v_init': -80.0, 'Neuron_model_array_prepared': 0, - 'Name_prepared_neuron_array': '0', - 'Global_rot': 1, - 'x_seed': 10.929, - 'y_seed': -12.117, - 'z_seed': -7.697, - 'x_steps': 10, # number of nodes around x - 'y_steps': 10, # number of nodes around y - 'z_steps': 10, # number of nodes around z - 'x_step': 1, # dx [mm] - 'y_step': 1, # dy [mm] - 'z_step': 1, # dz [mm] - 'alpha_array_glob': [0.0, 0.0, 0.0, 0.0, 45, 90, 135], - 'beta_array_glob': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], - 'gamma_array_glob': [0.0, 45.0, 90.0, 135.0, 0.0, 0.0, 0.0], - 'X_coord_old': 0, - 'Y_coord_old': 0, - 'Z_coord_old': 0, - 'YZ_angles': [0], - 'ZX_angles': [0], - 'XY_angles': [0], + 'Name_prepared_neuron_array': 'Allocated_axons.h5', 'EQS_core': 'QS', - 'Skip_mesh_refinement': 0, - 'refinement_frequency': [130.0], - 'num_ref_freqs': -1, - 'rel_div_CSF': 5.0, - 'CSF_frac_div': 5.0, - 'Min_Scaling': 1.0, - 'CSF_ref_reg': 10.0, - 'rel_div': 5.0, - 'Adaptive_frac_div': 5.0, - 'rel_div_current': 5.0, + 'Skip_mesh_refinement': 1, 'el_order': 2, - 'number_of_processors': 2, - 'current_control': 1, - 'Phi_vector': [0.0015, None, 0.0, None], + 'number_of_processors': 4, + 'FEniCS_MPI': 0, + 'current_control': 0, + 'Phi_vector': [None, 5, None, 0], + 'Solver_Type': 'GMRES', 'freq': 130.0, 'T': 60.0, 't_step': 1.0, @@ -90,23 +63,56 @@ 'Signal_type': 'Rectangle', 'Ampl_scale': 1.0, 'CPE_activ': 0, - 'beta': 0.0, + 'Full_Field_IFFT': 0, + 'spectrum_trunc_method': 'Octave Band Method', + 'trunc_param': 1000.0, + 'Truncate_the_obtained_full_solution': 0, + 'external_grounding': 0, + 'Stim_side': 0, + 'stretch': 1.0126193765016878, +'beta': 0.0, 'K_A': 0.0, 'beta_ground': 0.0, 'K_A_ground': 0.0, - 'Full_Field_IFFT': 0, + 'Global_rot': 1, + 'x_seed': 7.9418+5, # 2.5 mm away from the implantation site + 'y_seed': 24.0083+5, # 2.5 mm away from the implantation site + 'z_seed': 10.7399+5, # 2.5 mm away from the implantation site + 'x_steps': 20, + 'y_steps': 10, + 'z_steps': 20, + 'x_step': 0.5, + 'y_step': 1, + 'z_step': 1, + 'alpha_array_glob': [0, 20, 50, 0, 50, 20], + 'beta_array_glob': [0, 0, 0, 0, 0, 20], + 'gamma_array_glob': [0, 45, 50, 50, 0, 0], + 'X_coord_old': 0, + 'Y_coord_old': 0, + 'Z_coord_old': 0, + 'YZ_angles': [0], + 'ZX_angles': [0], + 'XY_angles': [0], 't_step_end': 1200, 'VTA_from_divE': 0, 'VTA_from_NEURON': 0, - 'VTA_from_E': 0, - 'Activation_threshold_VTA': 0, - 'spectrum_trunc_method': 'High Amplitude Method', - 'trunc_param': 5, - 'Truncate_the_obtained_full_solution': 0, - 'Show_paraview_screenshots': 0, - - 'Solver_Type': 'GMRES', - 'FEniCS_MPI': 0, - 'Axon_Model_Type': 'McIntyre2002', - 'Approximating_Dimensions': [40.0, 40.0, 40.0], - } + 'VTA_from_E': 1, + 'Activation_threshold_VTA': 0.2, + 'refinement_frequency': [0], + 'num_ref_freqs': 1, + 'rel_div_CSF': -1, + 'Adaptive_frac_div': 1.0, + 'Min_Scaling': 1.0, + 'CSF_ref_reg': 0.0, + 'rel_div': 5.0, + 'rel_div_current': 1.0, + 'VTA_type':'grid', # cylinder or grid + 'Nx': 40, # number of radial seeds + 'Nz_p': 100, # number of seeds above the implantation site + 'Nz_n':0, # number of seeds below the implantation site + 'dx':.1, # radial distance between seeds + 'dz':.1, # vertical distance between seeds + 'xmin':0.2, # minimum radial distance from the lead + 'z0': 1.5, # seeding starting level from the tip of the lead upward in S coordinates (mm) + 'n_dirs':5, # number of additional surfaces around the lead axis +} 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..af605f0 --- /dev/null +++ b/OSS_platform/VTA/visualize.py @@ -0,0 +1,260 @@ +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): + + + + args = {'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(**args) + else: + plot_grid_vta(**args) + + +def plot_cylindircal_vta(input_dict = '../GUI_inp_dict.py', + path_activation = '../Field_solutions/Activation/Network_status.npy', + path_target = '../Images/'): + """ + 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(figsize=(15,4)) + + 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) + plt.close() + + +def plot_grid_vta(path_target = '../Images/', + path_input_dict = '../GUI_input_dict.py', + path_nodes = '../Neuron_model_arrays/All_neuron_models.csv', + path_pattern = '../Neuron_model_arrays/default_pattern.csv', + path_activation = '../Field_solutions/Activation/Network_status.npy', + activate_close_axons=False): + """ + 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] = np.dot(R_tot, 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) + plt.close() + + From 3a4fa2dc0835cb1f67f484fc9d5fb2b5e737acf4 Mon Sep 17 00:00:00 2001 From: arashgmn Date: Tue, 29 Jun 2021 20:45:36 +0200 Subject: [PATCH 11/14] fixed the nifti affine transfromation for LeadDBS --- OSS_platform/VTA/visualize.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/OSS_platform/VTA/visualize.py b/OSS_platform/VTA/visualize.py index af605f0..cf4ed01 100644 --- a/OSS_platform/VTA/visualize.py +++ b/OSS_platform/VTA/visualize.py @@ -232,7 +232,7 @@ def plot_grid_vta(path_target = '../Images/', if Ny>1: aff_ = np.zeros((4,4)) scl = np.diag([dx, dy, dz]) - aff_[:3,:3] = np.dot(R_tot, scl) # ? + 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') From 2ceb0644e6d3587155251d5ae1b9da2ca109078a Mon Sep 17 00:00:00 2001 From: arashgmn Date: Sun, 11 Jul 2021 12:13:48 +0200 Subject: [PATCH 12/14] axisymmetric array generator and its visualzier --- OSS_platform/Dict_corrector.py | 4 +- OSS_platform/GUI_inp_dict.py | 30 +++++------ OSS_platform/Neural_array_processing.py | 71 ++++++++++++++++++------- 3 files changed, 69 insertions(+), 36 deletions(-) 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/GUI_inp_dict.py b/OSS_platform/GUI_inp_dict.py index 6a9191e..3a800c0 100644 --- a/OSS_platform/GUI_inp_dict.py +++ b/OSS_platform/GUI_inp_dict.py @@ -8,7 +8,7 @@ d = { 'voxel_arr_MRI': 1, - 'voxel_arr_DTI': 1, + 'voxel_arr_DTI': 0, 'Init_neuron_model_ready': 0, 'Init_mesh_ready': 0, 'Adjusted_neuron_model_ready': 0, @@ -18,7 +18,7 @@ 'Parallel_comp_ready': 0, 'Parallel_comp_interrupted':0, 'IFFT_ready': 0, - 'MRI_data_name': 'segmask.nii', + 'MRI_data_name': 'icbm_avg_152_segmented.nii.gz', 'MRI_in_m': 0, 'DTI_data_name': '', 'DTI_in_m': 0, @@ -28,14 +28,14 @@ 'default_material': 1, 'Electrode_type': 'Medtronic3389', 'Brain_shape_name': '0', - 'Aprox_geometry_center': [11.9418, 24.0083 ,13.7399], + 'Aprox_geometry_center': [11.7631, -13.9674 , -8.8594], 'Approximating_Dimensions': [80.0, 80.0, 80.0], - 'Implantation_coordinate_X': 7.9418, - 'Implantation_coordinate_Y': 24.0083, - 'Implantation_coordinate_Z': 10.7399, - 'Second_coordinate_X': 14.2483, - 'Second_coordinate_Y': 26.1642, - 'Second_coordinate_Z': 18.9309, + 'Implantation_coordinate_X': 11.7631, + 'Implantation_coordinate_Y': -13.9674 , + 'Implantation_coordinate_Z': -8.8594, + 'Second_coordinate_X': 13.7358, + 'Second_coordinate_Y': -10.0639, + 'Second_coordinate_Z': -3.7747, 'Rotation_Z': 0.0, 'encap_thickness': 0.1, 'encap_tissue_type': 2, @@ -51,7 +51,7 @@ 'EQS_core': 'QS', 'Skip_mesh_refinement': 1, 'el_order': 2, - 'number_of_processors': 4, + 'number_of_processors': 2, 'FEniCS_MPI': 0, 'current_control': 0, 'Phi_vector': [None, 5, None, 0], @@ -75,15 +75,15 @@ 'beta_ground': 0.0, 'K_A_ground': 0.0, 'Global_rot': 1, - 'x_seed': 7.9418+5, # 2.5 mm away from the implantation site - 'y_seed': 24.0083+5, # 2.5 mm away from the implantation site - 'z_seed': 10.7399+5, # 2.5 mm away from the implantation site + 'x_seed': 11.7631, #7.9418, # 2.5 mm away from the implantation site + 'y_seed': -13.9674,#24.0083, # 2.5 mm away from the implantation site + 'z_seed': -8.8594,#10.7399, # 2.5 mm away from the implantation site 'x_steps': 20, 'y_steps': 10, 'z_steps': 20, - 'x_step': 0.5, + 'x_step': 1, 'y_step': 1, - 'z_step': 1, + 'z_step': 0.5, 'alpha_array_glob': [0, 20, 50, 0, 50, 20], 'beta_array_glob': [0, 0, 0, 0, 0, 20], 'gamma_array_glob': [0, 45, 50, 50, 0, 0], diff --git a/OSS_platform/Neural_array_processing.py b/OSS_platform/Neural_array_processing.py index c261d24..3a3ac76 100644 --- a/OSS_platform/Neural_array_processing.py +++ b/OSS_platform/Neural_array_processing.py @@ -16,12 +16,14 @@ import h5py import numpy as np from dolfin import * +from numpy.core.fromnumeric import prod from pandas import read_csv from Axon_files.axon import Axon import time import os import pickle from scipy.spatial import distance +from VTA.cylindrical_array import generate_nodes_around_lead class Neuron_array(object): def __init__(self, input_dict, MRI_param): #for better readability, resaving the input dictionary to more explicit variables @@ -47,12 +49,23 @@ def __init__(self, input_dict, MRI_param): #for better readability, resaving t self.Type = 'Imported' elif input_dict['Global_rot'] == 1: # grid allocation of neurons for VTA according self.Type = 'VTA' + self.VTA_type = input_dict['VTA_type'] self.VTA_structure = {'Seed_coordinates' : [input_dict['x_seed'], input_dict['y_seed'], input_dict['z_seed']] , # usually at the electrode's tip or the last contact 'N_seeding_steps' : [input_dict['x_steps'], input_dict['y_steps'], input_dict['z_steps']], # how many additional axons will be placed along the axes (not counting the seeded one) 'Dist_seeding_steps': [input_dict['x_step'], input_dict['y_step'], input_dict['z_step']], # distance between axons (center nodes) along the axes 'X_angles_glob' : input_dict['alpha_array_glob'], # list containing rotational angles around X axis, already converted to radians in Dict_corrector.py 'Y_angles_glob' : input_dict['beta_array_glob'], - 'Z_angles_glob' : input_dict['gamma_array_glob']} + 'Z_angles_glob' : input_dict['gamma_array_glob'], + 'Nx' : input_dict['Nx'], + 'Nz_p' : input_dict['Nz_p'], + 'Nz_n' : input_dict['Nz_n'] , + 'dx' : input_dict['dx'], + 'dz' : input_dict['dz'], + 'xmin' : input_dict['xmin'], + 'z0' : input_dict['z0'], + 'n_dirs' : input_dict['n_dirs'] + } + else: self.Type = 'Custom' # custom allocation of neurons according to entries 'X_coord_old', ... and 'YZ_angles', ... in GUI_inp_dict.py self.custom_structure = {'Center_coordinates' : [input_dict['X_coord_old'], input_dict['Y_coord_old'], input_dict['Z_coord_old']] , # list of lists [[x1,x2,x2],[y1,y2,y3],...], usually at the electrode's tip or the last contact @@ -350,23 +363,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]) + 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'] - - 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 @@ -577,7 +608,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") @@ -597,6 +628,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]: @@ -618,8 +650,9 @@ 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))) + + #pdb.set_trace() + check3=False#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),)': points_csf=points_csf+1 From ab354ef447dadb188fdb665d1c01c9f0921fbffd Mon Sep 17 00:00:00 2001 From: arashgmn Date: Sun, 11 Jul 2021 12:19:10 +0200 Subject: [PATCH 13/14] fixed Neural array processing check3 flag --- OSS_platform/Neural_array_processing.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/OSS_platform/Neural_array_processing.py b/OSS_platform/Neural_array_processing.py index c9d50aa..3c63667 100644 --- a/OSS_platform/Neural_array_processing.py +++ b/OSS_platform/Neural_array_processing.py @@ -652,8 +652,7 @@ def adjust_neuron_models(self,Domains,MRI_param): check2_2=(voxel_array_CSF_shifted[:,1]>=Array_coord[inx,1]) check2_3=(voxel_array_CSF_shifted[:,2]>=Array_coord[inx,2]) - #pdb.set_trace() - check3=False#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))) + 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),)': points_csf=points_csf+1 From 6d1b3fec0e3c6a26ad9245c75434c76480d1e2cc Mon Sep 17 00:00:00 2001 From: arashgmn Date: Sun, 11 Jul 2021 17:53:34 +0200 Subject: [PATCH 14/14] fix figsize upon saving --- OSS_platform/VTA/visualize.py | 30 ++++++++++++------------------ 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/OSS_platform/VTA/visualize.py b/OSS_platform/VTA/visualize.py index cf4ed01..dc81e97 100644 --- a/OSS_platform/VTA/visualize.py +++ b/OSS_platform/VTA/visualize.py @@ -18,20 +18,18 @@ def plot_vta_figs(path_target = path+'/Images/', path_pattern = path+'/Neuron_model_arrays/default_pattern.csv', path_activation = path+'/Field_solutions/Activation/Network_status.npy', activate_close_axons=False): - - - args = {'path_target': path_target, 'path_nodes':path_nodes, 'path_pattern': path_pattern, 'path_activation': path_activation, 'activate_close_axons': activate_close_axons} + 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(**args) + plot_cylindircal_vta(**kwargs) else: - plot_grid_vta(**args) + plot_grid_vta(**kwargs) + -def plot_cylindircal_vta(input_dict = '../GUI_inp_dict.py', - path_activation = '../Field_solutions/Activation/Network_status.npy', - path_target = '../Images/'): +def plot_cylindircal_vta(path_activation, path_target, activate_close_axons, **kwargs): """ Visualizes the VTA cross-section for each surface. """ @@ -68,8 +66,7 @@ def plot_cylindircal_vta(input_dict = '../GUI_inp_dict.py', vta_plane = activations[n*axon_per_plane: (n+1)*axon_per_plane] vta_plane = vta_plane.reshape(len(xc), len(zc)) - plt.figure(figsize=(15,4)) - + 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) @@ -78,16 +75,12 @@ def plot_cylindircal_vta(input_dict = '../GUI_inp_dict.py', 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) + plt.savefig(path_target+'Cyl_VTA_dir'+str(n)+'.png', dpi=300, bbox_inches='tight') plt.close() -def plot_grid_vta(path_target = '../Images/', - path_input_dict = '../GUI_input_dict.py', - path_nodes = '../Neuron_model_arrays/All_neuron_models.csv', - path_pattern = '../Neuron_model_arrays/default_pattern.csv', - path_activation = '../Field_solutions/Activation/Network_status.npy', - activate_close_axons=False): +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 @@ -254,7 +247,8 @@ def plot_grid_vta(path_target = '../Images/', plt.colorbar(pcm, cax=cax) plt.tight_layout() - plt.savefig(path_target+'Grid_VTA_dir'+str(dir_id)+'.png', dpi=300) + plt.savefig(path_target+'Grid_VTA_dir'+str(dir_id)+'.png', dpi=300, + bbox_inches='tight') plt.close()