Brain-wide mapping of layer-specific functional connectivity in the human cortex at 3T using draining-vein-suppressed fMRI.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches

[1] § Materials and methods › Data preprocessing ↔ welton_toolbox_github/NIFTI_NORDIC_wt.m, lines 1–53 · score 0.74 · NIFTI_NORDIC, smoothing filter, fMRI, MATLAB, magnitude, width
[2] § Materials and methods › Common protocol in MR acquisitions ↔ welton_toolbox_github/gen_mb3d_encodemtx_i2k.m, lines 1–61 · score 0.65 · acceleration rate, phase encoding, imaging protocols, dimensional, SMS, slicing
[3] § Materials and methods › Common protocol in MR acquisitions ↔ welton_toolbox_github/mb3d_grappa_wt3.m, lines 1–28 · score 0.65 · acceleration rate, phase encoding, imaging protocols, dimensional, SMS, slicing
[4] § Materials and methods › Data preprocessing ↔ pproc2_s2025010601_run1_dn2.m, lines 108–176 · score 0.60 · motion corrected, preprocessing, band, MELODIC, ICA, smoothing

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

MATLAB · 1,148 lines · 36 KB · MIT · 1 match

function NIFTI_NORDIC_wt(fn_magn_in,fn_phase_in,fn_out,ARG)
% fMRI
% fn_magn_in='name.nii.gz';
% fn_phase_in='name2.nii.gz';
% fn_out=['NORDIC_' fn_magn_in(1:end-7)];
% ARG.temporal_phase=1;
% ARG.phase_filter_width=10;
% NIFTI_NORDIC(fn_magn_in,fn_phase_in,fn_out,ARG)
%
% dMRI
% fn_magn_in='name.nii.gz';
% fn_phase_in='name2.nii.gz';
% fn_out=['NORDIC_' fn_magn_in(1:end-7)];
% ARG.temporal_phase=3;
% ARG.phase_filter_width=3;
% NIFTI_NORDIC(fn_magn_in,fn_phase_in,fn_out,ARG)
%
%
% file_input assumes 4D data
%
%OPTIONS
% ARG.DIROUT VAL= string Default is empty
% ARG.noise_volume_last VAL = num specifiec volume from the end of the series
% 0 default
%
% ARG.factor_error val = num >1 use higher noisefloor <1 use lower noisefloor
% 1 default
%
% ARG.full_dynamic_range val = [ 0 1] 0 keep the input scale, output maximizes range.
% Default 0
% ARG.temporal_phase val = [1 2 3] 1 was default, 3 now in dMRI due tophase errors in some data
% ARG.NORDIC val = [0 1] 1 Default
% ARG.MP val = [0 1 2] 1 NORDIC gfactor with MP estimation.
% 2 MP without gfactor correction
% 0 default
% ARG.kernel_size_gfactor val = [val1 val2 val], defautl is [14 14 1]
% ARG.kernel_size_PCA val = [val1 val2 val], default is val1=val2=val3;
% ratio of 11:1 between spatial and temproal voxels
% ARG.magnitude_only val =[] or 1. Using complex or magntiude only. Default is []
% Function still needs two inputs but will ignore the second
%
% ARG.save_add_info val =[0 1]; If it is 1, then an additonal matlab file is being saved with degress removed etc.
% default is 0
% ARG.make_complex_nii if the field exist, then the phase is being saved in a similar format as the input phase
%
% ARG.phase_slice_average_for_kspace_centering val = [0 1]
% if val =0, not used, if val=1 the series average pr slice is first removed
% default is now 0
% ARG.phase_filter_width val = [1... 10] Specifiec the width of the smoothing filter for the phase
% default is now 3
%
% ARG.save_gfactor_map val = [1 2]. 1, saves the RELATIVE gfactor, 2 saves the
% gfactor and does not complete the NORDIC processing
% TODO
% Scaling relative to the width of the MP spectrum, if one wants to be
% conservative
%
% 4/15/21 swapped the uint16 and in16 for the phase
%
% VERSION 4/22/2021
if ~exist('ARG') % initialize ARG structure
ARG.DIROUT=[pwd '/'];
elseif ~isfield(ARG,'DIROUT') % Specify where to save data
ARG.DIROUT=[pwd '/'];
else
ARG.DIROUT=[ARG.DIROUT '/'];
end
if ~isfield(ARG,'noise_volume_last')
ARG.noise_volume_last=0; % there is no noise volume {0 1 2 ...}
end
if ~isfield(ARG,'factor_error')
ARG.factor_error=1.0; % error in gfactor estimatetion. >1 use higher noisefloor <1 use lower noisefloor
end
if ~isfield(ARG,'full_dynamic_range')
ARG.full_dynamic_range=0; % Format o
end
if ~isfield(ARG,'temporal_phase')
ARG.temporal_phase=1; % Correction for slice and time-specific phase
end
if ~isfield(ARG,'NORDIC') & ~isfield(ARG,'MP')
ARG.NORDIC=1; % threshold based on Noise
ARG.MP=0; % threshold based on Marchencko-Pastur
elseif ~isfield(ARG,'NORDIC') % MP selected
if ARG.MP==1
ARG.NORDIC=0;
else
ARG.NORDIC=1;
end
elseif ~isfield(ARG,'MP') % NORDIC selected
if ARG.NORDIC==1
ARG.MP=0;
else
ARG.MP=1;
end
end
if ~isfield(ARG,'phase_filter_width')
ARG.phase_filter_width=3; % default is [14 14 90]
end
if ~isfield(ARG,'NORDIC_patch_overlap')
ARG.NORDIC_patch_overlap=2; % default is [14 14 90]
end
if ~isfield(ARG,'gfactor_patch_overlap')
ARG.gfactor_patch_overlap=2; % default is [14 14 90]
end
if ~isfield(ARG,'kernel_size_gfactor')
ARG.kernel_size_gfactor=[]; % default is [14 14 90]
end
if ~isfield(ARG,'kernel_size_PCA')
ARG.kernel_size_PCA=[]; % default is 11:1 ratio
end
if ~isfield(ARG,'phase_slice_average_for_kspace_centering');
ARG.phase_slice_average_for_kspace_centering=0;
end
if ~isfield(ARG,'magnitude_only') % if legacy data
ARG.magnitude_only=0; %
end
if isfield(ARG,'save_add_info'); end % additional information is saved in matlab file
if isfield(ARG,'make_complex_nii'); end % two output NII files are saved
if ~isfield(ARG,'save_gfactor_map') % save out a map of a relative gfactor
ARG.save_gfactor_map=[]; %
end
if isfield(ARG,'use_generic_NII_read') % save out a map of a relative gfactor
if ARG.use_generic_NII_read==1
path(path,'/home/range6-raid1/moeller/matlab/ADD/NIFTI/');
end
else
ARG.use_generic_NII_read=0;
end
if ~isfield(ARG,'data_has_zero_elements') %
ARG.data_has_zero_elements=0; % % If there are pixels that are constant zero
end
ARG;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ARG.magnitude_only~=1
try
info_phase=niftiinfo(fn_phase_in);
info=niftiinfo(fn_magn_in);
catch; disp('The niftiinfo fails at reading the header') ;end
if ARG.use_generic_NII_read~=1
tmpvol=load_nii(fn_magn_in);
%I_M=abs(single(niftiread(fn_magn_in)));
I_M=abs(single(tmpvol.img));
tmpvol=load_nii(fn_phase_in);
%I_P=single(niftiread(fn_phase_in));
I_P=single(tmpvol.img*pi/180);
else
try
tmp=load_nii(fn_magn_in);
I_M=abs(single(tmp.img));
tmp=load_nii(fn_phase_in);
I_P=single(tmp.img);
catch
disp('Missing nfiti tool. Serach mathworks for load_nii fileexchange 8797')
end
end
phase_range=single(max(I_P(:)));
phase_range_min=single(min(I_P(:)));
if ~exist('info_phase')
info_phase.Datatype=class(I_P);
info.Datatype=class(I_M);
end
% Here, we combine magnitude and phase data into complex form
fprintf('Phase should be -pi to pi...\n')
% convert to single and then scale the phase
I_P = single(I_P);
range_norm=phase_range-phase_range_min;
range_center=(phase_range+phase_range_min)/range_norm*1/2;
I_P = (single(I_P)./range_norm -range_center)*2*pi;
II=single(I_M) .* exp(1i*I_P);
if 0
if strmatch(info_phase.Datatype,'uint16')
I_P = single(I_P)/phase_range*2*pi;
II=single(I_M) .* exp(1i*I_P);
elseif strmatch(info_phase.Datatype,'int16')
I_P = (single(I_P)+1-(phase_range+1)/2)/(phase_range+1)*2*pi;
II=single(I_M) .* exp(1i*I_P);
elseif strmatch(info_phase.Datatype,'single')
phase_range_min=min(I_P(:));
range_norm=phase_range-phase_range_min;
range_center=(phase_range+phase_range_min)/range_norm*1/2;
I_P = (single(I_P)./range_norm -range_center)*2*pi;
II=single(I_M) .* exp(1i*I_P);
end
end
fprintf('Phase data range is %.2f to %.2f\n', min(I_P(:)), max(I_P(:)))
else
try
info=niftiinfo(fn_magn_in);
catch; disp('The niftiinfo fails at reading the header') ;end
if ARG.use_generic_NII_read~=1
tmpvol=load_nii(fn_magn_in);
%I_M=abs(single(niftiread(fn_magn_in)));
I_M=abs(single(tmpvol.img));
else
tmp=load_nii(fn_magn_in);
I_M=abs(single(tmp.img));
end
if ~exist('info_phase')
info.Datatype=class(I_M);
end
end
if ~isempty(ARG.magnitude_only)
if ARG.magnitude_only==1
II=single(I_M);
ARG.temporal_phase=0;
end
end
TEMPVOL=abs(II(:,:,:,1));
ARG.ABSOLUTE_SCALE=min(TEMPVOL(TEMPVOL~=0));
II=II./ARG.ABSOLUTE_SCALE;
% load test_DATA
% II=II(:,:,:,1:95);
if size(II,4)<6
disp('Too few volumes')
% return
end
KSP2=II;
matdim=size(KSP2);
tt=mean(reshape(abs(KSP2),[],size(KSP2,4)));
[idx]=find(tt>0.95*max(tt));
meanphase=mean(KSP2(:,:,:,idx(1)),4);
if 1
disp('estimating slice-dependent phases ...')
meanphase=mean(KSP2(:,:,:,[1:end-ARG.noise_volume_last]),4);
for nsl=1:size(meanphase,3);
%meanphase2(:,:,nsl)=complex( medfilt2(squeeze(real(meanphase(:,:,nsl))),[7 7]), medfilt2(squeeze(imag(meanphase(:,:,nsl))),[7 7]) );
end
meanphase=meanphase*ARG.phase_slice_average_for_kspace_centering;
end
for slice=matdim(3):-1:1
for n=1:size(KSP2,4); % include the noise
KKSP2(:,:,slice,n)=KSP2(:,:,slice,n).*exp(-i*angle(meanphase(:,:,slice)));
end
end
DD_phase=0*KSP2;
if ARG.temporal_phase>0; % Standarad low-pass filtered map
for slice=matdim(3):-1:1
for n=1:size(KSP2,4);
tmp=KSP2(:,:,slice,n);
for ndim=[1:2]; tmp=ifftshift(ifft(ifftshift( tmp ,ndim),[],ndim),ndim+0); end
[nx, ny, nc, nb] = size(tmp(:,:,:,:,1,1));
tmp = bsxfun(@times,tmp,reshape(tukeywin(ny,1).^ARG.phase_filter_width,[1 ny]));
tmp = bsxfun(@times,tmp,reshape(tukeywin(nx,1).^ARG.phase_filter_width,[nx 1]));
for ndim=[1:2]; tmp=fftshift(fft(fftshift( tmp ,ndim),[],ndim),ndim+0); end
DD_phase(:,:,slice,n)=tmp;
end
end
end
if ARG.temporal_phase==2; % Secondary step for filtered phase with residual spikes
for slice=matdim(3):-1:1
for n=1:size(KSP2,4);
phase_diff=angle(KSP2(:,:,slice,n)./DD_phase(:,:,slice,n));
mask=abs(phase_diff)>1;
DD_phase2=DD_phase(:,:,slice,n);
tmp=(KSP2(:,:,slice,n));
DD_phase2(mask)= tmp(mask);
DD_phase(:,:,slice,n)=DD_phase2;
end
end
end
for slice=matdim(3):-1:1
for n=1:size(KSP2,4);
KSP2(:,:,slice,n)= KSP2(:,:,slice,n).*exp(-i*angle( DD_phase(:,:,slice,n) ));
end
end
disp('Completed estimating slice-dependent phases ...')
if isfield(ARG,'use_magn_for_gfactor')
if isempty(ARG.kernel_size_gfactor) | size(ARG.kernel_size_gfactor,2)<3
KSP2=abs(KSP2(:,:,1:end,1:min(90,end),1)); % should be at least 30 volumes
else
KSP2=abs(KSP2(:,:,1:end,1:min(ARG.kernel_size_gfactor(3),end),1));
end
else
if (isempty(ARG.kernel_size_gfactor) | size(ARG.kernel_size_gfactor,2)<3)
KSP2=(KSP2(:,:,1:end,1:min(90,end),1)); % should be at least 30 volumes
else
% KSP2=(KSP2(:,:,1:end,1:min(ARG.kernel_size_gfactor(3),end),1));
KSP2=(KSP2(:,:,1:end,1:min(ARG.kernel_size_gfactor(4),end),1));
end
end
KSP2(isnan(KSP2))=0;
KSP2(isinf(KSP2))=0;
master_fast=1;
KSP_recon=0*KSP2;
if isempty(ARG.kernel_size_gfactor)
ARG.kernel_size=[14 14 1];
else
ARG.kernel_size=[ARG.kernel_size_gfactor(1) ARG.kernel_size_gfactor(2) 1];
ARG.kernel_size=[ARG.kernel_size_gfactor(1) ARG.kernel_size_gfactor(2) ARG.kernel_size_gfactor(3)];
end
QQ.KSP_processed=zeros(1,size(KSP2,1)-ARG.kernel_size(1));
ARG.patch_average=0;
ARG.patch_average_sub= ARG.gfactor_patch_overlap;
ARG.LLR_scale=0;
ARG.NVR_threshold=1;
ARG.soft_thrs=10; % MPPCa (When Noise varies)
%ARG.soft_thrs=[]; % NORDIC (When noise is flat)
KSP_weight=KSP2(:,:,:,1)*0;
NOISE=KSP_weight;
Component_threshold=KSP_weight;
energy_removed=KSP_weight;
SNR_weight=KSP_weight;
QQ.KSP_processed=zeros(1,size(KSP2,1)-ARG.kernel_size(1));
if ARG.patch_average==0
KSP_processed=QQ.KSP_processed*0;
for nw1=2:max(1,floor(ARG.kernel_size(1)/ARG.patch_average_sub));
KSP_processed(1,nw1 : max(1,floor(ARG.kernel_size(1)/ARG.patch_average_sub)):end)=2;
end
KSP_processed(end)=0; % disp
QQ.KSP_processed=KSP_processed;
end
KSP_processed;
disp('estimating g-factor ...')
QQ.ARG=ARG;
for n1=1:size(QQ.KSP_processed,2)
% fprintf( [num2str(n1) ' '])
[KSP_recon,~,KSP_weight,NOISE,Component_threshold,energy_removed,SNR_weight]=sub_LLR_Processing(KSP_recon,KSP2,ARG,n1,QQ,master_fast,KSP_weight,NOISE,Component_threshold,energy_removed,SNR_weight) ; % save all files
% fprintf( [num2str(n1) ' '])
% if mod(n1,20)==0; disp(' ');end
end
KSP_recon=KSP_recon./repmat((KSP_weight),[1 1 1 size(KSP2,4)]);
ARG.NOISE= sqrt(NOISE./KSP_weight);
ARG.Component_threshold = Component_threshold./KSP_weight;
ARG.energy_removed = energy_removed./KSP_weight;
ARG.SNR_weight = SNR_weight./KSP_weight;
IMG2=KSP_recon;
ARG2=ARG;
disp('completed estimating g-factor')
gfactor=ARG.NOISE;
if size(KSP2,4)<6; % gfactor stimation most likely failed, replace with median estimated
gfactor(isnan(gfactor))=0;
gfactor(gfactor==0)=median(gfactor(gfactor~=0));
end
if sum(gfactor(:)==0)>0; % gfactor stimation most likely failed since it is zero
gfactor(isnan(gfactor))=0;
gfactor(gfactor<1)=median(gfactor(gfactor~=0));
ARG.data_has_zero_elements=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ARG.MP==2;
gfactor=ones(size(gfactor));
end
% gfactor=ones(size(gfactor));
if ( ARG.save_gfactor_map==2 ) | ( ARG.save_gfactor_map==1 )
g_IMG=abs(gfactor(:,:,:,1:end)); % remove g-factor and noise for DUAL 1
g_IMG(isnan(g_IMG))=0;
tmp=sort(abs(g_IMG(:))); sn_scale=2*tmp(round(0.99*end));%sn_scale=max();
gain_level=floor(log2(32000/sn_scale));
if ARG.full_dynamic_range==0; gain_level=0;end
if strmatch(info.Datatype,'uint16')
g_IMG= uint16(abs(g_IMG)*2^gain_level);
elseif strmatch(info.Datatype,'int16')
g_IMG= int16(abs(g_IMG)*2^gain_level);
else
g_IMG= single(abs(g_IMG)*2^gain_level);
end
%niftiwrite((g_IMG),[ARG.DIROUT 'gfactor_' fn_out(1:end) '.nii'])
nii=make_nii(g_IMG);
save_nii(nii,[ARG.DIROUT 'gfactor_' fn_out(1:end) '.nii']);
if ARG.save_gfactor_map==2
return
end
end
KSP2=II;
matdim=size(KSP2);
for slice=matdim(3):-1:1
for n=1:size(KSP2,4); % include the noise
KSP2(:,:,slice,n)=KSP2(:,:,slice,n).*exp(-i*angle(meanphase(:,:,slice)));
end
end
for n=1:size(KSP2,4);
KSP2(:,:,:,n)= KSP2(:,:,:,n)./ gfactor;
end
if ARG.noise_volume_last>0
KSP2_NOISE =KSP2(:,:,:,end+1-ARG.noise_volume_last);
end
if ARG.temporal_phase==3; % Secondary step for filtered phase with residual spikes
for slice=matdim(3):-1:1
for n=1:size(KSP2,4);
phase_diff=angle(KSP2(:,:,slice,n)./DD_phase(:,:,slice,n));
mask = abs(phase_diff)>1;
mask2 = abs(KSP2(:,:,slice,n))>sqrt(2);
DD_phase2=DD_phase(:,:,slice,n);
tmp=(KSP2(:,:,slice,n));
DD_phase2(mask.*mask2==1)= tmp(mask.*mask2==1);
DD_phase(:,:,slice,n)=DD_phase2;
end
end
end
for slice=matdim(3):-1:1
for n=1:size(KSP2,4);
KSP2(:,:,slice,n)= KSP2(:,:,slice,n).*exp(-i*angle( DD_phase(:,:,slice,n) ));
end
end
KSP2(isnan(KSP2))=0;
KSP2(isinf(KSP2))=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ARG.noise_volume_last>0
%tmp_noise=KSP2(:,:,:,end+1-ARG.noise_volume_last);
tmp_noise=KSP2_NOISE;
tmp_noise(isnan(tmp_noise))=0;
tmp_noise(isinf(tmp_noise))=0;
ARG.measured_noise=std(tmp_noise(tmp_noise~=0)); % sqrt(2) for real and complex
else
ARG.measured_noise=1; % IF COMPLEX DATA
end
if ~isfield(ARG,'use_magn_for_gfactor') & (isempty(ARG.magnitude_only) | ARG.magnitude_only==0) %% WOULD THIS BE THE ISSUE & replaced by |
ARG.measured_noise = ARG.measured_noise/sqrt(2)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ARG.data_has_zero_elements==1
MASK=(sum(abs(KSP2),4)==0);
Num_zero_elements=sum(MASK(:));
for nvol=1:size(KSP2,4)
tmp=KSP2(:,:,:,nvol);
tmp(MASK)=(randn(Num_zero_elements,1)+1i*randn(Num_zero_elements,1))/sqrt(2);
KSP2(:,:,:,nvol)=tmp;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for readout=1:size(KSP2,1)
TT(readout)= std(reshape(KSP2(readout,:,:,:),[],1));
TT1(readout)= mean(reshape(KSP2(readout,:,:,:),[],1));
end
% ARG.measured_noise=median(TT([2:8 end-8:end-1]))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
master_fast=1;
KSP_recon=0*KSP2;
ARG.kernel_size=repmat([ round((size(KSP2,4)*11)^(1/3)) ],1,3);
if isempty(ARG.kernel_size_PCA)
ARG.kernel_size=repmat([ round((size(KSP2,4)*11)^(1/3)) ],1,3);
else
ARG.kernel_size = ARG.kernel_size_PCA ;
end
if matdim(3) <= ARG.kernel_size(3) % Number of slices is less than cubic kernel
ARG.kernel_size = repmat([ round((size(KSP2,4)*11/matdim(3) )^(1/2)) ],1,2);
ARG.kernel_size(3)= matdim(3);
end
QQ.KSP_processed=zeros(1,size(KSP2,1)-ARG.kernel_size(1));
ARG.patch_average=0;
ARG.patch_average_sub= ARG.NORDIC_patch_overlap ;
% ARG.kernel_size=[7 7 7]; ARG.patch_average_sub=7; MPPCA
% ARG.soft_thrs=10; % MPPCa (When Noise varies)
ARG.LLR_scale=1;
ARG.NVR_threshold=0;
ARG.soft_thrs=[]; % NORDIC (When noise is flat)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ARG.NVR_threshold=0;
for ntmp=1:10
[~,S,~]=svd(randn(prod(ARG.kernel_size),size(KSP2,4)));
ARG.NVR_threshold=ARG.NVR_threshold+S(1,1);
end
if ARG.magnitude_only~=1 % 4/29/2021
ARG.NVR_threshold= ARG.NVR_threshold/10*sqrt(2)* ARG.measured_noise*ARG.factor_error; % sqrt(2) due to complex 1.20 due to understimate of g-factor
else
ARG.NVR_threshold= ARG.NVR_threshold/10*ARG.measured_noise*ARG.factor_error; % sqrt(2) due to complex 1.20 due to understimate of g-factor
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ARG.MP>0
ARG.soft_thrs=10;
end
if isfield(ARG,'soft_thrs_in')
ARG.soft_thrs=ARG.soft_thrs_in;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
KSP_weight=KSP2(:,:,:,1)*0;
NOISE=KSP_weight;
Component_threshold=KSP_weight;
energy_removed=KSP_weight;
SNR_weight=KSP_weight;
QQ.KSP_processed=zeros(1,size(KSP2,1)-ARG.kernel_size(1));
if ARG.patch_average==0
KSP_processed=QQ.KSP_processed*0;
for nw1=2:max(1,floor(ARG.kernel_size(1)/ARG.patch_average_sub));
KSP_processed(1,nw1 : max(1,floor(ARG.kernel_size(1)/ARG.patch_average_sub)):end)=2;
end
KSP_processed(end)=0; % disp
QQ.KSP_processed=KSP_processed;
end
disp('starting NORDIC ...')
KSP_processed;
QQ.ARG=ARG;
for n1=1:size(QQ.KSP_processed,2)
[KSP_recon,~,KSP_weight,NOISE,Component_threshold,energy_removed,SNR_weight]=sub_LLR_Processing(KSP_recon,KSP2,ARG,n1,QQ,master_fast,KSP_weight,NOISE,Component_threshold,energy_removed,SNR_weight) ; % save all files
end
KSP_recon=KSP_recon./repmat((KSP_weight),[1 1 1 size(KSP2,4)]); % Assumes that the combination is with N instead of sqrt(N). Works for NVR not MPPCA
ARG.NOISE= sqrt(NOISE./KSP_weight);
ARG.Component_threshold = Component_threshold./KSP_weight;
ARG.energy_removed = energy_removed./KSP_weight;
ARG.SNR_weight = SNR_weight./KSP_weight;
IMG2=KSP_recon;
disp('completing NORDIC ...')
if isfield(ARG,'save_residual_matlab')
if ARG.save_residual_matlab==1;
Residual=KSP2-KSP_recon;
save([ARG.DIROUT 'RESIDUAL' fn_out '.mat' ],'Residual','-v7.3')
end
end
for n=1:size(IMG2,4);
IMG2(:,:,:,n)= IMG2(:,:,:,n).* gfactor;
end
for slice=matdim(3):-1:1
for n=1:size(IMG2,4); % include the noise
IMG2(:,:,slice,n)=IMG2(:,:,slice,n).*exp(i*angle(meanphase(:,:,slice)));
end
end
for slice=matdim(3):-1:1
for n=1:size(IMG2,4);
IMG2(:,:,slice,n)= IMG2(:,:,slice,n).*exp(i*angle( DD_phase(:,:,slice,n) ));
end
end
IMG2=IMG2.*ARG.ABSOLUTE_SCALE;
IMG2(isnan(IMG2))=0;
if isfield(ARG,'make_complex_nii')
IMG2_tmp=abs(IMG2(:,:,:,1:end)); % remove g-factor and noise for DUAL 1
IMG2_tmp(isnan(IMG2_tmp))=0;
tmp=sort(abs(IMG2_tmp(:))); sn_scale=2*tmp(round(0.99*end));%sn_scale=max();
gain_level=floor(log2(32000/sn_scale));
%IMG2_tmp= int16(abs(IMG2_tmp)*2^gain_level);
if ARG.full_dynamic_range==0; gain_level=0;end
if strmatch(info.Datatype,'uint16')
IMG2_tmp= uint16(abs(IMG2_tmp)*2^gain_level);
elseif strmatch(info.Datatype,'int16')
IMG2_tmp= int16(abs(IMG2_tmp)*2^gain_level);
else
IMG2_tmp= single(abs(IMG2_tmp)*2^gain_level);
end
%niftiwrite((IMG2_tmp),[ARG.DIROUT fn_out 'magn.nii'],info)
nii=make_nii(IMG2_tmp);
save_nii(nii,[ARG.DIROUT fn_out 'magn.nii']);
IMG2_tmp=angle(IMG2(:,:,:,1:end));
if strmatch(info_phase.Datatype,'int16')
% IMG2_tmp=IMG2_tmp+pi;
end
IMG2_tmp= (IMG2_tmp/(2*pi)+range_center)*range_norm;
if strmatch(info_phase.Datatype,'uint16')
IMG2_tmp= uint16(IMG2_tmp);
elseif strmatch(info_phase.Datatype,'int16')
IMG2_tmp= int16(IMG2_tmp);
else
IMG2_tmp= single((IMG2_tmp));
end
if 0
if strmatch(info_phase.Datatype,'uint16')
IMG2_tmp=IMG2_tmp/(2*pi)*phase_range;
IMG2_tmp= uint16(abs(IMG2_tmp)*2^gain_level);
elseif strmatch(info_phase.Datatype,'int16')
IMG2_tmp=IMG2_tmp/(2*pi)*phase_range;
IMG2_tmp= int16((IMG2_tmp)*2^gain_level);
else
IMG2_tmp= single(abs(IMG2_tmp)*2^gain_level);
end
end
%niftiwrite((IMG2_tmp),[ARG.DIROUT fn_out 'phase.nii'],info_phase)
nii=make_nii(IMG2_tmp);
save_nii(nii,[ARG.DIROUT fn_out 'phase.nii']);
else
IMG2=abs(IMG2(:,:,:,1:end)); % remove g-factor and noise for DUAL 1
IMG2(isnan(IMG2))=0;
tmp=sort(abs(IMG2(:))); sn_scale=2*tmp(round(0.99*end));%sn_scale=max();
gain_level=floor(log2(32000/sn_scale));
if ARG.full_dynamic_range==0; gain_level=0;end
if strmatch(info.Datatype,'uint16')
IMG2= uint16(abs(IMG2)*2^gain_level);
elseif strmatch(info.Datatype,'int16')
IMG2= int16(abs(IMG2)*2^gain_level);
else
IMG2= single(abs(IMG2)*2^gain_level);
end
if ARG.use_generic_NII_read==0;
%niftiwrite((IMG2),[ARG.DIROUT fn_out(1:end) '.nii'],info)
nii=make_nii(IMG2);
save_nii(nii, [ARG.DIROUT fn_out(1:end) '.nii']);
else
nii=make_nii(IMG2);
save_nii(nii, [ARG.DIROUT fn_out(1:end) '.nii']);
end
end
if isfield(ARG,'save_add_info')
if ARG.save_add_info==1
disp('saving additional info')
save([ARG.DIROUT fn_out '.mat' ],'ARG2','ARG','-v7.3')
end
end
return
function [KSP_recon,KSP2,KSP2_weight,NOISE, Component_threshold,energy_removed,SNR_weight]=sub_LLR_Processing(KSP_recon,KSP2,ARG,n1,QQ,master,KSP2_weight,NOISE,Component_threshold,energy_removed,SNR_weight)
if ~exist('NOISE'); NOISE=[]; end
if ~exist('Component_threshold');Component_threshold=[]; end
if ~exist('energy_removed'); energy_removed=[]; end
if ~exist('SNR_weight'); SNR_weight=[]; end
% QQ.KSP_processed 0 nothing done, 1 running, 2 saved 3 completed and averaged
if master==0 && ARG.patch_average==0
OPTION='NO_master_NO_PA';
elseif master==0 && ARG.patch_average==1
OPTION='NO_master_PA' ;
elseif master==1 && ARG.patch_average==0
OPTION='master_NO_PA' ;
elseif master==1 && ARG.patch_average==1
OPTION='master_PA' ;
end
switch OPTION
case 'NO_master_NO_PA'
case 'NO_master_PA'
case 'master_NO_PA'
case 'master_PA'
end
if QQ.KSP_processed(1,n1)~=1 && QQ.KSP_processed(1,n1)~=3 % not being processed also not completed yet
if QQ.KSP_processed(1,n1)==2 && master==1% processed but not added.
% loading instead of processing
% load file as soon as save, if more than 10 sec, just do the recon
% instead.
try % try to load otherwise go to next slice
load([ARG.filename 'slice' num2str(n1) '.mat'],'DATA_full2')
catch;
QQ.KSP_processed(1,n1)=0; % identified as bad file and being identified for reprocessing
return ;end
end
if QQ.KSP_processed(1,n1)~=2
QQ.KSP_processed(1,n1)=1; % block for other processes
if ~exist('DATA_full2')
ARG2=QQ.ARG;
if master==0
QQ.KSP_processed(1,n1)=1 ; % STARTING
KSP2a=QQ.KSP2([1:ARG.kernel_size(1)]+(n1-1),:,:,:); lambda=ARG2.LLR_scale*ARG.NVR_threshold;
else
QQ.KSP_processed(1,n1)=1 ; % STARTING
KSP2a=KSP2([1:ARG.kernel_size(1)]+(n1-1),:,:,:); lambda=ARG2.LLR_scale*ARG.NVR_threshold;
end
if ARG.patch_average==1
% [DATA_full2, ~,NOISE, Component_threshold] =subfunction_loop_for_NVR_avg(KSP2a,ARG.kernel_size(3),ARG.kernel_size(2),ARG.kernel_size(1),lambda,1,ARG.soft_thrs);
[DATA_full2, KSP2_weight] =subfunction_loop_for_NVR_avg(KSP2a,ARG.kernel_size(3),ARG.kernel_size(2),ARG.kernel_size(1),lambda,1,ARG.soft_thrs,KSP2_weight);
else
KSP2_weight_tmp =KSP2_weight([1:ARG.kernel_size(1)]+(n1-1),:,:,:);
NOISE_tmp =NOISE([1:ARG.kernel_size(1)]+(n1-1),:,:,:);
Component_threshold_tmp =Component_threshold([1:ARG.kernel_size(1)]+(n1-1),:,:,:);
energy_removed_tmp =energy_removed([1:ARG.kernel_size(1)]+(n1-1),:,:,:);
SNR_weight_tmp =SNR_weight([1:ARG.kernel_size(1)]+(n1-1),:,:,:);
[DATA_full2,KSP2_weight_tmp,NOISE_tmp, Component_threshold_tmp,energy_removed_tmp,SNR_weight_tmp] =...
subfunction_loop_for_NVR_avg_update(KSP2a,ARG.kernel_size(3),ARG.kernel_size(2),ARG.kernel_size(1),lambda,1,ARG.soft_thrs,KSP2_weight_tmp,ARG,NOISE_tmp,Component_threshold_tmp,energy_removed_tmp,SNR_weight_tmp);
KSP2_weight([1:ARG.kernel_size(1)]+(n1-1),:,:,:)=KSP2_weight_tmp;
try; NOISE([1:ARG.kernel_size(1)]+(n1-1),:,:,:) =NOISE_tmp; catch;end
Component_threshold([1:ARG.kernel_size(1)]+(n1-1),:,:,:) = Component_threshold_tmp;
energy_removed([1:ARG.kernel_size(1)]+(n1-1),:,:,:) = energy_removed_tmp;
SNR_weight([1:ARG.kernel_size(1)]+(n1-1),:,:,:) = SNR_weight_tmp;
%DATA_full=subfunction_loop_for_NVR(KSP2a,ARG.kernel_size(3),ARG.kernel_size(2),ARG.kernel_size(1),lambda);
%DATA_full2(1, round(w2/2)+[1:size(DATA_full,1)],:,: )=DATA_full; % center plane only
end
end
end
if master==0
if QQ.KSP_processed(1,n1)~=2
save([ARG.filename 'slice' num2str(n1) '.mat'],'DATA_full2', '-v7.3' )
QQ.KSP_processed(1,n1)=2 ; % COMPLETED
end
else
if ARG.patch_average==1
tmp=KSP_recon([1:ARG.kernel_size(1)]+(n1-1),:,:,:) ;
KSP_recon([1:ARG.kernel_size(1)]+(n1-1),:,:,:)= tmp + DATA_full2;
else
KSP_recon([1:ARG.kernel_size(1)]+(n1-1) ,1:size(DATA_full2,2),:,:)= KSP_recon([1:ARG.kernel_size(1)]+(n1-1) ,1:size(DATA_full2,2),:,:) + DATA_full2;
end
QQ.KSP_processed(1,n1)=3 ;
end
end
return
function [KSP2_tmp_update, KSP2_weight]=subfunction_loop_for_NVR_avg(KSP2a,w3,w2,w1,lambda2,patch_avg, soft_thrs,KSP2_weight,ARG)
if ~exist('patch_avg'); patch_avg=1;end
if ~exist('soft_thrs'); soft_thrs=[]; end
if ~exist('KSP2_weight')
KSP2_weight=zeros(size(KSP2a(:,:,:,1)));
elseif isempty(KSP2_weight)
KSP2_weight=zeros(size(KSP2a(:,:,:,1)));
end
if ~exist('KSP2_tmp_update')
KSP2_tmp_update=zeros(size(KSP2a(:,:,:,:)));
elseif isempty(KSP2_tmp_update)
KSP2_tmp_update=zeros(size(KSP2a(:,:,:,:)));
end
% KSP2_weight=zeros(size(KSP2a(:,:,:,1)));
% KSP2_tmp_update=zeros(size(KSP2a));
% for n2=1:size(KSP2a,2)-w2+1;
% for n3=1:size(KSP2a,3)-w3+1;
for n2=[1: max(1,floor(w2/ARG.patch_average_sub)):size(KSP2a,2)*1-w2+1 size(KSP2a,2)-w2+1];
for n3=[1: max(1,floor(w3/ARG.patch_average_sub)):size(KSP2a,3)*1-w3+1 size(KSP2a,3)-w3+1 ];
KSP2_tmp=KSP2a(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:);
tmp1=reshape(KSP2_tmp,[],size(KSP2_tmp,4));
[U,S,V]=svd([(tmp1) ],'econ');
S=diag(S);
[idx]=sum(S<lambda2);
if isempty(soft_thrs)
S(S<lambda2)=0;
elseif soft_thrs==10 % USING MPPCA
% disp('test for zero entries')
Test_mat=sum(tmp1,2);
sum(Test_mat==0)
centering=0;
MM=size(tmp1,1);
NNN=size(tmp1,2);
R = min(MM, NNN);
scaling = (max(MM, NNN) - (0:R-centering-1)) / NNN;
scaling = scaling(:);
vals=S;
vals = (vals).^2 / NNN;
% First estimation of Sigma^2; Eq 1 from ISMRM presentation
csum = cumsum(vals(R-centering:-1:1)); cmean = csum(R-centering:-1:1)./(R-centering:-1:1)'; sigmasq_1 = cmean./scaling;
% Second estimation of Sigma^2; Eq 2 from ISMRM presentation
gamma = (MM - (0:R-centering-1)) / NNN;
rangeMP = 4*sqrt(gamma(:));
rangeData = vals(1:R-centering) - vals(R-centering);
sigmasq_2 = rangeData./rangeMP;
t = find(sigmasq_2 < sigmasq_1, 1);
S(t:end)=0;
else
S(max(1,end-floor(idx*soft_thrs)):end)=0;
end
tmp1=U*diag(S)*V';
tmp1=reshape(tmp1,size(KSP2_tmp));
if patch_avg==1
KSP2_tmp_update(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) =...
KSP2_tmp_update(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) +tmp1;
KSP2_weight(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) =...
KSP2_weight(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) +1;
else
KSP2_tmp_update(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) =...
KSP2_tmp_update(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) +tmp1(1,round(end/2),round(end/2),:);
KSP2_weight(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) =...
KSP2_weight(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) +1;
end
end
end
return
function [KSP2_tmp_update, KSP2_weight,NOISE,KSP2_tmp_update_threshold,energy_removed,SNR_weight]=subfunction_loop_for_NVR_avg_update(KSP2a,w3,w2,w1,lambda2,patch_avg, soft_thrs,KSP2_weight,ARG,NOISE,KSP2_tmp_update_threshold,energy_removed,SNR_weight)
if ~isfield(ARG,'patch_scale'); patch_scale=1; else; patch_scale=ARG.patch_scale;end
if ~exist('patch_avg'); patch_avg=1;end % patch_avg=0; means zero only
if ~exist('soft_thrs'); soft_thrs=[]; end
if ~exist('KSP2_weight')
KSP2_weight=zeros(size(KSP2a(:,:,:,1)));
elseif isempty(KSP2_weight)
KSP2_weight=zeros(size(KSP2a(:,:,:,1)));
end
if ~exist('NOISE_tmp')'% ~exist('NOISE_tmpKSP2_tmp_update')
NOISE_tmp=zeros(size(KSP2a(:,:,:,1)));
elseif isempty(KSP2_tmp_update)
NOISE_tmp=zeros(size(KSP2a(:,:,:,1)));
end
if ~exist('KSP2_tmp_update_threshold')
KSP2_tmp_update_threshold=zeros(size(KSP2a(:,:,:,1)));
elseif isempty(KSP2_tmp_update_threshold)
KSP2_tmp_update_threshold=zeros(size(KSP2a(:,:,:,1)));
end
if ~exist('energy_removed')
energy_removed=zeros(size(KSP2a(:,:,:,1)));
elseif isempty(energy_removed)
energy_removed=zeros(size(KSP2a(:,:,:,1)));
end
if ~exist('SNR_weight')
SNR_weight=zeros(size(KSP2a(:,:,:,1)));
elseif isempty(SNR_weight)
SNR_weight=zeros(size(KSP2a(:,:,:,1)));
end
KSP2_tmp_update=0*KSP2a;
%NOISE=[];
for n2=[1: max(1,floor(w2/ARG.patch_average_sub)):size(KSP2a,2)*1-w2+1 size(KSP2a,2)-w2+1];
for n3=[1: max(1,floor(w3/ARG.patch_average_sub)):size(KSP2a,3)*1-w3+1 size(KSP2a,3)-w3+1 ];
KSP2_tmp=KSP2a(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:);
tmp1=reshape(KSP2_tmp,[],size(KSP2_tmp,4));
[U,S,V]=svd([(tmp1) ],'econ');
S=diag(S);
[idx]=sum(S<lambda2);
if isempty(soft_thrs)
energy_scrub=sqrt(sum(S.^1)).\sqrt(sum(S(S<lambda2).^1));
S(S<lambda2)=0;
t=idx;
elseif soft_thrs~=10;
S=S-lambda2*soft_thrs;
S(S<0)=0;
energy_scrub=0;
t=1;
elseif soft_thrs==10 % USING MPPCA
% disp('test for zero entries')
Test_mat=sum(tmp1,2);
MM0=sum(Test_mat==0);
if MM0>1 & MM0<100
% 2
end
centering=0;
MM=size(tmp1,1)-MM0; % Correction for some zero entries
if MM>0
NNN=size(tmp1,2);
R = min(MM, NNN);
scaling = (max(MM, NNN) - (0:R-centering-1)) / NNN;
scaling = scaling(:);
vals=S;
vals = (vals).^2 / NNN;
% First estimation of Sigma^2; Eq 1 from ISMRM presentation
csum = cumsum(vals(R-centering:-1:1)); cmean = csum(R-centering:-1:1)./(R-centering:-1:1)'; sigmasq_1 = cmean./scaling;
% Second estimation of Sigma^2; Eq 2 from ISMRM presentation
gamma = (MM - (0:R-centering-1)) / NNN;
rangeMP = 4*sqrt(gamma(:));
rangeData = vals(1:R-centering) - vals(R-centering);
sigmasq_2 = rangeData./rangeMP;
t = find(sigmasq_2 < sigmasq_1, 1);
% NOISE(1:size(KSP2a,1),[1:w2]+(n2-1),[1:w3]+(n3-1),1) = sigmasq_2(t);
idx=size(S(t:end),1) ;
energy_scrub=sqrt(sum(S.^1)).\sqrt(sum(S(t:end).^1));
S(t:end)=0;
else % all zero entries
t=1;
energy_scrub=0;
sigmasq_2=0;
end
else
S(max(1,end-floor(idx*soft_thrs)):end)=0;
end
tmp1=U*diag(S)*V';
tmp1=reshape(tmp1,size(KSP2_tmp));
if patch_scale==1; else; patch_scale=size(S,1)-idx; end
if isempty(t); t=1; end % threshold removed all.
if patch_avg==1
KSP2_tmp_update(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) =...
KSP2_tmp_update(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) +patch_scale*tmp1;
KSP2_weight(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) =...
KSP2_weight(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) + patch_scale;
KSP2_tmp_update_threshold(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) =...
KSP2_tmp_update_threshold(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) +idx;
energy_removed(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) =...
energy_removed(:,[1:w2]+(n2-1),[1:w3]+(n3-1),:) +energy_scrub;
SNR_weight(:,[1:w2]+(n2-1),[1:w3]+(n3-1),1) =...
SNR_weight(:,[1:w2]+(n2-1),[1:w3]+(n3-1),1) + S(1)./S(max(1,t-1));
try
NOISE(1:size(KSP2a,1),[1:w2]+(n2-1),[1:w3]+(n3-1),1) = ...
NOISE(1:size(KSP2a,1),[1:w2]+(n2-1),[1:w3]+(n3-1),1) + sigmasq_2(t);; catch; end
else
KSP2_tmp_update(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) =...
KSP2_tmp_update(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) +patch_scale*tmp1(1,round(end/2),round(end/2),:);
KSP2_weight(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) =...
KSP2_weight(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) +patch_scale;
KSP2_tmp_update_threshold(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) =...
KSP2_tmp_update_threshold(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) +idx;
energy_removed(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) =...
energy_removed(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) +energy_scrub;
SNR_weight(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),1) =...
SNR_weight(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),1) + S(1)./S(max(1,t-1));
try
NOISE(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) = ...
NOISE(:,round(w2/2)+(n2-1),round(w3/2)+(n3-1),:) + sigmasq_2(t);; catch; end
end
% if MM0>1 & MM0<196
% [ sigmasq_2(t) MM0] % 2
% end
end
end
return

