import numpy as np
import matplotlib.pyplot as plt

from skimage import data, io, segmentation, color
from skimage.future import graph

%matplotlib inline

import requests
from PIL import Image
from io import StringIO

"""
Explanation: Objective
In this notebook I am testing a reduced workflow that will:
1) input a geophysical map as an image 
2) reduce the number of colours 
3) convert the image with reduced colours from RGB to HSL
4) extract unique values of H, then sort those and also the L values in the correponding pixels, by the same index
5) visualize the result as an H vs. L curve to QC of the method
6) run the perceptual test by checking for monotonicity
Preliminaries - import main libraries
End of explanation
"""


url = 'https://mycarta.files.wordpress.com/2015/04/jet_tight.png'
r = requests.get(url)
img = np.asarray(Image.open(StringIO(r.content)).convert('RGB'))

img = np.asarray(Image.open('data/cbar/test.png'))[...,:3]

# plot image
fig = plt.figure(figsize=(6,6))
ax1 = fig.add_subplot(111)
plt.imshow(img)
ax1.xaxis.set_ticks([])
ax1.yaxis.set_ticks([])
plt.show()

"""
Explanation: Import test image. The colormap is Matlab's Jet
End of explanation
"""


# parameters chosen by trial and error. Will have to find a way to automatically optimize them
labels1 = segmentation.slic(img, compactness=30, n_segments=32) 
out1 = color.label2rgb(labels1, img, kind='avg')

fig = plt.figure(figsize=(6,6))
ax1 = fig.add_subplot(111)
plt.imshow(out1)
ax1.xaxis.set_ticks([])
ax1.yaxis.set_ticks([])
plt.show()


"""
Explanation: Reduce number of colours
I use here Scikit-learn segmentation using k-means clustering in Color-(x,y,z) space:
http://scikit-image.org/docs/dev/api/skimage.segmentation.html#skimage.segmentation.slic
End of explanation
"""


width, height, dump = np.shape(out1)
print(width, height, dump)

# method 1
# extract lightness and hue, combine them into a 2D array

# extract
from skimage.color import rgb2lab, lab2lch, lch2lab, lab2rgb
lab = rgb2lab(out1)
lch = lab2lch(lab)
        
lab = np.asarray(lab)
lch = np.asarray(lch)

# reshape
pixels_lab = np.reshape(lab, (width*height, -1))
l1, a, b = np.split(pixels_lab, 3, axis=-1)
        
pixels_lch = np.reshape(lch, (width*height, -1))
l2, c, h = np.split(pixels_lch, 3, axis=-1)
 
# flatten    
import itertools
lM = list(itertools.chain.from_iterable(l2))
hM = list(itertools.chain.from_iterable(h))

# zip together to make 2D numpy array
lhM = np.asarray(list(zip(hM,lM)))

lhM

# Sorting unique rows
# Joe Kington's answer on Stackoverflow: http://stackoverflow.com/a/16971224
def unique_rows(data):
    print(data.shape)
    uniq = np.unique(data.view(data.dtype.descr * data.shape[1]))
    return uniq.view(data.dtype).reshape(-1, data.shape[1])

uniqLM = unique_rows(lhM)
uniqLM.shape

# method 2
# sorting both lightness and hue by hue separately

from skimage.color import rgb2lab, lab2lch, lch2lab, lab2rgb
lab = rgb2lab(out1)
lch = lab2lch(lab)
        
lab = np.asarray(lab)
lch = np.asarray(lch)

pixels_lab = np.reshape(lab, (width*height, -1))
l1, a, b = np.split(pixels_lab, 3, axis=-1)
        
pixels_lch = np.reshape(lch, (width*height, -1))
l2, c, h = np.split(pixels_lch, 3, axis=-1)
                
huniq, unidx = np.unique(h, return_index=True)
luniq = l2[unidx]
cuniq = c[unidx]

# flatten luniq, cuniq
import itertools 
luniqM = list(itertools.chain.from_iterable(luniq))
cuniqM = list(itertools.chain.from_iterable(cuniq))

# compare output of two methods
lhM2 = np.asarray(list(zip(huniq,luniqM)))
print('method 2')
print('         H            L')
print(lhM2[:4])
print(lhM2[-4:])

print('method 1')
print('         H            L')
print(uniqLM[:4])
print(uniqLM[-4:])

"""
Explanation: Convert from RGB to HSL, get  unique values of H, S, and L then sort both lightness L and hue H, by increasing values of H
I tried two methods that return the same result:
1) get H, L, combine into 2D array, sort unique rows
2) sort H, with index returned then sort L with index
End of explanation
"""


from matplotlib.collections import LineCollection
from matplotlib.colors import ListedColormap, BoundaryNorm


# Data manipulation:

def make_segments(x, y):
    '''
    Create list of line segments from x and y coordinates, in the correct format for LineCollection:
    an array of the form   numlines x (points per line) x 2 (x and y) array
    '''

    points = np.array([x, y]).T.reshape(-1, 1, 2)
    segments = np.concatenate([points[:-1], points[1:]], axis=1)
    
    return segments


# Interface to LineCollection:

