Add Sampsa's changes to NSE Tool and its dependecies in m/

- Apparently there was a mix-up with folder contents, so NSE Tool needed fixing.

m/barycentric/zef_surface_scalar_vector_Fn.m

@@ -38,7 +38,7 @@
 
 N = size(nodes,1);
 
-if nargin < 3
+if nargin < 4
     scalar_field = ones(size(tetra,1),1);
 end
 

m/forward_simulation/nse/zef_nse_poisson_dynamic.m

@@ -74,6 +74,11 @@
 area_arteries = sum(area(I));
 param_aux_integral_mean = param_aux_integral_mean/area_arteries;
 
+param_aux_ones = zeros(size(param_aux));
+param_aux_ones(t_ind) = 1;
+w_3 = zef_surface_scalar_vector_F(v_1_nodes, v_1_tetra, param_aux_ones);
+sum_w_3 = sum(w_3);
+
 gravity_amplitude = nse_field.gravity_amplitude;
 gravity_x = nse_field.gravity_x;
 gravity_y = nse_field.gravity_y;
@@ -110,6 +115,10 @@
 Q_2_v1 = spdiags(1./c_vec,0,size(v_1_nodes,1),size(v_1_nodes,1))*Q_2;
 Q_3_v1 = spdiags(1./c_vec,0,size(v_1_nodes,1),size(v_1_nodes,1))*Q_3;
 
+Q_1_v12 = Q_1_v1;
+Q_2_v12 = Q_2_v1;
+Q_3_v12 = Q_3_v1;
+
 Q_1_v2 = Q_1_v1(:,i_node_ind);
 Q_2_v2 = Q_2_v1(:,i_node_ind);
 Q_3_v2 = Q_3_v1(:,i_node_ind);
@@ -268,15 +277,15 @@
 end
 
 KDMD = @(x) zef_KDMD(x,K_1,M_1,D_1,nse_field.use_gpu);
-
+ 
 for i = 1 : n_time
 
     pressure_aux = (p + pressure_reference + p_hydrostatic)/hgmm_conversion;
-    total_flow_aux = sum(sqrt(u_1.^2 + u_2.^2 + u_3.^2).*w_1)/(ml_min_conversion);
+    total_flow_aux = (sum(((p + pressure_reference)/hgmm_conversion).*w_3)/sum_w_3)*nse_field.total_flow/nse_field.pressure;
     zef_waitbar(i,n_time,h_waitbar,['NSE solver: compute, total flow (ml/min): ' sprintf('%0.3g',total_flow_aux) ', 90 % pressure quantile (Hgmm): ' sprintf('%0.3g',quantile(pressure_aux,0.9)) ', 10 % pressure quantile (Hgmm): ' sprintf('%0.3g',quantile(pressure_aux,0.1)) '.']);
 
     y_0 = y(i);
-
+    
     for quadrature_step_ind = 1 : n_q_steps
 
         if nse_field.nse_type == 2
@@ -311,8 +320,10 @@
         if quadrature_step_ind == 1
             if i == 1
                 b = D_1 * (M_1 * (D_1 * (p_1  + source_vec * (y_0 - 2*y_1))));
+
             else
-                b = D_1 * (M_1 * (D_1 * (2*p_1 - p_2  + source_vec * (y_0 - 2*y_1 + y_2))));
+               b = D_1 * (M_1 * (D_1 * (2*p_1 - p_2  + source_vec * (y_0 - 2*y_1 + y_2))));
+
             end
 
             g_mu_1 = Q_1_v1*mu_vec(i_node_ind);
@@ -364,6 +375,8 @@
 
             p_2 = p_1;
             p_1 = p;
+            max(p_1);
+            min(p_1);
             y_2 = y_1;
             y_1 = y_0;
 
@@ -411,7 +424,7 @@
                 bp_vessels_aux(nse_field.bp_vessel_node_ind) =  abs(p + pressure_reference);
                 bf_vessels_to_capillaries = bp_vessels_aux(nse_field.bf_capillary_node_ind);
                 I_boundary = find(bf_vessels_to_capillaries);
-                total_flow_estimate= sum(sqrt(u_1.^2 + u_2.^2 + u_3.^2).*w_1);
+                total_flow_estimate = (sum(((p + pressure_reference)/hgmm_conversion).*w_3)/sum_w_3)*nse_field.total_flow/nse_field.pressure;
                 r = total_flow_estimate*bf_vessels_to_capillaries./sum(pressure_reference.*w_2(I_boundary));
 
                 if nse_field.use_gpu
@@ -487,36 +500,90 @@
         end
     end
 
+
+        u_1_2 = zeros(size(K_1,1),1);
+        u_1_2(i_node_ind,1) = u_1;
+        u_2_2 = zeros(size(K_1,1),1);
+        u_2_2(i_node_ind,1) = u_2;
+        u_3_2 = zeros(size(K_1,1),1);
+        u_3_2(i_node_ind,1) = u_3;
+%         collect1{i}=gather(u_1_2);
+%         collect2{i}=gather(u_2_2);
+%         collect3{i}=gather(u_3_2);
+       % trace_strain_rate2 = rateofshear(adj1,dist1,u_1_2,u_2_2,u_3_2);
+%        trace_strain_rate = Q_1_v12*u_1_2;
+        trace_strain_rate = 4.*(Q_1_v12*u_1_2).^2 + 4.*(Q_2_v12*u_2_2).^2 + 4.*(Q_3_v12*u_3_2).^2 + 2.*(Q_2_v12*u_1_2+Q_1_v12*u_2_2).^2 +...
+            2.*(Q_3_v12*u_1_2+Q_1_v12*u_3_2).^2+2.*(Q_3_v12*u_2_2+Q_2_v12*u_3_2).^2;
+        trace_strain_rate = sqrt((1/2).*trace_strain_rate);
+% if i >= n_time*0.5
     if nse_field.viscosity_model == 2
 
-        trace_strain_rate = abs(Q_1_v2*u_1) + abs(Q_2_v2*u_2) + abs(Q_3_v2*u_3);
+%         trace_strain_rate = abs(Q_1_v2*u_1) + abs(Q_2_v2*u_2) + abs(Q_3_v2*u_3);
         mu_vec = nse_field.mu*trace_strain_rate.^(nse_field.viscosity_exponent-1);
 
     elseif nse_field.viscosity_model == 3
 
-        trace_strain_rate = abs(Q_1_v2*u_1) + abs(Q_2_v2*u_2) + abs(Q_3_v2*u_3);
+%         trace_strain_rate = abs(Q_1_v2*u_1) + abs(Q_2_v2*u_2) + abs(Q_3_v2*u_3);
         mu_vec = nse_field.mu + nse_field.viscosity_delta.*(1 + (nse_field.viscosity_relaxation_time*trace_strain_rate).^nse_field.viscosity_transition).^((nse_field.viscosity_exponent-1)./nse_field.viscosity_transition);
 
+
+    elseif nse_field.viscosity_model == 4
+
+         C1 = 0.00797;
+         C2 = 0.0608;
+         C3 = 0.00499;
+         C4 = 14.585;
+         Hem = 40;
+         TPMA = 25.9;
+         mu_vec = 0.1*C1*exp(C2*Hem)*exp(C4*(TPMA/(Hem^2))).*trace_strain_rate.^-(C3*Hem);
+         mu_vec(mu_vec>0.056) = 0.056;
+         mu_vec(mu_vec<0.00345) = 0.00345;
+      elseif nse_field.viscosity_model == 5
+
+         lambda_g = zeros(size(K_1,1),1);
+         n_g = zeros(size(K_1,1),1);
+%          trace_strain_rate2 = abs(Q_1_v2*u_1) + abs(Q_2_v2*u_2) + abs(Q_3_v2*u_3);
+%          trace_strain_rate = rateofshear(kk,u_1_2,u_2_2,u_3_2);
+         lambda_g(trace_strain_rate>0.01) = 0.0345 + 0.25.*exp(-(1+trace_strain_rate(trace_strain_rate>0.01)./50).*exp(-3.*(trace_strain_rate(trace_strain_rate>0.01)).^-1));
+         n_g(trace_strain_rate>0.01) = 1 - 0.45.*exp(-(1+trace_strain_rate(trace_strain_rate>0.01)./50).*exp(-4.*(trace_strain_rate(trace_strain_rate>0.01)).^-1));
+         mu_vec(trace_strain_rate>0.01) = 0.1*lambda_g(trace_strain_rate>0.01).*trace_strain_rate(trace_strain_rate>0.01).^(n_g(trace_strain_rate>0.01)-1);
+         mu_vec(trace_strain_rate<=0.01)=0.056;
+         mu_vec(mu_vec>0.056) = 0.056;
+         mu_vec(mu_vec<0.00345)=0.00345;
+
+    elseif nse_field.viscosity_model == 6
+        mu_vec(trace_strain_rate>0.4e-3) = (sqrt(0.0012*((1-0.35)^(-5/2))) + sqrt(0.01*((0.625*0.35)^3)./trace_strain_rate(trace_strain_rate>0.4e-3))).^2;
+        mu_vec(mu_vec>0.056) = 0.056;
+        mu_vec(mu_vec<0.00345)=0.00345;
     end
