Merge pull request #49 from kschan0214/dev1.2.1

Dev1.2.1

README.md

@@ -41,7 +41,11 @@ If you have a more general question regarding the usage of SEPIA and/or other QS
 
 For full update log, please visit https://sepia-documentation.readthedocs.io/en/latest/getting_started/Release-note.html.
 
-### 1.2 (current master)
+### 1.2.1 (current master)
+* Fix bug for data with odd-number matrix size
+* Fix bug for missing file when using R2* mapping with NLLS algorithm
+
+### 1.2 (current d2f54a3)
 * Support several deep learning based methods (BFRnet, xQSM, QSMnet+ and LP-CNN) on Linux
 * Support atlas-based subcortical structure segmentation (CIT168 Reinforcement learning atlas, MuSus-100 and AHEAD) on Linux and Mac
 * Integrate R2* mapping toolbox into SEPIA

addons/qsm/MRI_susceptibility_calculation/Wrapper_QSM_iterTik.m

@@ -41,8 +41,12 @@
 weights = get_variable_from_headerAndExtraData(headerAndExtraData,'weights', matrixSize); % headerAndExtraData.weights;
 if isempty(weights)
     magn  = get_variable_from_headerAndExtraData(headerAndExtraData,'magnitude', matrixSize);  % you can access the magnitude and/or other data from the 'headerAndExtraData' variable
-    weights = sum(magn.^2,4);
-    clear magn
+    if ~isempty(magn)
+        weights = sum(magn.^2,4);
+        clear magn
+    else
+        weights = mask; % no weight and magnitude input
+    end
 end
 % masking weights
 weights = weights.*mask;

misc/r2s_mapping/lsq/Signal_mGRE.m

@@ -0,0 +1,79 @@
+%% S = Signal_mGRE(m0,x2s,te,varargin)
+%
+% Description: Multiple exponential decay model for multiecho GRE
+%
+% Input
+% ----------
+%   m0      : T1weighted signal
+%   x2s     : either T2* or R2*
+%   te      : echo times
+%   ***flag***
+%   ----------
+%   'r2s'   : input R2* instead of T2*
+%   'freq'  : frequency shift of all components (optional)
+%
+% Output
+% ----------
+%   s       : T2*-weighted signal
+%
+% Kwok-shing Chan @ DCCN
+% [email protected]
+% Date created: 23 February 2017
+% Date last modified: 13 October 2017
+%
+function s = Signal_mGRE(s0,x2s,te,varargin)
+%% check and parse input
+[mode,freq,verbose] = parse_varargin_Signal_mGRE(varargin);
+
+% number of exponential componets
+ncomp = length(s0);
+
+% check if the size of m0 and x2s match
+if length(x2s)~=ncomp
+    error('Mismatch of m0 size of relaxation constant size');
+end
+% check if frequency ok
+if length(freq)~=ncomp && freq(1)~=0
+    error('Mismatch of frequency size of relaxation constant size');
+end
+if freq(1)==0 && length(freq)==1
+    freq = repmat(freq(1),1,ncomp);
+end
+% convert t2s to r2s
+if strcmpi(mode,'r2s')
+    t2s = 1./x2s;
+else
+    t2s = x2s;
+end
+if verbose
+    fprintf('%i exponential component(s)\n',ncomp);
+end
+
+%% Core
+s = zeros(1,length(te));
+for kt=1:length(te)
+    s(kt) = sum(s0(:).*exp(-te(kt)./t2s(:)).*exp(1i*2*pi*freq(:)*te(kt))) ;
+end
+
+end
+
+%% Parsing varargin
+function [mode,freq,verbose] = parse_varargin_Signal_mGRE(arg)
+mode = 't2s';
+freq = 0;
+verbose=false;
+for kvar = 1:length(arg)
+    if strcmpi(arg{kvar},'r2s')
+        mode = 'r2s';
+        continue
+    end
+    if strcmpi(arg{kvar},'freq')
+        freq = arg{kvar+1};
+        continue
+    end
+    if strcmpi(arg{kvar},'-v')
+        verbose = true;
+        continue
+    end
+end
+end

sepia.m

@@ -262,6 +262,7 @@ function PushbuttonStart_Callback(source,eventdata)
 fprintf(fid,'sepia_addpath;\n\n');
 
 fprintf(fid,'%% Input/Output filenames\n');
+fprintf(fid,'input = struct();\n');
 % input data
 if isstruct(input)
     fprintf(fid,'input(1).name = ''%s'' ;\n',input(1).name);

sepia_universal_variables.m

@@ -14,7 +14,7 @@
 % DO NOT change the order of the entities, add a new one at the end instead
 %
 %% Version
-SEPIA_version = 'v1.2';
+SEPIA_version = 'v1.2.1';
 
 %% PATH
 SEPIA_HOME = fileparts(mfilename('fullpath'));

wrapper/estimateTotalField.m

@@ -49,8 +49,8 @@
 end
 
 %% ensure all variables have the same data type
-fieldMap	= double(fieldMap);
-mask       	= double(mask);
+fieldMap	= double(zeropad_odd_dimension(fieldMap,'pre'));
+mask       	= double(zeropad_odd_dimension(mask,'pre'));
 
 disp('--------------------');
 disp('Total field recovery');
@@ -99,4 +99,12 @@
 
 N_std = real(N_std);
 
+% remove zero padding 
+totalField              = double(zeropad_odd_dimension(totalField,'post',matrixSize));
+N_std                   = double(zeropad_odd_dimension(N_std,'post',matrixSize));
+mask                    = double(zeropad_odd_dimension(mask,'post',matrixSize));
+if ~isempty(fieldmapUnwrapAllEchoes)
+    fieldmapUnwrapAllEchoes = double(zeropad_odd_dimension(fieldmapUnwrapAllEchoes,'post',matrixSize));
+end
+
 end