NIFTI_NORDIC_wt.m at commit 29af031, under MIT · at the source

Overview

Authors: Wei-Tang Chang^1,2,3, Weili Lin^1,2, Kelly S Giovanello^1,4

ORCID iDs: Wei-Tang Chang

Biomedical Research Imaging Center, University of North Carolina at Chapel Hill Chapel Hill United States
Department of Radiology, University of North Carolina at Chapel Hill Chapel Hill United States
Department of Biomedical Engineering, University of North Carolina at Chapel Hill Chapel Hill United States
Department of Psychology and Neuroscience, University of North Carolina at Chapel Hill Chapel Hill United States

Institutions: University of North Carolina at Chapel Hill (United States)

Journal: eLife, volume 12, article RP92805

Dates: published online 21 April 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.7554/elife.92805 · PMID 42012994 · PMCID PMC13099135 · OpenAlex W4389922539

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: fMRI (modality), human (organism), systems (subfield)

Methods: Spectral & time-frequency, Connectivity, Statistics, Smoothing, state filtering, decompositions, Preprocessing, fMRI & imaging

Keywords: brain, layer fMRI, connectivity, Human

MeSH: Brain Mapping*, Cerebral Cortex*, Magnetic Resonance Imaging*, Adult, Female, Humans, Male (* major topic)

Journal subjects: Neuroscience

Topic: Advanced MRI Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex

Funding: National Institutes of Health (R21AG060324)

Citations: not cited yet (Europe PMC); 76 references in the paper

Abstract

Layer-dependent functional magnetic resonance imaging (fMRI) is a promising yet challenging approach for investigating layer-specific functional connectivity (FC). Achieving a brain-wide mapping of layer-specific FC requires several technical advancements, including sub-millimeter spatial resolution, sufficient temporal resolution, functional sensitivity, global brain coverage, and high spatial specificity. Although gradient echo (GE)-based echo planar imaging (EPI) is commonly used for rapid fMRI acquisition, it faces significant challenges due to the draining-vein contamination. In this study, we addressed these limitations by integrating velocity-nulling (VN) gradients into a GE-BOLD fMRI sequence to suppress vascular signals from the vessels with fast-flowing velocity. The extravascular contamination from pial veins was mitigated using a GE-EPI sequence at 3T rather than 7T, combined with phase regression methods. Additionally, we incorporated advanced techniques, including simultaneous multi-slice (SMS) acceleration and NOise Reduction with DIstribution Corrected principal component analysis (NORDIC PCA) denoising, to improve temporal resolution, spatial coverage, and signal sensitivity. This resulted in a VN fMRI sequence with 0.9 mm isotropic spatial resolution, a repetition time (TR) of 4 s, and brain-wide coverage. The VN gradient strength was determined based on results from a button-pressing task. Using resting-state data, we validated layer-specific FC through seed-based analyses, identifying distinct connectivity patterns in the superficial and deep layers of the primary motor cortex (M1), with significant inter-layer differences. Further analyses with a seed in the primary sensory cortex (S1) demonstrated the reliability of the method. Brain-wide layer-dependent FC analyses yielded results consistent with prior literature, reinforcing the efficacy of VN fMRI in resolving layer-specific functional connectivity. Given the widespread availability of 3T scanners, this technical advancement has the potential for significant impact across multiple domains of neuroscience research.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

