Phase Contrast Imaging

[drswgw]

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.

[drswgw]

Hi Korbinian, OK CLEARSWI is working, using MRItools 3.3.1 and being careful to use single quotes 'a' not "a". ClearSWI turns out to be no good for my application: asynchronous phase based encoding of motions/velocities etc. Currently SWI doesn't output the unwrapped phase data, just a magnitude image and the MIP, which for me the MIP is giving a different number of slices than the input data and the magnitude output just isn't my cup of tea. Some of the input slices come in blank, so that may be why the MIP is missing some slices, but it messes up the masking etc. to change the number of slices. So on to Julia to try and implement the Laplacian unwrapping.This is my version of Julia already installed, thought I would upgrade to version 1.7.2 to match the current version of the manual from Julia. Which version of Julia does ROMEO and clearSWI use?

[drswgw]

Hi Korbinian, Baby steps. Julia up and running, the import for "MriResearchTools" and "FFTW" went swimmingly. What about masking? Probably need to set a mask and call it from Julia, how to do that?

[drswgw]

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.

[korbinian90]

Unwrapping

As laplacian unwrapping produces arbitrary additions/subtractions and not +-2pi, best to use romeo. In the julia REPL, you can write ?romeo to get function hints on how to use it. For your data, the best option should be to use it on the 4D dataset (timepoints in the 4th dimension) like this:

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 TEs to know if the phase images are from multi-echo with different echo times or from epi with identical echo times for each echo.

Masking

romeo does not need a mask, however, if you want to create one:

• robustmask creates a mask from the magnitude (but I think you mentioned not having a magnitude?)
• a mask from the phase via romeo:

using MriResearchTools
qmap = romeovoxelquality(phase; TEs=ones(size(phase, 4)))
mask = mask_from_voxelquality(qmap, 0.5) # threshold between 0 and 1

Julia

You can type all these commands directly in the REPL or write them in a function and call with include like you did. If you get serious with julia, I would recommend to use VS code with the julia extension. There you can interactively work on writing scripts similar to matlab.

[korbinian90]

And your case of two disconnected regions might require some more tuning of romeo. You can have a look at ?romeo and try out the different options, otherwise I'm happy to help!

[drswgw]

Hi Korbinian, The Julia version is quickly reaching a critical mass, so that's nice. Slight correction, I do have magnitude images to draw from, they are mostly all the same, so typically we average those across the time domain and use those for mask construction e.g. a logical mask in matlab and cast those to double. Sometimes we will rearrange the matrix to take out the empty space, but that's kind of a lot of work for routine stuff. We do pig brains, human brains, phantoms, bioprints etc. so masking is an important part of things for us and one size most definitely does not fit all. I think I can supply the mask (as a .nii file) and put it in my julia_files folder along with the script(s). It's all about motion encoding, so the question is what does ROMEO do with the TE's. Our case is nominally close to the EPI case, all the echoes are acquired at a fixed interval, but the data in the phase domain is responding to an entirely different set of timings.....more like FMRI I suppose...with asynchronous tasks, in this case it's for heart driven motions...two clocks...one for acquisition...one for the heart....there are a variety of gating approaches...we are trying to do image processing INSTEAD of gating.....e.g. difference between breath hold (gated in a very draconian way) and free breathing (largely ungated). Since stopping the heart is not an option...here we are. One thing that kills us is offsets....the average phase for the original portions of the image that are not wrapped is in fact data for us...so moving that original average phase around is not ideal. We would rather preserve that then offset the wrapped portions either UP by 2 pi, 4 pi, 6 pi or DOWN by 2 pi, 4 pi, 6 pi. That's perfect for us. If we can do THAT, then we can move on to decoding what time it is in the acquisition frame and in the heart frame, which ultimately gives us something a bit like CINE, except with NUMBERS....as to just what the phases are, at a given time, in a QUANTITATIVE manner...it's a phase movie, with numerical values...that's the goal state. Will work on implementing the Julia side, and maybe look into VS. Another option is to pass to Julia from Matlab. For now, keeping it old school with gedit.

[drswgw]

Hi Korbinian, Here is the acquisition clock.....but this is NOT the heart clock....although the slice timing is directly related to the heart clock...kind of like a modulo thing e.g. moduloheart(slice timing )

