Poisson resampling of DICOM images

First we need to import some libraries, %matplotlib inline ensures that our images are displayed in the notebook.



In [1]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import dicom
import os
import shutil 
import time
import numpy as np

This function generates a new unique UID for our new image



In [2]:

    
def generate_new_UID():
    UID_prefix = "1.2.826.0.1.3680043.8.707" # change this if your not me!
    currentTime = (str(time.time()))
    time1, time2 = currentTime.split(".")
    UID = UID_prefix +"."+ time1 + ".1"+time2  # need to include the .1 to ensure UID element cannot start with 0
    time.sleep(2) # wait 1 second to ensure ech UID is different
    return UID

The next two functions perform a Poisson resampling of the image data. The first resamples each individual pixel. This was the way I did it in my python script that I wrote for the MSc project:



In [3]:

    
def half_time_values(oldImage):
	newImage = np.zeros_like(oldImage)
	for i in range(0, oldImage.shape[0]):
		for j in range(0, oldImage.shape[1]):
			origValue = oldImage[i,j]
			newValue=(origValue/2)
			Value = np.random.poisson(newValue)
			newImage[i,j] = Value		# write the value to newImage
	return newImage

Some time later I realised that you could perform the resampling on the complete array (which should be much faster)



In [4]:

    
def half_time_array(oldImage):
    newImage=(oldImage/2)
    newImage = np.random.poisson(newImage)
    return newImage

First we need to import our dataset using the dicom library



In [5]:

    
ds = dicom.read_file("testobject.dcm")

We can then display this using imshow on the pixel_array



In [6]:

    
image1 = plt.imshow(ds.pixel_array, cmap='bone')
image1.axes.xaxis.set_visible(False)
image1.axes.yaxis.set_visible(False);

max() gives us the maximum pixel value for this image:



In [7]:

    
ds.pixel_array.max()









    Out[7]:





16244

the ipython %timeit command allows use to time the execution of each of our poisson resamplaing routines:



In [8]:

    
%timeit newImage = half_time_values(ds.pixel_array)
%timeit newImage2 = half_time_array(ds.pixel_array)









    



1 loops, best of 3: 9.67 s per loop
10 loops, best of 3: 135 ms per loop

Doing the resampling on the entire array is much faster. We can then create a resampled copy of our original image



In [9]:

    
newImage2 = half_time_array(ds.pixel_array)
image2 = plt.imshow(newImage2, cmap='bone')



In [10]:

    
newImage2.max()









    Out[10]:





8441

The maximum value in our image is now approximately half the value of the original as expected



In [11]:

    
newImage2.min()









    Out[11]:





6556



In [ ]: