Hi,
I have previously posted a post about motion correction of complex images, however I have encountered a problem with the distortion correction process using TOPUP.
I have phase encode blips in opposite polarities:
AP.nii
PA.nii
topup_param.txt =
0 -1 0 0.072
0 1 0 0.072
fslmerge -t mAPPA mAP.nii.gz mPA.nii.gz
topup --imain=mAPPA --datain=topup_param.txt --config=b02b0.cnf --out=my_topup_results --fout=my_field --iout=uAPPA
% try correction on the magnitude image (4D image - 100 x 100 x 20 x 24)
applytopup --imain=m20 --inindex=1 --datain=topup_param.txt --topup=my_topup_results --method=jac --interp=spline --out=m20_DC
This looks good!
To correct the phase image, I know I need to convert the magnitude and phase images into a complex image, and then obtain the real and imaginary parts. I seem get the correct real and imaginary nifti files (preserving header info including voxel size), but the correction does not look right (not like the magnitude example above). The images look just as distorted, although they have definitely been altered:
%load phase and magnitude 4D images
p = load_untouch_nii('p20.nii');
m = load_untouch_nii('m20.nii');
p_img = double(p.img);
m_img = double(m.img);
cplx = m_img .* exp(1i.*p_img*2*pi/(4096*2));
re.img = real(cplx);
im.img = imag(cplx);
re.hdr.dime.datatype = 64;
im.hdr.dime.datatype = 64;
real = make_nii(re.img,[2.3 2.3 2.3], re.hdr.dime.datatype);
save_nii(real,'real.nii')
imag = make_nii(im.img,[2.3 2.3 2.3], im.hdr.dime.datatype);
save_nii(imag,'imag.nii')
%correction done here on real and imag.nii files (applytopup)
!$FSLDIR/bin/applytopup --imain=real --inindex=1 --datain=topup_param.txt --topup=my_topup_results --method=jac --interp=spline --out=re20_DC
!$FSLDIR/bin/applytopup --imain=imag --inindex=1 --datain=topup_param.txt --topup=my_topup_results --method=jac --interp=spline --out=im20_DC
im = load_untouch_nii('im20_DC.nii');
re = load_untouch_nii('re20_DC.nii');
im_img = double(im.img);
re_img = double(re.img);
DATA_ma = abs(re.img + 1i*im.img);
DATA_ph = angle(re.img + 1i*im.img)/pi*4096;
F_mag = make_nii(DATA_ma);
save_nii(F_mag,'fin_mag.nii')
F_phs = make_nii(DATA_ph);
save_nii(F_phs,'fin_phase.nii')
Any ideas of where I am going wrong? I just don't understand why the magnitude correction looks great, but why the same can't be done for the real and imaginary parts! I assume it is something to do with formatting?
Many thanks in advance for your help,
Lucy
|