dicom saved can't match the plt.show

[tanjia123456]

Hello, I want generated new dicom and to save, But what is saved is not the same as what is displayed with plt.show

# the orig dicom file
print("pixel_bytes=",type(ds.PixelData), len(ds.PixelData), "pix_array=",ds.pixel_array.shape)
# pixel_bytes= <class 'bytes'>  0\x00\x00\x00\x00\x00\x00\x00\x00\x00\x
# pixel_bytes= <class 'bytes'> 228160 pix_array= (368, 310)     368*310*2 = 228160 

# newimage is generated by Gernerator in GAN
ds.PixelData = newimg.tobytes()  # 替换矩阵
print("test1=", ds.pixel_array.shape, len(ds.PixelData))  # test1= (368, 310) 524288
ds.Rows, ds.Columns = newimg.shape
print("test2=", ds.pixel_array.shape, len(ds.PixelData)) # test2= (512, 512) 524288

ds.save_as(os.path.join(out_path2 + name))
# pydicom.dcmwrite(out_path2 + name, ds)

pyplot.imshow(ds.pixel_array,  cmap="gray")# cmap=pyplot.cm.bone
pyplot.show()

the plt.show is

, the saved dicom is

how to solve the problem?
pydicom 2.4.4
Python 3.7.16
@pdurbin @volrath @yarikoptic @cancan101

[scaramallion]

Please post the group 0x0028 element values: print(ds.group_dataset(0x0028))

Does your newimg ndarray's dtype match the dataset's Bits Allocated, Bits Stored and Pixel Representation values?

[tanjia123456]

yes, ds[0x0028, 0x0100].value is 16, So I converted pred to match it

if ds[0x0028, 0x0100].value == 16:  # 如果dicom文件矩阵是16位格式
    newimg = newimg.astype(np.int16)  # newimg 是图像矩阵 ds是dcm uint16
elif ds[0x0028, 0x0100].value == 8:
    newimg = newimg.astype(np.int8)
else:
    raise Exception("unknow Bits Allocated value in dicom header")

[scaramallion]

And Pixel Representation?

[scaramallion]

Can you post a minimum reproducible code example (including test data)? It's really difficult for anyone to figure out the problem without one.

Please also post the values of the group 0x0028 elements...

[tanjia123456]

he values of the group 0x0028 elements：

ds.group_dataset(0x0028)= (0028, 0002) Samples per Pixel US: 1
(0028, 0004) Photometric Interpretation CS: 'MONOCHROME2'
(0028, 0010) Rows US: 368
(0028, 0011) Columns US: 310
(0028, 0030) Pixel Spacing DS: [0.48913043737412, 0.48913043737412]
(0028, 0100) Bits Allocated US: 16
(0028, 0101) Bits Stored US: 12
(0028, 0102) High Bit US: 11
(0028, 0103) Pixel Representation US: 0
(0028, 0106) Smallest Image Pixel Value US: 0
(0028, 0107) Largest Image Pixel Value US: 1300
(0028, 1050) Window Center DS: '230.0'
(0028, 1051) Window Width DS: '565.0'
(0028, 1055) Window Center & Width Explanation LO: 'Algo1'
pixel_bytes= <class 'bytes'> 228160 pix_array= (368, 310)

the code example as follow:

import pydicom
import numpy as np
import os
from matplotlib import pyplot as plt
import lpips
from matplotlib import pyplot

ds = pydicom.dcmread("data/CTA-GAN/PA1/SE1/IM0.IMA", force=True) # 读取头文件
dsa = pydicom.uid.generate_uid()
ds.SeriesInstanceUID = dsa

The newimg generated by my network
newimg = np.ones((512, 512))

newimg[newimg == 0] = -2000
if ds[0x0028, 0x0100].value == 16: # 如果dicom文件矩阵是16位格式
    newimg = newimg.astype(np.int16)
    print("newimg 是图像矩阵 ds是dcm uint16")
elif ds[0x0028, 0x0100].value == 8:
    newimg = newimg.astype(np.int8)
else:
    raise Exception("unknow Bits Allocated value in dicom header")

ds.PixelData = newimg.tobytes() # 替换矩阵
print("test1=", ds.pixel_array.shape, len(ds.PixelData)) # test1= (368, 310) 524288
ds.Rows, ds.Columns = newimg.shape
print("test2=", ds.pixel_array.shape, len(ds.PixelData)) # test2= (512, 512) 524288

ds.save_as('test.IMA')

