Skip to content
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

Alternative approach to correcting gantry tilt #253

Closed
alund opened this issue Dec 18, 2018 · 12 comments

Comments

@alund
Copy link
Contributor

commented Dec 18, 2018

I would like to suggest an alternate method for determining when and how gantry tilt should be corrected.

My current understanding is that gantry tilt is being corrected now in dcm2niix like this: The tilt angle is extracted from the Gantry/Detector Tilt header (0018,1120). If this value is non-zero, and if the value of Derivation Description (0008,2111) is not "MEDCOM_RESAMPLED", the corrected image is formed by interpolating pixel values in an enlarged image, padding as necessary. The interpretation of the tilt angle varies by manufacturer: for some manufactures, the sign of the angle is reversed. Is that more or less correct?

In C.8.16.3.3 from the DICOM documentation, the description of Gantry/Detector Tilt says "Nominal angle of tilt in degrees of the scanning gantry. Not intended for mathematical computations." Further, ad hoc corrections for various manufacturers and other non-standard values found in headers make the current approach problematic. These should be acceptable only if there were no other way to perform the correction.

An alternative does exist, one which complies (as far as I understand) with the plain reading of the standard, and without the need for special handling of manufacturers and interpretation of other tags. Specifically, the absolute tilt angle can be computed using the angle between a vector normal to the two image orientation vectors and the vector difference between the image positions of two slices. Whether this value should be positive or negative can be computed (I think) using the same difference vector and the row or column orientation vectors. (In theory the tilt need not be aligned with the either of these, but I only propose for now to compute the tilt angle differently and leave the interpolation and padding as it currently exists.)

While not exactly the same use case as dcm2niix, this seems to be the approach taken by 3DSlicer. See for example the discussion here: https://discourse.slicer.org/t/how-to-adjust-the-distorted-image/817, especially one comment by Andras Lasso in October 2017: "The Gantry Tilt angle DICOM tag is for display purposes only. Accurate slice positions are defined by the slice position and axis directions DICOM tags. You may be able to compute an approximate transformation matrix from the Gantry Tilt tag, but it won’t be completely accurate (angle is rounded to one decimal, so you may get up to a few pixel error) and you may need some experimentation to get the rotation axis and direction right."

I suggest this in part because I have a CT where the tilt has already been corrected, the Gantry/Detector Tilt is populated but the Derivation Description is empty. dcm2niix therefore introduces shear in the "Tilt" image. Of course the other NIfTI is also available and correct, but choosing it requires manual inspection of the images.

If this seems like something worth pursuing, I can attempt at least a draft implementation. I do not have access to a large variety of examples on which to test, perhaps 6 or 7, several of which are probably from a single scanner.

Without this, I will at least try to detect cases in which the tilt correction is inappropriate and choose the appropriate NIfTI as produced by dcm2niix automatically. This I can do outside of dcm2niix, so a change to dcm2niix is not absolutely critical, but it does seem more widely beneficial than just me.

Thoughts?

@neurolabusc

This comment has been minimized.

Copy link
Collaborator

commented Dec 18, 2018

Feel free to develop your solution and submit a pull request. The current implementation works with all the samples I have. You seem to have run across a situation where some tool adjusted for the gantry tilt but left the gantry tilt tag unchanged in the DICOM header. I try not to second guess what is in the header - dcm2niix does not try to solve the Quis custodiet ipsos custodes? problem. It seems like another approach to your situation is to not apply the gantry tilt correction when the Derivation Description is empty. However, I am not sure how often this is populated. In any case, I would be happy to consider any alternative, but have a limited test set (locally we use CT systems that no longer apply gantry tilt).

@alund

This comment has been minimized.

Copy link
Contributor Author

commented Dec 18, 2018

Okay, I will see if I can make progress. I have just written a Python implementation of my proposed algorithm and, for those samples I have with a gantry tilt that actually requires correction, the computed value matches the Gantry/Detector Tilt value apart from sign. Since the value in the header sometimes needs to be reversed anyway, this may be correct or perhaps I need to reverse mine. I hope to start working on the dcm2niix implementation tomorrow, but it may take some time to figure out the correct way to do it there where I am less familiar with the available tools.

With regard to second guessing the header, I think this simply makes fewer assumptions. If the positions and orientations are not correct, the DICOM header has a serious problem.

neurolabusc added a commit that referenced this issue Dec 20, 2018
@neurolabusc

This comment has been minimized.

Copy link
Collaborator

