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

m32 = scipy.ndimage.imread("m32.jpg")[...,0]

m32.shape

plt.imshow(m32, interpolation="none", cmap="gray")

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

m32

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)

m32 = m31[5200:6000,3700:4500,0]

plt.imshow(m32, interpolation="none", cmap="gray")

plt.plot(m32[400,:])

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

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)

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)