From f41ad0cd1e8d53ce13d4390b89968a6b20b03f19 Mon Sep 17 00:00:00 2001 From: Ali Aghaeifar Date: Wed, 20 Mar 2024 11:14:04 +0100 Subject: [PATCH] tissue specific plot --- MATLAB/plot_results.m | 108 +++++++++++++++++++++++++---------------- MATLAB/read_spinwalk.m | 14 +++--- 2 files changed, 72 insertions(+), 50 deletions(-) diff --git a/MATLAB/plot_results.m b/MATLAB/plot_results.m index 6a8c987..576b1f9 100644 --- a/MATLAB/plot_results.m +++ b/MATLAB/plot_results.m @@ -1,53 +1,66 @@ clc clear -% fname = cell(3,1); -fname{1}{1} = '../../outputs/gre_m1_0.dat'; -fname{1}{2} = '../../outputs/gre_m1_1.dat'; - -fname{2}{1} = '../../outputs/se_m1_0.dat'; -fname{2}{2} = '../../outputs/se_m1_1.dat'; - -fname{3}{1} = '../../outputs/ssfp_m1_0.dat'; -fname{3}{2} = '../../outputs/ssfp_m1_1.dat'; - -% fname{4}{1} = '../../outputs/grase_m1_0.dat'; -% fname{4}{2} = '../../outputs/grase_m1_1.dat'; - +file_fieldmap = '../../field_maps/fieldmap_0.dat'; + +file_m1{1}{1} = '../../outputs/gre_m1_fieldmap_0.dat'; +file_m1{1}{2} = '../../outputs/gre_m1_fieldmap_1.dat'; +file_m1{2}{1} = '../../outputs/se_m1_fieldmap_0.dat'; +file_m1{2}{2} = '../../outputs/se_m1_fieldmap_1.dat'; +file_m1{3}{1} = '../../outputs/ssfp_m1_fieldmap_0.dat'; +file_m1{3}{2} = '../../outputs/ssfp_m1_fieldmap_1.dat'; + +file_xyz1{1}{1} = '../../outputs/gre_xyz1_fieldmap_0.dat'; +file_xyz1{1}{2} = '../../outputs/gre_xyz1_fieldmap_1.dat'; +file_xyz1{2}{1} = '../../outputs/se_xyz1_fieldmap_0.dat'; +file_xyz1{2}{2} = '../../outputs/se_xyz1_fieldmap_1.dat'; +file_xyz1{3}{1} = '../../outputs/ssfp_xyz1_fieldmap_0.dat'; +file_xyz1{3}{2} = '../../outputs/ssfp_xyz1_fieldmap_1.dat'; + +dim_xyz = 1; dim_echo = 2; dim_spin = 3; dim_vessel_size = 4; rad_ref_um = 53.367; -TR = 0.2; -EcoSpc = 0.005; +tissue_type = 0; % 0 = extra-vascular, 1 = intra-vascular, [0,1] = combined +[~, mask, fov] = read_fieldmap(file_fieldmap); -signal_magnitude = cell(numel(fname), 1); -for seq = 1:numel(fname) - spins_xy = []; - for i=1:numel(fname{seq}) - [m_xyz, dims, scales] = read_spinwalk(fname{seq}{i}); +signal_magnitude = cell(numel(file_m1), numel(file_m1{1})); +for seq = 1:numel(file_m1) + for i=1:numel(file_m1{seq}) + [m1, dims, scales] = read_spinwalk(file_m1{seq}{i}); + [xyz1, ~, ~] = read_spinwalk(file_xyz1{seq}{i}); if dims(4) ~= numel(scales) warning('Header info is confusing here?') end - spins_xy = cat(ndims(m_xyz)+1, spins_xy, m_xyz); + + m1_t = zeros(numel(scales), 1); + for s=1:numel(scales) + m1_f = filter_spinwalk(m1(:,end,:,s), xyz1(:,end,:,s), tissue_type, mask, fov * scales(s)); + m1_t(s) = abs(complex(sum(m1_f(1,:)), sum(m1_f(2,:)))); + end + signal_magnitude{seq, i} = m1_t; end - - spins_xy = complex(sum(spins_xy(1,:,:,:,:), dim_spin), sum(spins_xy(2,:,:,:,:), dim_spin) ); - signal_magnitude{seq} = abs(spins_xy); end clf vessel_radius = rad_ref_um * scales; -for seq = 1:numel(fname) -% if seq ~= numel(fname) -% subplot(2,1,1); - relative_signal = 100 * (1 - signal_magnitude{seq}(:,:,:,:,1)./ signal_magnitude{seq}(:,:,:,:,2)); - relative_signal = squeeze(relative_signal); - h = semilogx(vessel_radius, relative_signal); xlabel('Vessel radius (um)'); ylabel('BOLD Signal %'); - hold on; - ylim([0, 7]) -% else +for seq = 1:numel(file_m1) + relative_signal = 100 * (1 - signal_magnitude{seq, 1} ./ signal_magnitude{seq, 2}); + h = semilogx(vessel_radius, relative_signal); xlabel('Vessel radius (um)'); ylabel('BOLD Signal %'); + hold on; + ylim([0, 7]) +end +legend('GRE', 'SE', 'SSFP') + + +%% grase sequence + +% file_m1{4}{1} = '../../outputs/grase_m1_0.dat'; +% file_m1{4}{2} = '../../outputs/grase_m1_1.dat'; +TR = 0.2; +EcoSpc = 0.005; % subplot(2,1,2) % relative_signal = signal_magnitude{seq}(:,:,:,:,2) - signal_magnitude{seq}(:,:,:,:,1); % relative_signal = 100 * (1 - signal_magnitude{seq}(:,:,:,:,1)./ signal_magnitude{seq}(:,:,:,:,2)); @@ -67,15 +80,12 @@ % % title('BOLD Signal') % colormap(inferno); % colorbar -% end -end -legend('GRE', 'SE', 'SSFP') %% stimulated echo clc % clear -fname{1} = '/DATA/aaghaeifar/Nextcloud/Projects/microvascular/outputs/ste_m1_0.dat'; -fname{2} = '/DATA/aaghaeifar/Nextcloud/Projects/microvascular/outputs/ste_m1_1.dat'; +file_m1{1} = '/DATA/aaghaeifar/Nextcloud/Projects/microvascular/outputs/ste_m1_0.dat'; +file_m1{2} = '/DATA/aaghaeifar/Nextcloud/Projects/microvascular/outputs/ste_m1_1.dat'; dim_echo = 2; dim_spin = 3; @@ -84,8 +94,8 @@ spins_xyz = []; -for i=1:numel(fname) - [m_xyz, dims, scales] = read_spinwalk(fname{i}); +for i=1:numel(file_m1) + [m_xyz, dims, scales] = read_spinwalk(file_m1{i}); spins_xyz = cat(ndims(m_xyz)+1, spins_xyz, m_xyz); end vessel_radius = rad_ref_um * scales; @@ -124,7 +134,19 @@ hold off legend('STE Rest', 'STE Act', 'SE Rest', 'SE Act') +%% random-walk trajectory +clear +clc; +filename = '../../outputs/gre_trajectory_xyz1_fieldmap_0.dat'; +[xyz_all, dims, hdr_extra] = read_spinwalk(filename); +size(xyz_all) +xyz_all = xyz_all * 1e6; +n_spins = size(xyz_all, 3); +for s=round(linspace(1, n_spins, 10)) + xyz = xyz_all(:,:,s,2); + plot3(xyz(1,:), xyz(2,:), xyz(3,:)) + hold on; +end +hold off - - - +% nonzeros(xyz_all) diff --git a/MATLAB/read_spinwalk.m b/MATLAB/read_spinwalk.m index da20872..c92a00b 100644 --- a/MATLAB/read_spinwalk.m +++ b/MATLAB/read_spinwalk.m @@ -1,11 +1,11 @@ function [m_xyz, dims, hdr_extra] = read_spinwalk(filename) -fileID = fopen(filename); -hdr_size = fread(fileID, 1, 'int32=>int32'); -dims = fread(fileID, 4, 'int32=>int32'); % 3 * n_echo * n_spins * n_sample_length_scales -hdr_extra = fread(fileID, hdr_size/4-numel(dims), 'single=>single'); -m_xyz = fread(fileID, prod(dims), 'single=>single'); +fileID = fopen(filename); +hdr_size = fread(fileID, 1, 'int32=>int32'); +dims = fread(fileID, 4, 'int32=>int32'); % 3 * n_echo * n_spins * n_sample_length_scales +hdr_extra = fread(fileID, hdr_size/4-numel(dims), 'single=>single'); +m_xyz = fread(fileID, prod(dims), 'single=>single'); fclose(fileID); -m_xyz = reshape(m_xyz, dims(:)'); -m_xyz = reshape(m_xyz, [dims(1), dims(2), dims(3), dims(4)]); +m_xyz = reshape(m_xyz, dims(:)'); +m_xyz = reshape(m_xyz, [dims(1), dims(2), dims(3), dims(4)]);