commented Dec 20, 2018

@alund and @muschellij2 can you look at the latest commit in the development branch? It will report v1.0.20181220. The new method to determine gantry tilt is to read the row and column direction from Image Orientation (0020,0037): which I refer to as the read and phase direction, though note these are CT scans. I next read the Image Position Patient (0020,0032) for the first and last image in the volume. Subtracting these two gives the true slice direction. If the volume has no gantry tilt, the slice direction should be identical to the cross product of the read and phase directions. The angle between these two vectors is the gantry tilt. The cross product of these two vectors reveals whether the gantry tilt is positive or negative. The developmental release includes the compiler definition newTilt - this is defined by default but if undefined the old behavior is restored. I have tested this on samples from GE, Philips, Siemens and Toshiba. It is a very limited set, but the new code seems robust.

The Matlab code for my method is

Read =[1 0 0];
Phase =[0 0.939693 -0.34202];
CrossReadPhase = cross(Read,Phase)
Slice =[0 0 147.657];
Sign = cross(Slice,CrossReadPhase)
ThetaInDegrees = atan2d(norm(cross(Slice,CrossReadPhase)),dot(Slice,CrossReadPhase))

Benefits versus the old code:

  • The gantry tilt tag (0018,1120) is optional (type 3), and therefore a legal DICOM image might not include it.
  • @alund reports some tools correct gantry tilt in DICOM images but do not update the tag (e.g. the tag might refer to acquired gantry tilt, but not necessarily the stored gantry tilt).
  • Some have argued that 0018,1120 is often stored with low precision (though image position and orientation are often as well - if you want precision you should not be using DICOM).
  • It seems that different vendors specify the sign of 0018,1120 opposite of each other. So code needs a kludge to deal with this. On the other hand, the spatial mapping of the new approach is unambiguous.
@alund

This comment has been minimized.

Copy link
Contributor Author

commented Dec 21, 2018

I will take a look. I have also been working on this and was almost ready to submit a pull request.

@alund

This comment has been minimized.

Copy link
Contributor Author

commented Dec 21, 2018

A couple of comments. First, the tilt-corrected version looks correct for all of the examples I have, and it does not attempt to correct the tilt on the one image where the header says a tilt exists but the tilt was already corrected. In fact, I have a second image where I believe the tilt was corrected, but only imperfectly; there was a 0.02 degree residual. For this image, your code does not attempt to correct the tilt, which is acceptable. (This same image also has 0.01 degree tilt along the row orientation, which I suspect was an artifact of the original correction.)

Second, one of the things I noticed while working on my implementation of this is that the uncorrected output was affected by the tilt. Specifically, the slice spacing is being adjusted by the tilt angle, so that the spacing gets reduced even for the uncorrected output. I found this behavior surprising, and fixed it in my implementation. My reasoning was that either the tilt angle should be used (for the corrected version) or not at all (for the uncorrected).

