Play with some images. I want to make a point about the meaning of 'resolution'.
Let's say that these properties are the 3 horsemen of resolution.
import numpy as np
#from numpy.fft import fft, fftfreq
from scipy.fftpack import fft, rfft, irfft, fftfreq, rfftfreq
import scipy.signal
import matplotlib.pyplot as plt
%matplotlib inline
import scipy.ndimage
im = np.zeros((20, 20))
im[5:-5, 5:-5] = 1
im = scipy.ndimage.distance_transform_bf(im)
im_noise = im + 0.2 * np.random.randn(*im.shape)
im_med = scipy.ndimage.median_filter(im_noise, 3)
plt.imshow(im_med, interpolation="none", cmap="Greys")
plt.show()
Start by reading the red channel of an astronomical photograph from a file.
m32 = scipy.ndimage.imread("m32.jpg")[...,0]
m32.shape
(866, 866)
plt.imshow(m32, interpolation="none", cmap="gray")
<matplotlib.image.AxesImage at 0x111f19050>
m32_med = scipy.ndimage.median_filter(m32, 5)
plt.imshow(m32_med, interpolation="none", cmap="Greys")
plt.show()
m32_gss = scipy.ndimage.gaussian_filter(m32, 5)
plt.imshow(m32_gss, interpolation="none", cmap="Greys")
plt.show()
m32.shape
(866, 866)
m32
array([[ 7, 6, 5, ..., 24, 25, 24], [ 9, 5, 2, ..., 31, 27, 28], [10, 8, 4, ..., 28, 21, 26], ..., [16, 15, 17, ..., 40, 49, 49], [17, 23, 33, ..., 41, 46, 46], [17, 24, 34, ..., 43, 46, 46]], dtype=uint8)
Let's look for another image
import requests
from StringIO import StringIO
from PIL import Image
a = requests.get("http://svs.gsfc.nasa.gov/vis/a010000/a010400/a010485/M31_Wide.jpg")
m31 = np.array(Image.open(StringIO(a.content))) # ndimage can't work with this object
plt.imshow(m31)
<matplotlib.image.AxesImage at 0x11629afd0>
m32 = m31[5200:6000,3700:4500,0]
plt.imshow(m32, interpolation="none", cmap="gray")
<matplotlib.image.AxesImage at 0x1163ec390>
plt.plot(m32[400,:])
[<matplotlib.lines.Line2D at 0x116426810>]
def index(a, colours=8):
b = a.reshape(a.size)
b = float(colours) * b / np.amax(b)
bins = np.linspace(0, colours-1, colours)
c = np.digitize(b, bins)
return c.reshape(a.shape)
quantized = index(m32, 8)
plt.figure(figsize=(12,5))
plt.subplot(121)
plt.imshow(quantized)
plt.hlines(400,0,880,'red')
plt.subplot(122)
plt.plot(quantized[400,:], 'red')
plt.show()
We can attack it in the signal:noise space:
noise = np.random.random_integers(0, 255, m32.shape) # The image is 8-bit int.
m32_nse = m32 + noise
plt.figure(figsize=(15,8))
plt.subplot(121)
plt.imshow(m32, cmap="gray", interpolation="none")
plt.subplot(122)
plt.imshow(m32_nse, cmap="gray", interpolation="none")
plt.show()
def impair(image, zoom=0.25, gauss=4, colours=4, cmap="gray", plots=False, lines=False, offset=0):
im_low = scipy.ndimage.zoom(image, zoom)
im_gss = scipy.ndimage.gaussian_filter(image, gauss)
im_bit = index(image, colours)
im_nse = image * (0.5 + np.random.random(size=image.shape))
im_crop = image[image.shape[0]/3.:-image.shape[0]/3.,image.shape[1]/3.:-image.shape[1]/3.]
layout = 160
if plots:
plt.figure(figsize=(18,3))
plt.subplot(layout + 1)
plt.plot(image[image.shape[0]/2+offset,:])
plt.subplot(layout + 2)
plt.plot(im_low[im_low.shape[0]/2+offset*zoom,:])
plt.subplot(layout + 3)
plt.plot(im_gss[im_gss.shape[0]/2+offset,:])
ax = plt.subplot(layout + 4)
ax.plot(im_bit[im_bit.shape[0]/2+offset,:])
ax.set_ylim(1, colours+1)
plt.subplot(layout + 5)
plt.plot(im_nse[im_nse.shape[0]/2+offset,:])
plt.subplot(layout + 6)
plt.plot(im_crop[im_crop.shape[0]/2+offset,:])
plt.show()
else:
plt.figure(figsize=(18,3))
plt.subplot(layout + 1)
plt.imshow(image, cmap=cmap)
if lines: plt.hlines(image.shape[0]/2+offset, 0, image.shape[1], 'b')
plt.subplot(layout + 2)
plt.imshow(im_low, interpolation="none", cmap=cmap)
if lines: plt.hlines(im_low.shape[0]/2+offset*zoom, 0, im_low.shape[1], 'b')
plt.subplot(layout + 3)
plt.imshow(im_gss, interpolation="none", cmap=cmap)
if lines: plt.hlines(im_gss.shape[0]/2+offset, 0, im_gss.shape[1], 'b')
plt.subplot(layout + 4)
plt.imshow(im_bit, interpolation="none", cmap=cmap)
if lines: plt.hlines(im_bit.shape[0]/2+offset, 0, im_bit.shape[1], 'b')
plt.subplot(layout + 5)
plt.imshow(im_nse, interpolation="none", cmap=cmap)
if lines: plt.hlines(im_nse.shape[0]/2+offset, 0, im_nse.shape[1], 'b')
plt.subplot(layout + 6)
plt.imshow(im_crop, interpolation="none", cmap=cmap)
if lines: plt.hlines(im_crop.shape[0]/2+offset, 0, im_crop.shape[1], 'b')
plt.show()
# Bricks image flickr/Grzesiek CC-BY
bricks = scipy.ndimage.imread("bricks.jpg")[:400,:400,0]
bricks_200 = scipy.ndimage.zoom(bricks, 0.5)
impair(bricks_200)
# From my thesis
shell = scipy.ndimage.imread("photomic.jpg")[:,:1738,0] # Make it square
shell_400 = scipy.ndimage.zoom(shell, 400/1738.)
shell_200 = scipy.ndimage.zoom(shell_400, 0.5)
impair(shell_200)
Some things are not well localised in space.
m32_400 = scipy.ndimage.zoom(m32, 0.4618)
m32_200 = scipy.ndimage.zoom(m32, 0.2309)
impair(m32_200, gauss=2, lines=True, offset=4)
impair(m32_200, gauss=2, plots=True, offset=4)
# Image Flickr/FotoSleuth CC-BY
clouds = scipy.ndimage.imread("clouds.jpg")[:,:800,0] # Make it square
clouds_200 = scipy.ndimage.zoom(clouds, 0.25)
impair(clouds_200)