welton0411/matlab_sms_recon

License: MIT

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: 29af03114c6b85202e8a109e7107663ab2a89d3d, 20 April 2026

Languages: MATLAB (34)

Size: 36 files, 34 scripts

Software Heritage: archived

Found in: “Data availability”

Holds: README, license file

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Image Processing Toolbox (4 files), Tools for NIfTI and ANALYZE image (MATLAB) (3 files), Parallel Computing Toolbox (1 file), Signal Processing Toolbox (1 file), Statistics and Machine Learning Toolbox (1 file)

Availability: 1 check, the latest on 29 September 2026: the link answers

29 September 2026: the link answers

36 files

convertVE_s2025010601_run1.m, MATLAB, 189 lines
pproc2_s2025010601_run1_dn2.m, MATLAB, 1,399 lines, 1 match
recon_s2025010601_run1_dn2.m, MATLAB, 1,237 lines
welton_toolbox_github/NIFTI_NORDIC_wt.m, MATLAB, 1,148 lines, 1 match
welton_toolbox_github/NORDIC_wt.m, MATLAB, 1,220 lines
welton_toolbox_github/NORDIC_wt_ker11.m, MATLAB, 1,216 lines
welton_toolbox_github/NORDIC_wt_ker7.m, MATLAB, 1,216 lines
welton_toolbox_github/NORDIC_wt_ker9.m, MATLAB, 1,216 lines
welton_toolbox_github/coilCombine.m, MATLAB, 58 lines
welton_toolbox_github/coilCombine_wt2.m, MATLAB, 67 lines
welton_toolbox_github/datapack4_VE.m, MATLAB, 837 lines
welton_toolbox_github/datapack_VB.m, MATLAB, 396 lines
welton_toolbox_github/datapack_VE.m, MATLAB, 510 lines
welton_toolbox_github/deming.m, MATLAB, 188 lines
welton_toolbox_github/dice.m, MATLAB, 121 lines
welton_toolbox_github/ecalib_2d.m, MATLAB, 50 lines
welton_toolbox_github/ecalib_3d.m, MATLAB, 48 lines
welton_toolbox_github/fft2c.m, MATLAB, 7 lines
welton_toolbox_github/gen_mb3d_encodemtx_i2k.m, MATLAB, 203 lines, 1 match
welton_toolbox_github/gen_mussel_encodemtx.m, MATLAB, 166 lines
welton_toolbox_github/im2row.m, MATLAB, 13 lines
welton_toolbox_github/im2row_3d.m, MATLAB, 13 lines
welton_toolbox_github/jaccard.m, MATLAB, 170 lines
welton_toolbox_github/kernelEig.m, MATLAB, 67 lines
welton_toolbox_github/mapVBVD.m, MATLAB, 830 lines
welton_toolbox_github/mapVBVD_wt.m, MATLAB, 853 lines
welton_toolbox_github/mb3d_grappa_wt3.m, MATLAB, 732 lines, 1 match
welton_toolbox_github/mb3d_grappa_wt3p.m, MATLAB, 732 lines
welton_toolbox_github/noise_map_DWI_wt.m, MATLAB, 58 lines
welton_toolbox_github/pocs.m, MATLAB, 671 lines
welton_toolbox_github/pocs_wt.m, MATLAB, 693 lines
welton_toolbox_github/read_twix_hdr.m, MATLAB, 160 lines
welton_toolbox_github/twix_map_obj.m, MATLAB, 1,132 lines
welton_toolbox_github/zpad.m, MATLAB, 63 lines
LICENSE.txt, License, 21 lines
README.md, Text, 1 line

