Skip to content

Commit

Permalink
Working implementation with multiple contrast output. Need to test wh…
Browse files Browse the repository at this point in the history
…ether batch system works as expected. Note that input images are no longer denoising method specific.
  • Loading branch information
ledwards committed Jun 16, 2024
1 parent 5fa808e commit 0ac7eb5
Show file tree
Hide file tree
Showing 4 changed files with 73 additions and 74 deletions.
29 changes: 15 additions & 14 deletions hmri_denoising.m
Original file line number Diff line number Diff line change
@@ -1,15 +1,14 @@
function varargout = hmri_denoising(job)
function varargout = hmri_denoising(jobsubj)

% retrieve effective acquisition & processing parameters
jobsubj = job.subj;
denoisedout = get_denoising_params(jobsubj);
denoising_params = get_denoising_params(jobsubj);
protocolfield = fieldnames(jobsubj.denoisingtype);
denoising_protocol = protocolfield{1};

% execute the chosen denoising method and define output
switch denoising_protocol
case 'lcpca_denoise'
[output_mag, output_phase] = hmri_calc_lcpcadenoise(denoisedout);
[output_mag, output_phase] = hmri_calc_lcpcadenoise(denoising_params);
varargout{1} = output_mag;
varargout{2} = output_phase;
end
Expand Down Expand Up @@ -61,20 +60,22 @@
[v,r] = spm('Ver');
denoising_params.SPMver = sprintf('%s (%s)', v, r);

% Load input data
denoising_params.mag_img = cellstr(char(spm_file(jobsubj.mag_img,'number','')));
denoising_params.phase_img = cellstr(char(spm_file(jobsubj.phase_img,'number','')));
denoising_params.phase_bool = any(~cellfun(@isempty, denoising_params.phase_img));
denoising_params.mag_bool = any(~cellfun(@isempty, denoising_params.mag_img));

% processing can continue if only magnitude images were entered but
% only warn that optional phase img are missing
if ~denoising_params.phase_bool || ~isfield(jobsubj,'phase_img')
hmri_log('Warning: No (optional) phase images were entered, denoising will continue with only magnitude images', denoising_params.defflags);
end

% Denoising method-specific parameters
% Load all the batch entered and possibly user-modified parameters
switch denoising_protocol
case 'lcpca_denoise'
denoising_params.mag_img = cellstr(char(spm_file(jobsubj.denoisingtype.(denoising_protocol).mag_img,'number','')));
denoising_params.phase_img = cellstr(char(spm_file(jobsubj.denoisingtype.(denoising_protocol).phase_img,'number','')));
denoising_params.phase_bool = any(~cellfun(@isempty, denoising_params.phase_img));
denoising_params.mag_bool = any(~cellfun(@isempty, denoising_params.mag_img));

% processing can continue if only magnitude images were entered but
% only warn that optional phase img are missing
if ~denoising_params.phase_bool || ~isfield(jobsubj.denoisingtype.(denoising_protocol),'phase_img')
hmri_log('Warning: No (optional) phase images were entered, Lcpca-denoising continues with only magnitude images', denoising_params.defflags);
end

dnstruct = jobsubj.denoisingtype.lcpca_denoise;
denoising_params.ngbsize = dnstruct.ngbsize;
Expand Down
77 changes: 34 additions & 43 deletions hmri_run_denoising.m
Original file line number Diff line number Diff line change
Expand Up @@ -8,17 +8,24 @@

end

function out = hmri_denoising_local(job)
function out = hmri_denoising_local(jobsubj)

% Get denoising protocol
dntype=fieldnames(job.subj.denoisingtype);
% concatenate all contrasts into jobsubj.mag_img and phase_img
contrasts = {'pdw','t1w','mtw'};
jobsubj.mag_img = {};
jobsubj.phase_img = {};
for c = 1:length(contrasts)
con = contrasts{c};
jobsubj.mag_img = [jobsubj.mag_img; jobsubj.(con).mag_img];
jobsubj.phase_img = [jobsubj.phase_img; jobsubj.(con).phase_img];
end

