Multidimentional data - Matrices and Images



In [ ]:

    
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

from scipy import linalg



In [ ]:

    
plt.style.use('ggplot')
plt.rc('axes', grid=False)   # turn off the background grid for images

Let us work with the matrix: $ \left[ \begin{array}{cc} 1 & 2 \\ 1 & 1 \end{array} \right] $



In [ ]:

    
my_matrix = np.array([[1,2],[1,1]])

print(my_matrix.shape)



In [ ]:

    
print(my_matrix)



In [ ]:

    
my_matrix_transposed = np.transpose(my_matrix)

print(my_matrix_transposed)



In [ ]:

    
my_matrix_inverse = linalg.inv(my_matrix)

print(my_matrix_inverse)

numpy matrix multiply uses the `dot()` function:



In [ ]:

    
my_matrix_inverse.dot(my_matrix)

Caution the `*` will just multiply the matricies on an element-by-element basis:



In [ ]:

    
my_matrix_inverse * my_matrix_inverse

Solving system of linear equations

$$ \begin{array}{c} x + 2y = 4 \\ x + y = 3 \\ \end{array} \hspace{2cm} \left[ \begin{array}{cc} 1 & 2 \\ 1 & 1 \\ \end{array} \right] \left[ \begin{array}{c} x\\ y \end{array} \right] = \left[ \begin{array}{c} 4\\ 3\\ \end{array} \right] \hspace{2cm} {\bf A}x = {\bf b} \hspace{2cm} \left[ \begin{array}{c} x\\ y \end{array} \right] = \left[ \begin{array}{cc} 1 & 2 \\ 1 & 1 \\ \end{array} \right]^{-1} \left[ \begin{array}{c} 4\\ 3\\ \end{array} \right] = \left[ \begin{array}{c} 2\\ 1 \end{array} \right] $$



In [ ]:

    
A = np.array([[1,2],[1,1]])

print(A)



In [ ]:

    
b = np.array([[4],[3]])

print(b)



In [ ]:

    
# Solve by inverting A and then mulitply by b

linalg.inv(A).dot(b)



In [ ]:

    
# Cleaner looking

linalg.solve(A,b)

System of 3 equations

$$ \begin{array}{c} x + 3y + 5z = 10 \\ 2x + 5y + z = 8 \\ 2x + 3y + 8z = 3 \\ \end{array} \hspace{3cm} \left[ \begin{array}{ccc} 1 & 3 & 5 \\ 2 & 5 & 1 \\ 2 & 3 & 8 \end{array} \right] \left[ \begin{array}{c} x\\ y\\ z \end{array} \right] = \left[ \begin{array}{c} 10\\ 8\\ 3 \end{array} \right] $$



In [ ]:

    
A = np.array([[1,3,5],[2,5,1],[2,3,8]])
b = np.array([[10],[8],[3]])

print(linalg.inv(A))

print(linalg.solve(A,b))

Images are just 2-d arrays - `imshow` will display 2-d arrays as images



In [ ]:

    
print(A)

plt.imshow(A, interpolation='nearest', cmap=plt.cm.Blues);

Read in some data



In [ ]:

    
test_image = np.load("./MyData/test_data.npy")    # load in a saved numpy array



In [ ]:

    
test_image.ndim, test_image.shape, test_image.dtype



In [ ]:

    
print("The minimum value of the image is {0:.2f}".format(test_image.min()))
print("The maximum value of the image is {0:.2f}".format(test_image.max()))
print("The mean value of the image is {0:.2f}".format(test_image.mean()))
print("The standard deviation of the image is {0:.2f}".format(test_image.std()))



In [ ]:

    
#flatten() collapses n-dimentional data into 1-d

plt.hist(test_image.flatten(),bins=30);

Math on images applies to every value (pixel)



In [ ]:

    
another_test_image = test_image + 8

print("The minimum value of the other image is {0:.2f}".format(another_test_image.min()))
print("The maximum value of the other image is {0:.2f}".format(another_test_image.max()))
print("The mean value of the other image is {0:.2f}".format(another_test_image.mean()))
print("The standard deviation of the other image is {0:.2f}".format(another_test_image.std()))

Show the image represenation of `I` with a colorbar



In [ ]:

    
plt.imshow(test_image, cmap=plt.cm.gray)
plt.colorbar();

Colormap reference: http://matplotlib.org/examples/color/colormaps_reference.html



In [ ]:

    
fig, ax = plt.subplots(1,5,sharey=True)

fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].imshow(test_image, cmap=plt.cm.viridis)
ax[0].set_xlabel('viridis')

