1. timefreak
  2. images.ipynb

Images

Play with some images. I want to make a point about the meaning of 'resolution'.

  • In digital photography, resolution usually means sample interval (or, equivalently, sample rate). In reflection seismic, this is usually 250 to 1000 Hz, giving a half-Nyquist of 125 Hz to 500 Hz, which is usually plenty for our signals.
  • In signal analysis, we talk a lot about localization, which is the smearedness of things. This is akin to 'focus' in photography, but that's not quite it, because something can be fuzzy in the physical world not just because we took a poor photo. On other words, something canbe inherently unlocalized (see the M32 galaxy example, below).
  • In audio signals, 'resolution' usually means bit depth.

Let's say that these properties are the 3 horsemen of resolution.

Prelims

In [2]:
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

Image analogies

In [3]:
import scipy.ndimage
In [4]:
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()

Unlocalized object

Start by reading the red channel of an astronomical photograph from a file.

In [5]:
m32 = scipy.ndimage.imread("m32.jpg")[...,0]
In [7]:
m32.shape
Out[7]:
(866, 866)
In [8]:
plt.imshow(m32, interpolation="none", cmap="gray")
Out[8]:
<matplotlib.image.AxesImage at 0x111f19050>
In [9]:
m32_med = scipy.ndimage.median_filter(m32, 5)
plt.imshow(m32_med, interpolation="none", cmap="Greys")
plt.show()
In [10]:
m32_gss = scipy.ndimage.gaussian_filter(m32, 5)
plt.imshow(m32_gss, interpolation="none", cmap="Greys")
plt.show()
In [11]:
m32.shape
Out[11]:
(866, 866)
In [12]:
m32
Out[12]:
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

In [13]:
import requests
from StringIO import StringIO
from PIL import Image
In [14]:
a = requests.get("http://svs.gsfc.nasa.gov/vis/a010000/a010400/a010485/M31_Wide.jpg")
In [15]:
m31 = np.array(Image.open(StringIO(a.content))) # ndimage can't work with this object
In [16]:
plt.imshow(m31)
Out[16]:
<matplotlib.image.AxesImage at 0x11629afd0>
In [17]:
m32 = m31[5200:6000,3700:4500,0]
In [18]:
plt.imshow(m32, interpolation="none", cmap="gray")
Out[18]:
<matplotlib.image.AxesImage at 0x1163ec390>
In [19]:
plt.plot(m32[400,:])
Out[19]:
[<matplotlib.lines.Line2D at 0x116426810>]

Quantize image values

In [6]:
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)
In [7]:
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()

Add noise

We can attack it in the signal:noise space:

In [8]:
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()