+% end
 
-    if ismember(nse_field.viscosity_model,[2 3])
+
+    if ismember(nse_field.viscosity_model,[2 3 4 5 6])%2 3 4 5 6
 
         if nse_field.viscosity_smoothing > 0
             if nse_field.use_gpu
                 mu_vec = pcg_iteration_gpu(S_mu,I_mu*mu_vec,nse_field.pcg_tol,nse_field.pcg_maxit,DM_S_mu);
+                strain_rate2 = pcg_iteration_gpu(S_mu,I_mu*trace_strain_rate,nse_field.pcg_tol,nse_field.pcg_maxit,DM_S_mu);
             else
                 mu_vec = pcg_iteration(S_mu,I_mu*mu_vec,nse_field.pcg_tol,nse_field.pcg_maxit,DM_S_mu);
+                strain_rate2 = pcg_iteration(S_mu,I_mu*trace_strain_rate,nse_field.pcg_tol,nse_field.pcg_maxit,DM_S_mu);
             end
+        else
+            strain_rate2 = trace_strain_rate;
         end
         D_1 = spdiags(sqrt(mu_vec),0,size(M_1,1),size(M_1,2));
         KDMD = @(x) zef_KDMD(x,K_1,M_1,D_1,nse_field.use_gpu);
+    else
+        strain_rate2 = trace_strain_rate;
 
     end
 
     if ismember(i,time_frame_ind)
 
         i_aux = i_aux + 1;
 
+        nse_field.shearrate1{i_aux} = zeros(size(K_1,1),1);
         nse_field.bp_vessels{i_aux} = zeros(size(K_1,1),1);
         nse_field.mu_vessels{i_aux} = zeros(size(K_1,1),1);
         nse_field.bv_vessels_1{i_aux} = zeros(size(K_1,1),1);
@@ -526,6 +593,7 @@
         nse_field.bv_vessels_2{i_aux}(i_node_ind) = gather(u_2);
         nse_field.bv_vessels_3{i_aux}(i_node_ind) = gather(u_3);
         nse_field.bp_vessels{i_aux} = gather(p)/hgmm_conversion + gather(p_hydrostatic)/hgmm_conversion + gather(pressure_reference)/hgmm_conversion;
+        nse_field.shearrate1{i_aux} = gather(strain_rate2);
         nse_field.mu_vessels{i_aux} = gather(mu_vec);
         if nse_field.microcirculation_model
             nse_field.bf_capillaries{i_aux} = gather(mc_vec);

plugins/NSE_tool/m/zef_nse_calculate_perfusion.m