% Case indir versus outdir
try
outpath = job.subj.output.outdir{1}; % case outdir
outpath = jobsubj.output.outdir{1}; % case outdir
if ~exist(outpath,'dir'); mkdir(outpath); end
catch
Pin = char(job.subj.denoisingtype.(dntype{1}).mag_input);
Pin = char(jobsubj.mag_img);
outpath = fileparts(Pin(1,:)); % case indir
end

Expand Down Expand Up @@ -52,41 +59,41 @@
supplpath = fullfile(outpath, 'Results', 'Supplementary');
if ~exist(supplpath,'dir'); mkdir(supplpath); end

% save all these paths in the job.subj structure
job.subj.path.dnrespath = respath;
job.subj.path.respath = respath;
job.subj.path.supplpath = supplpath;
% save all these paths in the jobsubj structure
jobsubj.path.dnrespath = respath;
jobsubj.path.respath = respath;
jobsubj.path.supplpath = supplpath;

% save log file location
job.subj.log.logfile = fullfile(supplpath, 'hMRI_denoising_logfile.log');
job.subj.log.flags = struct('LogFile',struct('Enabled',true,'FileName','hMRI_denoising_logfile.log','LogDir',supplpath), ...
'PopUp',job.subj.popup,'ComWin',true);
flags = job.subj.log.flags;
jobsubj.log.logfile = fullfile(supplpath, 'hMRI_denoising_logfile.log');
jobsubj.log.flags = struct('LogFile',struct('Enabled',true,'FileName','hMRI_denoising_logfile.log','LogDir',supplpath), ...
'PopUp',jobsubj.popup,'ComWin',true);
flags = jobsubj.log.flags;
flags.PopUp = false;
hmri_log(sprintf('\t============ DENOISING MODULE - %s.m (%s) ============', mfilename, datetime('now')),flags);

if newrespath
hmri_log(sprintf(['WARNING: existing results from previous run(s) were found, \n' ...
'the output directory has been modified. It is now:\n%s\n'],outpath),job.subj.log.flags);
'the output directory has been modified. It is now:\n%s\n'],outpath),jobsubj.log.flags);
else
hmri_log(sprintf('INFO: the output directory is:\n%s\n',outpath),flags);
end

% check basic requirements and error if basic requirements fail
check_params.mag_input = cellstr(char(spm_file(job.subj.denoisingtype.(dntype{1}).mag_input,'number','')));
check_params.phase_input = cellstr(char(spm_file(job.subj.denoisingtype.(dntype{1}).phase_input,'number','')));
check_params.phase_bool = any(~cellfun(@isempty, check_params.phase_input));
check_params.mag_bool = any(~cellfun(@isempty, check_params.mag_input));
check_params.mag_img = cellstr(char(spm_file(jobsubj.mag_img,'number','')));
check_params.mag_bool = any(~cellfun(@isempty, check_params.mag_img));
check_params.phase_img = cellstr(char(spm_file(jobsubj.phase_img,'number','')));
check_params.phase_bool = any(~cellfun(@isempty, check_params.phase_img));

% Issue an error and abort in cases of non-existent magnitude image data and
% non-equal number of non-empty phase and magnitude images
if ~check_params.mag_bool || ~isfield(job.subj.denoisingtype.(dntype{1}),'mag_input')
if ~check_params.mag_bool || ~isfield(jobsubj,'mag_img')
msg = 'No magnitude images were entered, aborting! Please check your input data and try again!';
hmri_log(msg, flags);
error(msg)
end

if check_params.phase_bool && (length(check_params.mag_input) ~= length(check_params.phase_input))
if check_params.phase_bool && (length(check_params.mag_img) ~= length(check_params.phase_img))
msg = 'The number of phase and magnitude images are different, please check your input data and try again!';
hmri_log(msg, flags);
error(msg)
Expand All @@ -95,47 +102,31 @@
% save SPM version (slight differences may appear in the results depending
% on the SPM version!)
[v,r] = spm('Ver');
job.SPMver = sprintf('%s (%s)', v, r);
jobsubj.SPMver = sprintf('%s (%s)', v, r);

