Add files via upload

001_Preparing/RUN_THIS_FILE.m

@@ -0,0 +1,57 @@
+clear
+
+%{
+dir_root = {'Z:\Yongsoo_Kim_Lab\STP_processed\2019_optical\20191212_UC_U504_C57J_FITC-fill_M_p67_optical';
+            'Z:\Yongsoo_Kim_Lab_2\STP_processed\2020_optical\20200220_UC_U547_C57J_FITC-fill_M_p559_optical';
+            'Z:\Yongsoo_Kim_Lab_2\STP_processed\2020_optical\20200227_UC_U548_C57J_FITC-fill_M_p559_optical';
+            'Z:\Yongsoo_Kim_Lab_2\STP_processed\2020_optical\20200319_HB_U549_C57J_FITC-fill_M_p63_optical';
+            'Z:\Yongsoo_Kim_Lab_2\STP_processed\2020_optical\20200324_YK_U533_C57J_FITC-fill_F_p56_optical';
+            'Z:\Yongsoo_Kim_Lab_2\STP_processed\2020_optical\20200407_UC_U581_C57J_FITC-fill_F_p557_optical';
+            'Z:\Yongsoo_Kim_Lab_3\STP_processed\2020_optical\20200402_YK_U554_C57J_FITC-fill_F_p56_optical'};
+%}       
+%{
+dir_root = {'Z:\Yongsoo_Kim_Lab_3\STP_processed\2020_optical\20200328_YK_U573_QZ719_C57J_FITC-fill_F_24mo_optical';
+            'Z:\Yongsoo_Kim_Lab_3\STP_processed\2020_optical\20200412_YK_U550_C57J_FITC-fill_M_p63_optical';
+            'Z:\Yongsoo_Kim_Lab_3\STP_processed\2020_optical\20200424_UC_U582_C57J_FITC-fill_F_p557_optical';
+            'Z:\Yongsoo_Kim_Lab_3\STP_processed\2020_optical\20200501_UC_U598_C57J_FITC-fill_F_p69_optical'};
+
+dir_root = {
+'Z:\Yongsoo_Kim_Lab_3\STP_processed\2020_optical\20200515_UC_U583_C57J_FITC-fill_F_p557_optical'
+};
+%}        
+%{
+dir_root = {
+'Z:\Yongsoo_Kim_Lab_2\STP_processed\2020_optical\20200123_HB_HB106_PA5xFAD_WT_FITC-rat_optical';
+'Z:\Yongsoo_Kim_Lab_2\STP_processed\2020_optical\20200109_HB_HB107_PA5xFAD_MUT_FITC-rat_optical';
+'Z:\Yongsoo_Kim_Lab_2\STP_processed\2020_optical\20200312_HB_HB125_PA5xFAD_MUT_P72_FITC-fill_optical';
+'Z:\Yongsoo_Kim_Lab_3\STP_processed\2020_optical\20200522_HB_HB130_5xFAD_FITC-filll_F_MUT_p236_optical';
+'Z:\Yongsoo_Kim_Lab_3\STP_processed\2020_optical\20200529_HB_HB131_5xFAD_FITC-fill_F_WT_p236_optical';
+};
+%}
+
+%{
+% Batch #2
+dir_root = {
+'Y:/Yongsoo_Kim_Lab_3/STP_processed/2020_optical/20200725_HB_U584_C57J_FITC-fill_F_p558_optical';
+'Y:/Yongsoo_Kim_Lab_3/STP_processed/2020_optical/20200729_HB_U601_C57J_FITC-fill_F_p56_optical';
+'Y:/Yongsoo_Kim_Lab_3/STP_processed/2020_optical/20200803_HB_U602_C57J_FITC-fill_F_p56_optical';
+'Y:/Yongsoo_Kim_Lab_3/STP_processed/2020_optical/20200814_YK_U585_C57J_FITC-fill_F_18mo_optical';
+'Y:/Yongsoo_Kim_Lab_3/STP_processed/2020_optical/20200827_HB_U605_C57J_FITC-fill_M_p56_optical';
+%'Y:/Yongsoo_Kim_Lab_3/STP_processed/2020_optical/20201008_HB_HB292_C57J_FITC-fill_M_p64_optical';
+'Y:/Yongsoo_Kim_Lab_2/STP_processed/2020_optical/20200417_HB_HB123_5xFAD_WT_p120_FITC-rat_optical';
+'Y:/Yongsoo_Kim_Lab_3/STP_processed/2020_optical/20200522_HB_HB130_5xFAD_FITC-filll_F_MUT_p236_optical';
+'Y:/Yongsoo_Kim_Lab_3/STP_processed/2020_optical/20200529_HB_HB131_5xFAD_FITC-fill_F_WT_p236_optical';
+};
+%}
+
+
+% Batch #2
+dir_root = {
+'Y:\Yongsoo_Kim_Lab_3\STP_processed\2020_optical\20201013_HB_HB294_C57J-NIA_FITC-fill_F_18mo_optical';
+};
+
+
+
+for ii = 1:length(dir_root)
+    mat_clean_skel(dir_root{ii});
+end

