diff --git a/coreg/OspreyCoreg.m b/coreg/OspreyCoreg.m index 6f750352..80aefeed 100755 --- a/coreg/OspreyCoreg.m +++ b/coreg/OspreyCoreg.m @@ -73,78 +73,84 @@ gunzip(MRSCont.files_nii{kk}); MRSCont.files_nii{kk} = strrep(MRSCont.files_nii{kk},'.gz',''); end - % Call voxel mask generator depending on file type - switch MRSCont.vendor - case 'Siemens' - % Load the *.nii file provided in the job file - - vol_image = spm_vol(MRSCont.files_nii{kk}); - switch MRSCont.datatype - case 'TWIX' - [vol_mask, T1_max, voxel_ctr] = coreg_siemens(MRSCont.raw{kk}, vol_image, maskFile); - case 'RDA' - [vol_mask, T1_max, voxel_ctr] = coreg_siemens(MRSCont.raw{kk}, vol_image, maskFile); - case 'DICOM' - [vol_mask, T1_max, voxel_ctr] = coreg_siemens(MRSCont.raw{kk}, vol_image, maskFile); - otherwise - msg = 'Data type not supported. Please contact the Osprey team (gabamrs@gmail.com).'; - fprintf(msg); - error(msg); - end - case 'Philips' - % Load the *.nii file provided in the job file - vol_image = spm_vol(MRSCont.files_nii{kk}); - switch MRSCont.datatype - case 'SDAT' - if ~MRSCont.flags.isMRSI % SVS coregistration - [vol_mask, T1_max, voxel_ctr,~] = coreg_sdat(MRSCont.raw{kk}, vol_image, maskFile); - else - [vol_mask, T1_max, voxel_ctr,~] = coreg_sdat(MRSCont.raw{kk}, vol_image, maskFile,2); -% MRSCont.coreg.vol_mask_mrsi{kk} = vol_mask_mrsi; - end - case 'DATA' - if isfield(MRSCont.raw{kk}, 'geometry') - if ~MRSCont.flags.isPRIAM % SVS coregistration - [vol_mask, T1_max, voxel_ctr,~] = coreg_sdat(MRSCont.raw{kk}, vol_image, maskFile); - else - [vol_mask, T1_max, voxel_ctr,~] = coreg_sdat(MRSCont.raw{kk}, vol_image, maskFile, MRSCont.SENSE{kk}); - end - else - msg = 'Philips DATA files do not contain voxel geometry information.'; - fprintf(msg); - error(msg); - end - case 'RAW' - msg = 'Philips RAW files do not contain voxel geometry information.'; - fprintf(msg); - error(msg); - otherwise - msg = 'Data type not supported. Please contact the Osprey team (gabamrs@gmail.com).'; - fprintf(msg); - error(msg); - end - case 'GE' - [~,~,file_exten]=fileparts(MRSCont.files_nii{kk}); - if contains(file_exten,'.nii') + + if strcmp(MRSCont.datatype,'NIfTI-MRS') + vol_image = spm_vol(MRSCont.files_nii{kk}); + [vol_mask, T1_max, voxel_ctr] = coreg_nifti(MRSCont.raw{kk}, vol_image, maskFile); + else + % Call voxel mask generator depending on file type + switch MRSCont.vendor + case 'Siemens' % Load the *.nii file provided in the job file - vol_image = spm_vol(MRSCont.files_nii{kk}); - [vol_mask, T1_max, voxel_ctr] = coreg_ge_nifti(MRSCont.raw{kk}, vol_image, maskFile); - else + vol_image = spm_vol(MRSCont.files_nii{kk}); switch MRSCont.datatype - case 'P' - % Load the DICOM folder provided in the job file - [vol_mask, T1_max, vol_image, voxel_ctr] = coreg_p(MRSCont.raw{kk}, MRSCont.files_nii{kk}, maskFile); + case 'TWIX' + [vol_mask, T1_max, voxel_ctr] = coreg_siemens(MRSCont.raw{kk}, vol_image, maskFile); + case 'RDA' + [vol_mask, T1_max, voxel_ctr] = coreg_siemens(MRSCont.raw{kk}, vol_image, maskFile); + case 'DICOM' + [vol_mask, T1_max, voxel_ctr] = coreg_siemens(MRSCont.raw{kk}, vol_image, maskFile); otherwise msg = 'Data type not supported. Please contact the Osprey team (gabamrs@gmail.com).'; fprintf(msg); + error(msg); + end + case 'Philips' + % Load the *.nii file provided in the job file + vol_image = spm_vol(MRSCont.files_nii{kk}); + switch MRSCont.datatype + case 'SDAT' + if ~MRSCont.flags.isMRSI % SVS coregistration + [vol_mask, T1_max, voxel_ctr,~] = coreg_sdat(MRSCont.raw{kk}, vol_image, maskFile); + else + [vol_mask, T1_max, voxel_ctr,~] = coreg_sdat(MRSCont.raw{kk}, vol_image, maskFile,2); + % MRSCont.coreg.vol_mask_mrsi{kk} = vol_mask_mrsi; + end + case 'DATA' + if isfield(MRSCont.raw{kk}, 'geometry') + if ~MRSCont.flags.isPRIAM % SVS coregistration + [vol_mask, T1_max, voxel_ctr,~] = coreg_sdat(MRSCont.raw{kk}, vol_image, maskFile); + else + [vol_mask, T1_max, voxel_ctr,~] = coreg_sdat(MRSCont.raw{kk}, vol_image, maskFile, MRSCont.SENSE{kk}); + end + else + msg = 'Philips DATA files do not contain voxel geometry information.'; + fprintf(msg); error(msg); + end + case 'RAW' + msg = 'Philips RAW files do not contain voxel geometry information.'; + fprintf(msg); + error(msg); + otherwise + msg = 'Data type not supported. Please contact the Osprey team (gabamrs@gmail.com).'; + fprintf(msg); + error(msg); + end + case 'GE' + [~,~,file_exten]=fileparts(MRSCont.files_nii{kk}); + if contains(file_exten,'.nii') + % Load the *.nii file provided in the job file + vol_image = spm_vol(MRSCont.files_nii{kk}); + + [vol_mask, T1_max, voxel_ctr] = coreg_ge_nifti(MRSCont.raw{kk}, vol_image, maskFile); + else + switch MRSCont.datatype + case 'P' + % Load the DICOM folder provided in the job file + [vol_mask, T1_max, vol_image, voxel_ctr] = coreg_p(MRSCont.raw{kk}, MRSCont.files_nii{kk}, maskFile); + otherwise + msg = 'Data type not supported. Please contact the Osprey team (gabamrs@gmail.com).'; + fprintf(msg); + error(msg); + end end - end - otherwise - msg = 'Vendor not supported. Please contact the Osprey team (gabamrs@gmail.com).'; - fprintf(msg); - error(msg); + otherwise + msg = 'Vendor not supported. Please contact the Osprey team (gabamrs@gmail.com).'; + fprintf(msg); + error(msg); + end end % Save back the image and voxel mask volumes to MRSCont diff --git a/coreg/coreg_nifti.m b/coreg/coreg_nifti.m new file mode 100644 index 00000000..37078c79 --- /dev/null +++ b/coreg/coreg_nifti.m @@ -0,0 +1,66 @@ +% coreg_nifti.m +% Chris Davies-Jenkins, Johns Hopkins University 2022. +% +% USAGE: +% [vol_mask, T1_max, voxel_ctr] = coreg_nifti(in, vol_image, maskFile) +% +% DESCRIPTION: +% Creates a SPM volume containing a voxel mask with the same dimensions +% as the SPM volume containing a structural image. The voxel mask is +% created by transforming the NIfTI MRS voxel into T1 resolution, removing +% MRS deminsions, then writing this to maskFile using nii_tool. The +% resulting file is then loaded using SPM to remove NaN background and +% edit the header information. +% +% CREDITS: +% Thanks to Xiangrui Li, Ph.D. for his helpful suggestions using nii_tool. +% +% INPUTS: +% in = Input data structure. +% vol_image = SPM volume of the structural image that the voxel defined in +% the 'geometry' field of the input data structure 'in' is +% supposed to be co-registered to. +% maskFile = Filename under which the SPM volume of the co-registered +% voxel mask is supposed to be saved. +% +% OUTPUTS: +% vol_mask = SPM volume of the coregistered voxel mask. +% T1_max = maximum intensity of the image volume. + +function [vol_mask, T1_max, voxel_ctr] = coreg_nifti(in, vol_image, maskFile) + +image_fname = vol_image.fname; + +NiiStruct = nii_tool('load',image_fname);%load structural nifti +NiiVox = nii_tool('load',in.nii_mrs.hdr.file_name);%load voxel nifti + +% fh = nii_viewer(NiiStruct,NiiVox);%plot voxel over structural + +% Assume Nii voxel and structural are in same space: +NiiVox.hdr.sform_code = NiiStruct.hdr.sform_code; +NiiVox.hdr.qform_code = NiiStruct.hdr.qform_code; + +NiiVox.img = 1; % Overwrites image, so mask +NiiVox.hdr.dim(4:end) = 1; % remove additional MRS dimensions from header + + %transform voxel to image resolution and save under maskFile for now +nii_xform(NiiVox, NiiStruct.hdr, maskFile, 'linear', 0); + +% Load maskFile into spm to adapt some fields: +vol_mask = spm_vol(maskFile); +vol_mask.dt = vol_image.dt; + +%if ~isstruct(DualVoxel) + vol_mask.descrip = 'MRS_voxel_mask'; +%else % For PRIAM data create two voxel masks +% vol_mask.descrip = [ 'MRS_voxel_mask_' num2str(rr)]; +%end + +vol_mask = spm_write_vol(vol_mask,vol_mask.private.dat(:,:,:)); % write mask to vol + +[T1,~] = spm_read_vols(vol_image); +T1_max = max(T1(:)); + +voxel_ctr = [NiiVox.hdr.qoffset_x, NiiVox.hdr.qoffset_y, NiiVox.hdr.qoffset_z]; % CWDJ need to verify this + +end \ No newline at end of file diff --git a/fit/osp_fitInitialise.m b/fit/osp_fitInitialise.m index bb30ef42..ce8d66c3 100644 --- a/fit/osp_fitInitialise.m +++ b/fit/osp_fitInitialise.m @@ -38,6 +38,7 @@ end if ~strcmp(seq,'press') && ~strcmp(seq,'slaser') %Unable to find the localization type we will assume it is PRESS + warning('Unable to detect the localization type. We will assume it is PRESS') seq = 'press'; end diff --git a/libraries/FID-A/inputOutput/io_loadspec_niimrs.m b/libraries/FID-A/inputOutput/io_loadspec_niimrs.m index 01184d52..2be3df30 100644 --- a/libraries/FID-A/inputOutput/io_loadspec_niimrs.m +++ b/libraries/FID-A/inputOutput/io_loadspec_niimrs.m @@ -367,6 +367,13 @@ ppm = ppm + centerFreq; t = [0 : dt : (nPts-1)*dt]; +% Try and grab sequence name from header +if isfield(hdr_ext, 'SequenceName') + seq = hdr_ext.SequenceName; +else + seq = 'nii_spec'; +end +out.seq=seq; % MANDATORY FIELDS % Data & dimensions @@ -399,11 +406,6 @@ out.nii_mrs.hdr = hdr; % Save the header extension out.nii_mrs.hdr_ext = hdr_ext; -if isfield(hdr_ext, 'SequenceName') - out.seq = hdr_ext.SequenceName; -else - out.seq = 'nii_spec'; -end % Geometry geometry.size.dim1 = hdr.pixdim(2); % [mm] @@ -444,6 +446,19 @@ out.flags.isHERCULES = 0; out.flags.isPRIAM = 0; out.flags.isMRSI = 0; +if strcmp(seq,'PRESS') || strcmp(seq,'STEAM') || strcmp(seq,'SLASER') + out.flags.isUnEdited = 1; +end +if contains(seq,'MEGA') + out.flags.isMEGA = 1; +end +if strcmp(seq,'HERMES') + out.flags.isHERMES = 1; +end +if strcmp(seq,'HERCULES') + out.flags.isHERCULES = 1; +end + end diff --git a/libraries/FID-A/inputOutput/io_writeniimrs.m b/libraries/FID-A/inputOutput/io_writeniimrs.m index c624d552..bb1f345b 100644 --- a/libraries/FID-A/inputOutput/io_writeniimrs.m +++ b/libraries/FID-A/inputOutput/io_writeniimrs.m @@ -1,8 +1,10 @@ % io_writeniimrs.m % Georg Oeltzschner, Johns Hopkins University 2021. +% Adapted by +% Chris Davies-Jenkins, Johns Hopkins University 2022. % % USAGE: -% nii = io_writeniimrs(in, outfile); +% nii = io_writeniimrs(in, outfile, UserNames); % % DESCRIPTION: % Takes MRS data in matlab structure format and writes it to a NIfTI-MRS @@ -33,52 +35,102 @@ % OUTPUTS: % nii = Same as input. Not used. -function nii = io_writeniimrs(in, outfile, dim_names) -%function nii = io_writeniimrs(in, outfile) +function nii = io_writeniimrs(in, outfile, UserNames) +%function nii = io_writeniimrs(in, outfile, UserNames) % Bring FID array into NIfTI-MRS shape -% This will need work as we move to 2D and 3D arrays dims = in.dims; -% Find non-singleton dimensions -sqzDims = {}; -dimsFieldNames = fieldnames(dims); -for rr = 1:length(dimsFieldNames) - if dims.(dimsFieldNames{rr}) ~= 0 - sqzDims{end+1} = dimsFieldNames{rr}; + +% Check InNames supplied for additional dimensions +if ~exist('UserNames','var') + if dims.subSpecs || dims.extras + error('Need InNames for additional dimensions') + else + UserNames = {}; end end - -if ~isfield(dims, 'x') + +%% Manage essential dimensions +if ~isfield(dims, 'x') || dims.x==0 dims.x = 0; dim_1 = 1; +else + dim_1 = in.sz(dims.x); end -if ~isfield(dims, 'y') +if ~isfield(dims, 'y') || dims.y==0 dims.y = 0; dim_2 = 1; +else + dim_2 = in.sz(dims.y); end -if ~isfield(dims, 'z') +if ~isfield(dims, 'z')|| dims.z==0 dims.z = 0; dim_3 = 1; +else + dim_3 = in.sz(dims.z); end -if length(sqzDims) == 1 - fids = zeros(dim_1, dim_2, dim_3, size(in.fids, dims.t)); - fids(1,1,1,:) = in.fids; -elseif length(sqzDims) == 2 - if dims.extras==0 && dims.subSpecs==0 && dims.averages==0 - zeros(dim_1, dim_2, dim_3, size(in.fids, dims.t), size(in.fids, dims.coils)); - fids(1,1,1,:,:) = in.fids; - dim_5_name = {'DIM_COIL'}; - elseif dims.extras==0 && dims.subSpecs==0 && dims.coils==0 - zeros(dim_1, dim_2, dim_3, size(in.fids, dims.t), size(in.fids, dims.averages)); - fids(1,1,1,:,:) = in.fids; - dim_5_name = {'DIM_DYN'}; - elseif dims.extras==0 && dims.averages==0 && dims.coils==0 - zeros(dim_1, dim_2, dim_3, size(in.fids, dims.t), size(in.fids, dims.subSpecs)); - fids(1,1,1,:,:) = in.fids; - dim_5_name = dim_names{1}; - end - +dim_4 = in.sz(dims.t); + +ArraySize = [dim_1,dim_2,dim_3,dim_4]; +dimorder = [dims.x,dims.y,dims.z,in.dims.t]; +dimname = {}; + +%% Manage non-essential dimensions, starting with defaults +if dims.coils + ArraySize = [ArraySize, in.sz(dims.coils)]; + dimorder = [dimorder,dims.coils]; + dimname = [dimname,{'DIM_COIL'}]; +end +if dims.averages + ArraySize = [ArraySize, in.sz(in.dims.averages)]; + dimorder = [dimorder,dims.averages]; + dimname = [dimname,{'DIM_DYN'}]; +end + +%% Now manage user defined dimensions +Cnt = 1; +if in.dims.subSpecs + ArraySize = [ArraySize, in.sz(dims.subSpecs)]; + dimorder = [dimorder,dims.subSpecs]; + dimname = [dimname,UserNames(Cnt)]; + Cnt = Cnt+1; +end +if in.dims.extras + ArraySize = [ArraySize, in.sz(dims.extras)]; + dimorder = [dimorder,dims.extras]; + dimname = [dimname,UserNames(Cnt)]; + Cnt = Cnt+1; +end + +%% Do some error handling + +if length(dimorder)>7 + error('%i dimensions detected, but maximum allowed is 7',length(dimorder)) +elseif ~(length(UserNames)==(Cnt-1)) + error('%i user-defined dimensions, and %i dimension labels',Cnt-1,length(UserNames)) +elseif length(UserNames)>3 + error('%i user-defined names supplied. Maximum of 3 allowed',length(UserNames)) +end + +%% Format fid output and export to nii + +% If dimensions are ordered correctly, then just add singletons. Otherwise, +% reshape the array to correct for this. +% +% Note: Matlab will remove trailing singletons +if issorted(dimorder) + fids = reshape(in.fids,ArraySize); +else + %Permute fid non-singleton dimensions according to default BIDs precedence + NonSingletonOrdered = permute(in.fids,nonzeros(dimorder)); + + %Define overall array size including singleton dimensions + ArraySize_S = double(ArraySize==1); + ArraySize_S(~(ArraySize==1)) = size(NonSingletonOrdered); + + % Reshape output to include correct singletons + fids = reshape(NonSingletonOrdered,ArraySize_S); end % Initialize the NIfTI struct using the dicm2nii toolbox @@ -97,23 +149,16 @@ end end -% Update NIfTI-MRS header +%% Update NIfTI-MRS header newDim = nii.hdr.dim; newVoxOffset = nii.hdr.vox_offset; nii.hdr = in.nii_mrs.hdr; nii.img = conj(fids); -if exist('dim_5_name', 'var') - nii.hdr_ext.dim_5 = dim_5_name; - if exist('dim_6_name', 'var') - nii.hdr_ext.dim_6 = dim_6_name; - if exist('dim_7_name', 'var') - nii.hdr_ext.dim_7 = dim_7_name; - end - end +for JJ = 1:length(dimname) + in.nii_mrs.hdr_ext.(sprintf('dim_%i',JJ+4)) = dimname{JJ}; end - nii.ext.ecode = 44; nii.ext.edata_decoded = jsonencode(in.nii_mrs.hdr_ext); len = int32(numel(nii.ext.edata_decoded)); @@ -123,8 +168,5 @@ % nii.hdr.pixdim = in.nii_mrs.hdr.pixdim; % nii.hdr.vox_offset = in.nii_mrs.hdr.vox_offset; - - - % Save nii_tool('save', nii, outfile); \ No newline at end of file diff --git a/libraries/dicm2nii/nii_xform.m b/libraries/dicm2nii/nii_xform.m new file mode 100644 index 00000000..04cd8a36 --- /dev/null +++ b/libraries/dicm2nii/nii_xform.m @@ -0,0 +1,250 @@ +function varargout = nii_xform(src, target, rst, intrp, missVal) +% Transform a NIfTI into different resolution, or into a template space. +% +% NII_XFORM('source.nii', 'template.nii', 'result.nii') +% NII_XFORM(nii, 'template.nii', 'result.nii') +% NII_XFORM('source.nii', [1 1 1], 'result.nii') +% nii = NII_XFORM('source.nii', 'template.nii'); +% NII_XFORM('source.nii', {'template.nii' 'source2template.mat'}, 'result.nii') +% NII_XFORM('source.nii', {'template.nii' 'source2template_warp.nii.gz'}, 'result.nii') +% NII_XFORM('source.nii', 'template.nii', 'result.nii', 'nearest', 0) +% +% NII_XFORM transforms the source NIfTI, so it has the requested resolution or +% has the same dimension and resolution as the template NIfTI. +% +% Input (first two mandatory): +% 1. source file (nii, hdr/img or gz versions) or nii struct to be transformed. +% 2. The second input determines how to transform the source file: +% (1) If it is numeric and length is 1 or 3, [2 2 2] for example, it will be +% treated as requested resolution in millimeter. The result will be in +% the same coordinate system as the source file. +% (2) If it is a nii file name, a nii struct, or nii hdr struct, it will be +% used as the template. The result will have the same dimension and +% resolution as the template. The source file and the template must have +% at least one common coordinate system, otherwise the transformation +% doesn't make sense, and it will err out. With different coordinate +% systems, a transformation to align the two dataset is needed, which is +% the next case. +% (3) If the input is a cell containing two file names, it will be +% interpreted as a template nii file and a transformation. The +% transformation can be a FSL-style .mat file with 4x4 transformation +% matrix which aligns the source data to the template, in format of: +% 0.9983 -0.0432 -0.0385 -17.75 +% 0.0476 0.9914 0.1216 -14.84 +% 0.0329 -0.1232 0.9918 111.12 +% 0 0 0 1 +% The transformation can also be a FSL-style warp nii file incorporating +% both linear and no-linear transformation from the source to template. +% 3. result file name. If not provided or empty, nii struct will be returned. +% This allows to use the returned nii in script without saving to a file. +% 4. interpolation method, default 'linear'. It can also be one of 'nearest', +% 'cubic' and 'spline'. +% 5. value for missing data, default NaN. This is the value assigned to the +% location in template where no data is available in the source file. +% +% Output (optional): nii struct. +% NII_XFORM will return the struct if the output is requested or result file +% name is not provided. +% +% Please note that, once the transformation is applied to functional data, it is +% normally invalid to perform slice timing correction. Also the data type is +% changed to single unless the interpolation is 'nearest'. +% +% See also NII_VIEWER, NII_TOOL, DICM2NII + +% By Xiangrui Li (xiangrui.li@gmail.com) +% History(yymmdd): +% 151024 Write it. +% 160531 Remove narginchk so work for early matlab. +% 160907 allow src to be nii struct. +% 160923 allow target to be nii struct or hdr; Take care of logical src img. +% 161002 target can also be {tempFile warpFile}. +% 170119 resolution can be singular. +% 180219 treat formcode 3 and 4 the same. + +if nargin<2 || nargin>5, help('nii_xform'); error('Wrong number of input.'); end +if nargin<3, rst = []; end +if nargin<4 || isempty(intrp), intrp = 'linear'; end +if nargin<5 || isempty(missVal), missVal = nan; else, missVal = missVal(1); end +intrp = lower(intrp); +quat2R = nii_viewer('func_handle', 'quat2R'); + +if isstruct(src), nii = src; +else, nii = nii_tool('load', src); +end + +if isstruct(target) || ischar(target) || (iscell(target) && numel(target)==1) + hdr = get_hdr(target); +elseif iscell(target) + hdr = get_hdr(target{1}); + if hdr.sform_code>0, R0 = [hdr.srow_x; hdr.srow_y; hdr.srow_z; 0 0 0 1]; + elseif hdr.qform_code>0, R0 = quat2R(hdr); + end + + [~, ~, ext] = fileparts(target{2}); + if strcmpi(ext, '.mat') % template and xform file names + R = load(target{2}, '-ascii'); + if ~isequal(size(R), [4 4]), error('Invalid transformation file.'); end + else % template and warp file names + warp_img_fsl = nii_tool('img', target{2}); + if ~isequal(size(warp_img_fsl), [hdr.dim(2:4) 3]) + error('warp file and template file img size don''t match.'); + end + R = eye(4); + end + + if nii.hdr.sform_code>0 + R1 = [nii.hdr.srow_x; nii.hdr.srow_y; nii.hdr.srow_z; 0 0 0 1]; + elseif nii.hdr.qform_code>0 + R1 = quat2R(nii.hdr); + end + + % I thought it is something like R = R0 \ R * R1; but it is way off. It + % seems the transform info in src nii is irrevelant, but direction must be + % used: Left/right-handed storage of both src and target img won't affect + % FSL alignment R. Alignment R may not be diag-major, and can be negative + % for major axes (e.g. cor/sag slices). + + % Following works for tested FSL .mat and warp.nii files: Any better way? + % R0: target; R1: source; R: xform; result is also R + R = R0 / diag([hdr.pixdim(2:4) 1]) * R * diag([nii.hdr.pixdim(2:4) 1]); + [~, i1] = max(abs(R1(1:3,1:3))); + [~, i0] = max(abs(R(1:3,1:3))); + flp = sign(R(i0+[0 4 8])) ~= sign(R1(i1+[0 4 8])); + if any(flp) + rotM = diag([1-flp*2 1]); + rotM(1:3,4) = (nii.hdr.dim(2:4)-1) .* flp; + R = R / rotM; + end +elseif isnumeric(target) && any(numel(target)==[1 3]) % new resolution in mm + if numel(target)==1, target = target * [1 1 1]; end + hdr = nii.hdr; + ratio = target(:)' ./ hdr.pixdim(2:4); + hdr.pixdim(2:4) = target; + hdr.dim(2:4) = round(hdr.dim(2:4) ./ ratio); + if hdr.sform_code>0 + hdr.srow_x(1:3) = hdr.srow_x(1:3) .* ratio; + hdr.srow_y(1:3) = hdr.srow_y(1:3) .* ratio; + hdr.srow_z(1:3) = hdr.srow_z(1:3) .* ratio; + end +else + error('Invalid template or resolution input.'); +end + +if ~iscell(target) + s = hdr.sform_code; + q = hdr.qform_code; + sq = [nii.hdr.sform_code nii.hdr.qform_code]; + if s>0 && (any(s == sq) || (s>2 && (any(sq==3) || any(sq==4)))) + R0 = [hdr.srow_x; hdr.srow_y; hdr.srow_z; 0 0 0 1]; + frm = s; + elseif any(q == sq) || (q>2 && (any(sq==3) || any(sq==4))) + R0 = quat2R(hdr); + frm = q; + else + switch q + case 1 + targetqspace = 'Scanner space'; + case 2 + targetqspace = 'Coordinates aligned to another file''s, or to anatomical "truth". '; + case 3 + targetqspace = 'Coordinates aligned to Talairach-Tournoux Atlas'; + case 4 + targetqspace = 'MNI 152 normalized coordinates.'; + otherwise + targetqspace = 'Undefined coordinate system'; + end + switch sq(2) + case 1 + sourceqspace = 'Scanner space'; + case 2 + sourceqspace = 'Coordinates aligned to another file''s, or to anatomical "truth". '; + case 3 + sourceqspace = 'Coordinates aligned to Talairach-Tournoux Atlas'; + case 4 + sourceqspace = 'MNI 152 normalized coordinates.'; + otherwise + sourceqspace = 'Undefined coordinate system'; + end + error('%s\ntarget file: %s\nsource file: %s','No matching transformation between source and template:',targetqspace,sourceqspace); + end + + if sq(1) == frm || (sq(1)>2 && frm>2) || sq(2)<1 + R = [nii.hdr.srow_x; nii.hdr.srow_y; nii.hdr.srow_z; 0 0 0 1]; + else + R = quat2R(nii.hdr); + end +end + +d = single(hdr.dim(2:4)); +I = ones([4 d], 'single'); +[I(1,:,:,:), I(2,:,:,:), I(3,:,:,:)] = ndgrid(0:d(1)-1, 0:d(2)-1, 0:d(3)-1); +I = reshape(I, 4, []); % template ijk +if exist('warp_img_fsl', 'var') + warp_img_fsl = reshape(warp_img_fsl, [], 3)'; + if det(R0(1:3,1:3))<0, warp_img_fsl(1,:) = -warp_img_fsl(1,:); end % correct? + warp_img_fsl(4,:) = 0; + I = R \ (R0 * I + warp_img_fsl) + 1; % ijk+1 (fraction) in source +else + I = R \ (R0 * I) + 1; % ijk+1 (fraction) in source +end +I = reshape(I(1:3,:)', [d 3]); + +V = nii.img; isbin = islogical(V); +d48 = size(V); % in case of RGB +d48(numel(d48)+1:4) = 1; d48(1:3) = []; +if isbin + intrp = 'nearest'; missVal = false; + nii.img = zeros([d d48], 'uint8'); +elseif isinteger(V) + nii.img = zeros([d d48], 'single'); +else + nii.img = zeros([d d48], class(V)); +end +if ~isfloat(V), V = single(V); end +if strcmpi(intrp, 'nearest'), I = round(I); end % needed for edge voxels +if size(V,1)<2 + V = repmat(V,[3 1 1 1]); % replicate to help interp + I(:,:,:,1) = I(:,:,:,1)+1; % use middle slice + I(:,:,:,1) = I(:,:,:,1) + double(I(:,:,:,1)<1.5 | I(:,:,:,1)>2.5)*1e3; % needed for edge voxels +end +if size(V,2)<2, V = repmat(V,[1 3 1 1]); I(:,:,:,2) = I(:,:,:,2)+1; I(:,:,:,2) = I(:,:,:,2) + double(I(:,:,:,2)<1.5 | I(:,:,:,2)>2.5)*1e3; end +if size(V,3)<2, V = repmat(V,[1 1 3 1]); I(:,:,:,3) = I(:,:,:,3)+1; I(:,:,:,3) = I(:,:,:,3) + double(I(:,:,:,3)<1.5 | I(:,:,:,3)>2.5)*1e3; end + +try + F = griddedInterpolant(V(:,:,:,1), intrp, 'none'); % since 2014? + for i = 1:prod(d48) + F.Values = V(:,:,:,i); + nii.img(:,:,:,i) = F(I(:,:,:,1), I(:,:,:,2), I(:,:,:,3)); + end + if ~isnan(missVal), nii.img(isnan(nii.img)) = missVal; end +catch + for i = 1:prod(d48) + nii.img(:,:,:,i) = interp3(V(:,:,:,i), I(:,:,:,2), I(:,:,:,1), I(:,:,:,3), intrp, missVal); + end +end +if isbin, nii.img = logical(nii.img); end + +% copy xform info from template to rst nii +nii.hdr.pixdim(1:4) = hdr.pixdim(1:4); +flds = {'qform_code' 'sform_code' 'srow_x' 'srow_y' 'srow_z' ... + 'quatern_b' 'quatern_c' 'quatern_d' 'qoffset_x' 'qoffset_y' 'qoffset_z'}; +for i = 1:numel(flds), nii.hdr.(flds{i}) = hdr.(flds{i}); end + +if ~isempty(rst), nii_tool('save', nii, rst); end +if nargout || isempty(rst), varargout{1} = nii_tool('update', nii); end + +%% +function hdr = get_hdr(in) +if iscell(in), in = in{1}; end +if isstruct(in) + if isfield(in, 'hdr') % nii input + hdr = in.hdr; + elseif isfield(in, 'sform_code') % hdr input + hdr = in; + else + error('Invalid input: %s', inputname(1)); + end +else % template file name + hdr = nii_tool('hdr', in); +end \ No newline at end of file diff --git a/load/osp_LoadNII.m b/load/osp_LoadNII.m index edb1c73c..79b5e600 100644 --- a/load/osp_LoadNII.m +++ b/load/osp_LoadNII.m @@ -118,8 +118,51 @@ end end + end + + % Set flag and save data under appropriate name + if raw.dims.coils == 0 + MRSCont.flags.coilsCombined = 1; + MRSCont.raw{kk} = raw; + if MRSCont.flags.hasRef + MRSCont.raw_ref{kk} = raw_ref; + end + + if MRSCont.flags.hasWater + MRSCont.raw_w{kk} = raw_w; + end + + else + MRSCont.flags.coilsCombined = 0; + MRSCont.raw_uncomb{kk} = raw; + if MRSCont.flags.hasRef + MRSCont.raw_ref_uncomb{kk} = raw_ref; + end + + if MRSCont.flags.hasWater + MRSCont.raw_w_uncomb{kk} = raw_w; + end + end +end +% Try to get some params from NIFTI +if isfield(raw.nii_mrs.hdr_ext,'Manufacturer') + switch lower(raw.nii_mrs.hdr_ext.Manufacturer) + case 'siemens' + MRSCont.vendor = 'Siemens'; + case 'philips' + MRSCont.vendor = 'Philips'; + case 'ge' + MRSCont.vendor = 'GE'; end +else + MRSCont.vendor = 'NIFTI'; end +if isfield(raw.nii_mrs.hdr_ext,'ProtocolName') && ~isempty(raw.nii_mrs.hdr_ext.ProtocolName) + MRSCont.seq = raw.nii_mrs.hdr_ext.ProtocolName; +else + MRSCont.seq = 'NIFTI'; +end + time = toc(refLoadTime); [~] = printLog('done',time,MRSCont.nDatasets(1),MRSCont.nDatasets(2),progressText,MRSCont.flags.isGUI ,MRSCont.flags.isMRSI); diff --git a/plot/osp_plotSegment.m b/plot/osp_plotSegment.m index 31fe769d..279eb4ac 100755 --- a/plot/osp_plotSegment.m +++ b/plot/osp_plotSegment.m @@ -96,7 +96,11 @@ if strcmp(T1ext, '.gz') && ~exist(MRSCont.coreg.vol_image{kk}.fname,'file') gunzip(MRSCont.coreg.vol_image{kk}.fname) end - + [~, ~, Voxext] = fileparts(MRSCont.coreg.vol_mask{kk}.fname); + if ~strcmp(Voxext, '.gz') && ~exist(MRSCont.coreg.vol_mask{kk}.fname,'file') + gunzip([MRSCont.coreg.vol_mask{kk}.fname,'.gz']) + end + %%% 3. SET UP THREE PLANE IMAGE %%% if ~(isfield(MRSCont.flags,'isPRIAM') && (MRSCont.flags.isPRIAM == 1)) [img_montage,vox_t_size] = osp_extract_three_plane_image_seg(MRSCont.coreg.vol_image{kk}.fname, MRSCont.coreg.vol_mask{kk}.fname,vol_GM_mask,vol_WM_mask,vol_CSF_mask,MRSCont.coreg.voxel_ctr{kk},MRSCont.coreg.T1_max{kk}); diff --git a/process/OspreyProcess.m b/process/OspreyProcess.m index ab868ad1..4a5bacc0 100644 --- a/process/OspreyProcess.m +++ b/process/OspreyProcess.m @@ -787,7 +787,7 @@ [raw_C] = op_rmempty(raw_C); % Remove empty lines raw_D = op_takesubspec(raw,4); [raw_D] = op_rmempty(raw_D); % Remove empty lines - raw = op_concatAverages(raw_A,raw_B,raw_C,raw_D); + raw = op_concatAverages(op_concatAverages(raw_A,raw_B),op_concatAverages(raw_C,raw_D)); end % Align and verage the refernce data diff --git a/process/osp_combineCoils.m b/process/osp_combineCoils.m index 4e47ac77..d3e0c81f 100644 --- a/process/osp_combineCoils.m +++ b/process/osp_combineCoils.m @@ -72,7 +72,7 @@ raw_ref_B = op_takesubspec(MRSCont.raw_ref{ll,kk},2); raw_ref_C = op_takesubspec(MRSCont.raw_ref{ll,kk},3); raw_ref_D = op_takesubspec(MRSCont.raw_ref{ll,kk},4); - MRSCont.raw_ref{ll,kk} = op_concatAverages(raw_ref_A,raw_ref_B,raw_ref_C,raw_ref_D); + MRSCont.raw_ref{ll,kk} = op_concatAverages(op_concatAverages(raw_ref_A,raw_ref_B),op_concatAverages(raw_ref_C,raw_ref_D)); end end else