Skip to content

Commit

Permalink
Merge pull request #398 from CWDAVIESJENKINS/develop
Browse files Browse the repository at this point in the history 
Integrate NIfTI MRS read/write/coreg implementaion
  • Loading branch information
@HJZollner
HJZollner authored Apr 27, 2022
2 parents 2b713a6 + 8d022f8 commit 012cde4
Show file tree
Hide file tree
Showing 10 changed files with 544 additions and 117 deletions.
134 changes: 70 additions & 64 deletions coreg/OspreyCoreg.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ([email protected]).';
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 ([email protected]).';
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 ([email protected]).';
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 ([email protected]).';
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 ([email protected]).';
fprintf(msg);
error(msg);
end
end
end
otherwise
msg = 'Vendor not supported. Please contact the Osprey team ([email protected]).';
fprintf(msg);
error(msg);
otherwise
msg = 'Vendor not supported. Please contact the Osprey team ([email protected]).';
fprintf(msg);
error(msg);
end
end

% Save back the image and voxel mask volumes to MRSCont
Expand Down
66 changes: 66 additions & 0 deletions coreg/coreg_nifti.m
View file
Original file line number Diff line number Diff line change
@@ -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
1 change: 1 addition & 0 deletions fit/osp_fitInitialise.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
25 changes: 20 additions & 5 deletions libraries/FID-A/inputOutput/io_loadspec_niimrs.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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

Loading

0 comments on commit 012cde4

Please sign in to comment.