Determine Effective Echo Spacing for Philips

[neurolabusc]

The current stable version of dcm2niix does not populate the BIDS field EffectiveEchoSpacing for Philips data. While this is often not a crucial parameter, the sdcflows repository does attempt to estimate this. It would be great to validate this. I suspect the sdcflows formula does not account for partial Fourier acquisitions. Tag 0018,9081 and 2001,1019 only reveals if pF is on/off. This could be inferred by knowing the echo train length (from 2001,1013), the SENSE factor (2005,140f) and the un-interpolated matrix size. @drmclem do you have any hints on determining partial Fourier? Ideally, it would be nice to have a simple validation dataset from someone with access to Philips equipment (e.g. scans where interpolation, partial Fourier, and SENSE are varied).

[neurolabusc]

Closing this issue. I do not have access to Philips hardware to validate solutions and no Philips users have stepped up. The upcoming release will not populate BIDS EffectiveEchoSpacing and TotalReadoutTime for Philips Dicom images. If someone wants to resolve competing methods this can be added in a future release.

(f your readout time is identical for all acquisitions you don’t necessarily have to specify a valid value in this column (you can e.g. just set it to 1). So this is a 'want' not 'need' feature. Without being confident of the optimal solution, it is best not to report this tag.

[neurolabusc]

For Philips, (0018,0091) EchoTrainLength does not take partial Fourier into account. For example, one a Siemens scanner for an image with a base resolution of 90, the ETL will be 90 for full Fourier, 79 for 7/8th pF and 56 for 5/8th pF. In contrast, Philips will report 90 for all these values. Therefore, dcm2niix can not populate the tag EchoTrainLength for Philips data where PartialFourier (0018,9081) is [YES].

[drmclem]

Hi

I think there may be some historical changes here. I haven't tested exhaustively but the following seems to be consistent. I have checked these figures on a latest r5.6 system basic EPI scans

In the above table the first 4 columns are the information available to the user on the information tab - Sense factor, Halfscan factor Water fat shift and EPI Factor.
As you can see - Sense (parallel imaging) changes the displayed EPI factor but half scan does not.

Using some internal tooling the actual echo spacing and actual EPI repetition number are displayed.

Using the formula

EffectiveEchoSpacing = (((1000 * WaterFatShift)/(434.282 * (EPI Factor + 1))

produces echo spacings which closely match the internal values.

The logic for this is that the WFS value is our surrogate for distortion in combination with the EPI factor - half scan does not substantially change distortion but Sense does - so the combination of EPI Factor and WFS is sufficient to calculate the echo spacing.

I think there may be some historical confusion as in some very early releases the Sense factor did not change the EPIF on display and hence some formulas include a /(acceleration) factor but I'm not sure when this changed.

The magic number by the has some varaition in different sources but should the the fat water shift in Hz (3.4 * 0.000001 * 127.73e6 )

In terms of populating the echo train length value - thats no going to be easy to get the correct value with halfscan but if this is empty - I wouldnt want to prevent users from not being able to use distortion correction software if it needs a value here to calculate the actual number of interest - the echo spacing.

M

[neurolabusc]

@drmclem this is timely. @parekhpravesh and Paul Morgan have each sent me Philips validation datasets. I believe my formula is the same as yours, with 2 modifications:

1. Using the ImagingFrequency (0018,0084) allows the formula to work for scanners that are 1.5T or not precisely 3.0T (e.g. the Ingenia CX reports 127.775 MHz). I prefer this to Magnetic Field strength (0018,0087), as the latter is provided with low precision.
2. EffectiveEchoSpacing and TotalReadOutTime in seconds, as expected by eddy and BIDS.

EffectiveEchoSpacing = WaterFatShift/(ImagingFrequency*(EPI_Factor+1)*3.4))
TotalReadOutTime = EffectiveEchoSpacing * (ReconMatrixPE - 1)
WaterFatShift = (2001,1022)
ImagingFrequency = (0018,0084)
EPI_Factor =  (2001,1013)
ReconMatrixPE = Reconstructed Matrix Pixels in Phase Encoding Direction 

The TotalReadoutTime is in BIDS format (FSL definition) so this is not the actual readout time, for example for partial Fourier it is the time if all lines were acquired, and for interpolation it takes into account the number of lines in the output image.

@drmclem I would be grateful if you could look over the following Excel spreadsheet created by @ parekhpravesh and me. There are three worksheets, for no acceleration, partial Fourier, and SENSE. It would be great to have your thoughts on the red column (TotalReadoutTime):

Philips_IngeniaCX_EPI_Validation_Guide.xlsx

As a final check, I am running all images with Eddy, to see if the -fout image is calibrated across modifications of SENSE, pF, interpolation.

Thanks to a great team effort, it appears we are close to closing this issue.

[neurolabusc]

@drmclem @captainnova and other Philips users.

@parekhpravesh and Paul Morgan have kindly provided open source validation datasets. I have created a Open Science Framework resource named Calculating TotalReadoutTime for Philips MRI that includes the DICOM datasets, Excel spreadsheets (with critical variables clearly listed), Python scripts, generated field maps (in Hz) and a brief overview of the problem.

Based on this, my current working hypothesis is that the formula I describe above is robust. However, I would be grateful if others could either validate my formula or describe a better formula.

For those who are not members of an OSF institution, here is a public link

[neurolabusc]

I have updated the sample data, validation scripts and description. The values for TotalReadoutTIme and EffectiveEchoSpacing are our best approximation for the FSL definitions for these values.

ActualEchoSpacing = WaterFatShift / (ImagingFrequency * 3.4 * (EPI_Factor + 1)) 
TotalReadoutTIme = ActualEchoSpacing * EPI_Factor
EffectiveEchoSpacing = TotalReadoutTime / (ReconMatrixPE - 1)
WaterFatShift = 2001,1022
ImagingFrequency = 0018,0084
EPI_Factor = 0018,0091 or 2001,1013
ReconMatrixPE = 0028,0010 or 0028,0011 depending on 0018,1312

As @parekhpravesh notes:

Unfortunately, Philips does not allow to set echo spacing as a parameter; it gets determined based on other factors. For example, for TSE scans, the help manual states:
Echo spacing is the interval between two echoes
should be kept relatively short in order to prevent blurring
is no parameter, but is determined by the parameters TSE Factor, TE, Profile order, or can be defined with the asymmetric profile order (parameter TE spacing)
Another strange behaviour on the Philips console is that the acquisition matrix size in the phase direction is not exactly what you would want it to be. For example, the baseline scan has an acquisition matrix of 64x63 (MxP) rather than 64x64; this is a persistent behaviour. Only with a few "magic" numbers (for example a SENSE factor of 1.45 or 2.2) does the acquisition matrix become symmetric. This would explain, for example, why the EPI factor (and therefore ETL) do not exactly equal size of acquisition matrix in the PE direction. This is also revealed by the scan percentage property in the exam card which shows values like 98.4375%. There are rounding off issues at each stage (revealed by ExamCards). For example, the acquisition voxel sizes are 3.29787493 and 3.29999995 rather than strictly being the same size (noSENSE_noPF_noPO_acqrecon_AP_A).

Therefore, as far as I can see, the actual acquisition matrix would be 64scan_percentage. In case of phase oversampling, 50mm in A and P directions (noSENSE_noPF_50APPO_acqrecon_AP_A), the acquisition matrix in the PE direction would be 6498.9583359 = 63.333334976. Thus, the EPI factor would likely be 63.333334976 + ( 63.3333349760.5) = 95; similarly for 100mm phase overasampling, the acquisition matrix would be 6499.21875 = 63.5and EPI factor would be 63.5 + 63.5 = 127. These are the values reported in the ExamCard. Similarly, for SENSE factor of 2, the scan percentage is 96.875, therefore the PE matrix size is 62, giving an EPI factor of 31. The scan percentage does not change with PF and neither does EPI factor.

The fieldmaps generated by FSL's TOPUP look pretty similar across settings, but not identical. This could reflect a combination of an error in the TotalRedoutTIme formula or changes to the Fieldmap due to Philips sequences. The latter is clearly at least part of the issue, as the differences are not simply scaling issues. As a general rule, one should keep all parameters identical except phase encoding polarity for images provided to TOPUP. When this advice is obeyed, the TotalReadoutTime has no impact. Philips users should be particularly cautious if they ignore this guidance.

Unless someone provides a better solution, I consider this issue closed. Therefore, I plan to make a new stable release of dcm2niix next week.

[mharms]

If you look at the estimated field map images, this is clearly still not right for all situations, in my opinion. In which case, I would argue that it is dangerous to populate the sidecar json with these estimates for TotalReadoutTime and EffectiveEchoSpacing since users are naturally going to assume that they are correct.

On the OSF page you mention that "scanner shimming across these settings" may be a factor in the inconsistent maps, but I thought that the latest round of data was acquired on a phantom in a single imaging session with the same shim setting throughout the entire session...?

[parekhpravesh]

I just checked all the ExamCards. The preparation phases for all the scans are set to auto i.e. no reshimming would have happened between the protocols in the SE stack. And yes, they were acquired in the same session

[neurolabusc]

Pragmatically, the results seem in the same range of magnitudes, but the morphology looks different. This seems to suggest that adjusting the TotalReadoutTime is not sufficient to fix the variance. In contrast, the Siemens validation set does have strong morphological correspondence across the sequences, providing strong confidence in the scaling values.

It seems like dcm2niix can either populate the BIDS tags with our best guess, or we just provide no information reflecting our uncertainty in the solution.

Regardless, Philips users need to be especially careful when wanting to provide inputs to Eddy where any parameter beyond phase encoding polarity is changed.

[mharms]

My perspective is that an extremely small percentage of Philips users will ever be aware of these "warnings", in which case I would tend toward not providing the values, consistent with the uncertainty of the situation.

[neurolabusc]

Understood. Since it does not appear we have data to resolve this completely, and since it is typically not important but users expect these values, I propose we populate BIDS JSON fields with the Bogus prefix. Therefore, scripts will not automatically detect the value. Advanced users can see the Bogus title and find this issue. We can revisit this in the future if anyone makes a breakthrough. I really don't want to hold up a new stable release, as the development branch has many improvements over the current stable release.

BogusTotalReadoutTime
BogusEffectiveEchoSpacing

[parekhpravesh]

I agree with your suggestion. I think its a good idea to still populate these values (as Bogus fields) and let the user decide whether they think the method works for them or not. Perhaps populating this with Bogus field will also attract attention of other Philips users who might have additional insight.

[mharms]

Sounds good to me. Unproven or Questionable would be two possible other options, if you like either of those better.

[neurolabusc]

@effigies do you have a preference for a prefix for a BIDS tag when we are unable to confirm the accuracy? Bogus, Unproven or Questionable?

[effigies]

I don't and I can't think of any justification in BIDS that would lead me toward any particular one. None of them seems likely to ever become a prefix of some actual metadata field.

[drmclem]

Hi - I do feel these are a bit judgemental - how about Estimated ?

[neurolabusc]

Good point. This sort of makes sense, as FSL definitions for TotalReadoutTime and EffectiveEchoSpacing both often differ from the true absolute values (this explains why some Siemens JSON files show DerivedVendorReportedEchoSpacing that differs from EffectiveEchoSpacing). Further, it is possible that the variations we see between these sequences actually reflect clever tweaks Philips is applying to optimize each sequence. Therefore, the fact that we can not derive the perfect values expected by FSL could in theory reflect strengths rather weaknesses of these sequences.

[mharms]

I don't think that Estimated necessarily connotates the necessary caution. For all practical purposes, the values for TotalReadoutTime and EffectiveEchoSpacing are essentially defined in BIDS based on the values expected by FSL, so if we aren't generating those accurately per the FSL expectations, then we don't have them right.

[drmclem]

Hi - I'm not familiar with the calculations needed for the other manufacturers - but I thought the FSL definitions were pragmatic rather than actual (i.e. could either be based on the actual read out time, or a fictitious one based on the reconstructed matrix). The calculation in my post above seems to be valid for echo spacing and presumably total equivalent readout would be as good as FSL expects. Seems a bit harsh to single Philips data out as "bogus" - (disclaimer - I am a Philips Scientist :-) when others seem to be supplied on a "pragmatic approach to amibiguous data".

My suspicion is that there is high variation in the accuracy to which you can derive the field map from the associated distortion over the range of parameters in the Ingenia test set - between sense 4 and No sense, there is almost a factor of 5 difference in the distortion sensitivty.

Looking at Mharms data here https://github.com/neurolabusc/dcm_qa/tree/master/In/TotalReadoutTime the consistency has been based on a very different assesment.

A better comparision would be against an anatomical outline - i.e does the distortion correction based on the derived values actually work, and from my users this appears to be the case.

[mharms]

Well, as @neurolabusc has repeatedly pointed out to users, for most practical use cases in topup, the actual value for TotalReadoutTime isn't necessarily important for distortion correction (given certain, albeit important, assumptions detailed elsewhere, including on the FSL Wiki page for eddy). However, if you want the actual field map estimate (in Hz) to be a correct estimate of the field itself, then getting TotalReadoutTime correct is critical.

I'm not as familiar with fugue based distortion correction, but for that, I think that getting EffectiveEchoSpacing correct is important as well.

I agree that Bogus seems a bit harsh, which is why I had suggested Unproven or Questionable as alternative prefixes.

[neurolabusc]

@drmclem - I appreciate you are explicit about your potential conflict of interest, but once noted I think it is also appropriate for you be an advocate for your employer.

@mharms I am actually fine with Estimated as long as it is not precisely the BIDS key, someone is going to have to make an informed decision regarding how to proceed. For most users, it does not matter. The situations where it would matter reflect poorly designed sequences, and might be perilous even if the TotalReadoutTime is set correctly. In other words, I am not convinced TOPUP would do a good job correcting an image pair where the P->A was unaccelerated and the A->P was SENSE x2 (with much less distortion). I guess one could test with your Siemens sample, but I am not sure we could test all use cases.

From this perspective, a logical approach might be to just have the user (or sdcflows) use the EstimatedTotalReadoutTime but expect that subsequent software (e.g. sdcflows) will generate an error if this value differs between an image pair.

[mharms]

@drmclem, if you work for Philips, can you provide any insight as to what we might be missing and why the field maps that Chris calculated using topup don't all have the same morphology?

[drmclem]

Hi,

I think the biggest issue is that the phantom data is a bit extreme - coupled with the fact that it is higher SNR than brain data leads to interesting effects at the edges when the distortion is large due to bit reolution isssues.

To illustrate - I looked at the phantom set for the large WFS (series 4,5) -

The phantom images are pretty crazy and topup fails to correct them so the field maps are going to be off.

Compared to the head images with similar WFS (but different matrix)

Things are much more well behaved.

Also - I am not familiar with the FSL field output but I guess the field predicted outside the head is just some extraoplation of coefficients ? in which case for a good comparison shouldn't these be masked to head signal locations ?

I don't have access to scanners at the moment but a better validation set would be to collect on a head I feel.

M

[mharms]

I guess its possible that topup isn't well-behaved for some reason, although we didn't have that problem with the Siemens data, so the question would remain why it is having a particular challenge with the analog from Philips. Also, IIRC from the images that @neurolabusc distributed, the morphology of the field maps differed within the phantom, not just outside it.

The problem with using a person is that they are likely to move, esp. over the potentially longish session over which all the acquisition variants were collected. That movement will change the field, and thus add "error" to the metric that we are using here.

[drmclem]

Do you have an example Phantom image with the same distortion sensitivity from Siemens ? Without looking at this in detail I would suspect the phantom support - was this a perspex phantom holder? This can have very large susceptibility changes at the points of contact with the phantom.

[neurolabusc]

The Siemens TotalReadoutTime validation dataset is part of dcm_qa.

[oesteban]

Sorry for unearthing this issue again.

Coming back to it I noticed one more reason to get it right:

@mharms said: the actual value for TotalReadoutTime isn't necessarily important for distortion correction

It's not, so long:

i - The EPI up/down blips you supply have exactly the same total readout time (it is possible to have uneven EPI sequences for some reason) for the estimation of the fieldmap.
ii - The target EPI you want to correct (with applytopup) ALSO has the same total readout time. In this second step where the image is unwarped, knowing the effective echo spacing or total readout time is fundamental always (yes, fugue does need it).

This is to say that I agree with Michael that it is really important for Philips to help their users get accurate values. Has there been any relevant development on the matter? @drmclem @neurolabusc .

[neurolabusc]

@oesteban the Philips DICOM images do not contain the required information. This is a limitation of the Philips data, not dcm2niix. Therefore, I am unable to resolve this.

@sandeepganji is leading two independent paths to aid Philips users:

1. While the Philips DICOM images do not contain the required data, most of the information is preserved in the Philips DICOM exam cards. These non-image files are often discarded by PACS and store data in a proprietary format. However, they do provide a mechanism to extract required meta data.
2. Sandeep is developing a list of desired meta data to lobby for future Philips software releases.

I urge you to contact Sandeed directly. You can help him curate his wishlist for new meta data. Also, you can help him identify Philips users who can work together to lobby Philips.

[oesteban]

Thanks much!

[mharms]

Now that I'm actually working with some Philips (5.7) data, I went went back and looked at this issue a bit more closely. Something fundamentally doesn't make sense to me. Both in the examples shown previously by @drmclem (for "EffectiveEchoSpacing", under the Formula column: #377 (comment)), as well as in the spreadsheet Philips_IngeniaCX_EPI_Validation_Guide.xlsx (see #377 (comment)), the computed "EffectiveEchoSpacing" is not changing when some in-plane acceleration (SENSE) is added. It seems like that can't be correct. (The reason it is called "Effective" Echo Spacing is to reflect the fact that the value should at a minimum change when in-plane acceleration is used).

Am I missing something in reviewing this (complicated) issue?

[mharms]

Ahh, I just noticed that as of the March 23, 2020 post (#377 (comment)), you changed the calculation, such that what you were previously calling EffectiveEchoSpacing became ActualEchoSpacing, and now use the following formulation instead:

ActualEchoSpacing = WaterFatShift / (ImagingFrequency * 3.4 * (EPI_Factor + 1)) 
TotalReadoutTIme = ActualEchoSpacing * EPI_Factor
EffectiveEchoSpacing = TotalReadoutTime / (ReconMatrixPE - 1)

whereas the two comments and the spreadsheet that I cite above were created before that date.
So, I guess you never created an updated spreadsheet based on the presumably current formulas?

It's easy to miss that the stuff before March 23, 2020 in this Issue was superseded by revised formulas. Hopefully this will help others to avoid making the same mistake.

[mharms]

FWIW, it would be good to link to #377 on https://github.com/rordenlab/dcm2niix/tree/master/Philips, since it provides a more complete, and more recent discussion of the issue than the current link to nipreps/sdcflows#5

[mharms]

That said, taking another look at the numbers in that spreadsheet, and revising them for the March 23, 2020 equations, I still have some questions. For example, Partial Fourier should have no impact on the amount of distortion, and thus no impact on the computed values for TotalReadoutTime and EffectiveEchoSpacing. But that doesn't seem to necessarily be the case.

For example, Series 201 (noSENSE_noPF_noPO_acqrecon_AP_A) on the "PartI_noSENSE_NoPF" tab has EPI_Factor = 63 and WFS = 12.6507. The same series on the "PartIII_PF" tab has EPI_Factor = 63 still, but WFS = 16.9547, and PartialFourierFactor = 0.698. Why did WFS change upon introduction of partial fourier? Of note, the ratio of the two WFS values (16.9547/12.6507 = 1.34) is strangely close to the inverse of the PF factor (1/0.698 = 1.43). Which makes me wonder if Philips is adjusting WFS by the PF factor (@drmclem)? And if that's the case, then at least for the situation of partial fourier, the current equations wouldn't be accurate.

[neurolabusc]

@mharms my sense is that we close this issue. The current Philips DICOMs are underspecified to determine TotalReadoutTime as specified by BIDS and FSL. This is a limitation of Philips DICOM images, not dcm2niix and therefore is not actionable.

[mharms]

That's fine. Although I'm not sure if it is a matter of the Philips DICOMs being underspecified, confusion about what exactly WaterFatShift is intended to represent, or just that the current equations/logic would need to be further refined upon better understanding of exactly how Philips is populating the relevant fields.

For anyone wanting to take this up, I believe that the current calculations are fine regarding use of in-plane acceleration (SENSE), but I'm not convinced that TotalReadoutTime (and thus EffectiveEchoSpacing) is correct for the cases of Interpolation and Partial Fourier.

In particular, given that WaterFatShift doesn't change when you add Interpolation, the interpolation factor probably needs to be accounted for explicitly as part of the calculation. (e.g., I believe that TotalReadoutTime should double if Interpolation = 2).

And the interpretation of WaterFatShift in the presence of Partial Fourier is ambiguous to me based on the data in the aforementioned spreadsheet. Namely, if you compare its PartI_noSENSE_NoPF and PartIII_PF tabs, for the FatShiftDirection = L data, the WFS values are basically the same for a given set of "PhaseOversampling" values. But that's not the case for the FatShiftDirection = A data. In particular, the data in those tabs with FatShiftDirection = A have similar WFS for the EPI_Factor = 127 data, and again for the EPI_Factor = 95 data. But the WFS differs across tabs for the EPI_Factor = 111 and 63 data. Perhaps that discrepancy arises because the WFS wasn't explicitly matched when setting up the acquisition, when in principle it could have been matched...? Given the amount of data collected, I could easily see that happening. @parekhpravesh: If you generated that data, do you have any insight that you can recall about how you went about setting WFS for all the various series that you acquired?

[neurolabusc]

@mharms I wish to close this issue to push a new release. Do you want to suggest a modification to the Philips-specific handling by dcm2niix page (even if it is only a link to this issue)?

[mharms]

Maybe something like the following for the relevant paragraph in the Philips README?

Another value desirable for susceptibility distortion correction (e.g., TOPUP) is the "TotalReadoutTime". Be aware that the FSL definition is not intuitive for scans with parallel imaging, phase oversampling, partial Fourier, interpolation, etc. It has been challenging to establish and validate the correct equations to calculate "TotalReadoutTime" on Philips in the presence of all these possible factors, and the issue is further complicated by the fact that Philips have may changed the manner in which it reports the relevant information over time (#377 (comment)). For this reason, for Philips data, dcm2niix currently reports EstimatedTotalReadoutTime and EstimatedEffectiveEchoSpacing (note the inclusion of Estimated as part of those variable names). See #377 for a fuller discussion. However, it is also relevant to note that getting "TotalReadoutTime" correct is not critical to the TOPUP correction if the readout time is identical for all the relevant images (i.e., the spin-echo EPI images to estimate the field, as well as any subsequent images being corrected). If that is not the case, or if you want TOPUP to return a valid estimate of the field itself (in Hz; --fout flag), then getting "Total Readout Time" correct is important.

[parekhpravesh]

Hi! Regarding the water-fat shift, if I remember correctly, we did not set this parameter at any point...it was calculated automatically by Philips and displayed as part of the ExamCard summary. I do not recall seeing a way to set WFS to a given value. Digging through old emails, I see in one of our email exchanges ["Study Times", March 17, 2020] that we were picking up the WFS value from these ExamCards directly, which further convinces me that it was not a setting we could manipulate explicitly (at least not for the sequences we acquired)

[parekhpravesh]

Digging through the emails, I see that the last point we had discussed over emails (just before the pandemic lockdown) was that the Philips help documentation states that "echo spacing" is a paramater that is displayed on the console in case of a (multishot?) TSE sequence. However, despite that statement, the console summary did not show the parameter. Perhaps it is possible now with the newer software releases? I, unfortunately, do not have access to a Philips scanner anymore.

[mharms]

If I understand correctly, you should have been able to use Water-fat shift = user defined, and then set the value you want to use (subject to some minimum of course).

[neurolabusc]

@mharms I am closing this issue. The documentation is updated with your suggestions. Thanks!