001_Preparing/mat_clean_skel.m

@@ -0,0 +1,41 @@
+function mat_clean_skel(dir_root)
+
+
+data_name = regexp(dir_root, '[\\/]', 'split');
+
+while isempty(data_name{end}) & ~isempty(data_name)
+    data_name(end) = [];
+end
+if    isempty(data_name)
+    error('bad dir name')
+end
+
+
+dir_tif = [dir_root '/binarized'];
+dir_clean_skel = [dir_root '/clean_skel'];
+
+DirTif = dir([dir_tif '/*.tif']);
+
+FileTif=[DirTif(1).folder '/' DirTif(1).name];
+InfoImage=imfinfo(FileTif);
+mImage=InfoImage(1).Width;
+nImage=InfoImage(1).Height;
+lImage = length(InfoImage);
+numberFiles = length(DirTif);
+
+CropRange = [1, nImage; 1, mImage; 1, lImage.*numberFiles];
+CropSize = CropRange(:,2);
+
+
+
+fileIndicator = [dir_clean_skel '/clean_skel.bin' ];
+fileID = fopen(fileIndicator);
+S_skel = fread(fileID,'uint64');
+fclose(fileID);
+
+fileIndicator = [dir_clean_skel '/clean_radii.bin' ];
+fileID = fopen(fileIndicator);
+S_radii = fread(fileID,'double');
+fclose(fileID);
+
+save([ data_name{end}, '.mat'],'S_skel','S_radii', 'CropSize', '-v7.3');

002_Voxel_CFD/RUN_THIS_FILE_slurm_batch

