pydicom and numpy

[sebhaase]

Hi,

I just heard about pydicom package on the pythonxy site.

I'm wondering if/how one would get a 2D or 3D image data set (like MRI
data or so) into a numpy array for further analysis ?
Also, would it be (at all, technically, in principle) be possible to
open that data in a memmap'ed way.

Thanks,
Sebastian Haase

[Darcy Mason]

On Jun 4, 4:43 pm, sebhaase <seb.ha...@gmail.com> wrote:
> I'm wondering if/how one would get a 2D or 3D image data set (like MRI
> data or so) into a numpy array for further analysis ?
> Also, would it be (at all, technically, in principle) be possible to
> open that data in a memmap'ed way.
>
> Thanks,
> Sebastian Haase

Hi Sebastian,

First, some ideas on working with in-memory structures:
A numpy array is returned by the Dataset's pixel_array property (but
note to write changes see [1]). It works for 2d data or for 3d in a
single file (at least radiotherapy dose, the only one I have used it
for). At present pixel_array does not work with compressed image data
(see the recent post 'DICOM multiframe' on this group).

For 3D data which has separate DICOM files for each slice, you could
read each one and get pixel_array as above; Numpy should have a way to
concatenate the slices into a proper 3d array.

Now for memmap:
You can avoid reading pixel data into memory using pydicom's read_file
() with the defer_size argument (new in pydicom 0.9.3), then point the
numpy memmap to the right location in the file -- something like the
following (untested):

import dicom
import numpy
ds = dicom.read_file("huge_file.dcm", defer_size="10 KB") # don't read
value of items > 10 KB
# Calc start of pixel data in file
# If ImplicitVR, is 8 bytes after start of tag, +4 more if Expl VR
and OB/OW
pixel_offset = ds.data_element("PixelData").data_element_tell + 8
if ds.isExplicitVR:
pixel_offset += 4
fp = numpy.memmap("huge_file.dcm", <other options...>,
offset=pixel_offset)

You might need to think about byte order depending on the endian-ness
of the DICOM file and your computer -- ds.isLittleEndian can tell you
which the file is.

Hopefully one of these ideas can work for you; if not, please feel
free to ask again with more details.
Darcy

[1] http://code.google.com/p/pydicom/wiki/WorkingWithPixelData