pyplot.imshow(ds.pixel_array, cmap=plt.cm.gray) # cmap=pyplot.cm.bone
pyplot.show()

[scaramallion]

Issues

1. (0028, 0103) Pixel Representation US: 0: This indicates the original pixel data uses unsigned integers, but your newimg ndarray is signed, so you should change ds.PixelRepresentation to 1

2. You can't simply recast newimg using ndarray.astype() as you will get overflow issues (plus you haven't said what dtype newimg is originally):

>>> np.asarray([0, 32768], dtype="uint16").astype("uint8")
array([0, 0], dtype=uint8)

Other

It's much easier to access dataset elements using:

if ds.BitsAllocated == 16:  # instead of ds[0x0028, 0x0100].value
    # ...

As long as there are sufficient bytes in Pixel Data to cover the image size you won't get an error as the DICOM Standard allows the use of excess padding:

>>> ds.BitsAllocated = 8
>>> ds.PixelData = b"\x00" * 36
>>> ds.Rows, ds.Columns = 5, 5
# OK, sufficient bytes available for a 25 pixel image
>>> ds.pixel_array  # equivalent to np.frombuffer(ds.PixelData[:25]).reshape(5, 5)
# Raises exception, not enough bytes for a 40 pixel image
>>> ds.Rows = 8
>>> ds.pixel_array

[tanjia123456]

the values of the group 0x0028 elements：

ds.group_dataset(0x0028)= (0028, 0002) Samples per Pixel                   US: 1
(0028, 0004) Photometric Interpretation          CS: 'MONOCHROME2'
(0028, 0010) Rows                                US: 368
(0028, 0011) Columns                             US: 310
(0028, 0030) Pixel Spacing                       DS: [0.48913043737412, 0.48913043737412]
(0028, 0100) Bits Allocated                      US: 16
(0028, 0101) Bits Stored                         US: 12
(0028, 0102) High Bit                            US: 11
(0028, 0103) Pixel Representation                US: 0
(0028, 0106) Smallest Image Pixel Value          US: 0
(0028, 0107) Largest Image Pixel Value           US: 1300
(0028, 1050) Window Center                       DS: '230.0'
(0028, 1051) Window Width                        DS: '565.0'
(0028, 1055) Window Center & Width Explanation   LO: 'Algo1'
pixel_bytes= <class 'bytes'> 228160 pix_array= (368, 310)

the code example as follow:

import pydicom
import numpy as np
import os
from matplotlib import pyplot as plt
import lpips
from matplotlib import pyplot


ds = pydicom.dcmread("data/CTA-GAN/PA1/SE1/IM0.IMA", force=True)  # 读取头文件
dsa = pydicom.uid.generate_uid()
ds.SeriesInstanceUID = dsa

# The newimg generated by my network  
newimg = np.ones((512, 512))

# newimg[newimg == 0] = -2000
if ds[0x0028, 0x0100].value == 16:  # 如果dicom文件矩阵是16位格式
    newimg = newimg.astype(np.int16)
    print("newimg 是图像矩阵 ds是dcm uint16")
elif ds[0x0028, 0x0100].value == 8:
    newimg = newimg.astype(np.int8)
else:
    raise Exception("unknow Bits Allocated value in dicom header")


ds.PixelData = newimg.tobytes()  # 替换矩阵
print("test1=", ds.pixel_array.shape, len(ds.PixelData))  # test1= (368, 310) 524288
ds.Rows, ds.Columns = newimg.shape
print("test2=", ds.pixel_array.shape, len(ds.PixelData))  # test2= (512, 512) 524288

ds.save_as('test.IMA')