def colorline(x, y, z=None, cmap=plt.get_cmap('copper'), norm=plt.Normalize(0.0, 1.0), linewidth=3, alpha=1.0):
    '''
    Plot a colored line with coordinates x and y
    Optionally specify colors in the array z
    Optionally specify a colormap, a norm function and a line width
    '''
    
    # Default colors equally spaced on [0,1]:
    if z is None:
        z = np.linspace(0.0, 1.0, len(x))
           
    # Special case if a single number:
    if not hasattr(z, "__iter__"):  # to check for numerical input -- this is a hack
        z = np.array([z])
        
    z = np.asarray(z)
    
    segments = make_segments(x, y)
    lc = LineCollection(segments, array=z, cmap=cmap, norm=norm, linewidth=linewidth, alpha=alpha)
    
    ax = plt.gca()
    ax.add_collection(lc)
    
    return lc
        
    
def clear_frame(ax=None): 
    # Taken from a post by Tony S Yu
    if ax is None: 
        ax = plt.gca() 
    ax.xaxis.set_visible(False) 
    ax.yaxis.set_visible(False) 
    for spine in ax.spines.itervalues(): 
        spine.set_visible(False) 

"""
Explanation: Import a function to plot colored lines in the final plot using the colormap created above
From David Sanders:
https://github.com/dpsanders/matplotlib-examples/blob/master/colorline.ipynb
End of explanation
"""


# To color by L, it has to be normalized to [0 1]
luniqM_n=(luniqM-min(luniqM))/(max(luniqM)-min(luniqM))

fig = plt.figure(figsize=(16,4))
plt.xticks(np.arange(0, 2.25*np.pi,0.25*np.pi),[0.,   45.,   90.,  135.,  180.,  225.,  270.,  315., 360.]) 
# Hue as 0-360 angle

ax1 = fig.add_subplot(111)

# ax1.scatter(huniq, luniq)
colorline(huniq,luniqM, luniqM_n, linewidth=4,cmap='gray')
ax1.set_xlim(0, 2.25*np.pi)
ax1.set_ylim(-5, 105)
ax1.text(5, 95, 'H vs. L - colored by L', va='top')
plt.show()

# as of Nov 05, 2015 after reinstalling Anaconda, using Colorline gives a future warning:
# https://github.com/dpsanders/matplotlib-examples/issues/1
# rerun above to suppress warning

"""
Explanation: Make final plot of the sorted hue, H versus lightness, L, colored by L
End of explanation
"""


# Stackoverflow answer http://stackoverflow.com/a/4985520

def pairwise(seq):
    items = iter(seq)
    last = next(items)
    for item in items:
        yield last, item
        last = item

def strictly_increasing(L):
    return all(x<y for x, y in pairwise(L))

def strictly_decreasing(L):
    return all(x>y for x, y in pairwise(L))

def non_increasing(L):
    return all(x>=y for x, y in pairwise(L))

def non_decreasing(L):
    return all(x<=y for x, y in pairwise(L))

print(strictly_increasing(luniq))
print(non_decreasing(luniq))

print(strictly_decreasing(luniq))
print(non_increasing(luniq))

"""
Explanation: Run perceptual test checks for monotonicity
End of explanation
"""


# Originally from: http://bgfons.com/upload/rainbow_texture1761.jpg. Resized and saved as png
url = 'https://mycarta.files.wordpress.com/2015/11/rainbow_texture17611.png'
r = requests.get(url)
img = np.asarray(Image.open(StringIO(r.content)).convert('RGB'))

# plot image
fig = plt.figure(figsize=(6,6))
ax1 = fig.add_subplot(111)
plt.imshow(img)
ax1.xaxis.set_ticks([])
ax1.yaxis.set_ticks([])
plt.show()

labels1 = segmentation.slic(img, compactness=30, n_segments=32)
out1 = color.label2rgb(labels1, img, kind='avg')

fig = plt.figure(figsize=(6,6))
ax1 = fig.add_subplot(111)
plt.imshow(out1)
ax1.xaxis.set_ticks([])
ax1.yaxis.set_ticks([])
plt.show()

width, height, dump = np.shape(out1)
# method 2
# sorting both lightness and hue by hue separately

from skimage.color import rgb2lab, lab2lch, lch2lab, lab2rgb
lab = rgb2lab(out1)
lch = lab2lch(lab)
        
lab = np.asarray(lab)
lch = np.asarray(lch)


pixels_lab = np.reshape(lab, (width*height, -1))
l1, a, b = np.split(pixels_lab, 3, axis=-1)
        
pixels_lch = np.reshape(lch, (width*height, -1))
l2, c, h = np.split(pixels_lch, 3, axis=-1)
                
huniq, unidx = np.unique(h, return_index=True)
luniq = l2[unidx]
cuniq = c[unidx]

# flatten luniq, cuniq
import itertools 
luniqM = list(itertools.chain.from_iterable(luniq))
cuniqM = list(itertools.chain.from_iterable(cuniq))

# To color by L, it has to be normalized to [0 1]
luniqM_n=(luniqM-min(luniqM))/(max(luniqM)-min(luniqM))

