Resampling in wrong direction?

[ashgillman]

Hi,

I have used the nm_mrac utility to convert a u-map acquisition to interfile, and then ran STIR's calculate_attenuation_coefficients with the input sinogram as the template. I attach here a screenshot of projection views of the sinogram and of the ACFs.

As you can see, it appears the y-axis has been inverted (it is a little hard to tell w.r.t. x-axis also). Is this likely to be an issue with the resampling in nm_mrac do you think?

Cheers, Ash

[ashgillman]

Here is the same issue with a phantom acquisition (3 point sources and a water bottle). Sorry, the ACF and emission sinogram have swapped sides here.

It is harder to see, but the cursor is on the same point in both images, and erroneously shows one of the points to lie within the bottle, when it was above. I would have to confirm, but I should be able to share this phantom data if you'd like to compare the acquisition with that you acquire.

[ashgillman]

And a subset of the potentially relevant DICOM tags from the phantom MRAC:

(0002,0013) SH [MR_VE11P]                                         # 8,1 Implementation Version Name
(0008,0008) CS [DERIVED\PRIMARY\M\NORM\IN_PHASE\DIS3D\DIS2D\MRPET_UMAP3D]         # 56,2-n Image Type
(0008,0060) CS [MR]                                               # 2,1 Modality
(0008,0070) LO [SIEMENS ]                                         # 8,1 Manufacturer
(0008,1090) LO [Biograph_mMR]                                     # 12,1 Manufacturer's Model Name
(0018,0015) CS [BRAIN ]                                           # 6,1 Body Part Examined
(0018,0020) CS [GR]                                               # 2,1-n Scanning Sequence
(0018,0021) CS [SP]                                               # 2,1-n Sequence Variant
(0018,0022) CS [PFP ]                                             # 4,1-n Scan Options
(0018,0023) CS [3D]                                               # 2,1 MR Acquisition Type
(0018,0024) SH [*fl3d2]                                           # 6,1 Sequence Name
(0018,0025) CS [N ]                                               # 2,1 Angio Flag
(0018,0050) DS [2.0199999809265 ]                                 # 16,1 Slice Thickness
(0018,0080) DS [4.14]                                             # 4,1 Repetition Time
(0018,0081) DS [2.51]                                             # 4,1 Echo Time
(0018,0083) DS [1 ]                                               # 2,1 Number of Averages
(0018,0084) DS [123.260646]                                       # 10,1 Imaging Frequency
(0018,0085) SH [1H]                                               # 2,1 Imaged Nucleus
(0018,0086) IS [2 ]                                               # 2,1-n Echo Number(s)
(0018,0087) DS [3 ]                                               # 2,1 Magnetic Field Strength
(0018,0089) IS [134 ]                                             # 4,1 Number of Phase Encoding Steps
(0018,0091) IS [2 ]                                               # 2,1 Echo Train Length
(0018,0093) DS [75]                                               # 2,1 Percent Sampling
(0018,0094) DS [53.125]                                           # 6,1 Percent Phase Field of View
(0018,0095) DS [1300]                                             # 4,1 Pixel Bandwidth
(0018,1000) LO [51040 ]                                           # 6,1 Device Serial Number
(0018,1020) LO [syngo MR E11]                                     # 12,1-n Software Version(s)
(0018,1030) LO [Head_MRAC_Brain_HiRes ]                           # 22,1 Protocol Name
(0018,1251) SH [Body]                                             # 4,1 Transmit Coil Name
(0018,1314) DS [10]                                               # 2,1 Flip Angle
(0018,1315) CS [N ]                                               # 2,1 Variable Flip Angle Flag
(0018,1316) DS [0.07342813610811]                                 # 16,1 SAR
(0018,1318) DS [0 ]                                               # 2,1 dB/dt
(0018,5100) CS [HFS ]                                             # 4,1 Patient Position
(0019,0010) LO [SIEMENS MR HEADER ]                               # 18,1 Private Creator
(0019,1008) CS [IMAGE NUM 4 ]                                     # 12,1 CSA Image Header Type
(0019,1009) LO [1.0 ]                                             # 4,1 CSA Image Header Version ??
(0019,100b) DS [38827.5 ]                                         # 8,1 SliceMeasurementDuration
(0019,100f) SH [Fast* ]                                           # 6,1 GradientMode
(0019,1011) SH [No]                                               # 2,1 FlowCompensation
(0019,1012) SL 0\0\-1155                                          # 12,3 TablePositionOrigin
(0019,1013) SL 0\0\-1156                                          # 12,3 ImaAbsTablePosition
(0019,1014) IS [0\0\-1]                                           # 6,3 ImaRelTablePosition
(0019,1015) FD -250.003\-132.814\131.82                           # 24,3 SlicePosition_PCS
(0019,1017) DS [0.5 ]                                             # 4,1 SliceResolution
(0019,1018) IS [1000]                                             # 4,1 RealDwellTime
(0020,0032) DS [-250.0032\-132.8142\130.81996052645 ]             # 36,3 Image Position (Patient)
(0020,0037) DS [1\0\0\0\0\-1]                                     # 12,6 Image Orientation (Patient)
(0020,1040) LO (no value)                                         # 0,1 Position Reference Indicator
(0020,1041) DS [-132.8142 ]                                       # 10,1 Slice Location
(0020,4000) LT [uMap]                                             # 4,1 Image Comments
(0028,0002) US 1                                                  # 2,1 Samples per Pixel
(0028,0004) CS [MONOCHROME2 ]                                     # 12,1 Photometric Interpretation
(0028,0010) US 126                                                # 2,1 Rows
(0028,0011) US 240                                                # 2,1 Columns
(0028,0030) DS [2.08626\2.08626 ]                                 # 16,2 Pixel Spacing
(0028,0100) US 16                                                 # 2,1 Bits Allocated
(0028,0101) US 12                                                 # 2,1 Bits Stored
(0028,0102) US 11                                                 # 2,1 High Bit
(0028,0103) US 0                                                  # 2,1 Pixel Representation
(0028,0106) US 0                                                  # 2,1 Smallest Image Pixel Value
(0028,0107) US 0                                                  # 2,1 Largest Image Pixel Value
(0028,1050) DS [900 ]                                             # 4,1-n Window Center
(0028,1051) DS [900 ]                                             # 4,1-n Window Width
(0028,1055) LO [Algo1 ]                                           # 6,1-n Window Center & Width Explanation
(0040,0280) ST (no value)                                         # 0,1 Comments on the Performed Procedure Step
(0051,0010) LO [SIEMENS MR HEADER ]                               # 18,1 Private Creator
(0051,100a) LO (SH) [TA 38.82]                                    # 8,1 TimeOfAcquisition
(0051,100c) LO (SH) [FoV 500*262 ]                                # 12,1 Field of View
(0051,100d) SH [SP A132.8 ]                                       # 10,1 Slice Position
(0051,100e) LO (SH) [Cor ]                                        # 4,1 Slice Orientation
(0051,100f) LO [HEA;HEP;NEP ]                                     # 12,1 Coil String
(0051,1012) SH [TP F1 ]                                           # 6,1 ?Table Position?
(0051,1013) SH [+LPH]                                             # 4,1 Positive PCS Directions
(0051,1016) LO [M/NORM/IN_PHASE/DIS3D ]                           # 22,1 ?Image Type?
(0051,1017) SH [SL 2.02 ]                                         # 8,1 ?Slice Thickness?
(0051,1019) LO [A1/PFP]                                           # 6,1 ?Scan options?

[KrisThielemans]

Before going any further with testing orientation, can you use the STIR ITK_orientation branch? It'll only work for HFS (if it does work at all...).

[ashgillman]

Just checking the potentially obvious: I should specify the output extension of nm_mrac to be .nii.gz?

[KrisThielemans]

it shouldn't really matter I believe. as long as you use a format that supports image orientation (like NRRD, MetaIO, Nifti), then nm_mrac will use ITK to write, and STIR will use ITK to read, and everything should be fine...

it'd be fun (?) to compare the images that you get when using STIR to directly read the MRAC-dicom, the Nifti/MetaIO file output by nm_mrac and the Interfile output by nm_mrac. STIR should read these the same way. (ok, we didn't put in the offsets yet in the Interfile output of nm_mrac as we're not 100% sure what to do with those yet). Question is of course how you do that comparison. Recommended would be to read them via the Python/MATLAB interface to STIR and display there. STIR's manip_image is a terrible way to display. Or you could ask STIR to convert back to something else (e.g. via stir_math) and then use another display (but then you're never sure what conventions that display program uses).

[ashgillman]

Thanks, in the above examples I used Interfile as the output, which would bypass that reader.

The latter two should be identical within STIR, but the version read directly from STIR would have a different (the original) sampling. If you convert all three the same way (stir_math) it shouldn't matter the way the program displays if all you care is they are the same? I often use stir_math to convert to NIfTI then view (somewhat arbitrarily) with mrtrix's mrview.

[ashgillman]

Hmm, I have reached another issue:

I can't seem to run calculate_attenuation_coefficients on a .nii u-map derived from nm_mrac. (Neither master nor ITK_orientation branches).

Forward projector used:           
forward projector using matrix parameters :=                                                                        
matrix type := Ray Tracing        
ray tracing matrix parameters :=       
disable caching := 0              
store only basic bins in cache := 1                                                                                   
restrict to cylindrical fov := 1  
number of rays in tangential direction to trace for each bin := 1
use actual detector boundaries := 0
do symmetry 90degrees min phi := 1                                                                                       
do symmetry 180degrees min phi := 1
do symmetry swap segment := 1          
do symmetry swap s := 1           
do symmetry shift z := 1                                                                                   
end ray tracing matrix parameters :=

end forward projector using matrix parameters :=

WARNING:                           
WARNING: BinNormalisationFromAttenuationImage:
        attenuation image data are supposed to be in units cm^-1
        Reference: water has mu .096 cm^-1                                                                                                  Max in attenuation image: 0.244165

ERROR: DataSymmetriesForBins_PET_CartesianGrid: the shift in the z-direction of the origin (which is 98.084) should be a multiple of
 the plane separation (2.03125)              

terminate called after throwing an instance of 'std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >'    Aborted (core dumped)     
Return code: 134      

This wasn't the case when using Interfile output directly from nm_mrac2mu.

The image info is:

gil2a4 at 11:18:00 on daftpunk-rbh in /media/gil2a4/My_Passport/PISA_raw

♪ cat petmr-rd-tools/PISA_410141_01_HEAD_MRAC_BRAIN_HIRES_IN_UMAP_0010.hv
!INTERFILE:=
%comment:=created with nm_mrac2mu for mMR data
!originating system:=2008

!GENERAL DATA:=
!name of data file:=PISA_410141_01_HEAD_MRAC_BRAIN_HIRES_IN_UMAP_0010.raw
!GENERAL IMAGE DATA:=
!type of data := PET

%study date (yyyy:mm:dd):=2017:11:08
%study time (hh:mm:ss GMT+00:00):=15:15:12
imagedata byte order:=LITTLEENDIAN
%patient orientation:=HFS
!PET data type:=image
number format:=float
!number of bytes per pixel:=4
number of dimensions:=3
matrix axis label[1]:=x
matrix axis label[2]:=y
matrix axis label[3]:=z
matrix size[1]:=440
matrix size[2]:=328
matrix size[3]:=108
scaling factor (mm/pixel) [1]:=2.08626008
scaling factor (mm/pixel) [2]:=2.08626008
scaling factor (mm/pixel) [3]:=2.03125
start horizontal bed position (mm):=-10
end horizontal bed position (mm):=-10
start vertical bed position (mm):=0.0

!IMAGE DATA DESCRIPTION:=
!total number of data sets:=1
number of time frames:=1
!image duration (sec)[1]:=0
!image relative start time (sec)[1]:=0

%SUPPLEMENTARY ATTRIBUTES:=
quantification units:=1/cm
slice orientation:=Transverse
%image zoom:=1
%x-offset (mm):=0.0
%y-offset (mm):=0.0
%image slope:=1
%image intercept:=0.0
maximum pixel count:=0.244164839
minimum pixel count:=0
!END OF INTERFILE :=

gil2a4 at 11:21:02 on daftpunk-rbh in /media/gil2a4/My_Passport/PISA_raw

♪ mrinfo petmr-rd-tools/PISA_410141_01_HEAD_MRAC_BRAIN_HIRES_IN_UMAP_0010.nii
************************************************
Image:               "petmr-rd-tools/PISA_410141_01_HEAD_MRAC_BRAIN_HIRES_IN_UMAP_0010.nii"
************************************************
  Dimensions:        440 x 328 x 108
  Voxel size:        2.08626 x 2.08626 x 2.03125
  Data strides:      [ 1 2 -3 ]
  Format:            NIfTI-1.1
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1          -0           0      -457.2
                                0           1          -0      -340.4
                                0           0           1      -119.3

I am not exactly sure where the 98.084 number is coming from, or why that offset is necessary.

[KrisThielemans]

This happens because we currently don't write the offset in the Interfile header as I said. Arguably we should write the x,y offsets as I'm 99% confident that we'd handle those ok (but I suspect they're going to be mostly zero). However, the z offset is current problematic as we ignore bed position (i.e. z offset) for projection data. Current STIR assumes that the latter is always zero. So, when using the nii, we keep the original DICOM offset, and therefore once you start forward projecting it, STIR will think it's been shifted.

Summary: when using the Interfile header will (only) work when the MRAC was centred w.r.t. the scanned bed position. When using the nii/mhd/whatever output of nm_mrac it will always fail.

Clearly, this is problematic.The only proper way to fix this is to incorporate bed positions in the projection data, and the projectors. This has been done by others, but was never contributed. It'd probably be 1 day's work for me, but I don't have a spare day...

One potential other thing to do is to check in nm_mrac if the image is centred, or to resample it such that it is centred. That'd require knowing what the centre is though (which would depend on what bed position you want to use this for)

[KrisThielemans]

coming back to the previous point about comparing the interfile output, the "nii" and DICOM. The offsets of the first will be "wrong". You're right that the "nii" and DICOM will have different sampling, but otherwise they should be fine. If your display programme handles that, they should display the same. The Interfile-one should be displayed shifted.

To handle the resampling in STIR, you could use zoom_image (which will then rescale image values though as it uses count-preserving resampling, which isn't appropriate for mu-maps. hence the zoom_att_image.sh work-around).

All a bit of a mess...