Difference between revisions of "Image processing with Python and SciPy"

From AstroEd
Jump to navigation Jump to search
Line 49: Line 49:
  
 
  import numpy as np
 
  import numpy as np
 +
import matplotlib.pyplot as plt
 
  from scipy.misc import imread, imsave
 
  from scipy.misc import imread, imsave
  
  im1 = imread('test.jpg').astype(np.float32)
+
  image_data = imread('test.jpg').astype(np.float32)
  print 'Size: ', im1.size
+
  print 'Size: ', image_data.size
  print 'Shape: ', im1.shape
+
  print 'Shape: ', image_data.shape
 +
 +
plt.imshow(image_data)
 +
plt.show()
 +
 
 
  exit()
 
  exit()
  
Line 63: Line 68:
 
because the image is 512x512 pixels and has 3 planes -- red, green, and blue. When SciPy reads a jpg or png image it will separate the colors for you. The "image" is a data cube.
 
because the image is 512x512 pixels and has 3 planes -- red, green, and blue. When SciPy reads a jpg or png image it will separate the colors for you. The "image" is a data cube.
  
With that, we can slice the cube and see various planes out of the data.  Here's a nice example from the [http://stackoverflow.com/questions/14138899/how-can-an-almost-arbitrary-plane-in-a-3d-dataset-be-plotted-by-matplotlib web]
+
We can display images with matplotlib.pyplot using [http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.imshow imshow()] --
 +
 
 +
imshow(X, cmap=None, norm=None, aspect=None, interpolation=None,
 +
  alpha=None, vmin=None, vmax=None, origin=None, extent=None, **kwargs)
 +
 
 +
where "X" is the image array.  If X is 3-dimensional, imshow will display a color image though without limiting the ranges for red, green, and blue the image will look unrealistic.
 +
 
 +
 
 +
With that, we can slice the cube and see various planes out of the data.  Here's a nice example adapted from the [http://stackoverflow.com/questions/14138899/how-can-an-almost-arbitrary-plane-in-a-3d-dataset-be-plotted-by-matplotlib web].  Import numpy and pylab so that we can view the images, then create a small sample data set in a 64x64x3
  
 
  import numpy as np
 
  import numpy as np
  import pylab as plt
+
  import matplotlib.pyplot as plt
  
  data = np.arange((64**3))
+
  data = np.arange((64*64*3), dtype=np.float32)
  data.resize((64,64,64))
+
  data.resize((64,64,3))
  
 
  def get_slice(volume, orientation, index):
 
  def get_slice(volume, orientation, index):
Line 80: Line 93:
  
 
  plt.subplot(221)
 
  plt.subplot(221)
  plt.imshow(get_slice(data, "x", 10), vmin=0, vmax=64**3)
+
  plt.imshow(get_slice(data, "x", 10), vmin=0, vmax=64*64*3)
  
 
  plt.subplot(222)
 
  plt.subplot(222)
  plt.imshow(get_slice(data, "x", 39), vmin=0, vmax=64**3)
+
  plt.imshow(get_slice(data, "y", 10), vmin=0, vmax=64*64*3)
  
plt.subplot(223)
 
plt.imshow(get_slice(data, "y", 15), vmin=0, vmax=64**3)
 
 
  plt.subplot(224)
 
  plt.subplot(224)
  plt.imshow(get_slice(data, "z", 25), vmin=0, vmax=64**3)   
+
  plt.imshow(get_slice(data, "z", 2), vmin=0, vmax=64*64*3)   
  
 
  plt.show()
 
  plt.show()

Revision as of 21:12, 26 February 2013

Given that NumPy provides multidimensional arrays, and that there is core support through the Python Imaging Library to display images and manipulate images in the Python environment, it's easy to take the next step and combine these for scientific image processing. As part of our short course on Python for Physics and Astronomy we begin by exploring how Python handle image input and output


Python Imaging Library - PIL

Before we get into the broad area, there is a caveat for users of Python 3. The essential Python Imaging Library (PIL) is not yet completely compatible with new version of Python. Consequently the FITS tools we will need for astronomical image processing are also currently only supported in the mature versions of Python 2. The comments that follow are based on Python 2.7.

PIL provides functions to manipulate images, including reading, modifying and saving in various standard image formats. Its functions are documented in an on-line manual with a tutorial, and in this handy pdf guide.

