diff --git a/CHANGELOG.md b/CHANGELOG.md index 96c72fa..01fbc7f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ Most recent version numbers *should* follow the [Semantic Versioning](https://se - set default WM percent value in hmri_defaults. - spatial processing: add explicit mask creation and fix implicit mask (0 to NaN in float images) - update FIL seste seq parameters in get_metadata_val_classic +- denoising module-first part: Java-Matlab interface for LCPCA denoising ### Fixed - replace `datestr(now)` with `datetime('now')` in line with [MATLAB recommendation](https://mathworks.com/help/matlab/matlab_prog/replace-discouraged-instances-of-serial-date-numbers-and-date-strings.html) @@ -17,6 +18,7 @@ Most recent version numbers *should* follow the [Semantic Versioning](https://se - make B1-map creation using 3DEPI SE/STE and AFI methods fall back to defaults without sidecar files, rather than crash - Modify the filenames as files are copied to RFsensCalc to prevent overwriting in further processing - batch interface now enforces the number of B1 input images correctly for B1 mapping methods which only need two images. +- fix error if optimization toolbox not present during NLLS R2* calculation ## [v0.6.1] ### Fixed diff --git a/config/hmri_denoising_defaults.m b/config/hmri_denoising_defaults.m new file mode 100644 index 0000000..3196a23 --- /dev/null +++ b/config/hmri_denoising_defaults.m @@ -0,0 +1,19 @@ +function hmri_denoising_defaults + +%init the global variable which carries the params +global hmri_def; + +%Enter denoising default values as: hmri_def.denoising.(denoising-protocol).(default-value) + +%The default values for lcpca denoising protocol +%all optional parameters turned off +hmri_def.denoising.lcpca_denoise.mag_input = {}; %%required-null initialize here input with GUI +hmri_def.denoising.lcpca_denoise.phase_input = {}; %%optional-null initialize here input with GUI +hmri_def.denoising.lcpca_denoise.output_path=''; %%required-null initialize here input with GUI +hmri_def.denoising.lcpca_denoise.min_dimension= 0; %%required-initialize here +hmri_def.denoising.lcpca_denoise.max_dimension = -1; %%required-initialize here +hmri_def.denoising.lcpca_denoise.unwrap = false; %%optional +hmri_def.denoising.lcpca_denoise.rescale_phs = false; %%optional +hmri_def.denoising.lcpca_denoise.process_2d=false; %%optional +hmri_def.denoising.lcpca_denoise.use_rmt=false; %%optional +end \ No newline at end of file diff --git a/config/local/hmri_local_denoising_defaults.m b/config/local/hmri_local_denoising_defaults.m new file mode 100644 index 0000000..2b7c23f --- /dev/null +++ b/config/local/hmri_local_denoising_defaults.m @@ -0,0 +1,16 @@ +%init the global variable which carries the params +global hmri_def; + +%Enter denoising default values as: hmri_def.denoising.(denoising-protocol).(default-value) + +%The default values for lcpca denoising protocol +%all optional parameters turned off +hmri_def.denoising.lcpca_denoise.mag_input = {}; %%required-null initialize here input with GUI +hmri_def.denoising.lcpca_denoise.phase_input = {}; %%optional-null initialize here input with GUI +hmri_def.denoising.lcpca_denoise.output_path=''; %%required-null initialize here input with GUI +hmri_def.denoising.lcpca_denoise.min_dimension= 0; %%required-initialize here +hmri_def.denoising.lcpca_denoise.max_dimension = -1; %%required-initialize here +hmri_def.denoising.lcpca_denoise.unwrap = false; %%optional +hmri_def.denoising.lcpca_denoise.rescale_phs = false; %%optional +hmri_def.denoising.lcpca_denoise.process_2d=false; %%optional +hmri_def.denoising.lcpca_denoise.use_rmt=false; %%optional diff --git a/hmri_calc_R2s.m b/hmri_calc_R2s.m index 5acb6b2..25489f1 100644 --- a/hmri_calc_R2s.m +++ b/hmri_calc_R2s.m @@ -147,6 +147,18 @@ beta(2:end,:)=exp(beta(2:end,:)); case {'nlls_ols','nlls_wls1','nlls_wls2','nlls_wls3'} + %lsqcurvefit uses optimization toolbox + %Check both for the toolbox and/or an active license for it + %Error if one of them is missing + versionStruct = ver; + versionCell = {versionStruct.Name}; + ver_status = any(ismember(versionCell, 'Optimization Toolbox')); + [license_status,~] = license('checkout', 'Optimization_toolbox'); + if ver_status==0 || license_status==0 + error('hmri:NoOptimToolbox', "The methods 'nlls_ols','nlls_wls1','nlls_wls2','nlls_wls3' require Optimization Toolbox: either this toolbox and/or its license is missing." + ... + " Please use another method which does not need the Optimization Toolbox such as 'ols','wls1','wls2','wls3'. ") + end + % Check for NLLS case, where specification of the log-linear % initialisation method is in the method string following a hyphen r=regexp(lower(method),'^nlls_(.*)$','tokens'); @@ -167,7 +179,7 @@ end expDecay=@(x,D) (D(:,2:end)*x(2:end)).*exp(x(1)*D(:,1)); - + % Loop over voxels parfor n=1:size(y,2) beta(:,n)=lsqcurvefit(@(x,D)expDecay(x,D),beta0(:,n),D,y(:,n),[],[],opt); diff --git a/hmri_create_denoising.m b/hmri_create_denoising.m new file mode 100644 index 0000000..9be9928 --- /dev/null +++ b/hmri_create_denoising.m @@ -0,0 +1,371 @@ +function varargout = hmri_create_denoising(job) + +% retrieve effective acquisition & processing parameters +jobsubj = job.subj; +denoisedout = 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); +varargout{1} = output_mag; +varargout{2} = output_phase; +end +end + +%=========================================================================% +% Write/Arrange denoising parameters tp the denoising_params. +%=========================================================================% +function denoising_params = get_denoising_params(jobsubj) + +% get denoising method +dntypename = fieldnames(jobsubj.denoisingtype); +denoising_protocol = dntypename{1}; + +% init defaults filename and optional-defaults bool +deffnam = ''; +custom_def = false; + +% load customized defaults parameters from customized denoising defaults file if any +% (the customized defaults file must be run to overwrite the standard +% defaults parameters) +if isfield(jobsubj.denoisingtype.(denoising_protocol),'DNparameters') + % first reinitialise processing parameters to standard defaults: + hmri_denoising_defaults; + deffnam = fullfile(fileparts(mfilename('fullpath')),'config','hmri_denoising_defaults.m'); + custom_def = false; + + % then, if customized defaults file available, run it to overwrite + % standard defaults parameters. + if isfield(jobsubj.denoisingtype.(denoising_protocol).DNparameters,'DNdefaults') + deffnam = jobsubj.denoisingtype.(denoising_protocol).DNparameters.DNdefaults; + spm('Run',deffnam); + custom_def = true; + end +end + +% set all denoising defaults and set the defaults file to be true +denoising_params = hmri_get_defaults(['denoising.' denoising_protocol]); +denoising_params.defaults_file = deffnam; +denoising_params.custom_defaults = custom_def; + +% flags for logging information and warnings +denoising_params.defflags = jobsubj.log.flags; % default flags +denoising_params.nopuflags = jobsubj.log.flags; % force no Pop-Up +denoising_params.nopuflags.PopUp = false; + +% save SPM version (slight differences may appear in the results depending +% on the SPM version!) +[v,r] = spm('Ver'); +denoising_params.SPMver = sprintf('%s (%s)', v, r); + +%Denoising method-specific parameters +%Load all the batch entered and possibly user-modified parameters +switch denoising_protocol + case 'lcpca_denoise' + denoising_params.mag_input = cellstr(char(spm_file(jobsubj.denoisingtype.(denoising_protocol).mag_input,'number',''))); + denoising_params.phase_input = cellstr(char(spm_file(jobsubj.denoisingtype.(denoising_protocol).phase_input,'number',''))); + denoising_params.phase_bool = any(~cellfun(@isempty, denoising_params.phase_input)); + denoising_params.mag_bool = any(~cellfun(@isempty, denoising_params.mag_input)); + + + %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_input') + hmri_log('No (optional) phase images were entered, Lcpca-denoising continues with only magnitude images', denoising_params.nopuflags); + warning('No (optional) phase images were entered, Lcpca-denoising continues with only magnitude images') + end + + + dnstruct = jobsubj.denoisingtype.lcpca_denoise; + denoising_params.ngbsize =dnstruct.ngbsize; + denoising_params.std =dnstruct.std; + denoising_params.output_path = jobsubj.path.dnrespath; + denoising_params.supp_path = jobsubj.path.supplpath; + + %Print lcpca denoising parameters ngbsize and std cut off + print_lcpca_params.Neighborhood_Size = denoising_params.ngbsize; + print_lcpca_params.Standard_Deviation_Cuttoff = denoising_params.std; + printdnstruct = printstruct(print_lcpca_params); + hmri_log(sprintf('Lcpca Denoising Parameters:\n\n%s', ... + printdnstruct),denoising_params.nopuflags); + +end + +end + +%===============================================================================================% +% Calculate LCPCA-denoising: Java-Matlab interface for Lcpca-denoising +%=================================================================================================% +function [output_mag, output_phase] = hmri_calc_lcpcadenoise(lcpcadenoiseparams) + +%{ +Paper Reference for Lcpca denoising (Java module originally written by Pilou Bazin): + [1] Bazin, et al. (2019) "Denoising High-Field Multi-Dimensional MRI With Local +Complex PCA", Front. Neurosci. doi:10.3389/fnins.2019.01066 +%} + + +%get the path to the jar files +[mfilepath, ~,~] = fileparts(mfilename("fullpath")); +jarcommons = fullfile(fullfile(mfilepath,'jar'),'commons-math3-3.5.jar'); +jarmipav = fullfile(fullfile(mfilepath,'jar'),'Jama-mipav.jar'); +jarlcpca = fullfile(fullfile(mfilepath,'jar'),'lcpca2.jar'); + +%add required .jar files to the path +javaaddpath(jarcommons) +javaaddpath(jarmipav) +javaaddpath(jarlcpca) + +%Read from the input the processing parameters +image_list = cellstr(lcpcadenoiseparams.mag_input); +phase_list = cellstr(lcpcadenoiseparams.phase_input); +ngb_size = lcpcadenoiseparams.ngbsize; +stdev_cutoff = lcpcadenoiseparams.std; +min_dimension = lcpcadenoiseparams.min_dimension; +max_dimension = lcpcadenoiseparams.max_dimension; +unwrap = lcpcadenoiseparams.unwrap; +rescale_phs = lcpcadenoiseparams.rescale_phs; +process_2d = lcpcadenoiseparams.process_2d; +use_rmt = lcpcadenoiseparams.use_rmt; +output_path = cellstr(lcpcadenoiseparams.output_path); +supp_path = cellstr(lcpcadenoiseparams.supp_path); +lcpcaflags = lcpcadenoiseparams.defflags; + +%Get the full input file list +phase_bool = any(~cellfun(@isempty,phase_list)); +if phase_bool + fullim_list = [image_list; phase_list]; +else + fullim_list = image_list; +end + +%set the metadata mod +json = hmri_get_defaults('json'); + +%Create denoising object +noiseObj= javaObject('nl.uva.imcn.algorithms.LocalComplexPCADenoising'); + +%Get the image_number and dimensions from the first image of magnitude +%images +imdatavol = spm_vol(image_list{1}); +resolutions = sqrt(sum(imdatavol.mat(1:3,1:3).^2)); +imdatamatrix = spm_read_vols(imdatavol); +dimensions = size(imdatamatrix); +image_number = length(image_list)+1; + +%Set the image_number, dimensions and resolution for processing +noiseObj.setImageNumber(image_number); %number of echoes + +%Set dimensions based on whether it is 3D or 4D data +%Throw an error if neither +if length(dimensions) == 3 +noiseObj.setDimensions(dimensions(1),dimensions(2),dimensions(3)); +elseif length(dimensions) == 4 +noiseObj.setDimensions(dimensions(1),dimensions(2),dimensions(3), dimensions(4)); +else + hmri_log('The raw data does not have the correct dimensions (must be 3D or 4D data) for processing: please check your data!', lcpcaflags); + error('The raw data does not have the correct dimensions (must be 3D or 4D data) for processing: please check your data!') +end +noiseObj.setResolutions(resolutions(1),resolutions(2), resolutions(3)); + +%Loop through all reshaped echos +for echo = 1:length(image_list) + %read vol from filepath + datavol = spm_vol(image_list{echo}); + datamatrix = spm_read_vols(datavol); +%place vol for denoising +noiseObj.setMagnitudeImageAt(echo, reshape(datamatrix, [1 prod(size(datamatrix))])); + +%Add phase images if they exist +if phase_bool + %read vol from filepath + datavol = spm_vol(phase_list{echo}); + datamatrix = spm_read_vols(datavol); +%place vol for denoising +noiseObj.setPhaseImageAt(echo, reshape(datamatrix, [1 prod(size(datamatrix))])); +end +end + +%set all other params +noiseObj.setPatchSize(ngb_size) +noiseObj.setStdevCutoff(stdev_cutoff) +noiseObj.setMinimumDimension(min_dimension) +noiseObj.setMaximumDimension(max_dimension) +noiseObj.setUnwrapPhase(unwrap) +noiseObj.setRescalePhase(rescale_phs) +noiseObj.setProcessSlabIn2D(process_2d) +noiseObj.setRandomMatrixTheory(use_rmt) + +%With all the input parameters execute +%Catch error if the java script returns an error +lcpcaflags_nopopup = lcpcaflags; +lcpcaflags_nopopup.PopUp = false; +hmri_log(sprintf('Executing Lcpca-denoising (Java) module \n'), lcpcaflags_nopopup); +try +noiseObj.execute; +catch ME +hmri_log('There was an error in the Java code: could not execute Lcpca denosing!', lcpcaflags_nopopup); + disp(ME.message) + rethrow(ME) +end +hmri_log('Lcpca-denoising (Java) module succesfully executed', lcpcaflags_nopopup); + + +%initialize the cells to be populated with fullpaths of output +%magnitude and phase images +out_mag = {}; +idx_mag = 1; +out_phase = {}; +idx_phase =1; + +%Get the results for all echos and reshape +for echo = 1:length(image_list) + datamatrix = noiseObj.getDenoisedMagnitudeImageAt(echo); + volumedata = reshape(datamatrix, dimensions); + + %Write the volume to .nii with an update to standard .nii header + firstfile = image_list{echo}; + filehdr = spm_vol(image_list{echo}); + + [path,filename,ext] = fileparts(firstfile); + [~,mainfilename,~] = fileparts(firstfile); + filename = strcat('LcpcaDenoised_',filename,'.nii'); + + outfname = fullfile(output_path{1}, filename); + filehdr.fname = outfname; + filehdr.descrip = strcat(filehdr.descrip, ' + lcpca denoised'); + spm_write_vol(filehdr, volumedata); + + %write metadata as extended header and sidecar json + Output_hdr = init_dn_output_metadata(fullim_list, lcpcadenoiseparams); + Output_hdr.history.procstep.descrip = [Output_hdr.history.procstep.descrip ' (LCPCA)']; + Output_hdr.history.output.imtype = 'Denoising (LCPCA)'; + %add acquisition data if available (otherwise fields will be empty) + jsonfilename = fullfile(path,strcat(mainfilename,'.json')); + if exist(jsonfilename, 'file') ==2 + try + jsondata = spm_jsonread(jsonfilename); + data_RepetitionTime = get_metadata_val(jsondata,'RepetitionTime'); + data_EchoTime = get_metadata_val(jsondata,'EchoTime' ); + data_FlipAngle = get_metadata_val(jsondata, 'FlipAngle'); + Output_hdr.acqpar = struct('RepetitionTime',data_RepetitionTime, ... + 'EchoTime',data_EchoTime,'FlipAngle',data_FlipAngle); + catch + hmri_log('Although json sidecar file were found, the writing of acquisition metadata failed', lcpcaflags_nopopup); + + end + else + hmri_log('No json sidecar file were found, skipping the writing of acquisition metadata', lcpcaflags_nopopup); + end + + %Set all the metadata + set_metadata(outfname,Output_hdr,json); + + %add image to the output list + out_mag{idx_mag} = outfname; + idx_mag = idx_mag +1; + + + %Get also the phase images if applicable + + if phase_bool + datamatrix = noiseObj.getDenoisedPhaseImageAt(echo); + volumedata = reshape(datamatrix, dimensions); + + %Write the volume to .nii with an update to standard .nii header + firstfile = phase_list{echo}; + filehdr = spm_vol(phase_list{echo}); + + [path,filename,ext] = fileparts(firstfile); + [~,mainfilename,~] = fileparts(firstfile); + filename = strcat('LcpcaDenoised_',filename,'.nii'); + + outfname = fullfile(output_path{1}, filename); + filehdr.fname = outfname; + filehdr.descrip = strcat(filehdr.descrip, ' + lcpca denoised'); + spm_write_vol(filehdr, volumedata); + + %write metadata as extended header and sidecar json + Output_hdr = init_dn_output_metadata(fullim_list, lcpcadenoiseparams); + Output_hdr.history.procstep.descrip = [Output_hdr.history.procstep.descrip ' (LCPCA)']; + Output_hdr.history.output.imtype = 'Denoising (LCPCA)'; + %add acquisition data if available (otherwise fields will be empty) + jsonfilename = fullfile(path,strcat(mainfilename,'.json')); + if exist(jsonfilename, 'file') ==2 + try + jsondata = spm_jsonread(jsonfilename); + data_RepetitionTime = get_metadata_val(jsondata,'RepetitionTime'); + data_EchoTime = get_metadata_val(jsondata,'EchoTime' ); + data_FlipAngle = get_metadata_val(jsondata, 'FlipAngle'); + Output_hdr.acqpar = struct('RepetitionTime',data_RepetitionTime, ... + 'EchoTime',data_EchoTime,'FlipAngle',data_FlipAngle); + catch + hmri_log('Although json sidecar file were found, the writing of acquisition metadata failed', lcpcaflags_nopopup); + + end + else + hmri_log('No json sidecar file were found, skipping the writing of acquisition metadata', lcpcaflags_nopopup); + end + + %Set all the metadata + set_metadata(outfname,Output_hdr,json); + + %add image to the output list + out_phase{idx_phase} = outfname; + idx_phase = idx_phase +1; + + end + +end + +%Take out computed outputs +output_mag = out_mag; +output_phase = out_phase; + +%save estimated local dimensions and residuals (between input and denoised images) +dim_img = reshape(noiseObj.getLocalDimensionImage(), dimensions); +save(fullfile(supp_path{1}, 'dim_img.mat'), 'dim_img') +err_img = reshape(noiseObj.getNoiseFitImage(), dimensions); +save(fullfile(supp_path{1}, 'err_img.mat'), 'err_img') + +%Clear object and remove .jar from path properly +clear("noiseObj") +javarmpath(jarcommons) +javarmpath(jarmipav) +javarmpath(jarlcpca) + + +end + +%=========================================================================% +% To arrange the metadata structure for denoising output. +%=========================================================================% +function metastruc = init_dn_output_metadata(input_files, denoising_params) + +proc.descrip = ['hMRI toolbox - ' mfilename '.m - Denoising']; +proc.version = hmri_get_version; +proc.params = denoising_params; + +output.imtype = 'Denoised Image'; +output.units = ''; + +metastruc = init_output_metadata_structure(input_files, proc, output); + +end + +%=========================================================================% +% To print a structure into text - assumes simple structure (no +% sub-structure in it at this point). +%=========================================================================% +function s = printstruct(struc) + +s = ''; +fntmp = fieldnames(struc); +for cf = 1:length(fntmp) + s = sprintf('%s %16s: %s\n', s, fntmp{cf}, num2str(struc.(fntmp{cf}))); +end +end \ No newline at end of file diff --git a/hmri_run_denoise.m b/hmri_run_denoise.m new file mode 100644 index 0000000..4458f64 --- /dev/null +++ b/hmri_run_denoise.m @@ -0,0 +1,118 @@ +function out = hmri_run_denoise(job) + +%Get denoising protocol +dntype=fields(job.subj.denoisingtype); + +%Case indir versus outdir +try + outpath = job.subj.output.outdir{1}; % case outdir + if ~exist(outpath,'dir'); mkdir(outpath); end +catch + Pin = char(job.subj.denoisingtype.(dntype{1}).mag_input); + outpath = fileparts(Pin(1,:)); % case indir +end + +% save outpath as default for this job +hmri_get_defaults('outdir',outpath); + +% Directory structure for results: +% /Results +% /Results/Supplementary +% If repeated runs, is replaced by /Run_xx to avoid overwriting previous outputs. + +% define a directory for final results +% RESULTS contains the final denoised maps +% Supplementary containes additonal info and log +respath = fullfile(outpath, 'Results'); +newrespath = false; +if exist(respath,'dir') + index = 1; + tmpoutpath = outpath; + while exist(tmpoutpath,'dir') + index = index + 1; + tmpoutpath = fullfile(outpath,sprintf('Run_%0.2d',index)); + end + outpath = tmpoutpath; + mkdir(outpath); + respath = fullfile(outpath, 'Results'); + newrespath = true; +end +if ~exist(respath,'dir'); mkdir(respath); end +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 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; +flags.PopUp = false; +hmri_log(sprintf('\t============ DENOISING MODULE - %s.m (%s) ============', mfilename, datestr(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); +else + hmri_log(sprintf('INFO: the output directory is:\n%s\n',outpath),flags); +end + +%check basic requirements and error if basic reqirements fail +switch dntype{1} + case 'lcpca_denoise' + 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)); + + %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') + hmri_log('No magnitude images were entered, aborting! Please check your input data and try again!', flags); + error('Error: No magnitude images were entered, aborting! Please check your input data and try again!') + end + + if check_params.phase_bool && (length(check_params.mag_input) ~= length(check_params.phase_input)) + hmri_log('The number of phase and magnitude images are different, please check your input data and try again!', flags); + error("Error: The number of phase and magnitude images are different, please check your input data and try again!") + end +end + +% 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); + +% save original job to supplementary directory +spm_jsonwrite(fullfile(supplpath,'hMRI_denoising_job.json'),job,struct('indent','\t')); + +% run the denoising and collect the results +switch dntype{1} + case 'lcpca_denoise' +%Write to log and command window the specific denosing being applied +hmri_log(sprintf('\t============ %s - %s.m (%s) ============',"APPLYING LCPCA-DENOISING", mfilename, datestr(now)),flags); +[output_mag, output_phase] = hmri_create_denoising(job); + +%Assign output dependency to the denoised output images +phase_bool = any(~cellfun(@isempty,output_phase)); +arrayLength = length(output_mag); +for i = 1:2*arrayLength + if i<=arrayLength + idxstr = ['DenoisedMagnitude' int2str(i)]; + out.subj.(idxstr) = {output_mag{i}}; + elseif phase_bool && i>arrayLength + idxstr = ['DenoisedPhase' int2str(i-arrayLength)]; + out.subj.(idxstr) = {output_phase{i-arrayLength}}; + else + break + end +end +end +hmri_log(sprintf('\t============ DENOISING MODULE: completed (%s) ============', datestr(now)),flags); +end \ No newline at end of file diff --git a/jar/Jama-mipav.jar b/jar/Jama-mipav.jar new file mode 100644 index 0000000..76c5eb3 Binary files /dev/null and b/jar/Jama-mipav.jar differ diff --git a/jar/commons-math3-3.5.jar b/jar/commons-math3-3.5.jar new file mode 100644 index 0000000..db99f8c Binary files /dev/null and b/jar/commons-math3-3.5.jar differ diff --git a/jar/lcpca2.jar b/jar/lcpca2.jar new file mode 100644 index 0000000..591dd0e Binary files /dev/null and b/jar/lcpca2.jar differ diff --git a/tbx_cfg_hmri.m b/tbx_cfg_hmri.m index 6367727..ee84bc1 100644 --- a/tbx_cfg_hmri.m +++ b/tbx_cfg_hmri.m @@ -43,6 +43,6 @@ 'and will include a number of (as yet unspecified) extensions in ',... 'future updates. Please report any bugs or problems to lreninfo@gmail.com.'] }'; -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}; +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_create}; end diff --git a/tbx_scfg_hmri_denoising_create.m b/tbx_scfg_hmri_denoising_create.m new file mode 100644 index 0000000..e4272de --- /dev/null +++ b/tbx_scfg_hmri_denoising_create.m @@ -0,0 +1,148 @@ +function denoise = tbx_scfg_hmri_denoising_create + +%-------------------------------------------------------------------------- +% To enable/disable pop-up messages for all warnings - recommended when +% piloting the data processing. +%-------------------------------------------------------------------------- +popup = cfg_menu; +popup.tag = 'popup'; +popup.name = 'Pop-up warnings'; +popup.help = {['The user can review and keep track of all the information ' ... + 'collected, including warnings and other messages coming up during ' ... + 'the creation of the maps. By default, the information is logged in ' ... + 'the Matlab Command Window, in a log file saved in the "Results" ' ... + 'directory, and when more critical, displayed as a pop-up message.'], ... + ['The latter must be disabled for processing series of datasets (since it ' ... + 'blocks the execution of the code) but it is strongly recommended to ' ... + 'leave it enabled when piloting the data processing (single subject) ' ... + 'to read through and acknowledge every message and make sure ' ... + 'everything is set up properly before running the processing on a ' ... + 'whole group.'], ... + ['More information about the various messages and action to be taken ' ... + '(or not) accordingly can be found on the hMRI-Toolbox WIKI (http://hmri.info). ' ... + 'In particular, see the "Debug tips & tricks" section.']}; +popup.labels = {'Disable' 'Enable'}; +popup.values = {false true}; +popup.val = {true}; + +% --------------------------------------------------------------------- +% menu denoisingtype +% --------------------------------------------------------------------- +denoisingtype = tbx_scfg_hmri_denoising_menu; + +% --------------------------------------------------------------------- +% indir Input directory as output directory +% --------------------------------------------------------------------- +indir = cfg_entry; +indir.tag = 'indir'; +indir.name = 'Input directory'; +indir.help = {['Output files will be written to the same folder ' ... + 'as each corresponding input file.']}; +indir.strtype = 's'; +indir.num = [1 Inf]; +indir.val = {'yes'}; +% --------------------------------------------------------------------- +% outdir Output directory +% --------------------------------------------------------------------- +outdir = cfg_files; +outdir.tag = 'outdir'; +outdir.name = 'Output directory'; +outdir.help = {'Select a directory where output files will be written to.'}; +outdir.filter = 'dir'; +outdir.ufilter = '.*'; +outdir.num = [1 1]; +% --------------------------------------------------------------------- +% output Output choice +% --------------------------------------------------------------------- +output = cfg_choice; +output.tag = 'output'; +output.name = 'Output choice'; +output.help = {['Output directory can be the same as the input ' ... + 'directory for each input file or user selected']}; +output.values = {indir outdir}; +output.val = {indir}; + +% --------------------------------------------------------------------- +% subj Subject +% --------------------------------------------------------------------- +subj = cfg_branch; +subj.tag = 'subj'; +subj.name = 'Subject'; +subj.help = {'Specify a subject for denoising.'}; +subj.val = {output denoisingtype popup}; + +% --------------------------------------------------------------------- +% data Data +% --------------------------------------------------------------------- +sdata = cfg_repeat; +sdata.tag = 'data'; +sdata.name = 'Few Subjects'; +sdata.help = {'Specify the number of subjects.'}; +sdata.num = [1 Inf]; +sdata.val = {subj}; +sdata.values = {subj}; + +% --------------------------------------------------------------------- +% denoise images +% --------------------------------------------------------------------- +denoise = cfg_exbranch; +denoise.tag = 'denoise'; +denoise.name = 'Denoising'; +denoise.val = { sdata }; +denoise.help = {'Denoising of raw/processed images with different methods'}; +denoise.prog = @hmri_run_denoise; +denoise.vout = @vout_create; + +end + +% ======================================================================== +%% VOUT & OTHER SUBFUNCTIONS +% ======================================================================== +% The RUN function: +% - out = hmri_run_denoise(job) +% is defined separately. +%_______________________________________________________________________ + +function dep = vout_create(job) +% This depends on job contents, which may not be present when virtual +% outputs are calculated, depending on the denoising type. + +dnfield = fieldnames(job.subj.denoisingtype); +denoisingmethod = dnfield{1}; + +switch denoisingmethod + case 'lcpca_denoise' +%define variables and initialize cfg_dep based on availibility of phase images +arrayLength = numel(job.subj.denoisingtype.lcpca_denoise.mag_input); +%phase_bool= any(~cellfun(@isempty, job.subj.denoisingtype.lcpca_denoise.phase_input)); +phase_bool = isempty(job.subj.denoisingtype.lcpca_denoise.phase_input); +if phase_bool +cdep(1,2*arrayLength) = cfg_dep; +else + cdep(1,arrayLength) = cfg_dep; +end + +%iterate to generate dependency tags for outputs +for i=1:numel(job.subj) + for k =1:2*arrayLength + if k<=arrayLength + cdep(k) = cfg_dep; + cdep(k).sname = sprintf('lcpcaDenoised_magnitude%d',k); + idxstr = ['DenoisedMagnitude' int2str(k)]; + cdep(k).src_output = substruct('.','subj','()',{i},'.',idxstr,'()',{':'}); + cdep(k).tgt_spec = cfg_findspec({{'filter','image','strtype','e'}}); + + elseif k>arrayLength && ~phase_bool + cdep(k) = cfg_dep; + cdep(k).sname = sprintf('lcpcaDenoised_phase%d',k-arrayLength); + idxstr = ['DenoisedPhase' int2str(k-arrayLength)]; + cdep(k).src_output = substruct('.','subj','()',{i},'.',idxstr,'()',{':'}); + cdep(k).tgt_spec = cfg_findspec({{'filter','image','strtype','e'}}); + else + break + end + end +end +dep = cdep; +end +end \ No newline at end of file diff --git a/tbx_scfg_hmri_denoising_menu.m b/tbx_scfg_hmri_denoising_menu.m new file mode 100644 index 0000000..e11443e --- /dev/null +++ b/tbx_scfg_hmri_denoising_menu.m @@ -0,0 +1,130 @@ +function [denoisingtype,DNparameters] = tbx_scfg_hmri_denoising_menu + +%-------------------------------------------------------------------------- +% Denoising (DN) defaults file +%-------------------------------------------------------------------------- +DNdefaults = cfg_files; +DNdefaults.tag = 'DNdefaults'; +DNdefaults.name = 'Customised DN defaults file'; +DNdefaults.help = {['Select the [hmri_denoising*_defaults_*.m] file containing ' ... + 'the parameters to process the denoising.'], ... + ['Please make sure that the modified denoising (optional) parameters ' ... + 'are correct for your data. To create your own customised defaults file, ' ... + 'edit the distributed version and save it with a meaningful name such as ' ... + '[hmri_denoising_local_defaults.m].']}; +DNdefaults.filter = 'm'; +DNdefaults.dir = fullfile(fileparts(mfilename('fullpath')),'config','local'); +DNdefaults.ufilter = [... + '(?i)', ... ignore case + '^hmri_', ... filename starts with 'hmri_' + '.*denoising.*', ... filename includes 'denoising' somewhere in the middle + '\.m$' ... filename ends with '.m' + ]; +DNdefaults.num = [1 1]; + +% --------------------------------------------------------------------- +% Use metadata or standard defaults (no customization) +% --------------------------------------------------------------------- +DNmetadata = cfg_entry; +DNmetadata.tag = 'DNmetadata'; +DNmetadata.name = 'Use metadata or standard defaults'; +DNmetadata.help = {''}; +DNmetadata.strtype = 's'; +DNmetadata.num = [1 Inf]; +DNmetadata.val = {'yes'}; + +%-------------------------------------------------------------------------- +% DN processing parameters +%-------------------------------------------------------------------------- +DNparameters = cfg_choice; +DNparameters.tag = 'DNparameters'; +DNparameters.name = 'Processing parameters'; +DNparameters.help = {['You can either stick with metadata and standard ' ... + 'defaults parameters (recommended) or select your own customised defaults file. ' ... + 'For possible modification of default parameters, please be careful about which parameters are '... + 'required and which parameters are optional.']}; +DNparameters.values = {DNmetadata DNdefaults}; +DNparameters.val = {DNmetadata}; + +% --------------------------------------------------------------------- +% Magnitude input images +% --------------------------------------------------------------------- +mag_img = cfg_files; +mag_img.tag = 'mag_input'; +mag_img.name = 'Magnitude input (required)'; +mag_img.help = {'Select the (required) magnitude images to be denoised'}; +mag_img.filter = 'image'; +mag_img.ufilter = '.*'; +mag_img.num = [1 Inf]; + + +% --------------------------------------------------------------------- +% Phase input images +% --------------------------------------------------------------------- +phase_img = cfg_files; +phase_img.tag = 'phase_input'; +phase_img.name = 'Phase input (optional)'; +phase_img.help = {'Select the (optional) phase images to be denoised'}; +phase_img.filter = 'image'; +phase_img.ufilter = '.*'; +phase_img.num = [0 Inf]; + +% --------------------------------------------------------------------- +% Standard deviation parameter +% --------------------------------------------------------------------- +std = cfg_entry; +std.tag = 'std'; +std.name = 'Standard deviation cut-off'; +std.val = {[1.05]}; +std.strtype = 'e'; +std.num = [1 1]; +std.help = {['Specify the standard deviation cut off for denoising'], ['This parameter' ... + ' is the factor of the local noise level to remove PCA components. Higher values remove more components.' ... + ' The optimal cutoff can be increased when including more images. ' ... + ' We recommend the following cut-off sizes for various image numbers:' ... + ' 5-10: 1.05, 10-20: 1.10, 20+: 1.20']}; + + +% --------------------------------------------------------------------- +% Neighborhood size +% --------------------------------------------------------------------- +ngbsize = cfg_entry; +ngbsize.tag = 'ngbsize'; +ngbsize.name = 'Neighborhood size'; +ngbsize.val = {[4]}; +ngbsize.strtype = 'e'; +ngbsize.num = [1 1]; +ngbsize.help = {['Specify the neghborhood size'], ['This parameter' ... + ' sets the size of the local PCA neighborhood.' ... + ' Though the default size is 4, a size of 3 works better in lower resolutions.']}; + +% --------------------------------------------------------------------- +% LCPCA Denoising protocol +% --------------------------------------------------------------------- +denoisinginput_lcpca = cfg_branch; +denoisinginput_lcpca.tag = 'lcpca_denoise'; +denoisinginput_lcpca.name = 'LCPCA denoising'; +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_denoisinglocal_defaults_*.m] customised defaults file.']}; +denoisinginput_lcpca.val = {mag_img phase_img DNparameters std ngbsize}; + + +% --------------------------------------------------------------------- +% menu denoisingtype +% --------------------------------------------------------------------- +denoisingtype = cfg_choice; +denoisingtype.tag = 'denoisingtype'; +denoisingtype.name = 'Denoising method for raw/processed images'; +denoisingtype.help = {'Choose the methods for denoising.' + ['Various types of denoising methods can be handled by the hMRI ' ... + 'toolbox. See list below for a ' ... + 'brief description of each type.'] + ['- Lcpca denoising: Bazin, et al. (2019) Denoising High-Field Multi-Dimensional MRI With Local'... +'Complex PCA, Front. Neurosci. doi:10.3389/fnins.2019.01066'] + [' - No denoising:...']}; +denoisingtype.values = {denoisinginput_lcpca}; +denoisingtype.val = {denoisinginput_lcpca}; + + +end \ No newline at end of file