# pydicom_PIL.py """View DICOM images using Python image Library (PIL) Usage: >>> import pydicom >>> from pydicom.contrib.pydicom_PIL import show_PIL >>> ds = pydicom.dcmread("filename") >>> show_PIL(ds) Requires Numpy: http://numpy.scipy.org/ and Python Imaging Library: http://www.pythonware.com/products/pil/ """ # Copyright (c) 2009 Darcy Mason, Adit Panchal # This file is part of pydicom, relased under an MIT license. # See the file LICENSE included with this distribution, also # available at https://github.com/pydicom/pydicom # Based on image.py from pydicom version 0.9.3, # LUT code added by Adit Panchal # Tested on Python 2.5.4 (32-bit) on Mac OS X 10.6 # using numpy 1.3.0 and PIL 1.1.7b1 have_PIL = True try: import PIL.Image, PIL.ImageOps from pydicom.pixels import apply_voi_lut, apply_modality_lut, apply_windowing except ImportError: have_PIL = False have_numpy = True try: import numpy as np except ImportError: have_numpy = False def lin_stretch_img(img, low_prc, high_prc, do_ignore_minmax=True): """ Apply linear "stretch" - low_prc percentile goes to 0, and high_prc percentile goes to 255. The result is clipped to [0, 255] and converted to np.uint8 Additional feature: When computing high and low percentiles, ignore the minimum and maximum intensities (assumed to be outliers). """ # For ignoring the outliers, replace them with the median value if do_ignore_minmax: tmp_img = img.copy() med = np.median(img) # Compute median tmp_img[img == img.min()] = med tmp_img[img == img.max()] = med else: tmp_img = img lo, hi = np.percentile(tmp_img, (low_prc, high_prc)) # Example: 1% - Low percentile, 99% - High percentile if lo == hi: return np.full(img.shape, 128, np.uint8) # Protection: return gray image if lo = hi. stretch_img = (img.astype(float) - lo) * (255/(hi-lo)) # Linear stretch: lo goes to 0, hi to 255. stretch_img = stretch_img.clip(0, 255).astype(np.uint8) # Clip range to [0, 255] and convert to uint8 return stretch_img def get_LUT_value(data, window, level): """Apply the RGB Look-Up Table for the given data and window/level value.""" if not have_numpy: raise ImportError("Numpy is not available." "See http://numpy.scipy.org/" "to download and install") return np.piecewise(data, [data <= (level - 0.5 - (window - 1) / 2), data > (level - 0.5 + (window - 1) / 2)], [0, 255, lambda data: ((data - (level - 0.5)) / (window - 1) + 0.5) * (255 - 0)]) def get_PIL_image(dataset): """Get Image object from Python Imaging Library(PIL)""" if not have_PIL: raise ImportError("Python Imaging Library is not available. " "See http://www.pythonware.com/products/pil/ " "to download and install") if ('PixelData' not in dataset): raise TypeError("Cannot show image -- DICOM dataset does not have " "pixel data") # Handle multi-frame DICOM files by extracting the middle frame pixel_array = dataset.pixel_array if len(pixel_array.shape) == 3: pixel_array = pixel_array[pixel_array.shape[0] // 2] # Apply modality lut array = apply_modality_lut(pixel_array, dataset) # voi lut doesn't seem to help #array = apply_voi_lut(array, dataset) # can only apply LUT if these window info exists if ('WindowWidth' not in dataset) or ('WindowCenter' not in dataset): bits = dataset.BitsAllocated samples = dataset.SamplesPerPixel if bits == 8 and samples == 1: mode = "L" elif bits == 8 and samples == 3: mode = "RGB" elif bits == 16: # not sure about this -- PIL source says is 'experimental' # and no documentation. Also, should bytes swap depending # on endian of file and system?? mode = "I;16" else: raise TypeError("Don't know PIL mode for %d BitsAllocated " "and %d SamplesPerPixel" % (bits, samples)) # PIL size = (width, height) size = (dataset.Columns, dataset.Rows) # For multi-frame we already have the middle frame sliced as 2D in pixel_array array = pixel_array #array = apply_modality_lut(array, dataset) array = apply_voi_lut(array, dataset, index=0) #array = apply_windowing(array, dataset) array = lin_stretch_img(array, 0.1, 99.9) # Apply "linear stretching" (lower percentile 0.1 goes to 0, and percentile 99.9 to 255). if dataset["PhotometricInterpretation"].value == "MONOCHROME1": # We need to invert array = np.invert(array) return PIL.Image.fromarray(array).convert('L') else: ew = dataset['WindowWidth'] ec = dataset['WindowCenter'] ww = int(ew.value[0] if ew.VM > 1 else ew.value) wc = int(ec.value[0] if ec.VM > 1 else ec.value) ## Convert mode to L since LUT has only 256 values: ## http://www.pythonware.com/library/pil/handbook/image.htm image = get_LUT_value(array, ww, wc) im = PIL.Image.fromarray(image).convert('L') return im def get_PIL_image2(dataset): """Get Image object from Python Imaging Library(PIL)""" if not have_PIL: raise ImportError("Python Imaging Library is not available. " "See http://www.pythonware.com/products/pil/ " "to download and install") if ('PixelData' not in dataset): raise TypeError("Cannot show image -- DICOM dataset does not have " "pixel data") pixel_array = dataset.pixel_array if len(pixel_array.shape) == 3: pixel_array = pixel_array[pixel_array.shape[0] // 2] # can only apply LUT if these window info exists if ('WindowWidth' not in dataset) or ('WindowCenter' not in dataset): bits = dataset.BitsAllocated samples = dataset.SamplesPerPixel if bits == 8 and samples == 1: mode = "L" elif bits == 8 and samples == 3: mode = "RGB" elif bits == 16: # not sure about this -- PIL source says is 'experimental' # and no documentation. Also, should bytes swap depending # on endian of file and system?? mode = "I;16" else: raise TypeError("Don't know PIL mode for %d BitsAllocated " "and %d SamplesPerPixel" % (bits, samples)) im = PIL.Image.fromarray(pixel_array).convert('L') else: ew = dataset['WindowWidth'] ec = dataset['WindowCenter'] ww = int(ew.value[0] if ew.VM > 1 else ew.value) wc = int(ec.value[0] if ec.VM > 1 else ec.value) image = get_LUT_value(pixel_array, ww, wc) # Convert mode to L since LUT has only 256 values: # http://www.pythonware.com/library/pil/handbook/image.htm im = PIL.Image.fromarray(image).convert('L') return im def show_PIL(dataset): """Display an image using the Python Imaging Library (PIL)""" im = get_PIL_image(dataset) im.show()