% save original job to supplementary directory
spm_jsonwrite(fullfile(supplpath,'hMRI_denoising_job.json'),job,struct('indent','\t'));
spm_jsonwrite(fullfile(supplpath,'hMRI_denoising_job.json'),jobsubj,struct('indent','\t'));

% run the denoising and collect the results
hmri_log(sprintf('\t============ %s - %s.m (%s) ============',"APPLYING DENOISING", mfilename, datetime('now')),flags);

% concatenate all contrasts into jobsubj.denoisingtype.(denoising_protocol).mag_img and phase_img
contrasts = {'mtw','pdw','t1w'};
jobsubj.denoisingtype.(dntype{1}).mag_img = {};
jobsubj.denoisingtype.(dntype{1}).phase_img = {};
for c = 1:length(contrasts)
con = contrasts{c};

jobsubj.denoisingtype.(dntype{1}).mag_img =
[jobsubj.denoisingtype.(dntype{1}).mag_img;
jobsubj.denoisingtype.(dntype{1}).(con).mag_img];

jobsubj.denoisingtype.(dntype{1}).phase_img =
[jobsubj.denoisingtype.(dntype{1}).phase_img;
jobsubj.denoisingtype.(dntype{1}).(con).phase_img];
end

% run the denoising
[output_mag, output_phase] = hmri_denoising(job);
[output_mag, output_phase] = hmri_denoising(jobsubj);

% assign output dependencies to denoised output images of separate contrasts
iMag = 1;
iPhase = 1;
for c = 1:length(contrasts)
con = contrasts{c};

nMag = length(jobsubj.denoisingtype.(dntype{1}).(con).mag_img);
nMag = length(jobsubj.(con).mag_img);
idxstr = ['DenoisedMagnitude' con];
out.subj.(idxstr) = output_mag(iMag:iMag+nMag);
out.(idxstr) = output_mag(iMag:iMag+nMag-1);
iMag = iMag + nMag;

nPhase = length(jobsubj.denoisingtype.(dntype{1}).(con).phase_img);
nPhase = length(jobsubj.(con).phase_img);
idxstr = ['DenoisedPhase' con];
out.subj.(idxstr) = output_phase(iPhase:iPhase+nPhase);
out.(idxstr) = output_phase(iPhase:iPhase+nPhase-1);
iPhase = iPhase + nPhase;
end

Expand Down
4 changes: 2 additions & 2 deletions tbx_cfg_hmri.m
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
% Bogdan Draganski & Ferath Kherif, 2011
% ======================================================================

if ~isdeployed,
if ~isdeployed
hMRIpath = fileparts(mfilename('fullpath'));
addpath(hMRIpath);
addpath(fullfile(hMRIpath, 'config'));
Expand Down Expand Up @@ -43,6 +43,6 @@
'and will include a number of (as yet unspecified) extensions in ',...
'future updates. Please report any bugs or problems to [email protected].']
}';
hmri.values = {tbx_scfg_hmri_config tbx_scfg_hmri_dicom_import tbx_scfg_hmri_autoreorient tbx_scfg_hmri_imperf_spoil tbx_scfg_hmri_B1_create tbx_scfg_hmri_create tbx_scfg_hmri_quality tbx_scfg_hmri_proc tbx_scfg_hmri_QUIQI tbx_scfg_hmri_denoising};
hmri.values = {tbx_scfg_hmri_config tbx_scfg_hmri_dicom_import tbx_scfg_hmri_denoising tbx_scfg_hmri_autoreorient tbx_scfg_hmri_imperf_spoil tbx_scfg_hmri_B1_create tbx_scfg_hmri_create tbx_scfg_hmri_quality tbx_scfg_hmri_proc tbx_scfg_hmri_QUIQI};
end

