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]

Aaah...fascinating...and I was just about to launch into building some test 'checkerboards' in time and space. If you have good ideas for test images...now is the time....now is the space....I think it is important to understand the 'dynamic' range aspects of the whole thing....it is not a trivial undertaking to 'de-alias' the wrapped data. There is some control on the magnet side, but that is useless without knowing what romeo can handle or not.....

touching in time, or touching in space??? It's guilt by association, right?

P.S. ....Slick Julia syntax there... Float32.(regions) vs regions = broadcast(Float32, regions)

[drswgw]

Pretty basic test image.... 2 gaussians that run 0 to < 8192, find the region that is going to get folded in [-4095,4095], fold it.
Using Mayo Awave colormap (it's a pretty decent bipolar colormap) (sorry cheating here - this is in Matlab, could probably do it in Julia....but aiming for expediency here)

TEST, TEST, 1,2,3

Add in time variation

TIME VARYING (not folded, zero offset)

Temporal with Mask for folding

Temporal - Folded

Build a semi-decent noise map (inversely scaled by sqrt(signal strength)). Noise that is too high is just truncated back to +/- 4095

Noise Map - Semi-Decent

COMBINED SIGNAL + NOISE = TEST SPATIAL/TEMPORAL IMAGE

Ooops...will need more slices to get the snazzy breakdown into spatial and temporal (right now it thinks time is slices)

MRIcroGL VIEW of TEST IMAGE (with noise)

Repeated the test image out to 8 slices (took some matrix rearranging to switch slice and phase). So all slices have same repetitive phase profile.

REGION WITH NO WRAP (noise turned down)

REGION WITH WRAP (noise turned down)

Onwards....to unwrapping....

[drswgw]

Temporal Unwrapping Failing so far. Tried the test image, FAIL. Note the slices are the same, the variation (and wrapping) are in the time domain. Fmkmsk is a function that calls different masking options, in this case the msk type calls for a phase based mask.

msk4D = Fmkmsk(msktyp, 0, phase4D) # Fmkmsk(msk type, mag file, phs file) szmsk4D = size(msk4D) print("size msk4D:\n",szmsk4D,"\n")

% START UNWRAPPING print("Calculating ROMEO Masking/Unwrapping...\n") unwrapped = romeo(phase4D; TEs=ones(size(phase4D,4)), maxseeds=50, mask=msk4D)

NOT UNWRAPPED YET

[drswgw]

Spatial Unwrapping looking a LOT BETTER! Scaling looks decent too!

print("Calculating ROMEO Masking/Unwrapping...\n") unwrapped = romeo(phase4D; TEs=ones(size(phase4D,4)), maxseeds=50, mask=msk4D,individual=true,globalcorrect=true)

SPATIAL UNWRAPPING with global correct

[drswgw]

Some progress on the stronger gradients too! With strong gradients romeo really freaks out, the data is very un-smooth, and un-connected.

Stronger Gradients: Raw Spatial Unwrapping (poor)

Pre-sorting into pls/min gradients and even/odd slices (i.e. break into 4 submatrices), yields a much more amenable problem for unwrapping....temporal unwrapping works now!

EvenPls Submatrix (1 of 4): Temporal Unwrapping (good)

Then recombine the submatrices...yields a properly unwrapped whole matrix (crazy phase jumps intact...important since the phase jumps ARE the signal). Still a ways to go to parse the acquisition times, check for consistency etc., but pretty pleased with how it's working out so far.

FULL MATRIX: Unwrapped properly by combining unwrapped submatrices

[drswgw]

OK, unwrapping going pretty well...however, gradients are still too weak (I think?.... based on some behind the scenes shenanigans). So onwards towards stronger gradients. Thanks for the help so far, the divide and conquer approach is basically working....the trend is ongoing...fill up the bandwidth with signal, which pushes more and more towards LARGE phase oscillations, which are increasingly wrapped. CRITICALLY IMPORTANT to have robust phase unwrapping, with a very high level of control over the unwrapping routines.

The slope based unwrapping remains of interest....things ARE noisy in the real world though...so that is a concern. I like the regions back inspection approach...need to play with that in more detail.

[drswgw]

Still working on things...getting there. This is an error from trying to MASK single slice data. Routine works OK for multi-slice data.

I don't think it's on the NIFTI creation side, MRIcroGL uses dcn2niix which is pretty robust. Things look right for the NIFTI, so I think it's a ROMEO thing...admittedly, I may well be calling ROMEO wrong as I am trying to use the same code for multi-slice as single slice...am rewriting for JUST single slice, but nice if I can figure out what the error means?

[korbinian90]

The warning is from smoothing. There is somewhere smoothing applied along the slice direction, which has only size 1. However, the warning just says the kernel size is reduced in slice direction, so this shouldn't be a problem. In Julia, if you want to keep the dimension, Zou can index like this:

a = ones(2,3,4)
b = a[:,[1],:]
size(b) # (2,1,4)

If I get to it, I will see what happens when I try ROMEO on a single slice

[drswgw]

Hi Korbinian, No longer freaking out about a not very important robustmask warning....starting to get to the actual meat now (how can you have any pudding, if you don't eat your meat)...really liking this small tweak to Julia, the menu just calls slight variations on the same basic script to direct to different data files, maybe different unwrapping conditions etc.

I think the basic strategy of dividing up the different gradients (am working mostly on single slice now) is working pretty well...but the single slice unwrapping is languishing behind multi-slice unwrapping. From what I can ken, the problem is being able to identify relatively abrupt wraps...there's a spatial extent thing going on, these small imaging volumes are quite precise, but the unwrapping seems to be expecting LARGE homogeneous regions....which don't exist in this high gradient data...think the local hill versus K2. So I need to unwrap K2 (apparently). More pictures to follow....

[drswgw]

I'm going to dive back into the phase unwrap map/quality thing. This is for WEAK gradients, with stronger gradients, things don't get better. It's a spatial extent thing I think. Small regions aren't getting unwrapped, but they should be (or they are getting unwrapped wrong). Obviously the ultimate small region with large phase changes is pure noise, so I get that the algorithm must be weighted against chasing after every noisy voxel, but in this case the image features ARE small. I don't think I am getting control with maxseeds. What I need to set is the PATCH AREA directly....think unwrapping a leopards spots...versus a zebra's stripes....it's not necessarily about NUMBER, but it is about AREA. I have the leopards spots case, NOT the zebra stripes. Look at the kidney shaped region that is not getting unwrapped properly. I think it is too small to be getting 'recognized'. If I increase maxseeds, some micro-parts of the kidney region get unwrapped, but still not all of it. I think it's an area thing, and I don't think the algorithm is 'aware' of these small regions.

In a general way, the maxseeds 512 case seems to be performing a bit better, but I am NOT imbued with confidence that Romeo is 'listening' to the maxseeds command in a reliable way. Do the seeds compare notes? or not? How do I get the seeds to compare notes, or not? One side of this is the temporal aspect is more tightly linked than the individual spatial unwrapping (although with only one slice things are...limited, and as a general comment, the unwrapping is generally better with more slices). 1 slice is just a stop along the way...but it is an important milestone...

LEFT - 1slc, Spatial unwrapping, maxseeds = 2 RIGHT - 1slc, Temporal unwrapping, maxseeds = 2

LEFT - 1slc, Spatial Unwrap, maxseeds = 512 RIGHT - 1slc,Temporal Unwrap, maxseeds = 512

[drswgw]

The meat => Positive vs Negative gradients maxseeds = 4096, Temporal unwrapping (sorting into gradients/times hopefully correct). Need the unwrapping to be rock-solid in order to move forward to comparing different gradients.

[korbinian90]

Do you use the flag merge_regions=true? That is supposed to merge the individual patches from different seeds

[korbinian90]

There is some heuristic approach to adding a new seed when only low quality edges remain for unwrapping. So maxseeds is not always reached.
Before adjusting this you could inspect the ROMEO weights that guide the unwrapping with calculateweights(<same arguments as romeo>). This returns a 4D array, where weights[1,:,:,:] are edges in x-direction, weights[2,:,:,:] in y-direction, weights[3,:,:,:] in z-direction.

To see how the regions looked after detection, you can use:

regions=zeros(UInt8, size(phase)[1:3]);
uw = romeo(phase; ...., regions, merge_regions=true, correct_regions=true)
savenii(Float32.(regions), ....)

[drswgw]

Hi Korbinian, The current comments + the previous insights into making 'unwrapping maps'...should do it for now...here are the current unwrapping calls.... S => Spatial, T => Temporal, pls => Plus gradient, mns => Minus gradient, raw has the gradients combined (i.e. as acquired...so far unwrapping mostly fails when the gradients are combined..once they reach a certain strength....because...the changes are too abrupt). So the divide and conquer approach so far is hands down the best way to go. Notes: Spatial unwrapping uses a reduced dimension mask. header = header(julia array) generally fails with the lone exception of julia arrays made directly by readphase....how to debug this is beyond me right now. Header is handy though, because it allows the display/NIFTI ordering to be set to neurological or radiological, so I am sad to see it broken currently.

[drswgw]

OK, added merge_regions=true (better for sure). Here is the basic issue. About 95% sure this SHOULD be getting unwrapped, but it is not. Romeo just isn't 'seeing' the phase boundary. Slightly weaker gradient...maxseeds = 1024, merge_regions=true, 1slc, Temporal unwrapping

Apt3,PplsT,time A

Apt3,PplsT,timeB

[drswgw]

Here are some reasonably clear examples Apt3, SPATIAL UNWRAPPING, merge_regions=true, maxseeds = 1024.

I think what is going on is the actual data looks kind of boundary like, but the true phase boundaries REALLY look like boundaries - REAL BOUNDARIES vs FAKE BOUNDARIES...ROMEO is totally struggling to tell the difference.

Back to blood vessels....fast flow in the middle of the vessel is phase wrapped.

Follow the moving time line........|.........

[korbinian90]

e.g. header(phase) only works if it is a nifti type. You can check with typeof(phase) to see what it is. After reading with readphase or readmag it should be of type nifti. A julia array does not contain nifti information, so header does not work.

[korbinian90]

What is the phase difference of the boundaries in your example? If it is 2pi, it should be unwrapped. It is hard to see what the reason is that it fails. The "problem voxels" are still inside the mask? ROMEO does not stop before it encounters all voxels inside the mask. Problems could be that it takes a wrong path into the area, but that doesn't look like it. Otherwise, it has to be that region merging fails. I haven't used that feature much. It is experimental and might be not always managing to merge the regions. You can try to add a correct_regions=true flag. If there are multiple regions left after region merging, this shifts each region towards 0 median (by only adding subtracting 2pi).

[drswgw]

Hi Korbinian, The header thing makes perfect sense, the header is a NIFTI header, not part of Julia per se. If there was a way to see the NIFTI header directly? then I could add it to Julia e.g. a default or dummy header, or a template that is updated. Hmm, 2 pi or not 2 pi, that is the question, whether tis nobler to take arms against outrageous phase shifts...in this case yes...phase shifts must be conquered! I have correctglobal = true for spatial unwrapping.....but not for temporal.

I can add correct_regions=true, but on some level maybe that's the actual problem...the mean of the 'wrapped' region is not near zero, it is something more like ~3pi/2.

What about weights, my intuition says the weights figure in to the unwrapping... would that look like weights=romeo2?

[korbinian90]

The NIfTI header is from here: https://github.com/JuliaNeuroscience/NIfTI.jl/blob/master/src/headers.jl It has quite many parameters and I'm not sure how to create it. Maybe the NIfTI.jl package has some explanation.
For me it is no issue, since I always start with a NIfTI file on disk, read it in, store the header as e.g. hdr = header(phase) and use hdr throughout my code, when saving julia arrays so that I have the proper voxelsize and stuff.

[drswgw]

Strike That. Reverse It. Here is the MagIm data after separating into gradients and retiming (still single slice). Upshot => low Mag means noisy phase. So there is definitely motion going on, but the phase data is crap in those specific regions (what why how the MagIm data is being modified by motions...well that's interesting).... Tell me about the phase quality thing again. Choices are...mask the %#&^ out of everything...I'm actually OK with that....but would prefer not to throw out data if I don't have to. So need to combine the phase quality map AND the MagIm data.....to come up with something like a quality factor that I can build a proper mask from, i.e. instead of mask the ^$^*^ out of everything, mask like a surgeon with a Quality Factor scalpel. Good news...Romeo is doing as well as it can with crap data....like me, I think Romeo is having a tough time distinguishing a phase boundary from a crap phase region...certainly looked like a phase boundary...

[korbinian90]

I might have an idea what goes wrong. ROMEO has a weight that uses the magnitude ratio of two neighboring voxels to determine how well connected they are (magnitude_coherence). In this case the magnitude drops and the ratio is high -> low connectivity. That might be the reason the unwrapping doesn't want to cross the boundary. So in your case, the magnitude weight might be bad. You can manually activate each of the 6 implemented weights with:

weights = falses(6)
weights[[1,2]] .= true # only weight one and two active
romeo(phase; mag, TEs, weights)

Or you could run ROMEO without the magnitude and see how it looks. There are also two additional magnitude weights, 5 and 6, which are deactivated by default. Here are the 6 possible weights, optionally activated by flags: https://github.com/korbinian90/ROMEO.jl/blob/master/src/weights.jl

Might also be that it doesn't help, just a guess ;)

[drswgw]

Hi Korbinian, Those seem like solid suggestions (the divide and conquer idea was VERY good)...I'm going to resurrect the phase maps...play with/without Mag for unwrapping, and take a peek at weights. None of this juice would be worth the squeeze (i.e. the obvious fix is to get better data), but in this case - experimental sequence development - this is the best data I've got, so I need to make it work. Boundaries....I fell into a burning ring of fire...(two quotes in a row...Willy Wonka (Gene Wilder) and Johnny Cash). Turns out this is of more general utility because later on, reconstruction, a quality factor map would be directly useful. Define "romeovoxelquality" just what does that mean? I am intrigued by 'magcoherence' can you give a more general description....? Normally our Mag data is very spotty, but for this sequence the MagData is not half bad. P.S. Don't hesitate to say...you are using up too much of my time! I enjoy the help...but...

[drswgw]

Super Cool! Vortex rings for fluid mechanics (looks just like the MRI data!)

[drswgw]

Hi Korbinian, Slight sideline - diving a little deeper into slice timing, so using MATLAB to parse the raw dicom files from the magnet (prior to NIFTI ifying) in order to pull the 'AcquisitionTime' parameter for each subimage. NIFTI has some slice timing information but relatively limited. Anyway, seems like the simplest thing to do is to export from Matlab as HDF5. I can already read .jld files into MATLAB as if they were HDF5 (for the most part anyway). Which circles around to MriResearchTools.....reading Dicom, and HDF5 might be a nice addition in capability to consider. For now....Dicom => NIFTI (using dcm2nixx embedded in MRIcroGL) => Julia => Romeo Dicom => Matlab => HDF5 => Julia => Romeo You can probably see where I'm headed with this....MriResearchTools is more useful (to me) if it can handle a slightly wider variety of IO options e.g. Dicom => Julia is a lot more direct.....I'll play with DICOM.jl a bit....but that package seems reasonably primitive at first look. (a lot of the Julia stuff seems reasonably primitive...which can be charming...or just annoying...depending on my mood).

In the end, I am reasonably format agnostic, and semi platform independent, but cutting out the middle-apps seems useful.

In this case, what is emerging as the go to....sort/retime/resync...THEN...unwrap. That's backwards from normal (unwrap first), but the unwrapping is markedly BETTER when the sort/retime/resync options are done first. I think that is becoming more obvious...because the data has LARGE time effects, LARGE gradient effects, and Large signal changes as a result....all of which seem to impact unwrapping directly. So having more capability for sorting/retiming/resyncing is what I am developing now.

All is moving forward....although sometimes it seems like two steps forward, one step back.

[korbinian90]

The dicom reading is quite complicated and vendor specific. I once filed an issue on dicom.jl to read 3D images, but they were not too enthusiastic to reinvent the wheel, as dcm2nii already has all the complicated dicom parse logic and works well. I see that for your case dcm2nii cannot retain the information. So there is actually a use case where it would make sense. If reading hdf5 is fairly simple to add i might consider it.

[korbinian90]

And no worries about wasting my time. However, it might happen that I answer only parts or don't read it for some days, if there is a lot to do. But I think it is very interesting and I like following along with your story ;)

[drswgw]

A few days compared to the general pace of research is lightning fast! I also am enjoying the dialog, it is a pleasant break from chasing code errors, and your suggestions are mostly good.

.jld files are basically an enhanced version of hdf5, so it should be reasonably straightforward to read them in julia...

I will play with hdf5.jl and let you know how it goes....i.e. it doesn't need to be part of MriResearchTools if the hdf5.jl module is any good.

IO => there's a reason that every coding adventure begins with 'hello world'. dcm2niixx is super good.....and it may be that the slice timing in the .nii files is just an abbreviated version of what will come from the more detailed deep dive.....for 1 slice...there is no slice timing...only TR. I guess it's a question of whether TR is really the slice to slice to slice time interval or not. Typically, the true MR delays are not constant but tend to oscillate around some average value...whether that matters? Philosophically yes, practically?

I've always felt that MATLAB's approach to IO is pretty lame (I think because I come originally from basic and C where IO is very logical). Back to being a hacker not a coder.

[drswgw]

Hi Korbinian, OK, the hdf5.jl implementation is slightly awkward (they like to mix pointers and strings in a way that makes me just slightly squeamish, at least I think that's what happening) and the examples/manual could always be better....but it works without too much messing around...so I am happy.

.hdf5 is a convenient way to bring data in/out of julia that is more universal (admittedly a bit more restricted) than is .jld

Things are a bit more broken up now...I went more C style with MAIN, although admittedly I pass forward TO MAIN instead of calling FROM MAIN, but that's the great thing about being a hacker, I'm not constrained to follow rules or niceties.

This is working better & better...LOTS of MRI data in...short code snippets that set the path and specific romeo calls like maxseeds for instance (would like to learn more about romeo weights...what are they, what do they do, how to set them properly) then on to MAIN to do the heavy lifting of resysncing, retiming then unwrapping. In this case the .hdf5 file carries the raw dicom slice times (more detailed than the .nii file e.g. it's in ms since midnight and has a time to go with every slice, every gradient, every delay) comes in as well, along with an external mask if needed. I could conceivably bring that data in as .json, which is slightly more nifti-ish, but .hdf5 seems much more natural for julia. Why there is more than 1 data format is a lot like why there is more than 1 computer language. Note the trend...it's a tossed salad...the data is lettuce, romeo is salad dressing, but there are cucumbers and tomatoes too!

CALL GO.JL => PICK SCRIPT => PASS FORWARD TO => DATA SPECIFIC PREAMBLE

PASS FORWARD TO MAIN (just showing first part for the HDF5 IO)

1st part of the running JULIA output dialog

Another step forward! HDF5 IO.....check....

[drswgw]

Now that HDF5 is working looks like a NEW thing is Siemens => ISMRMRD format for raw Siemens (K-space data). Can you unwrap in k-space? Or is that an image space thing only?

[korbinian90]

I would say unwrapping is an image space thing, but it's only a Fourier transform away ;)

[korbinian90]

This package supports reading of the ISMRMD format: https://github.com/MagneticResonanceImaging/MRIReco.jl (I haven't used it yet though)

[drswgw]

Here's a cool opensource dataset....trying to come to grips with it...it's HUGE https://fastmri.med.nyu.edu/ Uses complex hdf5 (cool huh!)

3 steps...just give me 3 steps, give me 3 steps mister For me it seems more like 3 steps = 2 steps forward, 1 step back. This is mostly a backwards step...but seems like it's to the good...

I like the menu approach, but it was getting unwieldy, so broke into categories first, followed by a directory call into the individual code folders. Basically eliminates the manual entry for code names. Choose category, THEN find the right code snippet...which does the IO setup...then call MAIN.

All aimed at modularizing the code, removing inconsistencies etc. Sometimes it seems exhausting...just building a data pipeline...but in this case it is VERY specific in terms of what it does.

ONWARDS! More code to hack away at. This is when I wish my coding was of a slightly higher caliber...Einstein had his own personal mathematician (Levicivita). His quote was "while I slog through the muddy fields of mathematics, Levicivita rides by on a large white charger"

Right now I'm feeling that with the code...simple stuff...but ohhh, sooo, sloooww! Probably means great things are imminent.

Starting to crank the hdf5 acquire times from dicom via matlab. Getting spoiled...Julia is faster than matlab, especially for loops.

[drswgw]

If you are going to fish for 'AcquisitionTime' or any other single dicom attribute/field using MATLAB dicom this is worth knowing.....

%% clear subfldnm subfldrstrs; for jj = 1 : nsubflds; subfldnm = subfldr(jj).name; subfldrstrs(kk,jj) = {strcat(subcdr,'/',subfldnm)}; if kk == flocs(1) % This natural approach is VERY SLOW! % acqtimes(kk,jj) = str2num(dicominfo(subfldrstrs{kk,jj}).AcquisitionTime);

                            % This uses a MATLAB back door to create a storage object, 
                            % then fills it with just the one value
                                obj = images.internal.dicom.DICOMFile(subfldrstrs{kk,jj});
                                acqtimes(kk,jj) = str2num(obj.getAttributeByName('AcquisitionTime'));
                            end
                        end

%%

In summary, dicominfo is painfully slow, this matlab back door is a LOT faster (still slow but acceptably slow ~ 15 GB of dicom files in ~6 mins.....versus what seemed like an hour or more...don't know for sure...got tired of waiting around, soo...hacked up this fix..as usual largely from opc (other peoples code)...this was via a MATLAB issue resolution page...thank you internet!) The preamble serves to traverse the dicom file(s) tree(s).

[drswgw]

OK, forward. Almost complete chaos again. A total of six (6) gradients now. But the data is of better quality (phantom instead of free fluid flow). Just starting to pull in the super-detailed hdf5 timing data as opposed to the rougher .nii timing. A bunch of work/coding ahead, but romeo is behaving pretty much as expected...OK when the phase jumps are small...less well when the phase jumps are big, or the amplitudes are weak. Will look at tweaking the romeo weights as needed, just not there yet. Romeo is definitely happier now that things look more brain-like (i.e. single region). Look at that beautiful phase data....all six asynchronous gradients, sixty different acquire times worth!

[drswgw]

Small breakthrough/workaround. Hdf5 and .jld CAN handle high (dim > 4) dimensional data even if Romeo and MricroGL and by inference NIFTI cannot, and/or prefer 4D. There's nothing wrong with breaking things down to 4D for unwrapping, it's just 'extra'. The whole IO process is MUCH simplified for this data if the lingua fraca is 6D (x,y,slc,grdA,grdB,phs) Same argument for 5D (x,y,slc,grd,phs) It's a GIFT from hdf5 and.jld that high dimensionality is native. And a small criticism for ROMEO and NIFTI that those are capped at 4D.

At least that's the tack I'm taking currently to try and wrap all this up.

[drswgw]

Mixed bag for the unwrapping over multiple gradients..... Definitely divide and conquer is key to the whole thing, without that...nothing works. Some is definitely on my side, there is a residual delay that I am just realizing that I am not dealing with (shows up as residual stripes). I think that is also messing up the temporal unwrapping to some extent. This seems like a 'new' issue, where LARGE regions have LARGE changes in phase....comparatively un-brainlike where the norm is to have largish regions of ~constant phase offset by smaller regions with modest to substantial phase wrapping. Two specific cases seem difficult for Romeo to handle - very small regions with large changes in phase (looks like noise) or very large regions with large changes in phase (looks like it's already unwrapped...but it's not). Will play with fixing the residual delay(s), then I am thinking I will need to dive into the unwrapping weights. Romeo in general is rather timid about unwrapping large regions. I think I need it to be just a tad more adventuresome...back to setting a default 'area size' for unwrapping. I think that's different than maxseeds, because maxseeds is just a possibility. This is more like a constraint....a minimum or maximum area/volume where phase consistency is expected. For this phantom data, the consistent regions are relatively large, but the phase changes are BIG.

THIS ONE LOOKS OK (weaker gradients), still has some residual timing and delay issues that I am working on

THIS ONE IS STILL WRAPPED (stronger gradients) location as indicated or more obvious in other slices. Romeo seems to be struggling to get up to |phs| > 3*pi for large areas

And just for FUN! RAW DATA: VERY STRONG GRADIENTS Once I get the delays fixed would like to be able to unwrap this - a lot more work for unwrapping but the SNR is definitely better at higher gradient amplitudes.

[drswgw]

Progress is slow because errors in time look a lot like errors in phase/unwrapping, which is...confusing.

Trying a loopback approach now. Basically have time based signals going IN to the magnet data files (but those are not gates, just markers). So adding a second set of time based signals coming back FROM the magnet (those now are specific markers for the gates). Such redundancy should not in principle be needed at all (meaning the dicom times should? reflect where the gates are), but in practice am finding small delays that accumulate over time and those just &^$&^% everything. Once those timing issues are fixed, will be back knocking for more help with the unwrapping. I think it's a vulnerability of asynchronous acquisition (i.e. 'ungated' or 'semi-gated'), imagine a race where all the swimmers use the times logged by hand on a stopwatch by their specific coaches - that's kind of where I'm at - haven't figured out the TIMEX part yet, so it's close...but not quite good enough yet.

Give some thought to > 2 point constraints meaning....nearest neighbor in space (or time) is easy to code but not entirely useful (at least not for this). The poster child of such things is boxcar averaging, the minimum example is where each point is replaced by the average over that point and the 2 nearest neighbors (that's the 1D example) i.e. it's a 3 point (or 5 or 7 etc.) constraint. In general I have 10 time points (well more than that but 10 after deconvolution/resyncing/retiming), and an 80x80 or 96x96 base matrix, so nearest neighbor continuity is good, but so is 'boxcar' type continuity, at least that's the hope. Meaning I am willing to sacrifice some level of sharpness in favor of more robust unwrapping (it's not ideal, but as my very wise MRI mentor (Kim) once told me: EVERYTHING in MRI is a TRADEOFF). Perhaps there is a way to establish smoothed constraints, that map back to the raw data? (semi-sort of like an unsharp mask). It's always nerve-wracking because the perfect unwrapping is content neutral...BUT...that may be unrealistic. Instead some form of consensus reality is probably needed to unwrap 'effectively'.

[korbinian90]

Good luck with fixing the timings! If there are still weird wrapped slices afterwards, maybe there is a way to remove the "gradient artefacts" before unwrapping. ROMEO is only allowed to add/subtract multiples of 2pi, so it might be tricky if there are slice jumps from differing gradients that are not multiples of 2pi.

[drswgw]

Noted, not 2 pi. Just talking with a collaborator and their 1st level attempt uses derivatives to try and detect phase jumps. I think you have similar options?

Making progress with the timing, oh so slowly. In summary, Siemens uses multiple clocks, not all of which are synced to each other, or more precisely...how to sync them is not entirely obvious. Am more and more coming around to getting everything onto a single clock. Sadly that clock probably is not based on the dicom times, or at least not currently. Sad, because the Dicom clock is time stamped everywhere, but seems to be out of sync with what I need to do. Ultimately would like to get my hands on the FAST siemens clock (us) but it is buried VERY deep and so far is eluding me. But making progress with other clocks (tics of 2.5 ms) and that is finally starting to go OK.

More unwrapping soon!

[korbinian90]

I don't think I have an option in ROMEO that would help. There is linear unwrapping (--wrap-addition), but that still only adds/subtracts 2π. It can help to unwrap very steep phase (in 3D), but you have phase jumps that are not multiples of 2π.

You can look at the bipolar correction here: https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.26963 , search for "bipolar aspire". Maybe something similar can be applied for your data set? I have multiple (slightly differing) implementations of this, one is here, called bipolar correction: https://github.com/korbinian90/MriResearchTools.jl/blob/master/src/mcpc3ds.jl

I'm not sure if it helps, but it might spark an idea for correcting these phase artefacts.

[drswgw]

New Questions! Resurfacing from timing adventures (long and drawn out and not so exciting adventures). Spatial unwrapping addresses spatial continuity within a slice and presumably somewhat between slices? Temporal unwrapping addresses continuity along the time dimension IN ADDITION to the spatial unwrapping constraints? In my case 'da bear' is the slice to slice discontinuity. Some of that is real i.e. from timing and sampling issues, but I wonder if some of it is from unwrapping creating discontinuities between slices? Posed as a question +- 2pi, but 'who' is 'zero' or 'average' then? Back to the foothills of the mountain, is the height of the foothills preserved from slice to slice? e.g. slices through a mountain....do I get the whole mountain back...with fidelity? or are there artifacts from the unwrapping? Is some of the mountain lost?

Sagittal version of the same data - not so good slice to slice - divide and conquer approach: resysncing/reordering BEFORE unwrapping (temporal)

Axial - single slice - looks pretty good! divide and conquer approach: resysncing/reordering BEFORE unwrapping (temporal)

All of this is like an iceberg. A lot under the surface just to get to here...which is only about 2/3 of the way.

Here is some raw data (same data set but not same data as above - i.e. different gradients) to give a sense....

[drswgw]

Current romeo parameters....note that divide and conquer means...for the most part...that although unwrapping may span multiples x 2pi, the time to time changes should be less than 2 pi...what is probably needed though...may very well be the option to 'follow' say a pi/2 step across 10 phases which adds up to => 5* pi.....what happens then? (that's an exaggerated case, but it is possible, particularly at higher gradient strengths...actually that may be better posed as a question...is it possible, does it happen? don't know for sure....for phantoms at high gradient strength => probably, but for brain...not really sure).

pre-parameters:

runtime-parameters:

[korbinian90]

There are no artifacts from unwrapping. You can do rem2pi.(unwrapped_phase) and you get the original wrapped phase back. If you look at original_wrapped_phase .- rem2pi.(unwrapped_phase) this should be 0 (maybe some floating point errors < 1e-10) So you can be certain that these small changes (<2π) are from the acqisition and not from unwrapping.

The romeo temporal unwrapping should do exactly what you are asking for in the second part. It takes the first unwrapped 3D volume and uses it as a template to unwrap the next one by scaling it approprietly (using the echo times). This assumes that the phase relation is linear with echo time and there is no phase offset (phase at time 0). But it is still working if that assumption is slightly wrong, when there is less than π deviation from the linear phase evolution. So if your phase is off by π/2 for each subsequent volume that should still work. You see that problems arise, when individual noisy voxels are off by 2π or when wraps emerge in your image. This should be very clearly visible if it occurs.

[drswgw]

So TEs = 1's should give uniform scaling throughout? That is ideal!

The alternative, say for a 2nd echo, would beef up the 2nd echo according to? Wouldn't you need to know T2 (and or T1, T2* ) to scale it based on TE? Still using e^-kx ~ 1 - kx gives a reasonably good linear approximation (for small x). What happens if one sets the TE values to TEs = 0's?

[korbinian90]

ROMEO only scales the phase temporarily for unwrapping the next volume. The resulting unwrapped phases are not scaled. This is the basic temporal unwrapping code:

for ieco in 2:length(TEs)
    iref = ieco - 1
    refvalue = wrapped[:,:,:,iref] .* (TEs[ieco] / TEs[iref])
    w = view(wrapped,:,:,:,ieco)
    w .= unwrapvoxel.(w, refvalue) # temporal unwrapping
end
unwrapvoxel(new, old) = new - 2pi * round((new - old) / 2pi)

The refvalue is the temporarily scaled one, it's only used as a template and afterwards thrown away. The w is the unwrapped output (a view into the array wrapped).

The phase should evolve linearly with echo time (disregarding weird effects like partial volume): phase = B0 * TE * some_constant. There is no direct involvement of T1, T2 ...

[drswgw]

Hi Korbinian, Well, I'm sure you figured this project was long done. Quite the contrary, it's just getting interesting. Here is some poor man's phase unwrapping where the DC term in the phase domain is subtracted from the raw phase data. So those regions where the resting coil phase is preserved and presumably far enough from pi or -pi, are 'unwrapped' and registered relative to 'zero phase'. The regions where the data is crossing over pi/-pi are then thrown into total chaos. Why do I care? Because it gives a quick and easy check on Romeo. The power of Romeo is in unwrapping the 'pesky' regions where the phase is crossing over and back the phase boundaries. Nevertheless, my remaining questions/checks and balances center on things like residual 2 pi offsets, and or drift in the baseline phase. As it turns out, we actually don't care so much about the DC term (that's the baseline phase) but rather the time varying AC term...that's our bread and butter. These don't match up exactly but should give the general concept. Of course this poor man's approach is singularly bad where there are BIG phase changes, but our phase changes are mostly...smallish.

Mean data (resting phase)

Subtracted Data (raw phs - mean phase => poor man 'unwrapped')

Mainly, the 'right answer' is obtained, so now I should be able to find slices etc. where Romeo has a residual 2pi offset and so on just by minimizing against this poor man's unwrapping (needs some proper masking to do that). Anyway, should be able to sort Romeo outputs against poor man unwrapping and flag those slices/regions where chaos still reigns.

Have much more accurate timing info now (i.e. abandoned Siemen's logs, Dicom times etc. almost entirely and am doing trigger logging and heart logging on an A to D instead). Definitely better! Sometimes it's just better to do it yourself => no middle man.

I think you have a parameter that generates a local zero...the question now is whether your local zero is our local zero = voxel by voxel DC term in the time domain only. Admittedly head motion is a bummer, but I think head motion is always a bummer.

COLORMAP = MAYO_AWAVE

[drswgw]

Hi Korbinian, Not strictly belonging here (Julia side)...but playing with calling from Matlab (R2022) using MRItools 3.5.2....I get this bug...multiple additional parameters don't run only 1. As far as I can tell the compound string isn't being recognized as 'commands'. Everything else is working...well except for setting more than 1 additional parameter...kind of a major problem.

% ROMEO unwrap from matlab side clearvars % addpath('NIfTI_20140122') cd('E:\Matlab_Code_Working\INTRNSC_CRTS\MRE_PTNT4_SG_7_21_22_FLAT') tempdir = pwd;

%% Phase phase_fn = 'PHS_la1D_P.nii'; %'/path_to_data/Phase.nii"; phase = load_nii(phase_fn).img;

%% Optional Parameters mag_fn = 'MAG_la1D_P.nii'; %'/path_to_data/Magnitude.nii"; mask_fn = 'Mask_LoRes.nii'; %'/path_to_data/Mask.nii";

parameters.output_dir = fullfile(tempdir, 'romeo_tmp'); % 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(mag_fn).img; parameters.mask = 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 = [2,2,2]; %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 = '-i'; % '-q -i'; % settings are pasted directly to ROMEO cmd (see https://github.com/korbinian90/ROMEO for options)

%%% &&& CAN ONLY USE 1 parameter in the line above...or the matlab call crashes &&& %%%%

%% 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(phase_fn);

%%% Otherwise this code runs

The unwrapping is only so-so....ROMEO is struggling at the phase crossover region. Using the poor man's unwrapping I can fix the 2pi offset, so that's not a big deal...but I need the finer control...i.e. call more than one additional parameter at a time. This may mean back to Julia...which always works. The Matlab script is working better than before...but still not fully functional...at least not for me. All of the named parameters work OK, but the collected commands do not.

Raw phase...80 phases...

Unwrapped phases - note the artifact at the crossover region (and the 2 pi offset but that's less of a problem).

here it is from the other side....this is all spatial unwrapping I think...can't access higher functions from matlab without being able to call them....

[korbinian90]

Thanks for posting bugs. I'll try to reproduce and fix this next week. I was quite busy the last weeks so I didn't catch up yet with the rest of your posts

[korbinian90]

The matlab caller should now be fixed in version 3.5.3 (currently compiling)