As a simple starting example, suppose you have an image that was taken with the camera turned so that "up" is to the side when the image is displayed. Here's how you would rotate an image 90 degrees.

import Image as pil
parser= argparse.ArgumentParser(description = 'Rotate a png image 90 degrees')
if len(sys.argv) == 1:
  sys.exit("Usage: png_image_rotate file.png ")
  exit()
elif len(sys.argv) == 2:
  infilename = sys.argv[1]
else:
  sys.exit("Usage: png_image_rotate file.png ")
  exit() 
myimage = pil.open(infilename)

mirror = myimage.transpose(pil.ROTATE_90)
outfilename = os.path.splitext(os.path.basename(infilename))[0]+'_r90.png'
mirror.save(outfilename)


The first part of this is standard form to get the image name on the command line and make it available to the program. The PIL is imported with Image, and appears in the code as "pil". This is an amazingly short program, because in opening the image the library handles all the conversions in formatting and stores the image internally so that you refer to it only by the name assigned when it is loaded. We operate on the image with the transpose function, which has an argument that controls what it does. Here we rotate the image 90 degrees, and then save it to a file with a new name. The saving operation converts the internal data back to the file format set by the extension used in the file name.

You can transpose an image left-right with

mirror = myimage.transpose(pil.FLIP_LEFT_RIGHT)

or do both in one step with

mirror = myimage.transpose(Image.FLIP_LEFT_RIGHT).transpose(pil.ROTATE_90)

Processing is not limited to "PNG" files, though that file type is preferred because it is not a lossy storage option. Python reads and writes "jpg" files too. While PIL provides some essential functionality, for more practical uses in astronomy we need to read Flexible Image Transport or "FITS" files, and to enable numerical work on images in SciPy.

SciPy image processing

SciPy can read jpg and png images directly, without using PIL.

import numpy as np

import matplotlib.pyplot as plt

from scipy.misc import imread, imsave
image_data = imread('test.jpg').astype(np.float32)
print 'Size: ', image_data.size
print 'Shape: ', image_data.shape

plt.imshow(image_data)
plt.show()
exit()

For our test image, this returns

Size:  786432
Shape:  (512, 512, 3)

because the image is 512x512 pixels and has 3 planes -- red, green, and blue. When SciPy reads a jpg or png image it will separate the colors for you. The "image" is a data cube.

We can display images with matplotlib.pyplot using imshow() --

imshow(X, cmap=None, norm=None, aspect=None, interpolation=None,
  alpha=None, vmin=None, vmax=None, origin=None, extent=None, **kwargs)

where "X" is the image array. If X is 3-dimensional, imshow will display a color image though without limiting the ranges for red, green, and blue the image will look unrealistic.


With that, we can slice the cube and see various planes out of the data. Here's a nice example adapted from the web. Import numpy and pylab so that we can view the images, then create a small sample data set in a 64x64x3

import numpy as np
import matplotlib.pyplot as plt
data = np.arange((64*64*3), dtype=np.float32)
data.resize((64,64,3))
def get_slice(volume, orientation, index):
   orientation2slicefunc = {
       "x" : lambda ar:ar[index,:,:], 
       "y" : lambda ar:ar[:,index,:],  
       "z" : lambda ar:ar[:,:,index]
   }
   return orientation2slicefunc[orientation](volume)
plt.subplot(221)
plt.imshow(get_slice(data, "x", 10), vmin=0, vmax=64*64*3)
plt.subplot(222)
plt.imshow(get_slice(data, "y", 10), vmin=0, vmax=64*64*3)
plt.subplot(224)
plt.imshow(get_slice(data, "z", 2), vmin=0, vmax=64*64*3)  
plt.show()

Astronomical FITS files with PyFITS

PyFITS is available from the Space Telescope Science Institute, and can be added easily to a Python installation that already and NumPy and SciPy. As of January 2013, the current version 3.1.1 of PyFITS supports all the functions needed to manage image and table data in the standard Flexibile Image Transport System (FITS) files of astronomy.


FITS headers

Processing astronomical images

SciKit image processing

Examples

For examples of Python illustrating image processing, see the examples section.

Assignments

For the assigned homework to use these ideas, see the assignments section.