[drswgw]

Anyway, Not sure I understand what ROMEO is doing with TE = [1,1,1,1....] or TE = [1,2,3,4,....] the latter is what I have been using up till now. In terms of the right timing from the signal side...it probably looks something like (scale of 0-1 of the heart cycle.....[0.3, 0.45, 0.17, 0.22, 0.96, 0.11, 0.44, 0.34, 0.71.......]....not like a tango at all. On a certain level, it might make more sense to feed all the slices independently as one offs.....IF ROMEO is trying to match phases across 'times' e.g. phase at t1 matches phase at t2...that's a disaster....all the 'cuts' are jump cuts in time right now...so there is NO obvious relationship in time between adjacent acquisitions...i.e. asynchronous acquisition.

[drswgw]

Hi Korbinian, ROMEO does not need a mask? So what happens when it tries to unwrap the regions of the image that are random phases?

[korbinian90]

Not sure I understand what ROMEO is doing with TE

Two things:

1. To improve the weights that guide unwrapping with the first two volumes. It looks for phase consistency across time. If the second volume has a longer echo time, this needs to be considered. It usually greatly helps the unwrapping
2. To use temporal unwrapping. By default, only the second volume in spatially unwrapped and starting from that all other volumes are unwrapped. As long there is significantly less than π phase evolution, this works fine.

-> the output phase is still only changed by multiples of 2π. No weird scaling is happening in ROMEO

ROMEO does not need a mask? So what happens when it tries to unwrap the regions of the image that are random phases

No mask is required ROMEO starts with easy-to-unwrap regions with well behaved phase and only later goes into these random regions. It will mess up greatly there, but this doesn't influence any of the "good" areas. If the image consists of multiple parts seperated by noise, you probably have to set the seeds parameter to a higher value. Otherwise the unwrapping has to cross noise.

[korbinian90]

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 correctglobal flag

[drswgw]

Perfect, a lot to chew on, but enough to move forward...cool! Onwards from "hellojulia.jl" to "unwrap.jl". I think I will have to try a couple of things and compare those to make sure the results are consistent.

[korbinian90]

btw, I added documentation to my packages ROMEO and MriResearchTools, accesible via the little badge docs|dev (or here: https://korbinian90.github.io/MriResearchTools.jl/dev/)

[drswgw]

Hi Korbinian, First go round - only so-so. Need to force several things to move forward...first....apply the mask, and have romeo not unwrap outside the mask, and display the data MASKED. As near as I can tell without a mask...romeo is folding the phases the WRONG way....the flow is down then up....or up then down...relative to the gradient, so one side should be darker than the other throughout. Where the flow is concentrated, things should speed up i.e. get whiter (or blacker) depending on direction. seems like all the values near pi or -pi are being wrapped towards zero, or the base level is being shifted. Consider a mountain.....the data starts...partway up the mountain....a correct unwrapping should yield a topo map, everything BELOW the starting point should be LESS, and everything ABOVE the starting point should be MORE. This current unwrapping isn't doing that...or at least I don't think it is, instead it all seems to be getting compressed towards the starting point/average.

[drswgw]

Do I need to scale the initial spectra to -pi to pi, that is an acceptable, or acceptable enough, version of how the data really is.....the last frame above on the lower left looks right, darker on one side, lighter on the other, that's pre-unwrapping. And this one below looks OK too. If I was going to guess, I'd say some regions are unwrapping WELL, and other regions, not so well .

[drswgw]

Definitely promising, but not optimized yet, at least I don't think so.

[drswgw]

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]

[drswgw]

Hi Korbinian, Thought about it a little bit and what is needed is FIXED scaling. The data starts as 12 bit ([0,4095]) and typically is rescaled slightly to yield ([-4095,4095])...I think that's what the DICOM to NIFTI conversion is doing, and that's all copasetic. From there probably needs to be FIXED SCALED down to [-pi, pi] before unwrapping. Variable scaling ([min, max]) => [-pi,pi] won't work because some data is wrapped, while other data is not, the data is intrinsically scaled (by the gradients) so the original [0,4095] => [-4095,4095] => [-pi, pi] is all fine as long as it's a FIXED scaling. I can do that outside Julia, but it would be nice to be able to do fixed scaling inside Julia..... The cast to 16 bit and later on to double etc. doesn't change the original [0,4095] data density.

P.S. All the help is moving this along very nicely and rapidly. Liking the Julia interface a lot, and it runs FAST!

[drswgw]

MUCH BETTER! BAM!

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). And better control over the scaling - fixed not variable. Great progress!

[drswgw]

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. Still, this is increasingly workable...so can move on towards 'fmri' timing...yay! Getting happy!

[drswgw]

Can't import or export masks...... the BitArray{4} julia type is not apparently compatible with NIFTI, or at least not obviously. Can't turn off pre-scaling and replace it with my own absolute scaling scal*Phase4D Can't turn off zero shift

[drswgw]

Purely cosmetic touches, but moving towards something easier to run/play with... the JSON module is to pull slice timings from the NIFTI prequel/header via MRIcroGL (which is in JSON).

.

[korbinian90]

A lot to read and understand, I don't have a lot of time, so I will answer some parts Glad that it is moving along!

Masking

I 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: unwrapped = unwrapped .* mask or unwrapped[mask] .= 0 The datatype of the mask shouldn't matter for romeo. But you can just create a binary mask with: binary_mask = mask .== 0 (The . is like in matlab, but is generalized. It works nearly everywhere to make stuff run voxelwise on an array)

Scaling

The readphase() function should perform quite good scaling. You can compare it by just saving the image again. The output of phase = readphase(...) is a struct with header and image. You can access the unscaled values via phase.raw and scale them yourself with scaled_phase = phase.raw .* (pi / 2048). And best to save it back to disk and check if the range is 2pi afterwards.

Gradient

ROMEO 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.

[korbinian90]

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.

[drswgw]

Cool, all helpful...will play on!

[drswgw]

Hi Korbinian, Took a while to get Julia to plot anything, documentation is pretty bad/conflicting/confusing. For Ubuntu 16, and the python versions it uses + what I am running in this environment, matplotlib tops out at about vsn 3.0, so gave up on using the pyplot backend, the GR backend is completely hopeless (can only launch a single plot window), and plotlyjs is mediocre (can't resize plots after they are drawn, can probably pre-size...just a detail right now)....plotting is only an aid along the way...therefore plotlyjs is good enough. Working on slice timing, ultimately need to use the slice timing to reorder the unwrapped data in time...and thus check if the unwrapping has affected the temporal response. Since the unwrapping is using constant timing...it may be that the unwrapping will be better if I reorder first...we shall see.....trying to do that reordering in Julia, and it's going OK for now....I'm a hacker not a coder...so elegant...not....but getting pretty happy so far.

[drswgw]

Hi Korbinian, Made some real progress here. The julia code is pretty messy still, so I probably won't post it here - our application is pretty niche. If there's interest I'm not against posting it. Anyway...here is the update.... DATA AFTER TEMPORAL UNWRAPPING (note stripes left image coronal and only image sagittal - from asynchronous acquisition/slice timings, 2 sine waves per slice, alternating UP/DOWN phases). Data is acquired axial, but the slice timings manifest as phase offsets from slice to slice...this is a pretty common problem!

DATA AFTER RESORTING/RETIMING FOR OUR APPLICATION (note stripes are gone now, and only 1 sine wave per slice, continuous phase direction from slice to slice...SUCCESS!)

Residual issues are mostly to do with optmizing the ROMEO unwrapping. For instance comparing spatial unwrapping with temporal unwrapping, sort AFTER unwrapping, sort BEFORE unwrapping, play with other parameters etc. All very happy, Julia was the key....cutting out the middle man made a big difference here. Could probably port it all back to Matlab calls...but why...as a Julia "app" very happy...NIFTI IN => unwrapped/sorted OUT....BAM! issue closed

P.S. If you don't mind leaving your comments up, still working through/along those suggestions, very relevant to what I am doing. So long and thanks for all the fish!

[drswgw]

[drswgw]

WOW! BEAUTIFUL! (still some tweaking to do...but looks really good!)

[korbinian90]

Wow, looks very nice! Really cool video Great that you came back to report on the progress, I was wondering if you were getting good results.

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. I'm not sure if it is a problem for you, though.

P.S. I will just keep this open so all comments stay there ;)