37 changes: 22 additions & 15 deletions tbx_scfg_hmri_denoising.m
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,8 @@
mag_img.help = {'Select the (required) magnitude images to be denoised'};
mag_img.filter = 'image';
mag_img.ufilter = '.*';
mag_img.num = [1 Inf];
mag_img.num = [0 Inf];
mag_img.val = {''};

% ---------------------------------------------------------------------
% Phase input images
Expand All @@ -99,22 +100,28 @@
% ---------------------------------------------------------------------
mtw = cfg_branch;
mtw.tag = 'mtw';
mtw.name = 'MTw input';
mtw.name = 'MTw input (optional)';
mtw.help = {'Input Magnitude/Phase images from MTw data'};
mtw.val = {mag_img phase_img};

pdw = cfg_branch;
pdw.tag = 't1w';
pdw.name = 'T1w input';
pdw.help = {'Input Magnitude/Phase images from PDw data'};
pdw.val = {mag_img phase_img};

t1w = cfg_branch;
t1w.tag = 'pdw';
t1w.name = 'PDw input';
t1w.tag = 't1w';
t1w.name = 'T1w input (optional)';
t1w.help = {'Input Magnitude/Phase images from T1w data'};
t1w.val = {mag_img phase_img};

% PDw is required input
pdw_mag_img = mag_img;
pdw_mag_img.num = [1 Inf];
pdw_mag_img.val = {};

pdw = cfg_branch;
pdw.tag = 'pdw';
pdw.name = 'PDw input';
pdw.help = {'Input Magnitude/Phase images from PDw data', ...
'If you only have one kind of weighting, please put them here.'};
pdw.val = {pdw_mag_img phase_img};

% ---------------------------------------------------------------------
% Standard deviation parameter
% ---------------------------------------------------------------------
Expand Down Expand Up @@ -152,7 +159,7 @@
denoisinginput_lcpca.help = {'Input Magnitude/Phase images for Lcpca-denoising'
['Regarding processing parameters, you can either stick with metadata and standard ' ...
'defaults parameters (recommended) or select your own hmri_denoising_local_defaults_*.m customised defaults file.']};
denoisinginput_lcpca.val = {mtw pdw t1w DNparameters std ngbsize};
denoisinginput_lcpca.val = {DNparameters std ngbsize};

% ---------------------------------------------------------------------
% menu denoisingtype
Expand Down Expand Up @@ -211,7 +218,7 @@
subj.tag = 'subj';
subj.name = 'Subject';
subj.help = {'Specify a subject for denoising.'};
subj.val = {output denoisingtype popup};
subj.val = {output pdw t1w mtw denoisingtype popup};

% ---------------------------------------------------------------------
% data Data
Expand Down Expand Up @@ -251,7 +258,7 @@

% iterate to generate dependency tags for outputs
nsub = numel(job.subj);
contrasts = {'mtw','pdw','t1w'};
contrasts = {'pdw','t1w','mtw'};
ncon = length(contrasts);
dep = repmat(cfg_dep,1,2*ncon*nsub);
for i=1:nsub
Expand All @@ -261,15 +268,15 @@

% define variables and initialize cfg_dep for magnitude images
cdep = cfg_dep;
cdep.sname = sprintf('lcpcaDenoised_%s_magnitude',con);
cdep.sname = sprintf('Denoised_%s_magnitude',con);
idxstr = ['DenoisedMagnitude' con];
cdep.src_output = substruct('.','subj','()',{i},'.',idxstr,'()',{':'});
cdep.tgt_spec = cfg_findspec({{'filter','image','strtype','e'}});
dep((i-1)*ncon+k) = cdep;

% define variables and initialize cfg_dep for phase images
cdep = cfg_dep;
cdep.sname = sprintf('lcpcaDenoised_%s_phase',con);
cdep.sname = sprintf('Denoised_%s_phase',con);
idxstr = ['DenoisedPhase' con];
cdep.src_output = substruct('.','subj','()',{i},'.',idxstr,'()',{':'});
cdep.tgt_spec = cfg_findspec({{'filter','image','strtype','e'}});
Expand Down

0 comments on commit 0ac7eb5

Please sign in to comment.