diff --git a/ext_libs/EasyReg/ea_easyreg.m b/ext_libs/EasyReg/ea_easyreg.m index 8fec2cc69..caf2fdbdc 100644 --- a/ext_libs/EasyReg/ea_easyreg.m +++ b/ext_libs/EasyReg/ea_easyreg.m @@ -67,14 +67,14 @@ function [] = freesurfer_nii_to_itk_h5(warp_file_in, warp_file_out) -% substract mm coordinates for each voxel -n = load_nii(warp_file_in); +% Subtract millimeter coordinates for each voxel +n = load_untouch_nii(warp_file_in); s = n.hdr.dime.dim(2:4); index = 1:prod(s); [v1,v2,v3] = ind2sub(s,index); mm = ea_vox2mm([v1',v2',v3'], ea_get_affine(warp_file_in)); % Need to do vox2mm since EasyReg uses disp_crs format) mm = reshape(mm, [s,3]); -out = n.img - mm; +out = apply_nifti_scaling(n) - mm; % reshape output out_rows = [-reshape(out(:,:,:,1),1,[]); -reshape(out(:,:,:,2),1,[]); reshape(out(:,:,:,3),1,[])]; @@ -83,23 +83,27 @@ % copy template h5 file copyfile(fullfile(ea_getearoot, 'ext_libs', 'EasyReg', 'itk_h5_template.h5'), warp_file_out); -if ~strcmp(ea_getspace, 'MNI152NLin2009bAsym') - % calculate TransformFixedParameters - spacedef = ea_getspacedef; - primarytemplate = [ea_space, spacedef.templates{1}, '.nii']; - hdr = ea_fslhd(primarytemplate); - TransformFixedParameters = zeros(18,1); - TransformFixedParameters(1:3) = [hdr.dim1; hdr.dim2; hdr.dim3]; - TransformFixedParameters(4:6) = [-hdr.sto_xyz1(4); -hdr.sto_xyz2(4); hdr.sto_xyz3(4)]; % RAS to LPS applied - TransformFixedParameters(7:9) = [hdr.pixdim1; hdr.pixdim2; hdr.pixdim3]; - TransformFixedParameters(10:18) = [-hdr.sto_xyz1(1:3)'/hdr.pixdim1; -hdr.sto_xyz2(1:3)'/hdr.pixdim2; hdr.sto_xyz3(1:3)'/hdr.pixdim3]; % RAS to LPS applied - - % update TransformFixedParameters in h5 - h5write(warp_file_out, "/TransformGroup/0/TransformFixedParameters", TransformFixedParameters); -end +% update TransformFixedParameters in h5 +h5write(warp_file_out, "/TransformGroup/0/TransformFixedParameters", ea_field_ref2itk(warp_file_in)); % save TransformParameters in h5 h5create(warp_file_out, "/TransformGroup/0/TransformParameters", numel(out_column)); h5write(warp_file_out, "/TransformGroup/0/TransformParameters", out_column); end + + +function img = apply_nifti_scaling(nii) + +img = double(nii.img); +slope = double(nii.hdr.dime.scl_slope); +intercept = double(nii.hdr.dime.scl_inter); + +if isfinite(slope) && slope ~= 0 + if ~isfinite(intercept) + intercept = 0; + end + img = img .* slope + intercept; +end + +end diff --git a/ext_libs/SynthMorph/ea_synthmorph.m b/ext_libs/SynthMorph/ea_synthmorph.m index 67ca318b3..cb7c89a72 100644 --- a/ext_libs/SynthMorph/ea_synthmorph.m +++ b/ext_libs/SynthMorph/ea_synthmorph.m @@ -64,7 +64,7 @@ function [] = freesurfer_nii_to_itk_h5(warp_file_in, warp_file_out) -% substract mm coordinates for each voxel +% Subtract millimeter coordinates for each voxel n = load_nii(warp_file_in); out = n.img; @@ -75,20 +75,8 @@ % copy template h5 file copyfile(fullfile(ea_getearoot, 'ext_libs', 'EasyReg', 'itk_h5_template.h5'), warp_file_out); -if ~strcmp(ea_getspace, 'MNI152NLin2009bAsym') - % calculate TransformFixedParameters - spacedef = ea_getspacedef; - primarytemplate = [ea_space, spacedef.templates{1}, '.nii']; - hdr = ea_fslhd(primarytemplate); - TransformFixedParameters = zeros(18,1); - TransformFixedParameters(1:3) = [hdr.dim1; hdr.dim2; hdr.dim3]; - TransformFixedParameters(4:6) = [-hdr.sto_xyz1(4); -hdr.sto_xyz2(4); hdr.sto_xyz3(4)]; % RAS to LPS applied - TransformFixedParameters(7:9) = [hdr.pixdim1; hdr.pixdim2; hdr.pixdim3]; - TransformFixedParameters(10:18) = [-hdr.sto_xyz1(1:3)'/hdr.pixdim1; -hdr.sto_xyz2(1:3)'/hdr.pixdim2; hdr.sto_xyz3(1:3)'/hdr.pixdim3]; % RAS to LPS applied - - % update TransformFixedParameters in h5 - h5write(warp_file_out, "/TransformGroup/0/TransformFixedParameters", TransformFixedParameters); -end +% update TransformFixedParameters in h5 +h5write(warp_file_out, "/TransformGroup/0/TransformFixedParameters", ea_field_ref2itk(warp_file_in)); % save TransformParameters in h5 h5create(warp_file_out,"/TransformGroup/0/TransformParameters", numel(out_column)); diff --git a/helpers/ea_field_ref2itk.m b/helpers/ea_field_ref2itk.m new file mode 100644 index 000000000..79c2456a7 --- /dev/null +++ b/helpers/ea_field_ref2itk.m @@ -0,0 +1,25 @@ +function TransformFixedParameters = ea_field_ref2itk(reference_image) +%EA_FIELD_REF2ITK Build ITK displacement-field fixed parameters from a NIfTI grid. + +hdr = ea_fslhd(reference_image); + +spacing = [hdr.pixdim1, hdr.pixdim2, hdr.pixdim3]; +if any(~isfinite(spacing)) || any(spacing <= 0) + error('Invalid voxel spacing in reference image: %s', reference_image); +end + +stoX = hdr.sto_xyz1(1:3); +stoY = hdr.sto_xyz2(1:3); +stoZ = hdr.sto_xyz3(1:3); +voxToRAS = [stoX(:)'; stoY(:)'; stoZ(:)']; + +rasToLPS = diag([-1, -1, 1]); +direction = rasToLPS * voxToRAS * diag(1 ./ spacing); + +TransformFixedParameters = zeros(18, 1); +TransformFixedParameters(1:3) = [hdr.dim1; hdr.dim2; hdr.dim3]; +TransformFixedParameters(4:6) = [-hdr.sto_xyz1(4); -hdr.sto_xyz2(4); hdr.sto_xyz3(4)]; +TransformFixedParameters(7:9) = spacing(:); +TransformFixedParameters(10:18) = reshape(direction', [], 1); + +end