@@ -0,0 +1,56 @@
+function perfusion_estimate = zef_nse_calculate_perfusion(nse_field,nodes,tetra,domain_labels,mvd_length)
+
+mm_conversion = 0.001;
+ml_min_conversion = 1e-6/60;
+hgmm_conversion = 101325/760;
+
+mvd_length = 1E6.*mvd_length(:,1);
+
+c_ind_1_domain = find(ismember(domain_labels,nse_field.artery_domain_ind));
+[v_1_nodes, v_1_tetra, nse_field.bp_vessel_node_ind] = zef_get_submesh(nodes, tetra, c_ind_1_domain);
+v_1_nodes = mm_conversion*v_1_nodes;
+[~, det] = zef_volume_barycentric(v_1_nodes,v_1_tetra);
+volume = abs(det)/6;
+
+[~,~,t_ind,~,~,~,~,f_ind] = zef_surface_mesh(v_1_tetra);
+[t_ind_2, f_ind_2] = zef_find_adjacent_tetra(tetra, c_ind_1_domain(t_ind), f_ind);
+[b_vec,det] = zef_volume_barycentric(v_1_nodes,v_1_tetra(t_ind,:),f_ind);
+area = (abs(det)/2).*sqrt(sum(b_vec(:,1:3).^2,2));
+
+param_aux = zeros(length(c_ind_1_domain),1);
+param_aux(t_ind) = mvd_length(t_ind_2);
+param_aux_integral_mean = zef_surface_scalar_vector_F(v_1_nodes, v_1_tetra, param_aux);
+param_aux_integral_mean = sum(param_aux_integral_mean);
+I = find(param_aux(t_ind));
+area_arteries = sum(area(I));
+
+param_aux_ones = zeros(size(param_aux));
+param_aux_ones(t_ind) = 1;
+w_3 = zef_surface_scalar_vector_F(v_1_nodes, v_1_tetra,param_aux_ones);
+perfusion_estimate = zeros(size(nse_field.bv_vessels_1,2),1);
+
+gravity_amplitude = nse_field.gravity_amplitude;
+gravity_x = nse_field.gravity_x;
+gravity_y = nse_field.gravity_y;
+gravity_z = nse_field.gravity_z;
+
+gravity_vec = gravity_x.*v_1_nodes(:,1) + gravity_y.*v_1_nodes(:,2) + gravity_z.*v_1_nodes(:,3);
+gravity_vec = gravity_vec - min(gravity_vec);
+p_hydrostatic = nse_field.rho*gravity_vec;
+
+%n_1 = zef_surface_scalar_vector_Fn(v_1_nodes, v_1_tetra, 1);
+%n_2 = zef_surface_scalar_vector_Fn(v_1_nodes, v_1_tetra, 2);
+%n_3 = zef_surface_scalar_vector_Fn(v_1_nodes, v_1_tetra, 3);
+
+for i = 1 : size(nse_field.bv_vessels_1,2)
+
+perfusion_estimate(i) = (sum((nse_field.bp_vessels{i}-p_hydrostatic/hgmm_conversion).*w_3)/sum(w_3))*nse_field.total_flow/nse_field.pressure;
+
+end
+
+end
+
+
+
+
+

plugins/NSE_tool/m/zef_nse_dir_v_node.m

@@ -0,0 +1,16 @@
+h_axes = gca;
+if isempty(findobj(allchild(h_axes),'Type','DataTip'))~=1
+    h_datatip = findobj(allchild(h_axes),'Type','DataTip');
+    zef.nse_field.dir_v_x = [];
+    zef.nse_field.dir_v_y = [];
+    zef.nse_field.dir_v_z = [];
+    for i=1:size(h_datatip,1)
+    zef.nse_field.dir_v_x = [zef.nse_field.dir_v_x h_datatip(i).X];
+    zef.nse_field.dir_v_y = [zef.nse_field.dir_v_y h_datatip(i).Y];
+    zef.nse_field.dir_v_z = [zef.nse_field.dir_v_z h_datatip(i).Z];
+    end
+    zef.nse_field.h_dir_v_x.Value = num2str(zef.nse_field.dir_v_x);
+    zef.nse_field.h_dir_v_y.Value = num2str(zef.nse_field.dir_v_y);
+    zef.nse_field.h_dir_v_z.Value = num2str(zef.nse_field.dir_v_z);
+
+end

plugins/NSE_tool/m/zef_nse_mean_velocity_roi.m

@@ -10,7 +10,7 @@
 v_2 = v_2 + mean(zef.nse_field.bv_vessels_2{i}(roi_ind));
 v_3 = v_3 + mean(zef.nse_field.bv_vessels_3{i}(roi_ind));
 end
-v_vec = [v_1 v_2 v_3]./length(nse_field.bv_vessels_1);
+v_vec = [v_1 v_2 v_3]./length(zef.nse_field.bv_vessels_1);
 v_mag = sqrt(sum(v_vec.^2));
 v_dir = v_vec/v_mag;
 

plugins/NSE_tool/m/zef_nse_plot_epoched.m

@@ -1,34 +1,80 @@
 
-function  zef_nse_plot_epoched(zef, nse_field, plot_vec, y_label)
+function  zef_nse_plot_epoched(zef, nse_field, plot_vec, y_label, plot_vec_ref)
 
 h_axes = zef.h_axes1;
 axes(h_axes);
 
+perform_shift = true;
+
+if nargin < 5
+plot_vec_ref = [];
+elseif isempty(plot_vec_ref)
+    perform_shift = false;
+end
+
 plot_colors = {'c', 'c'};
-plot_quantiles = [ 0.25 0.5 0.75];
+plot_quantiles = [ 0.0 0.5 1]; 
 plot_linestyle = {'--','-','--'};
 n_time = length(plot_vec);
 time_step_length_aux = (nse_field.time_length - nse_field.start_time)/(n_time-1);
 time_vec_aux = [nse_field.start_time:time_step_length_aux:nse_field.time_length];
 n_epoch =  find(time_vec_aux-nse_field.start_time >= nse_field.cycle_length,1,'first')-1;
 
 epoch_data = zeros(n_epoch,ceil(n_time/n_epoch));
+if not(isempty(plot_vec_ref))
+    epoch_data_2 = epoch_data;
+end
 
 for i = 1 : ceil(n_time/n_epoch)
 start_ind = n_epoch*(i-1)+1;
 end_ind =  min(n_time,n_epoch*i);
 epoch_data(1:end_ind-start_ind+1,i) = plot_vec(start_ind:end_ind); 
+if not(isempty(plot_vec_ref))
+epoch_data_2(1:end_ind-start_ind+1,i) = plot_vec_ref(start_ind:end_ind); 
+end
 end
 
 plot_data = zeros(n_epoch,length(plot_quantiles));
+if not(isempty(plot_vec_ref))
+plot_data_2 = plot_data;
+end
+
+if end_ind-start_ind+1 < size(epoch_data,1)
+    epoch_data(end_ind-start_ind+2:end,end) = mean(epoch_data(end_ind-start_ind+2:end,1:end-1),2);
+if not(isempty(plot_vec_ref))
+epoch_data_2(end_ind-start_ind+2:end,end) = mean(epoch_data_2(end_ind-start_ind+2:end,1:end-1),2);
+end
+end
 
 for i = 1 : length(plot_quantiles)
+% if plot_quantiles(i)==0.5
+% plot_data(:,i) = 0.5*(quantile(epoch_data,0,2)+quantile(epoch_data,1,2));
+% else
+aux_plot_vec = repmat(quantile(epoch_data,plot_quantiles(i),2),3,1);
+aux_plot_vec = reshape(movmean(aux_plot_vec, ceil(0.05*size(epoch_data,1))),size(epoch_data,1),3);
+aux_plot_vec = aux_plot_vec(:,2);
+plot_data(:,i) = aux_plot_vec;
+if not(isempty(plot_vec_ref))
+aux_plot_vec = repmat(quantile(epoch_data_2,plot_quantiles(i),2),3,1);
+aux_plot_vec = reshape(movmean(aux_plot_vec, ceil(0.05*size(epoch_data,1))),size(epoch_data,1),3);
+aux_plot_vec = aux_plot_vec(:,2);
+plot_data_2(:,i) = aux_plot_vec;
+end
+
+%end
+end
 
-plot_data(:,i) = quantile(epoch_data,plot_quantiles(i),2);
 
+if not(isempty(plot_vec_ref))
+    [~, index_shift] = max(plot_data_2(:,length(plot_quantiles)-1)/2);
+else
+    if perform_shift 
+[~, index_shift] = max(plot_data(:,length(plot_quantiles)-1)/2);
+    else
+        index_shift = 0;
+    end
 end
 
-[~,index_shift] = max(plot_data(:,length(plot_quantiles)-1)/2);
 index_shift = index_shift + floor(size(plot_data,1)/2);
 p_vec_aux = round(mod(index_shift+1:n_epoch+index_shift,n_epoch+0.1));