ax[1].imshow(test_image, cmap=plt.cm.hot)
ax[1].set_xlabel('hot')

ax[2].imshow(test_image, cmap=plt.cm.magma)
ax[2].set_xlabel('magma')

ax[3].imshow(test_image, cmap=plt.cm.spectral)
ax[3].set_xlabel('spectral')

ax[4].imshow(test_image, cmap=plt.cm.gray)
ax[4].set_xlabel('gray')

WARNING! Common image formats DO NOT preserve dynamic range of original data!!

Common image formats: jpg, gif, png, tiff
Common image formats will re-scale your data values to [0:1]
Common image formats are NOT suitable for scientific data!



In [ ]:

    
plt.imsave('Splash.png', test_image, cmap=plt.cm.gray)     # Write the array I to a PNG file

my_png = plt.imread('Splash.png')                   # Read in the PNG file

print("The original data has a min = {0:.2f} and a max = {1:.2f}".format(test_image.min(), test_image.max()))

print("The PNG file has a min = {0:.2f} and a max = {1:.2f}".format(my_png.min(), my_png.max()))

Creating images from math



In [ ]:

    
X = np.linspace(-5, 5, 500)
Y = np.linspace(-5, 5, 500)

X, Y = np.meshgrid(X, Y)     # turns two 1-d arrays (X, Y) into one 2-d grid

Z = np.sqrt(X**2+Y**2)+np.sin(X**2+Y**2)

Z.min(), Z.max(), Z.mean()

Fancy Image Display



In [ ]:

    
from matplotlib.colors import LightSource

ls = LightSource(azdeg=0,altdeg=40)
shadedfig = ls.shade(Z,plt.cm.copper)

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

fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].imshow(shadedfig)

contlevels = [1,2,Z.mean()]

ax[1].axis('equal')
ax[1].contour(Z,contlevels)

ax[2].imshow(shadedfig)
ax[2].contour(Z,contlevels);

Reading in images (`imread`) - Common Formats



In [ ]:

    
my_doctor = plt.imread('./MyData/doctor5.png')

print("The image my_doctor has a shape [height,width] of {0}".format(my_doctor.shape))
print("The image my_doctor is made up of data of type {0}".format(my_doctor.dtype))
print("The image my_doctor has a maximum value of {0}".format(my_doctor.max()))
print("The image my_doctor has a minimum value of {0}".format(my_doctor.min()))



In [ ]:

    
plt.imshow(my_doctor,cmap=plt.cm.gray);

Images are just arrays that can be sliced.

### For common image formats the origin is the upper left hand corner



In [ ]:

    
fig, ax = plt.subplots(1,4)
fig.set_size_inches(12,6)

fig.tight_layout()

# You can show just slices of the image - Rememeber: The origin is the upper left corner

ax[0].imshow(my_doctor, cmap=plt.cm.gray)
ax[0].set_xlabel('Original')

ax[1].imshow(my_doctor[0:300,0:100], cmap=plt.cm.gray)
ax[1].set_xlabel('[0:300,0:100]')                 # 300 rows, 100 columns

ax[2].imshow(my_doctor[:,0:100], cmap=plt.cm.gray)       # ":" = whole range
ax[2].set_xlabel('[:,0:100]')                     # all rows, 100 columns

ax[3].imshow(my_doctor[:,::-1], cmap=plt.cm.gray);
ax[3].set_xlabel('[:,::-1]')                      # reverse the columns



In [ ]:

    
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)

fig.tight_layout()

CutLine = 300

ax[0].imshow(my_doctor, cmap=plt.cm.gray)
ax[0].hlines(CutLine, 0, 194, color='b', linewidth=3)

ax[1].plot(my_doctor[CutLine,:], color='b', linewidth=3)
ax[1].set_xlabel("X Value")
ax[1].set_ylabel("Pixel Value")

Simple image manipulation



In [ ]:

    
from scipy import ndimage



In [ ]:

    
fig, ax = plt.subplots(1,5)
fig.set_size_inches(14,6)

fig.tight_layout()

ax[0].imshow(my_doctor, cmap=plt.cm.gray)

my_doctor_2 = ndimage.rotate(my_doctor,45,cval=0.75)               # cval is the value to set pixels outside of image
ax[1].imshow(my_doctor_2, cmap=plt.cm.gray)                 # Rotate and reshape

my_doctor_3 = ndimage.rotate(my_doctor,45,reshape=False,cval=0.75) # Rotate and do not reshape
ax[2].imshow(my_doctor_3, cmap=plt.cm.gray)

