on obtaining phase image from rescaled DICOM file

[Alpay Özcan]

Good Afternoon,

I wanted to raise a concern regarding reading the phase values from the DICOM file

when using Read_Siemens_DICOM m-file, namely the lines:

ph = transpose(single(dicomread(filename)));

iPhase(:,:,rctr) = (ph*info2.RescaleSlope + info2.RescaleIntercept)/single(max(ph(:)))*pi;%phase

From what I can see Siemens puts the RescaleSlope=2 and RescaleIntercept=4096 for the unsigned integer values ranging from 0 to 4095 mapping the interval to [-4096 4094].

Although rare, but not improbable, if there is NO phase unwrapping in the data, the maximum of the unscaled values will remain below 4095 and the last m-file line above will stretch the interval thereby erroneously changing the phase values.

Notice also that the second m-file line, assuming phase wrapping occurs i.e. allowed maximum of 4095 was reached, produces an interval:

[-4096/4095 4094/4095] * pi = [-1.0002    0.9998] * pi

which is, in a sense adding  a 'bias' term to the phase values.

My suggestion would be to modify the scaling as:

( (pixelvalue-SmallestImagePixelValue)/(LargestImagePixelValue-SmallestImagePixelValue) *2 -1) * pi

without going through the rescaling procedure, directly from the DICOM pixel values.

I hope this helps in the future...

Thanks,

alpay

[Alpay Özcan]

Following up on the thread, when mapped to

[-4096/4095 4094/4095] * pi = [-1.0002    0.9998] * pi

the angle interval is in fact chopped.

Matlab outputs these values

angle(exp((4094/4095)*j*pi))=3.1408

angle(exp((-4096/4095)*j*pi))= 3.1408     -> or  -3.1424 < -pi

My recommendation of the previous message stands.

Please let me know if I am missing something. As it stands, this point looks important to me.

Thanks again,

alpay