-
Notifications
You must be signed in to change notification settings - Fork 0
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Phase Contrast Imaging #9
Comments
This movie (it's asynchronous, but there are beat frequencies) gives a reasonably good idea of the data, the phase values encode the motion(s). Need unwrapping that preserves the original phase values, just adds/subtracts 2 pi, 4 pi etc. as needed. Once that part is done can move on to the timing, gradients etc. |
UnwrappingAs laplacian unwrapping produces arbitrary additions/subtractions and not +-2pi, best to use romeo. using MriResearchTools
phase4D = readphase("path-to-phase/phase.nii")
unwrapped = romeo(phase4D; TEs=ones(size(phase,4)))
savenii(unwrapped, "unwrapped.nii"; header=header(phase4D)) romeo needs the Maskingromeo does not need a mask, however, if you want to create one:
using MriResearchTools
qmap = romeovoxelquality(phase; TEs=ones(size(phase, 4)))
mask = mask_from_voxelquality(qmap, 0.5) # threshold between 0 and 1 JuliaYou can type all these commands directly in the REPL or write them in a function and call with |
And your case of two disconnected regions might require some more tuning of romeo. You can have a look at |
Hi Korbinian, |
Anyway, |
Hi Korbinian, |
Two things:
-> the output phase is still only changed by multiples of 2π. No weird scaling is happening in ROMEO
No mask is required |
I didn't fully understand what your effective TE is. So what romeo assumes is that phase evolves linearly with time. For longer echo time, the phase changes are bigger. If the "fMRI" phase changes are small compared to π, it is probably fine using TEs=ones(...). If these changes are quite large, it would probably help if you add proportionality factors as TEs. I think best is to start with the default romeo configuration of "4D with temporal unwrapping" If that doesn't give satisfying results, you can try "4D without temporal unwrapping (unwrap-individual flag)" And if that fails you could try separate 3D unwrapping, but to avoid ending up with 2pi jumps between volumes, it is probably best to use the |
Perfect, a lot to chew on, but enough to move forward...cool! |
btw, I added documentation to my packages ROMEO and MriResearchTools, accesible via the little badge |
Definitely promising, but not optimized yet, at least I don't think so. |
This looks better, well it doesn't 'look' better, but it seems better, more consistent, more robust Took a little while to get the syntax, and still don't know how to include a mask to guide the unwrapping...and not sure if the mag image is being used or not......but...tentatively....better and/or moving forward! The numbers look wrong though...I think it is expecting the input in radians, so there is something funky going on with the scaling...should start [-pi,pi] and probably end something like [-5 pi, 5 pi] |
Hi Korbinian, P.S. All the help is moving this along very nicely and rapidly. Liking the Julia interface a lot, and it runs FAST! |
Still need to get better control of the masking (importing a NIFTI mask so far is not working, the default type for the mask in ROMEO is apparently Bool, while for NIFTI it is float/double...or more precisely from Matlab to NIFTI). |
Looks like it has shifted the original data away from the starting values...that's NOT ideal....but there are more gradients possible - this is a simple all plus gradient - but +/- are possible which allows for a subtraction later....still even in that case, the DC term is not meaningless but carries data, so would prefer not to mess with the DC component......the challenge is how to turn things OFF in ROMEO....turn it ALL OFF...then bring back ON states one at a time...that would be more useful. i.e. the initial data average was around pi/2 on one side and around pi on the other...the true mean is around 3 pi/4.....the unwrapped data has a mean near zero.....the 3pi/4 true mean has been REMOVED....the raw unwrapped data should look something like pi/2 - 5pi/4.....of course I can sacrifice the DC term if I must....but I'll be sad to see it go....back to the topo map...the top of Everest after shifting has zero elevation...so climbing it is...zero (not strictly true, because the bottom would be negative) but it's meant to be a metaphor about the importance of DC terms. |
A lot to read and understand, I don't have a lot of time, so I will answer some parts MaskingI think masking or not masking the data doesn't really influence the values inside your object after unwrapping, only in the background noise. I think this problem could be solved by adapting the scaling inside the viewer. Another option would be to mask after unwrapping with something like: ScalingThe GradientROMEO shouldn't be able to change a gradient in the phase as it adds/subtracts only 2pi from individual voxels. If there is an unwrapping error in the object, it should be clearly visible as a discrete phase jump of some voxels. |
As your data is not wrapped in most part, the unwrapped output should be identical to the input phase in large areas. You could confirm that ROMEO does what you want with: phase = readphase(...)
unwrapped = romeo(phase, ...)
difference = phase .- unwrapped
# difference .*= mask # optional masking of difference image to avoid confusing background
savenii(difference, "difference.nii") Viewing the difference.nii should show probably 0 in most of your object and 2pi / -2pi at wrapped parts. In noise the values get prabably large and random. |
Cool, all helpful...will play on! |
Hi Korbinian, |
WOW! BEAUTIFUL! (still some tweaking to do...but looks really good!) |
Wow, looks very nice! Really cool video In the video, there are some individual voxels at the edge of the top part of the object, which jump between black and white. Temporal unwrapping is supposed to keep those voxels from jumping around by always unwrapping into the same direction. P.S. I will just keep this open so all comments stay there ;) |
If borders of the brain or tissue interfaces ring that is expected, but if phase wraps ring, that means that the scanner reconstruction is messing up the data. I would love to help to fix it afterwards in romeo, but romeo is heavily relying on the fact that phase wraps are between -pi and pi in directly neighboring voxels. This is the central assumption. You could just smooth the unwrapped image, but that would destroy high frequency/small structure information. Maybe it is a setting in reconstruction on the scanner that you can change. There might be a filter that can be switched off or sometimes data is interpolated and written out in twice the resolution (for no obvious reason to me). |
Hi Kn, At this moment, I will probably do CANNY edge or something and try to repair the ringing at the phase crossover. It's NOT an ideal solution, but it moves the football. There's data in the can, then there's new data. There's a lot of data in the can, so I am all about identifying and patching as we move forward. Or said a different way, the perfect is the enemy of the good. I'm guessing here, but there is probably a filter that gets applied in order to generate the 12 bit export data from what nominally should be a higher bit count for the internal machine data. How about adaptive coil combine? That is a must for us, and this is 64 Ch data, that CERTAINLY looks like phase averaging??? |
adaptive coil combine should be fine. It computes coil combination in complex. At high field strengths (7T) it can produce open-ended fringe lines, but I haven't seen them in your data. So adaptive coil combine is for sure not the culprit. It must be something later in the chain. Usually you can look at the functor chain by opening "twix" on the Siemens scanner, clicking the dat file of the scan, opening it (I think edit or something). Then clicking the pipeline. It must be a functor after "extractmodecombine" -> "extractPHS" -> ... that messes up the phase data |
Dicom data (P dirn only) posted at https://github.com/drswgw/MRE_DATA.git it's all .7z (not bad after compression). First folder is MAG, 2nd folder is PHASE We are 3T Prisma |
Cool, so this is artifact and not waves running in/out from the beach right around pi? |
Here's a 1D look at the ringing artifact. |
hi Kn, |
no new code, but new data, I'll post it up in a little bit. |
new data.... P_9_29_22.7z, R_9_29_22.7z, S_9_29_22.7z posted at https://github.com/drswgw/MRE_DATA.git |
MRI_Tools 3.5.5 => matlab caller is indeed fixed. for instance using max seeds as an example..... from the matlab caller I can now call it from auxilary commands as '-s 32' usage: [-p PHASE] [-m MAGNITUDE] [-o OUTPUT] Now possibly this makes sense to you, but it makes little if any external sense, because the syntax is not demonstrated, elucidated, exampled. Manuals written by people who already know how to use the program...generally...are largely entirely incomprehensible to people who do not know how to use the program. Heres a code snippet in Haskell (which by the way I gave up on as a programming language, because it was so incomprehensible) to try and get the point across. It almost makes sense...but it's just not obvious how the syntax works. Now there is probably more in the help pages....but without examples the syntax is still pretty OPAQUE Which is just me bitching about the new data set still has the same phase overlap issues (I think), although some other aspects are looking better. Will have to dig in more detail to know for sure. |
Hi Kn, |
And yes, planning to add a bit better documentation and more examples to romeo, just didn't find the time yet |
Hi Kn, |
Hi Kn, |
It can do retro recon, and you can change a bunch of settings. However, only what was intended by the programmers who wrote the pipeline and you have to know which parameter to change.. |
From what you said, I am kind of wondering what the Maxwell correction is doing...seems weird that it's AFTER the conversion to PHS/MAG. |
Yes, exactly, it looks like only the Maxwell correction functor can be the culprit to mess up the phase. I think moving around in the functor chain is not possible, but it might be possible to switch it off. Two ideas:
|
Hi Kn, |
Hi Korbininan, perform unwrapping... Tried converting to single, adding pi...still no good. Maybe I need to use unmasked data? Or throw out the DC component? onephs = ones([1,szphs(4)]); |
OK, FIXED! |
Interesting.. I would have assumed that the robustmask algorithm wasn't robust enough and failed. |
Hello, Error using ROMEO could you please help me to fix it? |
Hi Korbinian & fans....
Here's my Matlab code....(2022B)...It won't necessarily make a lot of sense because my data uses a lot of gradients.
The important stuff though....
NIFTI works best for DIM 4 matrices
The new matlab niftiwrite & niftiread do work...but are VERY MATLAB CENTRIC....ROMEO uses it's own IO functions...don't mix and match.
You CAN & SHOULD ADDRESS KORBINIANS CODE using a mix of built in and external commands e.g.
"parameters.mag = load_nii(magfils{aa}).img;" uses ROMEO IO functions
IF you send YOUR code, it would make it a lot easier to see what's going on.
The errors are somewhat generic so seeing the error without the code...doesn't help much.
And types are not irrelevant...Matlab defaults to type double which is awkward in the extreme.
%%% need to beef up the discriminator(s) here
%%% now there will be 2 sets of unwrapping files...
%%% need to keep the code branches SEPARATE
% make sure the current path is at the OUT_DATA level
% this is in the 1 drn manifold
tmpdr = pwd;
if tmpdr(end-7:end) == 'OUT_DATA'
rgtdr = 1
else
rgtdr = 0
end
flds = Ffldlst();
lnflds = length(flds);
nifyes = 0;
% more discriminating trigger for THIS PARTICULAR Romeo Unwrap
for sa = 1:lnflds
if length(strfind(flds{sa},'NIFTI_OUT'))>0
nifyes = 1;
else
nifyes = 0;
end
end
if nifyes == 0 & rgtdr == 1
mkdir('NIFTI_OUT')
end
if rgtdr == 1
cd('NIFTI_OUT')
% make a proper sub-folder name
tmpnm = mltfils{Imltfno(bb)};
wrknm = strcat(tmpnm(1:end-4),'_nifti')
mkdir(wrknm)
cd(wrknm)
% drns = 1, 6 output files
%%%% PHASE
niftiwrite(squeeze(CprawSrtd3DphsDat(:,:,revslcordr,1,1,:)),'CpRS3DphsDat_1_1.nii')
niftiwrite(squeeze(CprawSrtd3DphsDat(:,:,revslcordr,2,1,:)),'CpRS3DphsDat_2_1.nii')
niftiwrite(squeeze(CprawSrtd3DphsDat(:,:,revslcordr,1,2,:)),'CpRS3DphsDat_1_2.nii')
niftiwrite(squeeze(CprawSrtd3DphsDat(:,:,revslcordr,2,2,:)),'CpRS3DphsDat_2_2.nii')
niftiwrite(squeeze(CprawSrtd3DphsDat(:,:,revslcordr,1,3,:)),'CpRS3DphsDat_1_3.nii')
niftiwrite(squeeze(CprawSrtd3DphsDat(:,:,revslcordr,2,3,:)),'CpRS3DphsDat_2_3.nii')
%%%%% MAGNITUDE
niftiwrite(squeeze(CpSrtd3DmagDat(:,:,revslcordr,1,1,:)),'CpS3DmagDat_1_1.nii')
niftiwrite(squeeze(CpSrtd3DmagDat(:,:,revslcordr,2,1,:)),'CpS3DmagDat_2_1.nii')
niftiwrite(squeeze(CpSrtd3DmagDat(:,:,revslcordr,1,2,:)),'CpS3DmagDat_1_2.nii')
niftiwrite(squeeze(CpSrtd3DmagDat(:,:,revslcordr,2,2,:)),'CpS3DmagDat_2_2.nii')
niftiwrite(squeeze(CpSrtd3DmagDat(:,:,revslcordr,1,3,:)),'CpS3DmagDat_1_3.nii')
niftiwrite(squeeze(CpSrtd3DmagDat(:,:,revslcordr,2,3,:)),'CpS3DmagDat_2_3.nii')
% ROMEO unwrap from matlab side
% clearvars
% addpath('NIfTI_20140122')
% cd('E:\Matlab_Code_Working\INTRNSC_CRTS\MRE_PTNT9_SG_9_29_22_FLAT\LOCAL_NIFTI_OUT')
tempdir = pwd;
%%% find the working files
wfils = Ffilst();
lnwfils = length(wfils);
phsfnd = strfind(wfils,'phs');
magfnd = strfind(wfils,'mag'); % bring mag files along too!!
%
pp = 1;
mm = 1;
for ab = 1:lnwfils
phslocns(ab) = length(phsfnd{ab});
if phslocns(ab) == 1
phsfils(pp) = wfils(ab);
subphsdrs(pp,:) = phsfils{pp}(end-10:end-4);
pp = pp + 1;
end
maglocns(ab) = length(magfnd{ab});
if maglocns(ab) == 1
magfils(mm) = wfils(ab);
mm = mm + 1;
end
end
%
lnphsfils = length(phsfils);
lnmagfils = length(magfils);
%%% check for ROMEO OUTPUT FILES
cflds = Ffldlst();
lncflds = length(cflds)
romflds = 0
for cc = 1:lncflds
if length(strfind(cflds{cc},'Romeo')) > 0
romflds = 1
end
end
% % phsfils{1:3}
% % magfils{1:3}
%
% %%%% Up to here...need to call the mag and phase files..using wfils + locns
%
%
% % %% Phase
% % phase_fn1 = 'PHS_R1D_pls_P.nii'; %'/path_to_data/Phase.nii";
% % phase_fn2 = 'PHS_R1D_mns_P.nii'; %'/path_to_data/Phase.nii";
% % phase_fn3 = 'PHS_R1D_pls_R.nii'; %'/path_to_data/Phase.nii";
% % phase_fn4 = 'PHS_R1D_mns_R.nii'; %'/path_to_data/Phase.nii";
% % phase_fn5 = 'PHS_R1D_pls_S.nii'; %'/path_to_data/Phase.nii";
% % phase_fn6 = 'PHS_R1D_mns_S.nii'; %'/path_to_data/Phase.nii";
% % % phase = load_nii(phase_fn).img;
% %
% % %% Optional Parameters
% % mag_fn1 = 'MAG_R1D_pls_P.nii'; %'/path_to_data/Magnitude.nii";
% % mag_fn2 = 'MAG_R1D_mns_P.nii'; %'/path_to_data/Magnitude.nii";
% % mag_fn3 = 'MAG_R1D_pls_R.nii'; %'/path_to_data/Magnitude.nii";
% % mag_fn4 = 'MAG_R1D_mns_R.nii'; %'/path_to_data/Magnitude.nii";
% % mag_fn5 = 'MAG_R1D_pls_S.nii'; %'/path_to_data/Magnitude.nii";
% % mag_fn6 = 'MAG_R1D_mns_S.nii'; %'/path_to_data/Magnitude.nii";
%
%
% mask_fn = 'Mask_LoRes.nii'; %'/path_to_data/Mask.nii";
if romflds == 0 % if no output romeo folders => unwrap, else...skip unwrapping
for aa = 1:lnphsfils % loop over all the individual 'beads on a string'
% % subdrs = ['pls_P'; 'mns_P'; 'pls_R';'mns_R'; 'pls_S'; 'mns_S'];
% % subphs = [phase_fn1;phase_fn2;phase_fn3;phase_fn4;phase_fn5;phase_fn6];
% % submag = [mag_fn1;mag_fn2;mag_fn3;mag_fn4;mag_fn5;mag_fn6];
outstr = strcat('Romeo_LOCAL_',subphsdrs(aa,:));
phase = load_nii(phsfils{aa}).img;
parameters.output_dir = fullfile(tempdir, outstr); % if not set pwd() is used
szphs = size(phase);
onephs = ones([1,szphs(4)]);
parameters.TE = onephs; %onephs; % [1,2,3]; % required for multi-echo
parameters.mag = load_nii(magfils{aa}).img;
parameters.mask = 'robustmask'; % load_nii(mask_fn).img; % can be an array or a string: 'nomask' | 'robustmask' | 'qualitymask'
parameters.calculate_B0 = true; % optional B0 calculation for multi-echo
parameters.phase_offset_correction = 'off'; % options are: 'off' | 'on' | 'bipolar'
parameters.voxel_size = [3,3,3]; %load_nii_hdr(phase_fn).dime.pixdim(2:4); % if set the written NIfTI files will have the matching voxelsize; is also used for optimal kernel size in MCPC-3D-S phase offset smoothing; can be given as [0.3, 0.3, 1.2]
%
parameters.additional_flags = '-Q -v -s 8 --template 5'; % '-q -i'; % settings are pasted directly to ROMEO cmd (see https://github.com/korbinian90/ROMEO for options)
%
% %% Suggested steps
mkdir(parameters.output_dir);
addpath(parameters.output_dir);
%
[unwrapped, B0] = ROMEO(phase, parameters);
% %
% % %rmdir(parameters.output_dir, 's') % remove the temporary ROMEO output folder
% %
unwrapped_nii = make_nii(unwrapped);
unwrapped_nii.hdr = load_nii_hdr(phsfils{aa});
end
end
%
% %%% go back up primary data directory
% cd ..
tmpdrA = pwd
% repeat this check, should have romeo folders now no matter what (somewhat duplicative)
cflds = Ffldlst();
lncflds = length(cflds)
r = 1
for cc = 1:lncflds
if length(strfind(cflds{cc},'Romeo')) > 0
romlst{r} = cflds{cc}
r = r+1;
end
end
lnromlst = length(romlst)
for dd = 1:lnromlst
% load the unwrapped data back in to Matlab
cd(romlst{dd})
% read the identifier at the end of the folder name
pmno = str2num(romlst{dd}(end-2));
drno = str2num(romlst{dd}(end));
uwPhs3Dat(:,:,:,:,pmno,drno) = niftiread('Unwrapped.nii');
Mag3Dat(:,:,:,:,pmno,drno) = niftiread('Mag.nii');
% sztmpdat = size(tmpdat)
cd ..
end
szuwPhs3Dat = size(uwPhs3Dat)
cd ..
cd ..
tmpdrB = pwd
end % ends if for rgt dr
%% DEFINE FUNCTIONS
%%% find the list of just folders in the current directory
function [wkflds] = Ffldlst();
wkdr = dir(); % files and folders together
wkflds = {wkdr([wkdr.isdir]).name}; % just folder names now
wkflds = wkflds(~ismember(wkflds,{'.','..'})); % no bogus folders/names
end
%%% find the list of just files in the current directory
function [wfils] = Ffilst();
wdr = dir(); % files and folders together
wfils = {wdr(~[wdr.isdir]).name}; % just file names now
end
…--
Scott Gordon, Ph.D.
Advanced Imaging Center at
Dartmouth Hitchcock Medical Center
& Thayer School of Engineering
Dartmouth College
802-355-3225
________________________________
From: tayebehebrahimi64 ***@***.***>
Sent: Friday, July 14, 2023 2:41 PM
To: korbinian90/ROMEO ***@***.***>
Cc: Scott W. Gordon ***@***.***>; Author ***@***.***>
Subject: Re: [korbinian90/ROMEO] Phase Contrast Imaging (Issue #9)
Hello,
I am completely new in julia and that's why I want to run romeo via matlab in linux. But I get an error says:
Error using ROMEO
ROMEO unwrapping failed! Check input files for corruption in
/N/project/dMRI1/CARE_CSI/ANALYSIS/QSM/temp/Tina/mritools_ubuntu-22.04_3.6.6/bin/romeo_tmp
could you please help me to fix it?
Thank you
—
Reply to this email directly, view it on GitHub<#9 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AO5OB72YCLPZGL7N3S4CKGTXQGHHPANCNFSM5PYJZW2Q>.
You are receiving this because you authored the thread.Message ID: ***@***.***>
|
Hi Korbinian,
Moving the dialog here....the clear SWI portion of things is all resolved. I've moved the newer ROMEO & Laplacian unwrapping stuff here, so you can close the SWI part at your leisure. Right now your answer to the question below is still in the SWI part.
The text was updated successfully, but these errors were encountered: