import numpy as np
import matplotlib.pylab as plt
from scipy import ndimage
import pandas as pd

%matplotlib inline

data = np.genfromtxt('data/signal.csv',  delimiter=',')

X1 = data[0:13]
X2 = data[13:26]
X3 = data[26:39]
X4 = data[39:52]

plt.figure(figsize=(8, 3))
plt.imshow(np.hstack([X1, X2, X3, X4]), interpolation='none')
plt.show()

th = 0.2
plt.figure(figsize=(8, 3))
plt.imshow(np.hstack([X1, X2, X3, X4]) >= th, interpolation='none')
plt.show()

img = X4

df = pd.DataFrame(img)
stacked = pd.DataFrame(df.stack())
stacked.reset_index(inplace=True)
stacked.columns = ('y', 'x', 'v')
filtered = stacked[stacked.v >= th]
X = filtered[['x', 'y']].values
X_cent = X - X.mean(axis=0)

plt.scatter(X_cent[:, 0], X_cent[:, 1], c='white')
plt.show()

_, _, V = np.linalg.svd(X_cent)
v = V[0, :]

plt.scatter(X_cent[:, 0], X_cent[:, 1], c='white')
plt.quiver(v[0], v[1], angles='xy', scale_units='xy', scale=0.2)
plt.show()

plt.subplot(1, 2, 1)

plt.title('original')
plt.imshow(img)
plt.quiver([6], [6], v[0], v[1],
           angles='xy',scale_units='xy',scale=0.4)

plt.subplot(1, 2, 2)
deg = 180 * np.arccos(v[0]) / np.pi
X_rot = ndimage.rotate(img, deg, reshape=False)

plt.title('rotated by %d C' % deg)
plt.imshow(X_rot)
plt.show()

plt.plot(X_rot.sum(axis=0))
plt.show()

plt.plot(X1.sum(axis=0))
plt.show()

base = X1.sum(axis=0)

plt.plot(X_rot.sum(axis=0) - base)
plt.show()

s = X_rot.sum(axis=0) - base
f = np.fft.fft(s)
fabs = np.abs(np.fft.fftshift(f))
fabs = fabs / fabs.max()

plt.plot(fabs[7:])
plt.show()

for idx, img in enumerate([X1, X2, X3, X4]):
    print 'data set #', idx
    
    df = pd.DataFrame(img)
    stacked = pd.DataFrame(df.stack())
    stacked.reset_index(inplace=True)
    stacked.columns = ('y', 'x', 'v')
    filtered = stacked[stacked.v >= th]
    X = filtered[['x', 'y']].values
    X_cent = X - X.mean(axis=0)
    _, _, V = np.linalg.svd(X_cent)
    v = V[0, :]
    
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.quiver([6], [6], v[0], v[1],
               angles='xy',scale_units='xy',scale=0.4)

    plt.subplot(1, 2, 2)
    deg = 180 * np.arccos(v[0]) / np.pi

    if abs(deg) > 1e-6:
        X_rot = ndimage.rotate(img, deg, reshape=False)
    else:
        X_rot = img

    plt.title('rotated by %d C' % deg)
    plt.imshow(X_rot)
    plt.show()

    
    plt.subplot(2, 1, 1)
    s = X_rot.sum(axis=0) - base
    plt.title('1d signal, time')
    plt.plot(s)
    plt.xticks([])
    
    plt.subplot(2, 1, 2)
    f = np.fft.fft(s)
    fabs = np.abs(np.fft.fftshift(f))
    fabs = fabs / fabs.max()
    plt.title('1d signal, frequency')
    plt.plot(fabs)
    plt.xticks([])
    plt.show()

def extract_features(img, th, base):
    df = pd.DataFrame(img)
    stacked = pd.DataFrame(df.stack())
    stacked.reset_index(inplace=True)
    stacked.columns = ('y', 'x', 'v')
    filtered = stacked[stacked.v >= th]
    X = filtered[['x', 'y']].values
    X_cent = X - X.mean(axis=0)
    _, _, V = np.linalg.svd(X_cent)
    v = V[0, :]

    deg = 180 * np.arccos(v[0]) / np.pi

    if abs(deg) > 1e-6:
        X_rot = ndimage.rotate(img, deg, reshape=False)
    else:
        X_rot = img

    s = X_rot.sum(axis=0) - base
    f = np.fft.fft(s)
    fabs = np.abs(np.fft.fftshift(f))
    fabs = fabs / fabs.max()
    return fabs[[6, 8]]

f2 = extract_features(X2, th, base)
f3 = extract_features(X3, th, base)
f4 = extract_features(X4, th, base)

F = np.array([f2, f3, f4])
print F

plt.scatter(F[:, 0], F[:, 1], c=['white', 'black', 'black'])
plt.show()