my_doctor_4 = ndimage.shift(my_doctor,(10,30),cval=0.75)           # Shift image      
ax[3].imshow(my_doctor_4, cmap=plt.cm.gray)

my_doctor_5 = ndimage.gaussian_filter(my_doctor,5)                # Blur image
ax[4].imshow(my_doctor_5, cmap=plt.cm.gray);

`ndimage` can do much more: http://scipy-lectures.github.io/advanced/image_processing/

FITS file (Flexible Image Transport System) - Standard Astro File Format

FITS format preserves dynamic range of data
FITS format can include lists, tables, images, and combunations of different types of data
FITS fiels consiste of at least two parts - A Header and Data



In [ ]:

    
import astropy.io.fits as fits



In [ ]:

    
my_image_file = "./MyData/bsg01.fits"

image_data = fits.getdata(my_image_file)
image_header = fits.getheader(my_image_file)



In [ ]:

    
image_header



In [ ]:

    
print("The image has a shape [height,width] of {0}".format(image_data.shape))
print("The image is made up of data of type {0}".format(image_data.dtype))
print("The image has a maximum value of {0}".format(image_data.max()))
print("The image has a minimum value of {0}".format(image_data.min()))



In [ ]:

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

fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].imshow(image_data,cmap=plt.cm.gray)

ax[1].hist(image_data.flatten(),bins=20);

You can use masks on images



In [ ]:

    
CopyData = np.copy(image_data)

CutOff = 40

mask = np.where(CopyData > CutOff)
CopyData[mask] = 50                   # You can not just throw data away, you have to set it to something.

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

fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].imshow(CopyData,cmap=plt.cm.gray)

ax[1].hist(CopyData.flatten(),bins=20);

You can add and subtract images



In [ ]:

    
another_image_file = "./MyData/bsg02.fits"

another_image_data = fits.getdata(another_image_file)



In [ ]:

    
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].imshow(image_data, cmap=plt.cm.gray)
ax[1].imshow(another_image_data, cmap=plt.cm.gray);

The two images above may look the same but they are not! Subtracting the two images reveals the truth.



In [ ]:

    
fig, ax = plt.subplots(1,3)
fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].imshow(image_data, cmap=plt.cm.gray)
ax[1].imshow(another_image_data, cmap=plt.cm.gray);

real_image = image_data - another_image_data                 # Subtract the images pixel by pixel

ax[2].imshow(real_image, cmap=plt.cm.gray);

FITS Tables - An astronomical example

Stellar spectra data from the ESO Library of Stellar Spectra



In [ ]:

    
my_spectra_file = './MyData/Star_G5.fits'

spectra_data = fits.getdata(my_spectra_file)
spectra_header = fits.getheader(my_spectra_file)



In [ ]:

    
spectra_header



In [ ]:

    
# The FITS header has the information to make an array of wavelengths

Start = spectra_header['CRVAL1']
Number = spectra_header['NAXIS1']
Delta  = spectra_header['CDELT1']

End = Start + (Number * Delta)

Wavelength = np.arange(Start,End,Delta)



In [ ]:

    
fig, ax = plt.subplots(2,1)
fig.set_size_inches(11,8.5)

fig.tight_layout()

# Full spectra

ax[0].plot(Wavelength, spectra_data, color='b')
ax[0].set_ylabel("Flux")
ax[0].set_xlabel("Wavelength [angstroms]")

# Just the visible range with the hydrogen Balmer lines

ax[1].set_xlim(4000,7000)
ax[1].set_ylim(0.6,1.2)
ax[1].plot(Wavelength, spectra_data, color='b')
ax[1].set_ylabel("Flux")
ax[1].set_xlabel("Wavelength [angstroms]")

H_Balmer = [6563,4861,4341,4102,3970,3889,3835,3646]

ax[1].vlines(H_Balmer,0,2, color='r', linewidth=3, alpha = 0.25)

Stellar spectral classes



In [ ]:

    
import glob

star_list = glob.glob('./MyData/Star_*.fits')

star_list



In [ ]:

    
fig, ax = plt.subplots(1,1)
fig.set_size_inches(9,5)

fig.tight_layout()

# Full spectra

ax.set_ylabel("Flux")
ax.set_xlabel("Wavelength [angstroms]")
ax.set_ylim(0.0, 1.05)

for file in star_list:
        
    spectra = fits.getdata(file)
    
    spectra_normalized = spectra / spectra.max()
    
    ax.plot(Wavelength, spectra_normalized, label=file)

ax.legend(loc=0,shadow=True);

FITS Images - An astronomical example

World Coordinate System wcs



In [ ]:

    
from astropy.wcs import WCS

from matplotlib.colors import LogNorm



In [ ]:

    
my_image_file = "./MyData/m51.fits"

image_data = fits.getdata(my_image_file)
image_header = fits.getheader(my_image_file)



In [ ]:

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

fig.set_size_inches(12,4)

fig.tight_layout()

ax[0].set_title("Upside down!")
ax[1].set_title("Right side up!")
ax[2].set_title("LOG image intensity")

ax[0].imshow(image_data, cmap=plt.cm.gray)
ax[1].imshow(image_data, origin='lower', cmap=plt.cm.gray)
ax[2].imshow(image_data, origin='lower', cmap=plt.cm.gray, norm=LogNorm());



In [ ]:

    
image_header



In [ ]:

    
my_wcs = WCS(image_header)



In [ ]:

    
my_wcs



In [ ]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection=my_wcs)

fig.set_size_inches(6,6)

fig.tight_layout()

ax.grid(color='r', ls='--')
ax.set_xlabel('Right Ascension (J2000)')
ax.set_ylabel('Declination (J2000)')

ax.imshow(image_data, origin='lower', cmap=plt.cm.gray);



In [ ]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection=my_wcs)

fig.set_size_inches(6,6)

fig.tight_layout()

ax.set_xlabel('Right Ascension (J2000)')
ax.set_ylabel('Declination (J2000)')

ax.grid(color='r', ls='--')

overlay = ax.get_coords_overlay('galactic')

overlay.grid(color='y', ls='dotted')

overlay[0].set_axislabel('Galactic Longitude')
overlay[1].set_axislabel('Galactic Latitude')

ax.imshow(image_data, origin='lower', cmap=plt.cm.gray);

Pseudocolor - All color astronomy images are fake.

Color images are composed of three 2-d images:

JPG images are 3-d, even grayscale images



In [ ]:

    
redfilter = plt.imread("./MyData/sphereR.jpg")

redfilter.shape,redfilter.dtype

We just want to read in one of the three channels



In [ ]:

    
redfilter = plt.imread("./MyData/sphereR.jpg")[:,:,0]

redfilter.shape,redfilter.dtype



In [ ]:

    
plt.imshow(redfilter,cmap=plt.cm.gray);



In [ ]:

    
greenfilter = plt.imread("./MyData/sphereG.jpg")[:,:,0]
bluefilter = plt.imread("./MyData/sphereB.jpg")[:,:,0]



In [ ]:

    
fig, ax = plt.subplots(1,3)
fig.set_size_inches(12,3)

fig.tight_layout()

ax[0].set_title("Red Filter")
ax[1].set_title("Green Filter")
ax[2].set_title("Blue Filter")

ax[0].imshow(redfilter,cmap=plt.cm.gray)
ax[1].imshow(greenfilter,cmap=plt.cm.gray)
ax[2].imshow(bluefilter,cmap=plt.cm.gray);

Need to create a blank 3-d array to hold all of the images



In [ ]:

    
rgb = np.zeros((480,640,3),dtype='uint8')

print(rgb.shape, rgb.dtype)

plt.imshow(rgb,cmap=plt.cm.gray);

Fill the array with the filtered images



In [ ]:

    
rgb[:,:,0] = redfilter
rgb[:,:,1] = greenfilter
rgb[:,:,2] = bluefilter



In [ ]:

    
fig, ax = plt.subplots(1,4)
fig.set_size_inches(14,3)

fig.tight_layout()

ax[0].set_title("Red Filter")
ax[1].set_title("Green Filter")
ax[2].set_title("Blue Filter")
ax[3].set_title("All Filters Stacked")

ax[0].imshow(redfilter,cmap=plt.cm.gray)
ax[1].imshow(greenfilter,cmap=plt.cm.gray)
ax[2].imshow(bluefilter,cmap=plt.cm.gray)
ax[3].imshow(rgb,cmap=plt.cm.gray);



In [ ]:

    
print("The image rgb has a shape [height,width] of {0}".format(rgb.shape))
print("The image rgb is made up of data of type {0}".format(rgb.dtype))
print("The image rgb has a maximum value of {0}".format(rgb.max()))
print("The image rgb has a minimum value of {0}".format(rgb.min()))



In [ ]:

    
rgb[:,:,0] = redfilter * 1.5

plt.imshow(rgb)



In [ ]: