Why is MATLAB rounding to whole number after multiplying and summing DICOM structure fields?

Working with DICOM files, thinking a meshgrid the best way to work with DICOM coordinates between structure and dose (please tell me if there's a better way), I'm trying to prepare a meshgrid to represent the DICOM coordinates, given the RTDOSE matrix with Attributes specifying the first element's DICOM coordinate:

>>                   doseinfo.ImagePositionPatient(1)
        ans = -191.1930
>> gridxmin = double(doseinfo.ImagePositionPatient(1))
   gridxmin = -191.1930
>> doseinfo.Width  * doseinfo.PixelSpacing(1)
        ans =  380
>> gridxmax=double(gridxmin+ doseinfo.Width *doseinfo.PixelSpacing(1))
   gridxmax =  189
>>                gridxmin + doseinfo.Width * doseinfo.PixelSpacing(1)
        ans =  189
>> gridxmax - gridxmin
        ans =  380.1930

As you can see, double appears not to be the problem, but rather, somehow multiplying causes the result to round. Why is it summing to 189 instead of the correct value 188.8070?

It may be worth noting that if I try inserting a factor of 1000 (to convert from millimeter to micrometer), the max value is 0 instead of 188807 (gridxmin = -191193, gridxmax = 0) -- although such a meshgrid can't be used anyway, as it throws the error

Error using repmat
Maximum variable size allowed by the program is exceeded.
Error in meshgrid (line 77)
      xx = repmat(xx, ny, 1, nz);

Something weird seems to be happening with this use of structures created from dicominfo. The problem is not resolved if I extract the data from the structure before calculation:

>> xstart = doseinfo.ImagePositionPatient(1);
xwidth = doseinfo.Width; xspacing = doseinfo.PixelSpacing(1);
gridxmin = xstart;
gridxmax = xstart + xwidth * xspacing
gridxmax =
  189

Answer 1

[Walter Roberson]

When you multiply a double by an integer data type, the result is calculated as double and then converted to the integer data type.

[Guillaume]

To clarify, one (or both) of your doseinfo.Width or doseinfo.PixelSpacing is an integer type. Type

class(doseinfo.Width)
class(doseinfo.PixelSpacing)

Matlab is one of the rare language where narrowing conversions are the default. That is when combining two classes of numbers, the one with more information gets converted to the one with less. So combine double with integer, the double gets converted to integer.

To avoid this, convert the integer to double before combining it with the double:

gridxmax = gridxmin + double(doseinfo.Width * doseinfo.PixelSpacing(1))

Note that, depending on the type of integer that is Width and pixelSpacing that operation could also overflow and be clamped to whatever is the maximum value of the smaller integer type in the expression. So, for safety, this would be better:

gridxmax = gridxmin + double(doseinfo.Width) * double(doseinfo.PixelSpacing(1))

[Daniel Bridges]

You mean the same thing happens when a uint16 is added to a double?

>> class(doseinfo.ImagePositionPatient(1))
ans = double

>> class(doseinfo.Width)
ans = uint16

>> class(doseinfo.PixelSpacing)
ans = double

Managing data types irks me. So DICOM image width is stored as an unsigned 16-bit integer while DICOM coordinates and Pixel Spacing are 64 bit floating point numbers.

" Arithmetic Operations on Integer Classes" demonstrates what you say about multiplication, and I suppose addition and multiplication are regarded as the same -- i.e. that MATLAB treats multiplication as addition. But what is the good in converting it back to the Integer class rather than leaving it as double?

[Walter Roberson]

Daniel, consider

K = uint8(0);
K = K + 1;

What data type is K afterwards? Your proposal would say that it should be double() afterwards because you suggest the operation of an integer and a double should result in a double. MATLAB's current answer is that it will be uint8 afterwards, because although the 0 + 1 might have been calculated in double, it will be narrowed to uint8 afterwards.

The automatic narrowing saves a lot of nuisance in writing code, allowing you to be sloppy about calculations like

K = K * 2 + 1

instead of having to write

K = K * uint8(2) + uint8(1)

[Guillaume]

"The automatic narrowing saves a lot of nuisance in writing code". I completely disagree with that! Particularly, for a weakly typed language, where you should not have to care about the type of your inputs.

Mathworks answer is that if you're using a narrower type it's on purpose and thus you want the result to be narrower but in practice, as is the case here, you're using a narrow type because that's what you've been given by whatever function you're calling.

The default conversion is to lose information. That does not sound like a good idea to me. Admittedly, I come from languages where widening conversion is the default, so I may be biased.

[Walter Roberson]

In C, does K * 2 + 1 mean K * (short int)(2) + (short int)1 with the short int widened if necessary depending on the type of K? Or does it mean K * (int)(2) + (int)(1) with the int widened as needed? Or does it mean K * (long int)(2) + (long int)(1) with the int converted to float or double if needed? If you have a long int, float is not necessarily widening...

It has been a while for me; I would need to re-check the C standards to find out exactly how the expression would be interpreted. As (int) is what my memory is suggesting.

[Guillaume]

In C (and C++), an integer literal with no suffix is of type int.

In C#, an integer literal with no suffix is of type int, uint, long, or ulong (whichever smallest type can fit the constant).

If K is narrower than int, then the result will be int. If K is wider, the result will be the type of K.