@@ -0,0 +1,27 @@
+#!/bin/bash -l
+#SBATCH --account=lab_kim
+#SBATCH --partition=compute
+#SBATCH --job-name=YTW_CFD
+#SBATCH --output=GPU_XXX.%j.%N.out
+#SBATCH --error=GPU_XXX.%j.%N.err
+#SBATCH --nodes=1
+#SBATCH --ntasks=1
+#SBATCH --cpus-per-task=44
+#SBATCH --array=0-219%35 # job array index
+
+
+#SBATCH --time=100000:00:00
+#SBATCH --mem=300000
+
+module load matlab
+
+names=($(cat job_file_list))
+brain_name=`expr ${SLURM_ARRAY_TASK_ID} / 20`
+brain_part=`expr ${SLURM_ARRAY_TASK_ID} % 20`
+brain_part=`expr ${brain_part} + 1`
+
+echo ${names[${brain_name}]}
+echo ${brain_part}
+
+rm ~/.matlab/* -rf
+matlab -nodisplay -batch "mat_name = ${names[${brain_name}]}; HPC_poking_point = ${brain_part} ; cubical_CFD_micro_fff"

002_Voxel_CFD/cubical_CFD_micro_fff.m

@@ -0,0 +1,90 @@
+%function cubical_CFD_micro_function(mat_name)
+%clear 
+%mat_name = '20191212_UC_U504_C57J_FITC-fill_M_p67_optical.mat';
+
+
+data_name = regexp(mat_name, '[\\/.]', 'split');
+coreCount =22;
+aa = parcluster
+temp_folder = ['./matlab_cluster_', datestr(now,'yyyy-mm-dd-HH-MM-SS-FFF')];
+mkdir(temp_folder);
+aa.JobStorageLocation = temp_folder;
+parpool(aa,coreCount ,'IdleTimeout',inf)
+
+%[S_skel, S_radii, CropSize] = read_clean_skel(dir_root);
+load(mat_name);
+poking_size = [500, 500, 500];
+poking_distance = 100;
+
+[xxx,yyy,zzz] = ind2sub(CropSize,S_skel);
+
+
+[X,Y,Z] = meshgrid(1:poking_distance:CropSize(1),1:poking_distance:CropSize(2),1:poking_distance:CropSize(3));
+
+[Xi,Yi,Zi] = meshgrid(1:1:size(X,1),1:1:size(X,2),1:1:size(X,3));
+
+total_poking = length(X(:));
+%poking_point_list = zeros(total_poking,3);
+poking_point_list = [X(:),Y(:),Z(:)];
+easy_list = [Xi(:),Yi(:),Zi(:)];
+
+
+the_tensor = zeros(3,3,total_poking);
+the_tensor = num2cell(the_tensor,[1 2]);
+
+the_vector = zeros(3,1,total_poking);
+the_vector = num2cell(the_vector,[1 2]);
+
+
+
+poking_group_size = 100;
+poking_HPC_devider = 20;
+total_HPC_poking_group = ceil(total_poking./poking_group_size);
+HPC_poking_group_size = ceil(total_HPC_poking_group./poking_HPC_devider);
+
+%HPC_poking_point = 1~20
+
+starting_HPC = (HPC_poking_point-1).* HPC_poking_group_size ; 
+if HPC_poking_point == 20
+    ending_HPC = total_HPC_poking_group-1;
+else
+    ending_HPC = (HPC_poking_point).* HPC_poking_group_size -1; 
+end
+
+
+
+for jj = starting_HPC:ending_HPC
+    tic
+    starting_shift = jj.*poking_group_size ;
+    ending = (jj+1).*poking_group_size;
+    if jj == ceil(total_poking./poking_group_size)-1
+        ending = total_poking;
+    end
+
+    for ii = 1:(ending - starting_shift)
+
+        poking_point = poking_point_list(ii+starting_shift,:);
+
+        flag = (floor((xxx-poking_point(1))./poking_size(1)) == 0) & ...
+            (floor((yyy-poking_point(2))./poking_size(2)) == 0) & ...
+            (floor((zzz-poking_point(3))./poking_size(3)) == 0);
+
+        P_sk{ii} = S_skel(flag);
+        P_radii{ii} = S_radii(flag);
+        xyz_lower_upper{ii} = [poking_point(1), poking_point(1)+poking_size(1)-1;
+            poking_point(2), poking_point(2)+poking_size(2)-1;
+            poking_point(3), poking_point(3)+poking_size(3)-1]';
+
+
+    end
+
+    parfor ii = 1:(ending - starting_shift )
+
+        [the_tensor{ii+starting_shift}, the_vector{ii+starting_shift}] = making_tensor_micro(P_sk{ii}, P_radii{ii}, CropSize, xyz_lower_upper{ii});
+
+    end
+    toc
+end
+
+
+save( [data_name{end-1},'_',num2str(HPC_poking_point,'%03d'), '_out.mat'],'the_tensor','the_vector', 'poking_point_list', 'easy_list', '-v7.3');

002_Voxel_CFD/job_file_list

@@ -0,0 +1,11 @@
+'20191212_UC_U504_C57J_FITC-fill_M_p67_optical.mat'
+'20200220_UC_U547_C57J_FITC-fill_M_p559_optical.mat'
+'20200227_UC_U548_C57J_FITC-fill_M_p559_optical.mat'
+'20200319_HB_U549_C57J_FITC-fill_M_p63_optical.mat'
+'20200324_YK_U533_C57J_FITC-fill_F_p56_optical.mat'
+'20200328_YK_U573_QZ719_C57J_FITC-fill_F_24mo_optical.mat'
+'20200402_YK_U554_C57J_FITC-fill_F_p56_optical.mat'
+'20200407_UC_U581_C57J_FITC-fill_F_p557_optical.mat'
+'20200412_YK_U550_C57J_FITC-fill_M_p63_optical.mat'
+'20200424_UC_U582_C57J_FITC-fill_F_p557_optical.mat'
+'20200501_UC_U598_C57J_FITC-fill_F_p69_optical.mat'

002_Voxel_CFD/pre_statistics.m

@@ -0,0 +1,122 @@
+function [S_link, S_link_group, wierdNode, wierdNode_group] = pre_statistics(S_skel, CropSize)
+
+
+tic
+[~,nh] = get_nh_save_memory(S_skel,CropSize);
+toc
+distance_table = zeros(3,3,3);
+
+distance_table(:,:,1) = [1.732, 1.414, 1.732; 1.414, 1, 1.414; 1.732, 1.414, 1.732];
+distance_table(:,:,2) = [1.414, 1, 1.414; 1, 0, 1; 1.414, 1, 1.414];
+distance_table(:,:,3) = [1.732, 1.414, 1.732; 1.414, 1, 1.414; 1.732, 1.414, 1.732];
+
+distance_table = distance_table(:);
+
+S_link.length = nh.*distance_table';
+S_link.length = sum(S_link.length,2);
+sum_nh_full = sum(nh,2);
+
+
+S_link.name = find(sum_nh_full==3);
+S_link.length = S_link.length(sum_nh_full==3);
+
+tic
+[S_link_link_table] = get_all_link_par(S_skel(S_link.name),CropSize);
+
+S_link.group = grouping_remake(length(S_link.name), S_link_link_table);
+toc
+tic
+clear S_link_link_table
+
+S_link.group = S_link.name(S_link.group);
+S_link_group.name = unique(S_link.group);
+[~,loc] = ismember(S_link.group,S_link_group.name);
+
+S_link_group.length = accumarray(loc,S_link.length(S_link.name));
+
+badLinks = S_link_group.name(S_link_group.length<10);
+
+
+
+removal = ismember( S_link.group,badLinks);
+
+wierdNode.name = S_link.name(removal);
+wierdNode.name = [find(sum_nh_full>3 | sum_nh_full <= 2 ); wierdNode.name ];
+clear sum_nh_full
+
+S_link.name = S_link.name(~removal);
+S_link.group = S_link.group(~removal);
+S_link.length = S_link.length(~removal);
+
+removal = ismember( S_link_group.name,badLinks);
+
+S_link_group.name = S_link_group.name(~removal);
+S_link_group.length = S_link_group.length(~removal);
+
+wierdNode.name = sort(wierdNode.name);
+
+
+%fprintf('Group and map the connectivity of the problematic nodes \n');
+
+
+%%% group and map the connectivity of the problematic nodes
+
+
+[link_table_wierd] = get_all_link_par(S_skel(wierdNode.name),CropSize);
+wierdNode.group = grouping_remake(length(wierdNode.name), link_table_wierd);
+wierdNode.group = wierdNode.name(wierdNode.group);
+clear link_table_wierd
+
+
+[wierdNode_group.name] = unique(wierdNode.group);
+
+[link_table_full] = get_all_link_par(S_skel,CropSize);
+
+link_table_wierd_out = link_table_full(xor(ismember(link_table_full(:,1),wierdNode.name) , ismember(link_table_full(:,2),wierdNode.name)),:);
+link_table_wierd_full = link_table_full(ismember(link_table_full(:,1),wierdNode.name) | ismember(link_table_full(:,2),wierdNode.name),:);
+clear link_table_full
+
+
+
+lia1 = ismember(link_table_wierd_out(:,1),wierdNode.name);
+lia2 = ismember(link_table_wierd_out(:,2),wierdNode.name);
+
+
+node_wierd_out_in = zeros(length(link_table_wierd_out),1);
+node_wierd_out_out = zeros(length(link_table_wierd_out),1);
+node_wierd_out_out(lia1) = link_table_wierd_out(lia1,2);
+node_wierd_out_out(lia2) = link_table_wierd_out(lia2,1);
+node_wierd_out_in(lia1) = link_table_wierd_out(lia1,1);
+node_wierd_out_in(lia2) = link_table_wierd_out(lia2,2);
+
+[~, loc] = ismember(node_wierd_out_in,wierdNode.name);
+
+node_wierd_out_group = wierdNode.group(loc);
+
+clear lia1 lia2 loc link_table_wierd_out node_wierd_out_in
+
+
+
+temp_w = [node_wierd_out_out, node_wierd_out_group];
+temp_w = sortrows(temp_w,2);
+[~,loc] = ismember(temp_w(:,2), wierdNode_group.name);
+binsize = histcounts(loc,[0; find(wierdNode_group.name)]+0.5);
+
+wierdNode_out_name_par = mat2cell(temp_w(:,1),binsize);
+
+
+wierdNode_group.connection_count = zeros(size(wierdNode_group.name));
+
+for ii = 1:length(wierdNode_group.name)
+    wierdNode_group.connection_count(ii) = length(wierdNode_out_name_par{ii});
+end
+fprintf(['Isolated short segments count: ' num2str(sum(wierdNode_group.connection_count(:) == 0))  ' \n' ]);
+fprintf(['Single connected short segments count: ' num2str(sum(wierdNode_group.connection_count(:) == 1))  ' \n' ]);
+fprintf(['2x connected short segments count: ' num2str(sum(wierdNode_group.connection_count(:) == 2))  ' \n' ]);
+fprintf(['3x connected pile count: ' num2str(sum(wierdNode_group.connection_count(:) == 3))  ' \n' ]);
+fprintf(['4x connected pile count: ' num2str(sum(wierdNode_group.connection_count(:) == 4))  ' \n' ]);
+fprintf(['>=5x connected pile count: ' num2str(sum(wierdNode_group.connection_count(:) >=5 ))  ' \n' ]);
+
+
+
+toc

002_Voxel_CFD/private/Pre_cal_length.m

@@ -0,0 +1,28 @@
+function [sphere_mask, cubic_path,cubic_length] = Pre_cal_length(gap_length)
+
+
+cubic_size = 2.*gap_length+1;
+cubic_ind = ones(cubic_size,cubic_size,cubic_size) ;
+
+cubic_ind(:) = find(cubic_ind);
+[cubic_xxx, cubic_yyy, cubic_zzz] = ind2sub(size(cubic_ind),cubic_ind(:));
+
+cubic_length = zeros(cubic_size,cubic_size,cubic_size) ;
+cubic_length(:) = sqrt((cubic_xxx-(gap_length+1)).^2+(cubic_yyy-(gap_length+1)).^2+(cubic_zzz-(gap_length+1)).^2);
+
+sphere_mask = (cubic_length <=gap_length);
+
+cubic_CropSize = size(cubic_ind);
+goodLinkTable = getAllLinks_v3(cubic_ind(:),cubic_CropSize);
+
+distance_table = [3 2 3 2 1 2 3 2 3 2 1 2 1 0 1 2 1 2 3 2 3 2 1 2 3 2 3 ];
+distance_table = sqrt(distance_table)';
+goodLinkTable(:,3) = distance_table(goodLinkTable(:,3));
+
+cubic_graph = graph(goodLinkTable(:,1),goodLinkTable(:,2),goodLinkTable(:,3));
+center_node = ceil(length(cubic_ind(:))./2);
+
+parfor ii = 1:length(cubic_ind(:))
+    cubic_path{ii} = shortestpath(cubic_graph,center_node,ii);
+    cubic_path{ii} = cubic_path{ii}(cubic_path{ii}~= center_node & cubic_path{ii}~= ii);
+end

002_Voxel_CFD/private/Pre_cal_mask.m

@@ -0,0 +1,13 @@
+function [sphere_mask, cubic_length] = Pre_cal_mask(gap_length)
+
+
+cubic_size = 2.*gap_length+1;
+cubic_ind = ones(cubic_size,cubic_size,cubic_size) ;
+
+cubic_ind(:) = find(cubic_ind);
+[cubic_xxx, cubic_yyy, cubic_zzz] = ind2sub(size(cubic_ind),cubic_ind(:));
+
+cubic_length = zeros(cubic_size,cubic_size,cubic_size) ;
+cubic_length(:) = sqrt((cubic_xxx-(gap_length+1)).^2+(cubic_yyy-(gap_length+1)).^2+(cubic_zzz-(gap_length+1)).^2);
+
+sphere_mask = (cubic_length <=gap_length);