From 3a99c0bf466e282c03f48421d12e61f310adb5c2 Mon Sep 17 00:00:00 2001 From: arashgmn Date: Sat, 15 May 2021 16:09:58 +0200 Subject: [PATCH 1/9] 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 2/9] 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 3/9] 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 4/9] 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 5/9] 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 6/9] 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 7/9] 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 8/9] 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 9/9] 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=" ")