Third, there may be a bug on line 3062:

    if ((nConvert > 1) && (dcmList[indx0].modality == kMODALITY_CT) || (dcmList[indx0].isXRay) || (dcmList[indx0].gantryTilt > 0.0)) {

I suspect you want additional parentheses around the three conditions joined by ||.

Fourth, and unrelated to these changes, please double check the logic on line 1048 (the second line of this excerpt):

    if ((h->dim[2] > 0) && (h->dim[1] > 0)) {
		if  (h->dim[2] == h->dim[2]) //phase encoding does not matter
			reconMatrixPE = h->dim[2];
		else if (d.phaseEncodingRC =='R')
			reconMatrixPE = h->dim[2];
		else if (d.phaseEncodingRC =='C')
			reconMatrixPE = h->dim[1];
    }

Specifically, comparing h->dim[2] to itself will always succeed unless it is NAN, in which case an explicit test for NAN would be preferred.

neurolabusc added a commit that referenced this issue Dec 21, 2018
alund added a commit to alund/dcm2niix that referenced this issue Dec 21, 2018
@alund

This comment has been minimized.

Copy link
Contributor Author

commented Dec 21, 2018

If you are interested in looking at my implementation, you can find it here:
https://github.com/alund/dcm2niix/tree/ComputeTilt

I find my computation of the tilt angle easier to understand, but it may be a matter of taste. As I mentioned above, I also changed the interslice distance calculation, which now takes a parameter to determine whether the tilt should be corrected or not.

I ignored dcmList[indx0].isResampled and used only the computed tilt angle to determine whether there is a tilt to be corrected.

Separate from the tilt angle, I also made some adjustments to the slice spacing consistency check: rather than comparing all spacings to the first one, I compare the min and max. And rather than using the first spacing when the consistency check passes, I switched to the mean spacing. For images where limited floating point precision in the distance calculations generates slight variations in interslice distance, this should better limit the impact on the overall dimensions of the image.

These changes were, of course, based of a prior revision so that a merge would now be more difficult. Please feel free to take whatever part of these changes you like.

@neurolabusc

This comment has been minimized.

Copy link
Collaborator

commented Dec 21, 2018

Sorry for any redundancy in our work. I have implemented all of your suggestions except for your second comment.

As you note, intersliceDistance() accounts adjusts slice distance. The image below describes the situation. The input slices are shown in black, and the Euclidean distance between two successive slices is shown by the blue line. The function intersliceDistance computes the distance between planes, as shown by the red line. The tilt correction accurately pads the missing voxels from the rectangular volume, reconstructing the dotted red volume. For this corrected volume, it is obvious that the red line is the appropriate measure: consider measuring the distance between the purple and green voxels.

Your question regards the value stored in the uncorrected volume. I believe you are correct that one wants to store the scalar value of the blue line in pixdim(3). However, the crucial thing would be to make sure the sform accurately reports the shear matrix. If computed correctly, a tool like SPM should be able to reconstruct accurate space from the sheared volume. Have you tested whether this is the case? During my development, I noted that most tools could not handle shears (since most tools are optimized for MRI where these are not done). Therefore, I focused my efforts on creating the tilt-corrected images. But it is probably worth ensuring the images with tilt are correct.

dx

@alund

This comment has been minimized.

Copy link
Contributor Author

commented Dec 21, 2018

I agree that the best solution for the uncorrected image, at least in principal, is to make the sform correct. Then you would have an image with no resampling which still accurately represents the data. Without that, we end up discussing the right way to do something wrong. I picked one wrong way, but I understand the rationale for your choice as well.

neurolabusc added a commit that referenced this issue Dec 21, 2018
@neurolabusc

This comment has been minimized.

Copy link
Collaborator

commented Dec 21, 2018

The latest commit applies the shear to the sform of scans with gantry tilt where the image is not corrected. I now feel convinced that pixdim(3) should be the inter-plane distance: this preserves the voxel volume (and is indeed what is reported in the DICOM header for slice thickness). The image below shows the same image viewed with SPM, where one image saves the original rhomboidal image data but inserts the shear amount into the sform while the other pads the volume and saves a rectangular volume. In my experience, with the exception of SPM few neuroimaging tools cope with NIfTI images where a shear is present, so use with caution.

If this passes your tests, please close this issues. Thanks for your suggestions and improvements.

dx2

@alund

This comment has been minimized.

Copy link
Contributor Author

commented Dec 27, 2018

Sorry for the delay, been off for holidays.

My tests pass with this latest change. (They do not check the sform changes though.)

I think you made an inadvertent change to README.md, substituting git for pigz where pigz was correct. I am not sure if you want to change that back before I close this?

neurolabusc added a commit that referenced this issue Dec 28, 2018
neurolabusc added a commit that referenced this issue Dec 29, 2018
@neurolabusc

This comment has been minimized.

Copy link
Collaborator

commented Dec 29, 2018

@alund there was a flaw in my solution. The cross product of the image position (0020,0032) row and column of vectors from a single slice provide the red arrow in my schematic above, and comparing the image position of the first and last slice provide the blue arrow. However, the polarity for these two arrows is not well defined: the image instance numbers may begin at the top or bottom of the table (influencing the polarity of the blue arrow) and the order of the cross product vectors defines the polarity of the red arrow. My new solution ensures that the magnitude of each of these two arrows is ascending, not descending.

All my examples with gantry tilt use ascending slices (e.g. instance number positively correlated with Z-position). However, I do have a sample with descending slices without gantry tilt, and my prior solution reported a 180-degree gantry tilt. It would be great to have an example of descending slices and a tilt applied - I am not 100% sure I have the polarity of the gantry tilt calculated correctly for this case.

@neurolabusc

This comment has been minimized.

Copy link
Collaborator

commented Dec 29, 2018

I just tested my solution where I reversed the instance numbers of an tilted dataset (this fools dcm2niix sorting algorithm to assume the data was acquired in descending rather than ascending order) and the current solution seems robust.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
2 participants
You can’t perform that action at this time.