rordenlab / dcm2niix

dcm2nii DICOM to NIfTI converter: compiled versions available from NITRC
https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage
Other
872 stars 229 forks source link

Philips scaling with when RW units present not including scale slope ? #493

Closed drmclem closed 3 years ago

drmclem commented 3 years ago

I have had a query from one of our users - It looks as if there is a problem with the scaling logic when the real world tags are present for the rescale slope and rescale intercept. For the "precise" conversion (FP in Philips terminology) the scale slope needs to also be taken into consideration but I think even though displayed in the logs only the RS and RI are being used leading to a conversion to the displayed value (matching the scanner console display) rather than the FP value.

Using RWVSlope:RWVIntercept = 0.596581:0 Philips Scaling Values RS:RI:SS = 0.596581:0:1.6765 (see PMC3998685) Convert 20 DICOM as DICOM\DICOM_pCASL_rest_20210309125250_201a (128x128x20x1)

To Reproduce Compare the nii output from data with and without the real world scaling tags but with identical conventional rescale tags - the scaling is different

Expected behavior For FP output the results should be the same.

Chris Rorden's dcm2niiX version v1.0.20190902 (JP2:OpenJPEG) (JP-LS:CharLS) MSC1900 (64-bit Windows) Found 4 DICOM file(s) Philips Scaling Values RS:RI:SS = 0.596581:0:1.6765 (see PMC3998685) Convert 1 DICOM as Z:\mjf_head\PCASL_ENHANCED_DICOM\DICOM\20210309125250_pCASL_rest_201_real (128x128x20x1)

Further info - It looks like for this system that the real-world tags are present when the export is in classic dicom but are not present when the export is enhanced (no idea why this decision was made) but the conventional rescale tags are present in both versions.

Matthew

neurolabusc commented 3 years ago

@drmclem I propose changing the Image Scaling section of the Philips conversion notes. The new text is now explicit that precedence is given to the real world scaling. Further, the alternative scaling factors are available in the BIDS sidecar. The variation described by the user appears to describe differences in image export between Classic and Enhanced DICOM data. As you are one of the experts in this topic, I would be grateful if you could proof the text below and suggest edits.

Image Scaling

dcm2niix losslessy copies the raw data from DICOM to NIfTI format. These values are typically stored as 16-bit integers in the range -32768..32767. Both the DICOM and NIfTI formats describe how scaling intercept and slope values can be used to convert these raw values into calibrated values. For example, with an intercept of 0 and slope of 0.01 the raw value of 50 would be converted to 0.5.

The NIfTI header provides the scl_slope and scl_inter fields so each voxel value in the dataset is scaled as:

I = scl_slope * R + scl_inter

where R is the voxel value stored and I is the "true" transformed voxel intensity.

For all manufacturers except Philips, the scl_slope and scl_inter correspond to the DICOM tags 0028,1053 and 0028,1052, respectively.

Philips has three possible intensity transforms for their DICOM images (world (W), display (D), precise (P)). All of these transforms might be provided in a single DICOM image, while the NIfTI standard only designates a single scl_slope and scl_inter for each image. dcm2niix will retain the raw value (R) and sets the NIfTI scl_inter and scl_slope values for the desired intensity transform. dcm2niix will use W if possible. If this is not possible it will use P unless the user specifies -p n to disable the precise calculation resulting in D.

The formulas are provided below. The DICOM tags are in brackets (e.g. (0040,9225)) and the BIDS tag is in double quotes (e.g. "PhilipsRWVSlope"). Since all the scaling values are stored in the BIDS sidecar, you can always use these to generate your preferred intensity transform.

Inputs:
 R = raw stored value of voxel in DICOM without scaling
 WS = RealWorldValue slope (0040,9225) "PhilipsRWVSlope"
 WI = RealWorldValue intercept (0040,9224) "PhilipsRWVIntercept"
 RS = rescale slope (0028,1053) "PhilipsRescaleSlope"
 RI = rescale intercept (0028,1052) "PhilipsRescaleIntercept"
 SS = scale slope (2005,100E) "PhilipsScaleSlope"
Outputs:
 W = real world value
 P = precise value
 D = displayed value
Formulas:
 W = R * WS + WI
 D = R * RS + RI
 P = D / (RS * SS)
drmclem commented 3 years ago

Hi Chris,

I have made some changes below which you are free to use or not use as you see fit - has a slightly different take I think but please feel free to adapt as necessary.

I think the yellow bits may need revisiting depending on how you feel about the earlier text.

Matthew

From: Chris Rorden @.> Sent: Thursday, March 11, 2021 11:55 AM To: rordenlab/dcm2niix @.> Cc: Clemence, Matthew @.>; Mention @.> Subject: Re: [rordenlab/dcm2niix] Philips scaling with when RW units present not including scale slope ? (#493)

Caution: This e-mail originated from outside of Philips, be careful for phishing.

@drmclemhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdrmclem&data=04%7C01%7C%7C515789d6dcf94de5906208d8e4849037%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510605290613238%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=QJIFpC4WTyJzk2w1zIbxy7ElMDW%2FiM0dOEnRiAOFsNA%3D&reserved=0 I propose changing the Image Scaling section of the Philips conversion noteshttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Frordenlab%2Fdcm2niix%2Ftree%2Fmaster%2FPhilips&data=04%7C01%7C%7C515789d6dcf94de5906208d8e4849037%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510605290623233%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=NKvVMonglfpJ1wlVXQQ%2FZ%2Bl%2FprY5CKQ9FGQ6H21MsDI%3D&reserved=0. The new text is now explicit that precedence is given to the real world scaling. Further, the alternative scaling factors are available in the BIDS sidecar. The variation described by the user appears to describe differences in image export between Classic and Enhanced DICOM data. As you are one of the experts in this topic, I would be grateful if you could proof the text below and suggest edits.

Image Scaling

How data is represented in DICOM for MR has several challenges and the technology and standard has evolved over the years to accommodate new uses. Unlike CT, where the signal is naturally displayed in Hounsfield units, MR has no natural signal units and the magnitude is influenced by the electronics and the software processing required to bring this to the final image. Secondly most of the original DICOM implementations used small bit number integers to store the underlying images for economy of storage. As a result it is necessary to apply scaling from the internal DICOM storage to a form suitable for radiographic display or quantitative measurement. There remain several challenges with this process, ensuring that the mapping to the integer values makes best use of the available bit depth for images with large dynamic range, or large changes between images, without clipping the data while also preserving the appearance of the noise field which is demanded by the needs of radiographic visual review.

At its simplest this requires a rescale and the DICOM standard tags 0028,1053https://eur01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fdicomlookup.com%2Flookup.asp%3Fsw%3DTnumber%26q%3D(0028%2C1053)&data=04%7C01%7C%7C515789d6dcf94de5906208d8e4849037%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510605290633228%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=KCbY5f9kpSQy67TaU%2BEfZD06R%2BPMn0ZfXiSdDoi4yXg%3D&reserved=0 and 0028,1052https://eur01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fdicomlookup.com%2Flookup.asp%3Fsw%3DTnumber%26q%3D(0028%2C1053)&data=04%7C01%7C%7C515789d6dcf94de5906208d8e4849037%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510605290643227%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=mRHsqOhrNjTOA0N8Y5JV1AH7yNm2V41VakG7FrfUS7I%3D&reserved=0. Whether these values are the same for all images, or image specific depends upon the implementation and potentially the location of these tags withing the DICOM tag structure.

In addition, the DICOM standard introduced the concept of "real world units" see http://dicom.nema.org/dicom/2013/output/chtml/part03/sect_A.46.html . This allows the storage of 1 or more mappings to allow selective viewing of the data mapped into different value ranges (which may also be non-linear mappings).

Clearly, scaling is a crucial aspect to be aware of for quantitative methods and which representation is used depends upon your needs.

For Philips, we think in terms of 3 different representations (using the terminology of the documentation available to Philips collaborators)

Stored value SV The Integer representation which is contained in the DICOM file tag PIXEL DATA

Displayed Value DV The value which is shown to the user when using scanner interface, ROIS, measurements etc.

Floating Point FP An internal value at a point earlier in the reconstruction chain before the conversion to DICOM/integer for image presentation.

Usage cases

In general SV should not be used for quantitative measurements as it is an integer format. In practice, if the Rescale values are the same for all images (the typical case, but not guaranteed) SV can be used to compare signal intensities between images from the same scan.

DV can be used for quantitative comparison of signal intensities between images in the same scan as long as the relevant rescale values are taken into account. These rescale values may come from the tags standard tags 0028,1053https://eur01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fdicomlookup.com%2Flookup.asp%3Fsw%3DTnumber%26q%3D(0028%2C1053)&data=04%7C01%7C%7C515789d6dcf94de5906208d8e4849037%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510605290633228%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=KCbY5f9kpSQy67TaU%2BEfZD06R%2BPMn0ZfXiSdDoi4yXg%3D&reserved=0 and 0028,1052https://eur01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fdicomlookup.com%2Flookup.asp%3Fsw%3DTnumber%26q%3D(0028%2C1053)&data=04%7C01%7C%7C515789d6dcf94de5906208d8e4849037%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510605290643227%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=mRHsqOhrNjTOA0N8Y5JV1AH7yNm2V41VakG7FrfUS7I%3D&reserved=0 or from a relevant RealWorld block if present.

If the DV is derived from a RealWorld block with defined units (tag (0008,0104) such as Hz or ms rather than "no units") or a RescaleType (0028,1054) with a non-US type (not defined by the standard), then the DV is already quantitative and cross scan comparison may be done.

However, in general DV is not sufficient to compare images from different scans, especially if the signal intensity varies a lot (eg multiple inversion recovery scans) in which case the FP value may be used as this may be compared (with some caveats) across scans and across timescales. This scaling requires an additional scale factor on top of the DV value, the Scale Slope (private tag (2005,100E))

dcm2niix duplicates (is this correct? Rather than losslessly) the raw data from DICOM to NIfTI format. These values are typically stored as 16-bit integers in the range -32768..32767. Both the DICOM and NIfTI formats describe how scaling intercept and slope values can be used to convert these raw values into calibrated values. For example, with an intercept of 0 and slope of 0.01 the raw value of 50 would be converted to 0.5.

The NIfTIhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fnifti.nimh.nih.gov%2Fpub%2Fdist%2Fsrc%2Fniftilib%2Fnifti1.h&data=04%7C01%7C%7C515789d6dcf94de5906208d8e4849037%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510605290623233%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=fFmkjcSmhy9HkDgwslsTftuoBr3DJa88byM3zMz%2BHEY%3D&reserved=0 header provides the scl_slope and scl_inter fields so each voxel value in the dataset is scaled as:

I = scl_slope * R + scl_inter lessy

where R is the voxel value stored and I is the "true" transformed voxel intensity.

Philips has three possible intensity transforms for their DICOM images (world (W), display (D), precise (P)). All of these transforms might be provided in a single DICOM image, while the NIfTI standard only designates a single scl_slope and scl_inter for each image. dcm2niix will retain the raw value (R) and sets the NIfTI scl_inter and scl_slope values for the desired intensity transform. dcm2niix will use W if possible. If this is not possible it will use P unless the user specifies -p n to disable the precise calculation resulting in D.

The formulas are provided below. The DICOM tags are in brackets (e.g. (0040,9225)) and the BIDS tag is in double quotes (e.g. "PhilipsRWVSlope"). Since all the scaling values are stored in the BIDS sidecar, you can always use these to generate your preferred intensity transform.

Inputs:

R = raw stored value of voxel in DICOM without scaling

WS = RealWorldValue slope (0040,9225) "PhilipsRWVSlope"

WI = RealWorldValue intercept (0040,9224) "PhilipsRWVIntercept"

RS = rescale slope (0028,1053) "PhilipsRescaleSlope"

RI = rescale intercept (0028,1052) "PhilipsRescaleIntercept"

SS = scale slope (2005,100E) "PhilipsScaleSlope"

Outputs:

W = real world value

P = precise value

D = displayed value

Formulas:

W = R * WS + WI

D = R * RS + RI

P = D / (RS * SS)

- You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Frordenlab%2Fdcm2niix%2Fissues%2F493%23issuecomment-796683536&data=04%7C01%7C%7C515789d6dcf94de5906208d8e4849037%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510605290643227%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=CaFtlV8dqGLmHtYO6uz%2B%2FtiTT9%2F7RaUoN1wcdoJRsQM%3D&reserved=0, or unsubscribehttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAIRGATYZ4J7TKMJTMPDVZOTTDCOR5ANCNFSM4ZAALTPA&data=04%7C01%7C%7C515789d6dcf94de5906208d8e4849037%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510605290653221%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=6WIxcYcDKj1jCsKlpGJdivOPrYbTnYxal9TJ2OkYbkY%3D&reserved=0.


The information contained in this message may be confidential and legally protected under applicable law. The message is intended solely for the addressee(s). If you are not the intended recipient, you are hereby notified that any use, forwarding, dissemination, or reproduction of this message is strictly prohibited and may be unlawful. If you are not the intended recipient, please contact the sender by return e-mail and destroy all copies of the original message.

drmclem commented 3 years ago

HI

I still think there may be a bug, or at least surprising behaviour. The problem is that the real world scaling with no units is not sufficient for cross scan comparison - either the scale slope needs to be included as well, or the user needs a third option - just using real world produces scans which can't be compared.

Essentially - the expected behaviour for users is to use all three values (RI,RS,SS) with RI and RS either coming from the traditional values or the RW values. I'll check on the edge case where RW has units .

Matthew

neurolabusc commented 3 years ago

So do you recommend the default behavior be to provide the transforms for P, and save the D if the user explicitly requests -p n? In other words, the W values are not desired? Alternatively, should W transform be applied if the units are known, e.g. Fieldmap in Hz?

drmclem commented 3 years ago

Hi

I need to look at an edge case test which will take me a couple of days but I think when using the P option it should always be the floating point value which always needs SS value.

Without it then I'm not sure the best version - I wasn't sure if -p = n gets you the raw pixels or the displayed values - my preference would be for display (so rescale only).

The edge case is what happens for, for example b0 maps, where the value is already in Hz - adding the SS may lose that - that I need to check - in which case not sure how dependently you want to handle it.

If "precise" is a synonym for "most comparable" then it should only add SS when the units are not specified- but let me check if this is case in the current software (it may be when we have units SS=1 and so this becomes a null coice)

Also bear in mind (we don't do this at the moment) but DICOM allows for multiple RealWorld blocks to allow the viewing system to switch between representations. In an ideal world this is how all scaling would be handled.

Matthew

From: Chris Rorden @.> Sent: Thursday, March 11, 2021 4:44 PM To: rordenlab/dcm2niix @.> Cc: Clemence, Matthew @.>; Mention @.> Subject: Re: [rordenlab/dcm2niix] Philips scaling with when RW units present not including scale slope ? (#493)

Caution: This e-mail originated from outside of Philips, be careful for phishing.

So do you recommend the default behavior be to provide the transforms for P, and save the D if the user explicitly requests -p n? In other words, the W values are not desired?

- You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Frordenlab%2Fdcm2niix%2Fissues%2F493%23issuecomment-796876218&data=04%7C01%7C%7C792885098aa94e542c2c08d8e4acd95c%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510778327652625%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=6Qj%2F%2FiahvwuftNz5N5F1cmXKo%2BlZ4%2Brkew%2BWYOrhJC8%3D&reserved=0, or unsubscribehttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAIRGAT3ZZBC4YBMXVBV3VRDTDDQLJANCNFSM4ZAALTPA&data=04%7C01%7C%7C792885098aa94e542c2c08d8e4acd95c%7C1a407a2d76754d178692b3ac285306e4%7C0%7C0%7C637510778327662621%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=t6B%2Bv%2F%2BzK7%2B4btRpc%2FFlCDQiTHWSk9ZAxQOjnro1F5E%3D&reserved=0.


The information contained in this message may be confidential and legally protected under applicable law. The message is intended solely for the addressee(s). If you are not the intended recipient, you are hereby notified that any use, forwarding, dissemination, or reproduction of this message is strictly prohibited and may be unlawful. If you are not the intended recipient, please contact the sender by return e-mail and destroy all copies of the original message.

neurolabusc commented 3 years ago

You can find a Philips Fieldmap here. The stored values range from 18..4084 (~12 bit). For this example, you get the same rescale regardless of if you choose W,D,P intensity transforms:

BIDS JSON:

    "PhilipsRWVSlope": 0.2442,
    "PhilipsRWVIntercept": -500,
    "PhilipsRescaleSlope": 0.2442,
    "PhilipsRescaleIntercept": -500,
    "PhilipsScaleSlope": 4.095,

Derived NIfTI rescales:

W = scl_slope 0.24420, scl_intercept -500
D = scl_slope 0.24420, scl_intercept -500
P = scl_slope 0.24420, scl_intercept -499.99994
drmclem commented 3 years ago

Hi

Think I have got to the bottom of this with some limited testings (so caveats etc.). To summarise the summary - if the image derives from a quantitaive method, or inline subtractions (ASL) then DV and FP scalings are identical - in otherwords applying the Scale Slope modification to the scaling is a NULL operation (give or take rounding error).

If however the dataset is an unaltered stream the applying the SS correction is needed to enable cross scan comparsion. So as of the moment (probably since release 5.6 and current 5.7)- adding the SS correction to the real world scaling values if present is required for normal scans for cross scan comparison (DV!=FP) and has no effect if the images are quantitive (DV==FP).

If you wanted to avoid the unecessary additional correction for on quantitive images (and not require the private tag 2005,100e in that case) it may be possible to identify them through the tag (0028,1054) RescaleType. Based on very limited checking so far this has the value "normalized" for images where DV!=FP and somthing else (US,Hz etc) when DV=FP.

This is only applicable for the most recent releases (5.6+) and applies whether the Realworld block is present or not.

Hope this is clear

Matthew

neurolabusc commented 3 years ago

@drmclem thanks. So are you happy for dcm2niix to provide the Precise scaling values (P, aka FP) by default and the Display values (D, aka DV) if the user selected -p n. In no case will dcm2niix provide the RealWorldValue scaling. Regardless of selection, the operands will be provided in the BIDS JSON file so a user can select any scaling factor they prefer later.

 D = R * RS + RI
 P = D / (RS * SS)
drmclem commented 3 years ago

Hi - as an aside I think I assumed that that -p n did no conversion at all and returned the raw integers from the file.

I think also conceptually this isn't an issue with RealWorldValue per se as it it is currently an alternative location to the standard rescale values - in the cases I have seen if RWorld is present so have been the rescale tags - but I know cases in the wild where rescale has been missing (due to other software) but the RWorld has been the alternative.

I think the issue is whether it is desirable to always apply the SS value - at the moment it deosn't matter as it is either a null operation (SS=1/RS) when units/subtractions have been applied or it is isn't when the FP conversion is necessary. Just thinking to the future proof it - ideally you wouldn't apply the SS if not required but detecting the units or the equivalence of the reciprocal is error prone. The problem may be in the future if multiple RW blocks are added, or the private tag being set to the reciprocal is lost in the RW case.

M

neurolabusc commented 3 years ago

The raw pixel data values is always losslessly copied from the DICOM to the NIfTI. The only thing that changes is the NIfTI scl_slope and scl_inter values. With -p n' these are set to the rescale slope (0028,1053)and rescale intercept (0028,1052) , identical to behavior to the other manufacturers. One the other hand, the precise values adjustscl_slopeandscl_inter` based on the Philips precise scale slope (2005,100E). This does not matter for most MRI data that uses relative signal, but can matter for sequences like ASL, see PMC3998685.