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()

Impair images in various ways



In [107]:

    
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()



In [108]:

    
# Bricks image flickr/Grzesiek CC-BY
bricks = scipy.ndimage.imread("bricks.jpg")[:400,:400,0]
bricks_200 = scipy.ndimage.zoom(bricks, 0.5)



In [109]:

    
impair(bricks_200)



In [110]:

    
# From my thesis
shell = scipy.ndimage.imread("photomic.jpg")[:,:1738,0] # Make it square



In [111]:

    
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.



In [116]:

    
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)



In [113]:

    
impair(m32_200, gauss=2, plots=True, offset=4)



In [114]:

    
# Image Flickr/FotoSleuth CC-BY
clouds = scipy.ndimage.imread("clouds.jpg")[:,:800,0] # Make it square



In [115]:

    
clouds_200 = scipy.ndimage.zoom(clouds, 0.25)
impair(clouds_200)



In [106]:



In [ ]: