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>