Biomedical Images_1

 Imageio

# Import ImageIO

import imageio

import matplotlib.pyplot as plt 

# Load "chest-220.dcm"

im = imageio.imread("chest-220.dcm")

# Print image attributes

print('Image type:', type(im))

print('Shape of image array:', im.shape)

# Print the available metadata fields

print(im.meta.keys())

print(im.meta['PatientSex']) 

# Draw the image with greater contrast

plt.imshow(im, cmap ='gray', vmin = -200, vmax = 200)

# Remove ticks and axises

plt.axis('off')

# Render the image

plt.show()

# Import ImageIO and NumPy

import imageio

import numpy as np

# Read in each 2D image

im1 = imageio.imread('chest-220.dcm')

im2 = imageio.imread('chest-221.dcm')

im3 = imageio.imread('chest-222.dcm')

# Stack images into a volume

vol = np.stack([im1, im2, im3])

print('Volume dimensions:', vol.shape)

 

# Load the "tcia-chest-ct" directory

vol = imageio.volread("tcia-chest-ct")

# Print image attributes

print('Available metadata:', vol.meta.keys())

print('Shape of image array:', vol.shape)

 

In [5]:

vol.meta['sampling']

Out[5]:

(3.270000000000001, 0.976562, 0.976562)

In [6]:

vol.shape

Out[6]:

(25, 512, 512)

 The field of view = (82, 500, 500)

 

To select a 2D frame, pick a frame for the first axis and select all data 

from the remaining two: vol[0, :, :]

For this exercise, use for loop to plot every 40th slice of vol on a 

separate subplot. matplotlib.pyplot (as plt) has been imported for you.

# Plot the images on a subplots array

fig, axes = plt.subplots(nrows=1, ncols=4)

# Loop through subplots and draw image

for ii in range(4):

im = vol[ii*40, :, :]

axes[ii].imshow(im, cmap='gray')

axes[ii].axis('off')

 

# Select frame from "vol"

im1 = vol[:, 256, :]

im2 = vol[:, :, 256]

# Compute aspect ratios

d0, d1, d2 = vol.meta['sampling']

asp1 = d0 / d2

asp2 = d0/d1

# Plot the images on a subplots array

fig, axes = plt.subplots(nrows=2, ncols=1)

axes[0].imshow(im1, cmap='gray',aspect=asp1)

axes[1].imshow(im2, cmap='gray', aspect=asp2)

plt.show()

 

Intensity: 

# Import SciPy's "ndimage" module

import scipy.ndimage as ndi

# Create a histogram, binned at each possible value

# Change data type could save spaces.

# im_int8 = im.astype(np.unit8)

hist = ndi.histogram(im, min=0, max=256, bins=256)

# Create a cumulative distribution function

cdf = hist.cumsum() / hist.sum()

# Plot the histogram and CDF

fig, axes = plt.subplots(2, 1, sharex=True)

axes[0].plot(hist, label='Histogram')

axes[1].plot(cdf, label='CDF')

format_and_render_plot()

Mask:

# Create skin and bone masks

mask_bone = im>=145

mask_skin = (im>=45)&(im<145)

# Plot the skin (0) and bone (1) masks

fig, axes = plt.subplots(1,2)

axes[0].imshow(mask_skin, cmap='gray')

axes[1].imshow(mask_bone, cmap='gray')

format_and_render_plot()

 

# Import SciPy's "ndimage" module

import scipy.ndimage as ndi

# Screen out non-bone pixels from "im"

mask_bone = im >= 145

im_bone = np.where(mask_bone, im, 0)

# Get the histogram of bone intensities

hist = ndi.histogram(im_bone, min=1, max=255,bins=255)

# Plot masked image and histogram

fig, axes = plt.subplots(2,1)

axes[0].imshow(im_bone)

axes[1].plot(hist)

format_and_render_plot()

 

 

 

 

 

 

 

# Create and tune bone mask

mask_bone = im >=145

mask_dilate = ndi.binary_dilation(mask_bone, iterations=5)

mask_closed = ndi.binary_closing(mask_bone, iterations=5)

# Plot masked images

fig, axes = plt.subplots(1,3)

axes[0].imshow(mask_bone)

axes[1].imshow(mask_dilate)

axes[2].imshow(mask_closed)

format_and_render_plot()

Filters:

# Set filter weights

weights = [[0.11, 0.11, 0.11],

[0.11, 0.11, 0.11],

[0.11, 0.11, 0.11]]

# Convolve the image with the filter

im_filt = ndi.convolve(im, weights)

# Plot the images

fig, axes = plt.subplots(1,2)

axes[0].imshow(im)

axes[1].imshow(im_filt)

format_and_render_plot()

 

 

 

 

 

 

 

 

 

# Smooth "im" with Gaussian filters

im_s1 = ndi.gaussian_filter(im, sigma=1)

im_s3 = ndi.gaussian_filter(im, sigma=3)

# Draw bone masks of each image

fig, axes = plt.subplots(1,3)

axes[0].imshow(im >= 145)

axes[1].imshow(im_s1 >= 145)

axes[2].imshow(im_s3 >= 145)

format_and_render_plot()

 

 

 

 

Edges:
 

# Set weights to detect vertical edges

weights = [[1,0,-1], [1,0,-1], [1,0,-1]]

# Convolve "im" with filter weights

edges = ndi.convolve(im, weights)

# Draw the image in color

plt.imshow(edges, cmap = 'seismic', vmin=-150, vmax=150)

plt.colorbar()

format_and_render_plot()

 

 

 

 

 

 

 

 

 

 

# Apply Sobel filter along both axes

sobel_ax0 = ndi.sobel(im, axis=0)

sobel_ax1 = ndi.sobel(im, axis=1)

# Calculate edge magnitude

edges = np.sqrt(np.square(sobel_ax0) + np.square(sobel_ax1))

# Plot edge magnitude

plt.imshow(edges, cmap='gray', vmax=75)

format_and_render_plot()