Lossless scaling of 16-bit images

[neurolabusc]

A recent discussion on neurostars got me thinking about how we convert DICOM data to NIfTI. Most raw DICOM data is saved as 16-bit integers in the range -1024..1024 (CT), 0..4095 (12-bit MRI, e.g. most Siemens pre D13) or 0..65535 (16-bit MRI). For both 12 and 16-bit MRI, the actual data range is often a small part of the possible range as FFT scaling is set conservatively to prevent clipping. Therefore, in practice raw DICOM rarely uses much of the available 16-bit range. When the image intensity range is less than 32767 we can losslessly scale the data by an integer multiple. While it is true that the true SNR/CNR is probably limited, this is free precision.

The reason this matters is that many imaging tools will save images using the same data type as the input, so normalization, gaussian smooths, etc are often limited by this input precision. For example, a HCP T1 sequence which is at 0.8mm isotropic may be normalized to the space of SPM's TPM.nii that has a resolution of 1.5mm with better SNR than the source image. Saving these images as 32-bit float would preserve precision, but impact disk space and file IO.

The latest commit of dcm2niix contains a new integer scaling feature (opts.isMaximize16BitRange) which can be enabled by the hidden argument -is y (e.g. dcm2niix -is y -f %p_%s ~/myDICOMs). Consider an image where the darkest voxel is mn and the brightest is mx. If mn >= 0, we compute scale = trunc(64000/mx) and we save as UINT16. If mn < 0, scale = trunc(32000/max(mx,abs(mn)) and we save as INT16. In both cases, as long as scale > 1, all raw data is multiplied by this factor and the scl_slope is divided by this factor. Since the image data and scale factor are changed reciprocally, the calibrated voxel value does not change. Since scale is an integer, this is lossless (to the precision of scl_slope which is a 32-bit float).

In addition, the user gets a message such as Maximizing 16-bit range: raw 0..1764. This output might be useful: a user that sees this regularly for a particular sequence may decide that their FFT scaling is too conservative allowing them to capture more SNR from the reconstructor.

I realize this may have only marginal benefits, but I can think of no downside to this approach. Does anyone have any qualms with this? I am happy for comments from anyone, but off the top of my head I think this might interest @mharms @chrisfilo @yarikoptic @xiangruili @oesteban @hjmjohnson @ericearl

[yarikoptic]

I do not see major flaws (yet), just a few questions

• what if mx and mn <= 0? ;-) (not sure if realistically possible but some evil could make it happen)
• why 64000 and not just 65535? (similarly for half value)

and comments

• probably shouldn't ever become a default behavior happen someone relies on reading/analyzing raw values (although probably shouldn't without consulting scl_ fields)

[neurolabusc]

@yarikoptic

• My current solution works fine with mx and mn less than zero. Whichever is furthest from zero sets the scaling factor. In theory, we could just compute the range of mx-mn and adjust both the scl_slope and scl_inter appropriately.
• We could use 65535 for UINT16 and 32767 for INT16. I think in theory it is possible for windowed sinc functions to generate ringing where the output values exceed the input. In general one would hope that any algorithm would clip this, and for UINT16 I think zeros are so common this is probably done. In any case, I chose 64000 and 32000 for just a tiny bit of padding, but it is probably hard to justify.
• I will be interested what others say. Unlike you, I am tempted to make this the default behavior: in most imaging modalities the intensity is relative and software that ignores the scl_slope is likely cause other problems. It just seems like free precision should be the default behavior.

[coalsont]

I second the sentiment that it should probably not be made a default behavior. If someone simply adds two scans together, expecting a sane result, and instead it saturates the datatype because the converter did something more than just convert things, then that is surprising. I don't think it would affect us (connectome-workbench uses float32 no matter what, not sure about FSL and freesurfer), but it doesn't seem friendly to change the default output in a nontrivial way between versions of the tool.

Sinc and spline interpolations do in fact generate ringing, including negatives near the edge of the skull, so some built in padding is a reasonable idea. In an information-theory way, if you clip the overshoots, negative or positive, you are losing precision (in that interpolating back to the original grid will match better if you don't clip them). However, there are also valid concerns that the ringing artifacts can cause problems.

I replied via a lengthy email, mostly positive, but I will copy a specific concern here:

I am apprehensive about it changing the output datatype depending on the data values in the file, since apparently other tools write their output with the input datatype, as I assume by default you use the same signedness as the dicom encoding (disclaimer: I know very little about dicom, and have barely used it). If your scanner usually produces a few negative values, just because of noise and calibration, and then it drifts positive over time, resulting in subjects starting to have some scans with no negatives, and your tools suddenly can't put negative values into some processing intermediate file that you don't QC (because the tools aren't smart enough to change the type when needed, or to use a type that is always safe), that would be a rather cryptic error to chase down, and worse, it could be very subtle, too.

[coalsont]

A detail I missed earlier: If you don't constrain your expansion factor to be a power of 2, then you will have floating point rounding error in the scl_slope field (remember, binary floating point can't even exactly represent 0.2). Thus, the interpreted values from such a converted file will no longer be exact integers (and the rounding error from reconstruction will vary). Furthermore, you should not attempt to use the full range with an output datatype of INT32 or UINT32, because people may be computing the effective value in float32 (since that is the type of scl_slope), which can't represent odd integers greater than 2^24.

A second argument against making this a default: tools that are already smart enough to make use of scl_slope when the output values are non-integers (which is a simple, sensible approach) could be hurt by this behavior, from a different direction. Consider:

You convert an integer-valued dicom into a nifti volume with a scl_slope of 0.25 (the raw values in the data section of the file are now 0, 4, 8, 12, etc). Someone uses a tool, which likes to write things as the input datatype, to multiply the image by 2. Their software looks at the output values (which it calculates in a higher precision, as it should) and sees that they are integers, so it tries to reuse the input header. However, it then sees that it would saturate the output datatype if it used the same scaling as the input, and is now in a previously rarely-used code path. It comes up with a new scaling factor from scratch (because this is general i/o code, so it shouldn't care what computation was done), and decides the output should be scaled by 0.35, because it does the obvious thing, spreading the range of output values as wide over the output integers as possible, assuming the common case would be writing continuous-valued data, not discrete. So, now, a value that was 1 in the dicom, got multiplied to an exact 2 inside the tool, is now rounded to 2.1. A value that was 2 in the dicom, and got multiplied to an exact 4, is now rounded to 3.85 (since that is closer than 4.2).

Multiplying by 2 may be somewhat contrived, but adding two volumes together would trip the same circumstances (with the complication of 2 input headers, which may have different scl_slope).

I can also imagine that some tools might check whether scl_slope and scl_inter are anything other than identity scaling, and if so, switch their output to floating-point (and maybe even take a different computation code path).

[chrisgorgo]

@matthew-brett and @effigies might have some input

[xiangruili]

I don't see any downside, except someone makes mistake, like saturating the default datatype @coalsont mentioned, or ignoring the scale parameters. So it is indeed free precision.

If I remember correctly, any FSL image operation will result in float32 NIfTI, probably for simplicity of implementation. When a tool modifies image, I think it is the tool's duty to decide if it is necessary to switch to float32, or to apply this trick to keep int16/uint16 to reduce file size by a factor or 2.

[neurolabusc]

@coalsont thanks for your comment.

• I take your concern of converting a INT16 with positive range to a UINT16. In theory this could cause clipping of values below zero when functions like windowed sinc are applied. I did think about this in my code, and if you change the line bool isConvertToUint16 = true; to read bool isConvertToUint16 = false; than INT16 always remain INT16 (while an input UINT16 remains UINT16). I want to say no one mentioned any issues when Siemens output moved from INT16 to UINT16 with many magnitude images having lots of near-zero values that necessarily must elicit such behavior. Therefore, I assume this concern is largely academic and all tools handle clipping (presumable by zeroing negative values). Regardless, I do think you this is a valid point and in future commits I will set the default to isConvertToUint16 = false.
• Regarding your concern of scales other than power of 2, I did note that this is lossless to the precision of scl_slope - which is a 32-bit float in NIfTI-1 (7-8 decimal places precision). However, it is worth noting that this is far beyond the SNR of CT or MRI, and a veritable rounding error relative to the extra precision offered by providing an integer slope scaling.
• It is also worth pointing out that while the DICOM Rescale Slope tag (0028,1053) is often 1.0 (which can be handled losslessly as a FLOAT32), it is reported as a decimal string and in many instances has very low precision. While this does not matter for most MRI scans (where intensity scaling is relative), one should not have any illusions regarding the precision of this value for calibrated modalities.
• I do not think the concern of adding two images together exceeding the image range should be a legitimate concern - this reflects an error in the software that is doing the computation: if it fails adding two INT16 images that had the integer scaling, it would also fail in adding two Siemens 16-bit images where the FFT scale was set optimally. Further, tools that have such issues with scl_slope will probably have really issues with Philips data where the user may want the DV or FP. These faulty tools should be identified and updated. I think this is why tools like fslmaths default to saving any computation as FLOAT32. The tools that tend to preserve the input format are generally those that do not intend to increase or reduce the gain of the signal.
• Your concern with regards to UINT32 and INT32 is well taken. In practice, I have never seen this datatype in DICOM (high SNR data is saved as FLOAT32), and my integer scaling is specifically for the common case of INT16 and UINT16 data.
• As an example, consider an integer image with a raw value of 0..3200 and a scl_slope of 1.0 and scl_inter of 0.0, with integer scaling (and isConvertToUint16 = false) this will be saved as 0..32000 with a scl_slope of 0.1. The calibrated value of both is 0..3200, however any tool that smooths/interpolates/modifies this data and saves in the original format will be required to save in discrete bins of 1.0 (e.g. either 301 or 302 with no value in between) for the former, while in the scaled case the data can be saved in bins of 0.1 (e.g. 301.0, 301.1, 301.2...302).
• I agree that it is possible that tools could decide that images without a scl_slope of 1.0 should be saved as float, and this might create larger files and slow IO. However, I have tested the major popular tools (e.g. FSL's FLIRT, SPM's normalization, etc) and they those that have the default behavior of saving images in their input format do so regardless of integer scaling. In brief, while it is possible to think of edge cases where integer scaling could create larger files and cause slower IO, in the typical cases it offers free precision.
• I really do not have an agenda here, and appreciate the frank discussion. If the majority of people think the default behavior should be to retain the raw data range, I am happy to oblige.

[neurolabusc]

Just a clarification to @xiangruili's comment, FSL tools that can combine data like fslmaths do convert data to FLOAT32, but many (flirt, mcflirt, etc) use the input datatype. Likewise, this is the default for all SPM operations.

[neurolabusc]

Since I am eager to generate a new stable release (very overdue), and so far there is not a clear consensus, here is my idea for the imminent stable release:

• Integer scaling is turned off by default, and is switched on with the -l y argument
• INT16 input uses scale = trunc(32000/max(mx,abs(mn)) and is always saved as INT16
• UINT16 input uses scale = trunc(64000/mx) and is always saved as UINT16
• Over the next few months people can weigh on on whether this should be the default behavior. My sense is most users just run the software out of the box, so I like the idea of providing as much free precision as possible. However, I realize I might be in a minority.
• Integer scaling can be stored as a user preference with the -g y or -g o arguments (generate defaults file). This allows a user to make this the default behavior. This is shown below: initially, the default is to not use integer scaling, but the user can modify this behavior by generating a default file. The user can always override the default behavior by explicitly including either -l y or -l n. This default will remain until the user generates new defaults:

$ ./dcm2niix
Chris Rorden's dcm2niiX version v1.0.20180614 GCC6.1.0 (64-bit MacOS)
...
  -l : losslessly scale 16-bit integers to use dynamic range (y/n, default n)
...
$ ./dcm2niix -l y -g y
Chris Rorden's dcm2niiX version v1.0.20180614 GCC6.1.0 (64-bit MacOS)
Saving defaults file /Users/rorden/.dcm2nii.ini
Example output filename: '/myFolder_MPRAGE_19770703150928_1.nii'
$ ./dcm2niix
Chris Rorden's dcm2niiX version v1.0.20180614 GCC6.1.0 (64-bit MacOS)
...
  -l : losslessly scale 16-bit integers to use dynamic range (y/n, default y)
...

[captainnova]

Hi Chris, I have no objection, but will probably stick with the current behavior, at least for diffusion MRI from Siemens. There we often see clipping at 4095 (the 12 bit max) for CSF in b=0 volumes; with CSF being ~2x as bright as tissue and us caring more about the precision of tissue at b=2000, we do hit that ceiling. But we don't want to clip too much, and I like keeping an eye on how many voxels are at the magic 4095 I know to watch for. Apparently there's a way to get Siemens to store DICOM with different scaling for each volume, but it was a bit too fiddly for ADNI. Rob

…

On Thu, Jun 14, 2018, 1:31 PM Chris Rorden ***@***.***> wrote: Since I am eager to generate a new stable release (very overdue), and so far there is not a clear consensus, here is my idea for the imminent stable release: - Integer scaling is turned off by default, and is switched on with the -l y argument - INT16 input uses scale = trunc(32000/max(mx,abs(mn)) and is always saved as INT16 - UINT16 input uses scale = trunc(64000/mx) and is always saved as UINT16 - Over the next few months people can weigh on on whether this should be the default behavior. My sense is most users just run the software out of the box, so I like the idea of providing as much free precision as possible. However, I realize I might be in a minority. - Integer scaling can be stored as a user preference with the -g y or -g o arguments (generate defaults file). This allows a user to make this the default behavior. This is shown below: initially, the default is to not use integer scaling, but the user can modify this behavior by generating a default file. The user can always override the default behavior by explicitly including either -l y or -l n. This default will remain until the user generates new defaults: $ ./dcm2niix Chris Rorden's dcm2niiX version v1.0.20180614 GCC6.1.0 (64-bit MacOS) ... -l : losslessly scale 16-bit integers to use dynamic range (y/n, default n) ... $ ./dcm2niix -l y -g y Chris Rorden's dcm2niiX version v1.0.20180614 GCC6.1.0 (64-bit MacOS) Saving defaults file /Users/rorden/.dcm2nii.ini Example output filename: '/myFolder_MPRAGE_19770703150928_1.nii' $ ./dcm2niix Chris Rorden's dcm2niiX version v1.0.20180614 GCC6.1.0 (64-bit MacOS) ... -l : losslessly scale 16-bit integers to use dynamic range (y/n, default y) ... — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#198 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/APHk6Q_HYmzp8eusTux2G-R4Rk_iFSoVks5t8qvkgaJpZM4Um6oB> .

[coalsont]

That plan looks safe to me.

I had no idea that the command line tool would read a preferences file by default, changing how it interpreted any not-provided options, this makes it more difficult to write scripts with reproducible output across even different user accounts on the same machine. I have commented on #152 to ask for a preferences mode that can be friendly to both scripting and user preferences.

I also have a nitpick for the phrasing/framing used: you aren't providing more actual precision (it is a converter, after all, not reacquisition). The way I think of what this is doing, is trying to trick suboptimally-written tools into writing their outputs with more precision than they otherwise would have, by pretending that the data has more precision than it actually does. This can absolutely be a positive result for people using such tools, as long as it doesn't come as a surprise (investigators really don't like getting different results from rerunning the same scripts on the same data, regardless of the shortcomings of the tools used).

[xiangruili]

@coalsont That is exactly right. This has no way to increase resolution, but to help other tools to potentially give better resolution. That is why I think it should be other tool's duty. But without apparent downside, it is nice to take care of things for others :)
I am giving this option for dicm2nii.m, default without scaling.

[neurolabusc]

@coalsont thanks for your comments on issue 152. I have added your suggestion to allow scripts to ignore user defaults without destroying them (-g i). I agree that using the full dynamic range could change results by reducing rounding errors. Therefore, in theory rerunning a pipeline on an archival dataset might give a slightly different result. However, this is true whenever anyone improves any component of your pipeline, and this seems to be a tiny but unqualified improvement. I will certainly keep the old option available for those who want to replicate work, but as I have noted my personal preference is to make improvements the default behavior.

Also, it is worth pointing out that dcm2niix only saves three values in its default file, and I think each is justified:

• isGZ : save .nii (0) or .nii.gz (1) mirrors FSL's ~/.fslconf/fsl.sh FSLOUTPUTTYPE and allows users of SPM (which does not by default handle gz images) to ensure regular compatibility
• isMaximize16BitRange : save raw data (0) or lossless scaling (1): the topic of this issue.
• filename : default is "%f_%p_%t_%s" - I would suspect most scripts would explicitly set their file name preferences. The reason I made this an option is that for many centers storing the time in the filename is a key attribute to associate the image with a participant, but some users told me that they worried that the study time might be identifiable information. Therefore, for these users they can change this default.

[neurolabusc]

@captainnova - understood. Though do note that 4095 will be transformed predictably, so you could either have your QA scripts count either the occurrences of 28865 or the number of values 4095/scl_slope (since all the Siemens raw DTI data I have seen use a raw scale slope of 1).

[neurolabusc]

To help people like Rob out, if integer scaling is used this will be reported in the NIfTI description text field using the " isN" format. For example, if an integer scaling of 14 was applied, this field might look like "Time=181433.000 is14".

[neurolabusc]

I am closing this issue. This is an optional feature of the current release.