At the source: github.com/welton0411/matlab_sms_recon

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
34 scripts, each with its path and the digest of its content;
4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

openneuro:ds007543, at OpenNeuro; found in “Data availability”

Data availability

All imaging data has been deposited at OpenNeuro (https://openneuro.org/datasets/ds007543). Image reconstruction code is available at GitHub (https://github.com/welton0411/matlab_sms_recon) (copy archived at welton0411, 2026). Please note that the code provided is in its raw form; it has not been optimized, cleaned, or commented. While this may affect the ease of use or adaptation, we believe it remains a valuable resource for those interested in understanding or extending our analytical methods.

The following dataset was generated:

ChangWT LinW GiovanelloK OpenNeuro2026Open data of VNfMRIds007543

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 3 authors, 4 keywords, 7 MeSH terms, 1 funder, 74 references.

Cite

This paper

Chang, W.-T., Lin, W., & Giovanello, K. S. (2026). Brain-wide mapping of layer-specific functional connectivity in the human cortex at 3T using draining-vein-suppressed fMRI. eLife, 12, RP92805. https://doi.org/10.7554/elife.92805

BibTeX

@article{chang2026brain, author = {Chang, Wei-Tang and Lin, Weili and Giovanello, Kelly S}, title = {{Brain-wide mapping of layer-specific functional connectivity in the human cortex at 3T using draining-vein-suppressed fMRI}}, journal = {eLife}, year = {2026}, month = apr, volume = {12}, pages = {RP92805}, publisher = {eLife Sciences Publications, Ltd}, issn = {2050-084X}, doi = {10.7554/elife.92805}, url = {https://doi.org/10.7554/elife.92805}, pmid = {42012994}, pmcid = {PMC13099135} }

RIS

TY - JOUR AU - Chang, Wei-Tang AU - Lin, Weili AU - Giovanello, Kelly S TI - Brain-wide mapping of layer-specific functional connectivity in the human cortex at 3T using draining-vein-suppressed fMRI T2 - eLife J2 - eLife PY - 2026 DA - 2026/04/21 VL - 12 SP - RP92805 SN - 2050-084X PB - eLife Sciences Publications, Ltd DO - 10.7554/elife.92805 UR - https://doi.org/10.7554/elife.92805 LA - en ER -

CSL-JSON

{ "id": "10.7554/elife.92805", "type": "article-journal", "title": "Brain-wide mapping of layer-specific functional connectivity in the human cortex at 3T using draining-vein-suppressed fMRI", "container-title": "eLife", "author": [ { "family": "Chang", "given": "Wei-Tang" }, { "family": "Lin", "given": "Weili" }, { "family": "Giovanello", "given": "Kelly S" } ], "container-title-short": "eLife", "volume": "12", "page": "RP92805", "DOI": "10.7554/elife.92805", "PMID": "42012994", "PMCID": "PMC13099135", "ISSN": "2050-084X", "publisher": "eLife Sciences Publications, Ltd", "URL": "https://doi.org/10.7554/elife.92805", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 21 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1371/journal.pcbi.1014576 [code]: A synthetic 3D human cerebrovascular model informed by histology for simulating the cortical depth-dependent BOLD fMRI signal.
Journal: PLoS computational biology
In common: fMRI, 11 references
[2] doi:10.7554/elife.108408 [code]: Frequency and laminar profile of feature-specific visual activity revealed by interleaved EEG-fMRI.
Journal: eLife
In common: Tools for NIfTI and ANALYZE image (MATLAB), Image Processing Toolbox, Signal Processing Toolbox, 1 other tool, fMRI, systems, 7 references
[3] doi:10.1162/imag.a.1197 [code]: Blood volume-sensitive laminar fMRI with VASO in human hippocampus: Capabilities and biophysical challenges at clinical 7T scanners.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Tools for NIfTI and ANALYZE image (MATLAB), Image Processing Toolbox, fMRI, systems, 7 references
[4] doi:10.1162/imag.a.1340: Biophysical simulations of fMRI responses using realistic microvascular models: Insights into distinct hemodynamics in humans and mice.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: fMRI, 8 references
[5] doi:10.1162/imag.a.1247 [code]: Evaluating BOLD functional MRI biophysical simulation approaches: Impact of vascular geometry, magnetic field calculations, and water diffusion models.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: fMRI, 8 references
[6] doi:10.1038/s41467-026-74215-5 [code]: Multi-metric evaluations of acute psychedelic effects on fMRI brain entropy.
Journal: Nature communications
In common: Tools for NIfTI and ANALYZE image (MATLAB), Parallel Computing Toolbox, Image Processing Toolbox, 2 other tools, fMRI, 2 references
[7] doi:10.1038/s41467-026-71842-w [code]: Layer-specific attentional modulation in the human primary somatosensory cortex.
Journal: Nature communications
In common: Image Processing Toolbox, Statistics and Machine Learning Toolbox, 5 references
[8] doi:10.1038/s41467-026-73540-z [code]: Predictive acoustical processing in human cortical layers.
Journal: Nature communications
In common: Image Processing Toolbox, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 4 references
[9] doi:10.1038/s41598-026-58377-2 [code]: A novel distortion-matched anatomical imaging sequence for high-fidelity functional mapping in submillimeter-resolution fMRI.
Journal: Scientific reports
In common: Image Processing Toolbox, fMRI, 5 references
[10] doi:10.1162/imag.a.1279 [code]: Multimodal laminar characterization of visual areas along the cortical hierarchy.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: systems, 6 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.