import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from skimage.color import convert_colorspace
from skimage.io import imread
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1 import ImageGrid
%matplotlib inline

bluesky = imread("http://i.stack.imgur.com/fbIyO.jpg")
greysky = imread("http://i.stack.imgur.com/hBmW9.jpg")

bluesky_xyz = convert_colorspace(bluesky,"RGB","XYZ")
greysky_xyz = convert_colorspace(greysky,"RGB","XYZ")

print bluesky_xyz.min(), bluesky_xyz.max()
print greysky_xyz.min(), greysky_xyz.max()

#better colorscheme
plt.rc("image", cmap="cubehelix")

fig,axs = plt.subplots(1,2, subplot_kw={"xticks":[], "yticks":[]})
axs[0].imshow(bluesky_xyz[:,:,1],vmin=0,vmax=1);
axs[0].set_title("Blue");
axs[1].imshow(greysky_xyz[:,:,1],vmin=0,vmax=1);
axs[1].set_title("Gray");

fig,axs = plt.subplots(1,2, subplot_kw={"xticks":[], "yticks":[]})
axs[0].imshow(bluesky_xyz[:,:,2],vmin=0,vmax=1);
axs[0].set_title("Blue");
axs[1].imshow(greysky_xyz[:,:,2],vmin=0,vmax=1);
axs[1].set_title("Gray");

fig = plt.figure(1, (27,8))
grid = ImageGrid(fig, 111, nrows_ncols=(3,2), axes_pad=0.1, 
                 share_all=True, label_mode="L", cbar_location="bottom", cbar_mode="single" )

grid[0].imshow(bluesky);
grid[1].imshow(greysky);
grid[2].imshow(bluesky_xyz[:,:,1],vmin=0,vmax=1);
grid[3].imshow(greysky_xyz[:,:,1],vmin=0,vmax=1);
grid[4].imshow(bluesky_xyz[:,:,2],vmin=0,vmax=1);
im = grid[5].imshow(greysky_xyz[:,:,2],vmin=0,vmax=1);
grid.cbar_axes[0].colorbar(im)
grid[0].set_title("Blue")
grid[1].set_title("Gray")
grid[0].set_ylabel("original")
grid[2].set_ylabel("Y = luminance")
grid[4].set_ylabel("Z ~ blue")
grid.axes_llc.set_xticks([]);
grid.axes_llc.set_yticks([]);

fig,axs = plt.subplots(2, figsize=(10,3), sharex=True, sharey=True);
bins = np.linspace(0,1,200)
x = bins[1:] - 0.5*(bins[1]-bins[0])
axs[0].plot(x,np.histogram(bluesky_xyz[:,:,1].ravel(),bins,normed=True)[0],c='gray',lw=2, alpha=0.9, label="Y");
axs[0].plot(x,np.histogram(bluesky_xyz[:,:,2].ravel(),bins,normed=True)[0],c='blue',lw=2, alpha=0.7, label="Z");
axs[1].plot(x,np.histogram(greysky_xyz[:,:,1].ravel(),bins,normed=True)[0],c='gray',lw=2, alpha=0.9);
axs[1].plot(x,np.histogram(greysky_xyz[:,:,2].ravel(),bins,normed=True)[0],c='blue',lw=2, alpha=0.7);
plt.ylim((0,5));
plt.xlim((0,1.0));
axs[0].legend();
plt.yticks([]);
axs[0].vlines([0.3,0.6],0,5,linestyle="--",alpha=0.3);
axs[1].vlines([0.35,0.6],0,5,linestyle="--",alpha=0.3);
axs[0].set_ylabel("Blue Sky");
axs[1].set_ylabel("Gray Sky");

fig,axs = plt.subplots(1, 2, figsize=(10,5), subplot_kw={"xticks":[], "yticks":[]});
skyfilt_blue = (bluesky_xyz[:,:,2] > 0.3)*(bluesky_xyz[:,:,2] < 0.6)
skyfilt_grey = (greysky_xyz[:,:,2] > 0.35)*(greysky_xyz[:,:,2] < 0.6)
axs[0].imshow( skyfilt_blue, cmap=plt.bone());
axs[1].imshow( skyfilt_grey, cmap=plt.bone());

fig,axs = plt.subplots(1, 2, figsize=(10,5), subplot_kw={"xticks":[], "yticks":[]});
blue_emph = bluesky.copy()
blue_emph[~skyfilt_blue] *= 0
grey_emph = greysky.copy()
grey_emph[~skyfilt_grey] *= 0
axs[0].imshow( blue_emph, cmap=plt.bone());
axs[0].set_title("Blue Sky");
axs[1].set_title("Gray Sky");
axs[1].imshow( grey_emph, cmap=plt.bone());

fig,axs = plt.subplots(2, figsize=(10,3), sharex=True, sharey=True);
bins = np.linspace(0,1,200)
x = bins[1:] - 0.5*(bins[1]-bins[0])
axs[0].plot(x,np.histogram(bluesky_xyz[skyfilt_blue][:,1],bins,normed=True)[0],c='gray',lw=2, alpha=0.9, label="Y");
axs[0].plot(x,np.histogram(bluesky_xyz[skyfilt_blue][:,2],bins,normed=True)[0],c='blue',lw=2, alpha=0.7, label="Z");
axs[1].plot(x,np.histogram(greysky_xyz[skyfilt_grey][:,1],bins,normed=True)[0],c='gray',lw=2, alpha=0.9);
axs[1].plot(x,np.histogram(greysky_xyz[skyfilt_grey][:,2],bins,normed=True)[0],c='blue',lw=2, alpha=0.7);
plt.ylim((0,15));
plt.xlim((0.2,0.6));
plt.yticks([]);
axs[0].legend();
axs[0].vlines([0.3,0.6],0,15,linestyle="--",alpha=0.3);
axs[1].vlines([0.35,0.6],0,15,linestyle="--",alpha=0.3);
axs[0].set_ylabel("Blue Sky");
axs[1].set_ylabel("Gray Sky");

fig,axs = plt.subplots(2, figsize=(10,3), sharex=True, sharey=True);
bins = np.arange(256)
x = bins[:-1]
axs[0].plot(x,np.histogram(bluesky[skyfilt_blue][:,0],bins,normed=True)[0],c='red',lw=2, alpha=0.7, label="R");
axs[0].plot(x,np.histogram(bluesky[skyfilt_blue][:,1],bins,normed=True)[0],c='green',lw=2, alpha=0.7, label="G");
axs[0].plot(x,np.histogram(bluesky[skyfilt_blue][:,2],bins,normed=True)[0],c='blue',lw=2, alpha=0.7, label="B");
axs[1].plot(x,np.histogram(greysky[skyfilt_grey][:,0],bins,normed=True)[0],c='red',lw=2, alpha=0.7);
axs[1].plot(x,np.histogram(greysky[skyfilt_grey][:,1],bins,normed=True)[0],c='green',lw=2, alpha=0.7);
axs[1].plot(x,np.histogram(greysky[skyfilt_grey][:,2],bins,normed=True)[0],c='blue',lw=2, alpha=0.7);
# plt.ylim((0,15));
plt.xlim((100,210));
plt.yticks([]);
axs[0].legend();
axs[0].set_ylabel("Blue Sky");
axs[1].set_ylabel("Gray Sky");