fig = plt.figure(figsize=(8,4))
plt.xticks(np.arange(0, 2.25*np.pi,0.25*np.pi),[0.,   45.,   90.,  135.,  180.,  225.,  270.,  315., 360.]) 
# Hue as 0-360 angle

ax1 = fig.add_subplot(111)

# ax1.scatter(huniq, luniq)
colorline(huniq,luniqM, luniqM_n, linewidth=4,cmap='gray')
ax1.set_xlim(0, 2.25*np.pi)
ax1.set_ylim(-5, 105)
ax1.text(5, 95, 'H vs. L - colored by L', va='top')
plt.show()

print strictly_increasing(luniq)
print non_decreasing(luniq)

print strictly_decreasing(luniq)
print non_increasing(luniq)

"""
Explanation: Now we try it on an abstract rainbow image
End of explanation
"""


url = 'https://mycarta.files.wordpress.com/2015/04/cubic_no_red_tight.png'
r = requests.get(url)
img = np.asarray(Image.open(StringIO(r.content)).convert('RGB'))

# plot image
fig = plt.figure(figsize=(6,6))
ax1 = fig.add_subplot(111)
plt.imshow(img)
ax1.xaxis.set_ticks([])
ax1.yaxis.set_ticks([])
plt.show()

labels1 = segmentation.slic(img, compactness=30, n_segments=32)
out1 = color.label2rgb(labels1, img, kind='avg')

fig = plt.figure(figsize=(6,6))
ax1 = fig.add_subplot(111)
plt.imshow(out1)
ax1.xaxis.set_ticks([])
ax1.yaxis.set_ticks([])
plt.show()


width, height, dump = np.shape(out1)
# method 2
# sorting both lightness and hue by hue separately

from skimage.color import rgb2lab, lab2lch, lch2lab, lab2rgb
lab = rgb2lab(out1)
lch = lab2lch(lab)
        
lab = np.asarray(lab)
lch = np.asarray(lch)


pixels_lab = np.reshape(lab, (width*height, -1))
l1, a, b = np.split(pixels_lab, 3, axis=-1)
        
pixels_lch = np.reshape(lch, (width*height, -1))
l2, c, h = np.split(pixels_lch, 3, axis=-1)
                
huniq, unidx = np.unique(h, return_index=True)
luniq = l2[unidx]
cuniq = c[unidx]

# flatten luniq, cuniq
import itertools 
luniqM = list(itertools.chain.from_iterable(luniq))
cuniqM = list(itertools.chain.from_iterable(cuniq))

# To color by L, it has to be normalized to [0 1]
luniqM_n=(luniqM-min(luniqM))/(max(luniqM)-min(luniqM))

fig = plt.figure(figsize=(8,4))
plt.xticks(np.arange(0, 2.25*np.pi,0.25*np.pi),[0.,   45.,   90.,  135.,  180.,  225.,  270.,  315., 360.]) 
# Hue as 0-360 angle

ax1 = fig.add_subplot(111)

# ax1.scatter(huniq, luniq)
colorline(huniq,luniqM, luniqM_n, linewidth=4,cmap='gray')
ax1.set_xlim(0, 2.25*np.pi)
ax1.set_ylim(-5, 105)
ax1.text(5, 95, 'H vs. L - colored by L', va='top')
plt.show()

print strictly_increasing(luniq)
print strictly_decreasing(luniq)

print non_increasing(luniq)
print non_decreasing(luniq)

"""
Explanation: Try it on mycarta perceptual rainbow
End of explanation
"""


print luniqM[:15]
plt.plot(luniqM[:15])

"""
Explanation: The test should have worked but it did not.
We need to include some smoothing, or despiking to deal with small non monotonic value pairs. See below tests
End of explanation
"""


def moving_average(a, length, mode='valid'):
    
    #pad = np.floor(length/2)
    pad = int(np.floor(length/2))   # replace to prevent a deprecation warning 
                                    # due to passing a float as an index
    
    if mode == 'full': 
        pad *= 2
    
    # Make a padded version, paddding with first and last values
    r = np.empty(a.shape[0] + 2*pad)
    r[:pad] = a[0]
    r[pad:-pad] = a
    r[-pad:] = a[-1]
    
    # Cumsum with shifting trick
    s = np.cumsum(r, dtype=float)
    s[length:] = s[length:] - s[:-length]
    out = s[length-1:]/length
    
    # Decide what to return
    if mode == 'same':
        if out.shape[0] != a.shape[0]: 
            # If size doesn't match, then interpolate.
            out = (out[:-1,...] + out[1:,...]) / 2
        return out
    elif mode == 'valid':
        return out[pad:-pad]
    else: # mode=='full' and we used a double pad
        return out

avg = moving_average(np.asarray(luniqM), 7, mode='same')
plt.plot(avg[:15])


print strictly_increasing(avg)
print non_decreasing(avg)

print strictly_decreasing(avg)
print non_increasing(avg)

"""
Explanation: From Matt Hall
https://github.com/kwinkunks/notebooks/blob/master/Backus.ipynb
End of explanation
"""