pyplot.imshow(ds.pixel_array, cmap=plt.cm.gray)  # cmap=pyplot.cm.bone
pyplot.show()`

[scaramallion]

Bits Allocated is 16, so its 2 bytes per pixel: 512 x 512 x 2 = 524288

[tanjia123456]

I don't quite understand your second point, which means that directly using newimg = newimg.astype(np.int8) is too rough and will cause memory overflow.

So I changed it to the following code, and the problem still occurs that the images saved and displayed do not match

` ds.SeriesInstanceUID = dsa
# print("ds.group_dataset(0x0028)=", ds.group_dataset(0x0028))
shutil.copy(file_path, file_path0)
shutil.copy(file_path.replace('SE0', 'SE1'), file_path1)

                # newimg[newimg == 0] = -2000
                if ds.BitsAllocated == 16:  # 如果dicom文件矩阵是16位格式
                    **# newimg = newimg.astype(np.int16)  # newimg 是图像矩阵 ds是dcm uint16
                    newimg = np.asarray(newimg, dtype="uint16").astype(np.int16)**
                elif ds.BitsAllocated == 8:
                    # newimg = newimg.astype(np.int8)
                    newimg = np.asarray(newimg, dtype="uint8").astype(np.int8)
                else:
                    raise Exception("unknow Bits Allocated value in dicom header")

                ds.PixelData = newimg.tobytes()  # 替换矩阵
                ds.PixelRepresentation = 1
                ds.Rows, ds.Columns = newimg.shape
                ds.save_as(os.path.join(out_path2 + name)) # pydicom.dcmwrite(out_path2 + name, ds)
                pyplot.imshow(ds.pixel_array,   cmap=plt.cm.gray)# cmap=pyplot.cm.bone
                pyplot.show()`

[scaramallion]

Add casting="safe" to your original code and see what happens.

# Still no mention of what dtype `newimg` actually is
# Is it bool? float? uint? int? What's the itemsize?
print("Before casting", newimg.dtype, newimg.min(), newimg.max())
if ds.BitsAllocated == 16:
    newimg = newimg.astype(np.int16, casting="safe")
elif ds.BitsAllocated == 8:
    newimg = newimg.astype(np.int8, casting="safe")

print("After casting", newimg.dtype, newimg.min(), newimg.max())

[tanjia123456]

It does seem to be a data type error, so what should i do to solve the problem？

Before casting float32 0.0 2076.7463 Traceback (most recent call last): File "E:\CT2CTA\CTA-GAN-main1\train.py", line 53, in <module> main() File "E:\CT2CTA\CTA-GAN-main1\train.py", line 47, in main trainer.test() File "E:\CT2CTA\CTA-GAN-main1\trainer\Hd_Trainer_x.py", line 491, in test newimg = newimg.astype(np.int16, casting="safe") TypeError: Cannot cast array data from dtype('float32') to dtype('int16') according to the rule 'safe'

[scaramallion]

I have no idea what the appropriate conversion process is, it depends on how newimg is being created and what it represents.

Generally though, you'll want to rescale newimg to a suitable range such as (0.0 to 1.0), then convert that to the value range of uint8 or int16 or whatever it is you want, then use numpy.clip() and finally ndarray.astype().

[tanjia123456]

Also , I want know how the resample the dicom file 0.67 × 0.67 × 1.25-mm3 volumes and consisted of 75–490 sections with a matrix size of 512 × 512 pixels, Or a fixed space and size

`
pixelspace1= []
slicethick1= []
rows1 = []
column1 = []

for i in range(1, 11):

ds = pydicom.dcmread("data/CTA-GAN/PA"+str(i)+"/SE1/IM0.IMA", force=True)  # 读取头文件

# (0028, 0030) Pixel Spacing                       DS: [0.48913043737412, 0.48913043737412]
pixelspace1.append(ds.PixelSpacing)
slicethick1.append(ds.SliceThickness)
rows1.append(ds.Rows) #图像的宽
column1.append(ds.Columns) #图像的长

print(pixelspace1)
print(slicethick1)
print(rows1)
print(column1)
`

the orig info is:

[[0.48913043737412, 0.48913043737412], [0.48913043737412, 0.48913043737412], [0.48913043737412, 0.48913043737412], [0.48913046717644, 0.48913046717644], [0.48913043737412, 0.48913043737412], [0.48913043737412, 0.48913043737412], [0.48913043737412, 0.48913043737412], [0.48913043737412, 0.48913043737412], [0.48913046717644, 0.48913046717644], [0.48913043737412, 0.48913043737412]]
['0.5', '0.5', '0.5', '0.5', '0.5', '0.5', '0.5', '0.5', '0.5', '0.5']
[368, 368, 368, 368, 368, 368, 368, 368, 368, 368]
[310, 320, 310, 300, 310, 310, 310, 320, 300, 310]

[tanjia123456]

okk thanks, i will change it according to my need~hhh

…

---Original--- From: ***@***.***> Date: Mon, Feb 26, 2024 16:18 PM To: ***@***.***>; Cc: ***@***.******@***.***>; Subject: Re: [pydicom/pydicom] dicom saved can't match the plt.show(Discussion #2017) I have no idea what the appropriate conversion process is, it depends on how newimg is being created and what it represents. Generally though, you'll want to rescale newimg to a suitable range such as (0.0 to 1.0), then convert that to the value range of uint8 or int16 or whatever it is you want, then use numpy.clip() and finally ndarray.astype(). — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: ***@***.***>