0

I'm reading DICOM gray image file as

gray = dicom.dcmread(file).pixel_array

There I've got (x,y) shape but I need RGB (x,y,3) shape

I'm trying to convert using CV

img = cv2.cvtColor(gray, cv2.COLOR_GRAY2RGB)

And for testing I'm writing it to file cv2.imwrite('dcm.png', img)

I've got extremely dark image on output which is wrong, what is correct way to convert pydicom image to RGB?

cnd
  • 32,616
  • 62
  • 183
  • 313

2 Answers2

4

To answer your question, you need to provide a bit more info, and be a bit clearer.

First what are you trying to do? Are you trying to only get an (x,y,3) array in memory? or are you trying to convert the dicom file to a .png file? ...they are very different things.

Secondly, what modality is your dicom image?

It's likely (unless its ultrasound or perhaps nuc med) a 16 bit greyscale image, meaning the data is 16 bit, meaning your gray array above is 16 bit data.

So the first thing to understand is window levelling and how to display a 16-bit image in 8 bits. have a look here: http://www.upstate.edu/radiology/education/rsna/intro/display.php.

If it's a 16-bit image, if you want to view it as a greyscale image in rgb format, then you need to know what window level you're using or need, and adjust appropriately before saving.

Thirdly, like lenik mention above, you need to apply the dicom slope/intercept values to your pixel data prior to using.


If your problem is just making a new array with extra dimension for rgb (so sizes (r,c) to (r,c,3)), then it's easy

# orig is your read in dcmread 2D array:
r, c = orig.shape

new = np.empty((w, h, 3), dtype=orig.dtype)

new[:,:,2] = new[:,:,1] = new[:,:,0] = orig
# or with broadcasting
new[:,:,:] = orig[:,:, np.newaxis]

That will give you the 3rd dimension. BUT the values will still all be 16-bit, not 8 bit as needed if you want it to be RGB. (Assuming your image you read with dcmread is CT, MR or equivalent 16-bit dicom - then the dtype is likely uint16).

If you want it to be RGB, then you need to convert the values to 8-bit from 16-bit. For that you'll need to decide on a window/level and apply it to select the 8-bit values from the full 16-bit data range.


Likely your problem above - I've got extremely dark image on output which is wrong - is actually correct, but it's dark because the window/level cv is using by default makes it 'look' dark, or it's correct but you didn't apply the slope/intercept.


If what you want to do is convert the dicom to png (or jpg), then you should probably use PIL or matplotlib rather than cv. Both of those offer easy ways to save a 16 bit 2D array (which is what you 'gray' is in your code above), both which allow you to specify window and level when saving to png or jpg. CV is complete overkill (meaning much bigger/slower to load, and much higher learning curve).

Some psueudo code using matplotlib. The vmin/vmax values you need to adjust - the ones here would be approximately ok for a CT image.


    import matplotlib.pyplot as plt

    df = dcmread(file)
    slope = float(df.RescaleSlope)
    intercept = float(df.RescaleIntercept)
    df_data = intercept + df.pixel_array * slope

    # tell matplotlib to 'plot' the image, with 'gray' colormap and set the
    # min/max values (ie 'black' and 'white') to correspond to 
    # values of -100 and 300 in your array
    plt.imshow(df_data, cmap='gray', vmin=-100, vmax=300)

    # save as a png file
    plt.savefig('png-copy.png')

that will save a png version, but with axes drawn as well. To save as just an image, without axes and no whitespace, use this:

    inches = (3,3)
    dpi = 150
    fig, ax = plt.subplots(figsize=inches, dpi=dpi)
    fig.subplots_adjust(left=0, right=1, top=1, bottom=0, wspace=0, hspace=0)
    ax.imshow(df_data, cmap='gray', vmin=-100, vmax=300)

    fig.save('copy-without-whitespace.png')

Richard
  • 3,024
  • 2
  • 17
  • 40
  • actually I need image in (x,y,3) shape, I've used to save just to test to see how it looks like, maybe I will need it anyways but not sure for now, even in case of saving I'd like to avoid using matplot, I didn't need any visualization. – cnd Nov 14 '19 at 14:44
  • Gotcha. I've added a section above on how to do that. Don't need matplotlib of course, but you will need to 'window/level' the array if what you want is true RGB (which is 8-bit values, not 16-bit dicom values) – Richard Nov 15 '19 at 06:11
  • yes I needed some RGB image which looks not as dark as with my OpenCV method – cnd Nov 15 '19 at 06:16
  • Then you need to adjust the vmin/vmax levels yourself to make it "look" brighter. Again this isn't a problem or limitation in the data, but simply an aspect of the visualisation parameters. Note this job of tweaking vmin/vmax is common/typical/needed when actually visualising dicom data - don't think the need to do this is something related to you not doing something correctly etc. – Richard Nov 15 '19 at 06:54
  • ...probably the easiest way for you to both see how this works, and also figure out which values you want exactly, are to download any of the many dicom file viewing apps, and interactively play with the WL level settings using that to open your file. – Richard Nov 15 '19 at 06:56
  • a bit weird, image looks correct if I do `df_data = intercept + df.pixel_array * slope / 20` (And 20 is random number here) in addition I need to invert it alike `img = 255 - img` and to be able to use `cvtColor` I need to use `img = df_data.astype(np.float32)` so in general that works however that 20 number is way too magical... – cnd Nov 15 '19 at 08:01
  • I'm not sure why that would change things. You definately shouldn't being change the slope/intercept calculation. Instead, when you display it change vmin/vmax to change the brightness/contrast. Honestly, try using a dicom file viewer that lets you interactively change window/level and it will help you understand what's going on before you try in python. – Richard Nov 15 '19 at 09:34
  • It's not a viewer, it's classification algorithm which works with x,y,3 rgb images, usually pngs or jpegs... so I'm trying to get DICOM into same state... – cnd Nov 15 '19 at 12:42
1

The full tutorial on reading DICOM files is here: https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial

Basically, you have to extract parameters slope and interception from the DICOM file and do the math for every pixel: hu = pixel_value * slope + intercept -- all this explained in the tutorial with the code samples and pictures.

lenik
  • 23,228
  • 4
  • 34
  • 43