def plot_batch(color_model, q_ab, X_batch_black, X_batch_color, batch_size, h, w, nb_q, epoch):

    # Format X_colorized
    X_colorized = color_model.predict(X_batch_black / 100.)[:, :, :, :-1]
    X_colorized = X_colorized.reshape((batch_size * h * w, nb_q))
    X_colorized = q_ab[np.argmax(X_colorized, 1)]
    X_a = X_colorized[:, 0].reshape((batch_size, 1, h, w))
    X_b = X_colorized[:, 1].reshape((batch_size, 1, h, w))
    X_colorized = np.concatenate((X_batch_black, X_a, X_b), axis=1).transpose(0, 2, 3, 1)
    X_colorized = [np.expand_dims(color.lab2rgb(im), 0) for im in X_colorized]
    X_colorized = np.concatenate(X_colorized, 0).transpose(0, 3, 1, 2)

    X_batch_color = [np.expand_dims(color.lab2rgb(im.transpose(1, 2, 0)), 0) for im in X_batch_color]
    X_batch_color = np.concatenate(X_batch_color, 0).transpose(0, 3, 1, 2)

    list_img = []
    for i, img in enumerate(X_colorized[:min(32, batch_size)]):
        arr = np.concatenate([X_batch_color[i], np.repeat(X_batch_black[i] / 100., 3, axis=0), img], axis=2)
        list_img.append(arr)

    plt.figure(figsize=(20,20))
    list_img = [np.concatenate(list_img[4 * i: 4 * (i + 1)], axis=2) for i in range(len(list_img) / 4)]
    arr = np.concatenate(list_img, axis=1)
    plt.imshow(arr.transpose(1,2,0))
    ax = plt.gca()
    ax.get_xaxis().set_ticks([])
    ax.get_yaxis().set_ticks([])
    plt.tight_layout()
    plt.savefig("../../figures/fig_epoch%s.png" % epoch)
    plt.clf()
    plt.close()

 def test_rgb_lch_roundtrip(self):
     rgb = img_as_float(self.img_rgb)
     lab = rgb2lab(rgb)
     lch = lab2lch(lab)
     lab2 = lch2lab(lch)
     rgb2 = lab2rgb(lab2)
     assert_array_almost_equal(rgb, rgb2)

def LAB(img, k, filename):
    # print 'lab'
    # restructure image pixel values into range from 0 to 1 - needed for library
    img = img * 1.0 / MAX_COLOR_VAL

    # convert rgb to LAB
    pixels_lab = color.rgb2lab(img)
    # remove the L channel
    L = pixels_lab[:, :, 0]

    # reshape, cluster, and retrieve quantized values
    pixels_l = np.reshape(L, (L.shape[0] * L.shape[1], 1))
    clustered = cluster_pixels(pixels_l, k, (L.shape[0], L.shape[1]))
    pixels_lab[:, :, 0] = clustered[:, :, 0]

    # convert result to 255 RGB space
    quanted_img = color.lab2rgb(pixels_lab) * MAX_COLOR_VAL
    quanted_img = quanted_img.astype('uint8')

    fig = plt.figure(1)
    plt.imshow(quanted_img)
    plt.title("LAB quantization where k is " + str(k))
    plt.savefig('Q2/' + filename + '_LAB.png')
    plt.close(fig)
    return quanted_img

def imshow_rand(im, axis=None, labrandom=True):
    """Show a segmentation using a random colormap.

    Parameters
    ----------
    im : np.ndarray of int, shape (M, N)
        The segmentation to be displayed.
    labrandom : bool, optional
        Use random points in the Lab colorspace instead of RGB.

    Returns
    -------
    fig : plt.Figure
        The image shown.
    """
    if axis is None:
        fig, axis = plt.subplots()
    rand_colors = np.random.random(size=(ceil(np.max(im)), 3))
    if labrandom:
        rand_colors[:, 0] = rand_colors[:, 0] * 81 + 39
        rand_colors[:, 1] = rand_colors[:, 1] * 185 - 86
        rand_colors[:, 2] = rand_colors[:, 2] * 198 - 108
        rand_colors = color.lab2rgb(rand_colors[np.newaxis, ...])[0]
        rand_colors[rand_colors < 0] = 0
        rand_colors[rand_colors > 1] = 1
    rcmap = cm.colors.ListedColormap(np.concatenate((np.zeros((1, 3)),
                                                     rand_colors)))
    return axis.imshow(im, cmap=rcmap)

def colorizeTest(test_img_lum, trainData):
	n = (sampling_side-1)/2
	x_min = n
	x_max = test_img_lum.shape[0] - n - 1
	y_min = n
	y_max = test_img_lum.shape[1] - n - 1

	# classify
	output = np.zeros(test_img_lum.shape + (3,))
	count = 0
	total_pixels = (x_max - x_min + 1)*(y_max - y_min + 1) / 100.0
	for (x, y), lum in np.ndenumerate(test_img_lum):

		# make sure we don't look at pixels without enough neighbors
		if x < x_min or x > x_max or y < y_min or y > y_max :
			continue

		count += 1
		if count%100 == 0:
			print 1.0 * count / total_pixels

		output[x,y,0] = lum
		stddev = computeStdDevLuminance(test_img_lum, x, y)
		
		closestIndex = getClosestTraining([lum, stddev], trainData[:,[0,3]])
		output[x, y, 1:3] = trainData[closestIndex, 1:3]

	return color.lab2rgb(output[x_min : x_max + 1, y_min : y_max + 1])

def check_HDF5(size=64):
    """
    Plot images with landmarks to check the processing
    """

    # Get hdf5 file
    hdf5_file = os.path.join(data_dir, "CelebA_%s_data.h5" % size)

    with h5py.File(hdf5_file, "r") as hf:
        data_color = hf["training_color_data"]
        data_lab = hf["training_lab_data"]
        data_black = hf["training_black_data"]
        for i in range(data_color.shape[0]):
            fig = plt.figure()
            gs = gridspec.GridSpec(3, 1)
            for k in range(3):
                ax = plt.subplot(gs[k])
                if k == 0:
                    img = data_color[i, :, :, :].transpose(1,2,0)
                    ax.imshow(img)
                elif k == 1:
                    img = data_lab[i, :, :, :].transpose(1,2,0)
                    img = color.lab2rgb(img)
                    ax.imshow(img)
                elif k == 2:
                    img = data_black[i, 0, :, :] / 255.
                    ax.imshow(img, cmap="gray")
            gs.tight_layout(fig)
            plt.show()
            plt.clf()
            plt.close()

def save_img(img_path, lab):   
    lab = clamp_lab_img(lab)    
    rgb = color.lab2rgb(np.float64(lab))
    print 'rgb min max', np.min(rgb[:, :, 0]), np.max(rgb[:, :, 0]), \
    np.min(rgb[:, :, 1]), np.max(rgb[:, :, 1]), np.min(rgb[:, :, 2]), np.max(rgb[:, :, 2])
    
    scipy.misc.imsave(img_path, rgb)

def imshow_rand(im, labrandom=True):
    """Show a segmentation using a random colormap.

    Parameters
    ----------
    im : np.ndarray of int, shape (M, N)
        The segmentation to be displayed.
    labrandom : bool, optional
        Use random points in the Lab colorspace instead of RGB.

    Returns
    -------
    fig : plt.Figure
        The image shown.
    """
    rand_colors = np.random.rand(ceil(im.max()), 3)
    if labrandom:
        rand_colors[:, 0] = rand_colors[:, 0] * 60 + 20
        rand_colors[:, 1] = rand_colors[:, 1] * 185 - 85
        rand_colors[:, 2] = rand_colors[:, 2] * 198 - 106
        rand_colors = color.lab2rgb(rand_colors[np.newaxis, ...])[0]
        rand_colors[rand_colors < 0] = 0
        rand_colors[rand_colors > 1] = 1
    rcmap = matplotlib.colors.ListedColormap(np.concatenate(
        (np.zeros((1,3)), rand_colors)
    ))
    return plt.imshow(im, cmap=rcmap, interpolation='nearest')

 def sync_buffers_from_lab(self):
     lab = self.lab.value
     rgb = color.lab2rgb(lab[:,:,:3])
     rgb *= 255
     rgb_ = np.zeros((rgb.shape[0], rgb.shape[1], 4), dtype=np.uint8)
     rgb_[:,:,:3] = rgb[:,:,:3]
     self.rgb.value = rgb_   

 def test_lab_rgb_outlier(self):
     lab_array = np.ones((3, 1, 3))
     lab_array[0] = [50, -12, 85]
     lab_array[1] = [50, 12, -85]
     lab_array[2] = [90, -4, -47]
     rgb_array = np.array([[[0.501, 0.481, 0]], [[0, 0.482, 1.0]], [[0.578, 0.914, 1.0]]])
     assert_almost_equal(lab2rgb(lab_array), rgb_array, decimal=3)

def JPEG_decompression(data, channels=3):
	#Meta Data
	height = data[-1]
	width = data[-2]
	quality = data[-3]
	d_height = data[-4]
	d_width = data[-5]

	#Remove Meta Data
	data = data[:-5]
	#Unzigzag
	data_z = zigzag_decode(data, d_height, d_width, channels)
	#Unquantize
	im_q = unquantize(data_z,quality)
	#IDCT
	im_idct = idct_2d(im_q)
	#Unblock and Unpad
	im = unblock_image(im_idct,d_height,d_width)
	#Upsample

	#Undo offset and return to RGB
	im[:,:,[1,2]] -= 128
	im = lab2rgb(im) * 255 # lab2rgb converts to float64
	#Undo Padding
	im = im[:d_height,:d_width]

	#Upsample
	im = utils.upsample(im,(height,width))

	return im.astype(np.uint8)

def sshow(im, labrandom=True):
    """Show a segmentation (or cross-section) using a random colormap.

    Parameters
    ----------
    im : np.ndarray of int
        The segmentation to be displayed.
    labrandom : bool, optional
        Use random points in the Lab colorspace instead of RGB.

    Returns
    -------
    ax : matplotlib AxesImage object
        The figure axes.
    """
    if im.ndim > 2:
        mid = im.shape[0] // 2
        ax = sshow(im[mid], labrandom)
    else:
        rand_colors = np.random.rand(np.ceil(im.max()), 3)
        if labrandom:
            rand_colors[:, 0] = rand_colors[:, 0] * 60 + 20
            rand_colors[:, 1] = rand_colors[:, 1] * 185 - 85
            rand_colors[:, 2] = rand_colors[:, 2] * 198 - 106
            rand_colors = color.lab2rgb(rand_colors[np.newaxis, ...])[0]
            rand_colors[rand_colors < 0] = 0
            rand_colors[rand_colors > 1] = 1
        rcmap = colors.ListedColormap(np.concatenate((np.zeros((1, 3)),
                                                      rand_colors)))
        ax = plt.imshow(im, cmap=rcmap, interpolation='nearest')
    return ax

def run_color(image, image_out):
    caffe.set_mode_cpu()
    net = caffe.Net('colorization_deploy_v0.prototxt', 'colorization_release_v0.caffemodel', caffe.TEST)

    (H_in,W_in) = net.blobs['data_l'].data.shape[2:] # get input shape
    (H_out,W_out) = net.blobs['class8_ab'].data.shape[2:] # get output shape
    net.blobs['Trecip'].data[...] = 6/np.log(10) # 1/T, set annealing temperature
    
    img_rgb = caffe.io.load_image(image)
    img_lab = color.rgb2lab(img_rgb) # convert image to lab color space
    img_l = img_lab[:,:,0] # pull out L channel
    (H_orig,W_orig) = img_rgb.shape[:2] # original image size

    # resize image to network input size
    img_rs = caffe.io.resize_image(img_rgb,(H_in,W_in)) # resize image to network input size
    img_lab_rs = color.rgb2lab(img_rs)
    img_l_rs = img_lab_rs[:,:,0]

    net.blobs['data_l'].data[0,0,:,:] = img_l_rs-50 # subtract 50 for mean-centering
    net.forward() # run network

    ab_dec = net.blobs['class8_ab'].data[0,:,:,:].transpose((1,2,0)) # this is our result
    ab_dec_us = sni.zoom(ab_dec,(1.*H_orig/H_out,1.*W_orig/W_out,1)) # upsample to match size of original image L
    img_lab_out = np.concatenate((img_l[:,:,np.newaxis],ab_dec_us),axis=2) # concatenate with original image L
    img_rgb_out = np.clip(color.lab2rgb(img_lab_out),0,1) # convert back to rgb

    scipy.misc.imsave(image_out, img_rgb_out)

def dominant_colors(image, num_colors, mask=None):
    """Reduce image colors to a representative set of a given size.

    Args:
        image (ndarray): BGR image of shape n x m x 3.
        num_colors (int): Number of colors to reduce to.
        mask (array_like, optional): Foreground mask. Defaults to None.

    Returns:
        list: The list of Color objects representing the most dominant colors in the image.

    """
    image = rgb2lab(image / 255.0)

    if mask is not None:
        data = image[mask > 250]
    else:
        data = np.reshape(image, (-1, 3))

    # kmeans algorithm has inherent randomness - result will not be exactly the same
    # every time. Fairly consistent with >= 30 iterations
    centroids, labels = kmeans2(data, num_colors, iter=30)
    counts = np.histogram(labels, bins=range(0, num_colors + 1), normed=True)[0]

    centroids_RGB = lab2rgb(centroids.reshape(-1, 1, 3))[:, 0, :] * 255.0
    colors = [Color(centroid, count) for centroid, count in zip(centroids_RGB, counts)]
    colors.sort(key=lambda color: np.mean(color.BGR))

    return colors

def snap_ab(input_l, input_rgb, return_type='rgb'):
    ''' given an input lightness and rgb, snap the color into a region where l,a,b is in-gamut
    '''
    T = 20
    warnings.filterwarnings("ignore")
    input_lab = rgb2lab_1d(np.array(input_rgb))  # convert input to lab
    conv_lab = input_lab.copy()  # keep ab from input
    for t in range(T):
        conv_lab[0] = input_l  # overwrite input l with input ab
        old_lab = conv_lab
        tmp_rgb = color.lab2rgb(conv_lab[np.newaxis, np.newaxis, :]).flatten()
        tmp_rgb = np.clip(tmp_rgb, 0, 1)
        conv_lab = color.rgb2lab(tmp_rgb[np.newaxis, np.newaxis, :]).flatten()
        dif_lab = np.sum(np.abs(conv_lab-old_lab))
        if dif_lab < 1:
            break
        # print(conv_lab)

    conv_rgb_ingamut = lab2rgb_1d(conv_lab, clip=True, dtype='uint8')
    if (return_type == 'rgb'):
        return conv_rgb_ingamut

    elif(return_type == 'lab'):
        conv_lab_ingamut = rgb2lab_1d(conv_rgb_ingamut)
        return conv_lab_ingamut

def save_image_LAB_into_sRGB(img_path, lab):
    lab = np.float64(lab)
    srgb = color.lab2rgb(lab)
    # clamping to [0,1]
    idx1, idx2 = np.nonzero(srgb < 0), np.nonzero(srgb > 1)
    srgb[idx1[0], idx1[1], idx1[2]] = 0
    srgb[idx2[0], idx2[1], idx2[2]] = 1
    scipy.misc.imsave(img_path, srgb)

def lab2rgb_1d(in_lab, clip=True, dtype='uint8'):
    warnings.filterwarnings("ignore")
    tmp_rgb = color.lab2rgb(in_lab[np.newaxis, np.newaxis, :]).flatten()
    if clip:
        tmp_rgb = np.clip(tmp_rgb, 0, 1)
    if dtype == 'uint8':
        tmp_rgb = np.round(tmp_rgb * 255).astype('uint8')
    return tmp_rgb

def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01,
                       nbins=256):
    """Contrast Limited Adaptive Histogram Equalization.

    Parameters
    ----------
    image : array-like
        Input image.
    ntiles_x : int, optional
        Number of tile regions in the X direction.  Ranges between 2 and 16.
    ntiles_y : int, optional
        Number of tile regions in the Y direction.  Ranges between 2 and 16.
    clip_limit : float: optional
        Clipping limit, normalized between 0 and 1 (higher values give more
        contrast).
    nbins : int, optional
        Number of gray bins for histogram ("dynamic range").

    Returns
    -------
    out : ndarray
        Equalized image.

    Notes
    -----
    * The algorithm relies on an image whose rows and columns are even
      multiples of the number of tiles, so the extra rows and columns are left
      at their original values, thus  preserving the input image shape.
    * For color images, the following steps are performed:
       - The image is converted to LAB color space
       - The CLAHE algorithm is run on the L channel
       - The image is converted back to RGB space and returned
    * For RGBA images, the original alpha channel is removed.

    References
    ----------
    .. [1] http://tog.acm.org/resources/GraphicsGems/gems.html#gemsvi
    .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
    """
    args = [None, ntiles_x, ntiles_y, clip_limit * nbins, nbins]
    if image.ndim > 2:
        lab_img = color.rgb2lab(skimage.img_as_float(image))
        l_chan = lab_img[:, :, 0]
        l_chan /= np.max(np.abs(l_chan))
        l_chan = skimage.img_as_uint(l_chan)
        args[0] = rescale_intensity(l_chan, out_range=(0, NR_OF_GREY - 1))
        new_l = _clahe(*args).astype(float)
        new_l = rescale_intensity(new_l, out_range=(0, 100))
        lab_img[:new_l.shape[0], :new_l.shape[1], 0] = new_l
        image = color.lab2rgb(lab_img)
        image = rescale_intensity(image, out_range=(0, 1))
    else:
        image = skimage.img_as_uint(image)
        args[0] = rescale_intensity(image, out_range=(0, NR_OF_GREY - 1))
        out = _clahe(*args)
        image[:out.shape[0], :out.shape[1]] = out
        image = rescale_intensity(image)
    return image

def applyNailPolish(x , y , r = Rg, g = Gg, b = Bg):
	val = color.rgb2lab((im[x, y]/255.).reshape(len(x), 1, 3)).reshape(len(x), 3)
	L, A, B = mean(val[:,0]), mean(val[:,1]), mean(val[:,2])
	L1, A1, B1 = color.rgb2lab(np.array((r/255., g/255., b/255.)).reshape(1, 1, 3)).reshape(3,)
	ll, aa, bb = L1 - L, A1 - A, B1 - B
	val[:, 0] = np.clip(val[:, 0] + ll, 0, 100)
	val[:, 1] = np.clip(val[:, 1] + aa, -127, 128)
	val[:, 2] = np.clip(val[:, 2] + bb, -127, 128)
	im[x, y] = color.lab2rgb(val.reshape(len(x), 1, 3)).reshape(len(x), 3)*255

def fuckwith(image, p, colors=None):
    
    if isinstance(image, Image.Image):
        pil = True
        #image = np.array(image.convert('RGB'))
    else:
        pil = False

    a, b, t = p
    h, w = image.size
    
    lab = color.rgb2lab(image)
    
    for y in range(h):
        lmax = np.max(lab[y,:,0])
        lavg = np.mean(lab[y,:,0])
        saw = ( (y + 25*t) % 25 )*0.15
        shift = int( (lavg*0.1 + lmax*0.15)*a + rand()*5 + saw )
    
        lab[y,:,:1] = np.roll(lab[y,:,:1], shift, 0)
        lab[y,:,2] = np.roll(lab[y,:,2], -shift//2, 0)
        
        lab[y,:,0] += lmax * ( (1.2 ** a) - 1 )
        lab[y,:,0] *= ((a * 0.8) + 1.0)
        lab[y,:,0] += (a*5) - 1
        
        lab[y,:,1] += lmax * 0.2 * a
        lab[y,:,2] += lmax * 0.5 * a
    
    for x in range(w):
        lmax = np.max(lab[:,x,0])
        #shift = int(lmax*0.5*b + rand()*2)
        #lab[:,x] = np.roll(lab[:,x], shift, 0)
        lab[:,x,0] -= lmax * ( (1.2 ** b) - 1 )
        lab[:,x,1] -= lmax * 0.2 * b
        lab[:,x,2] -= lmax * 0.3 * b
    
    
    rgb = color.lab2rgb(lab)
    rgb = np.clip(rgb, 0, 1)
       
    if colors:
        rgb_scaled = (rgb*255).astype('uint8')  
        mod = Image.fromarray(rgb_scaled)
        mod = mod.convert('P', palette=Image.ADAPTIVE, colors=colors)
        if pil == True:
            return mod
        else:
            mod = mod.convert('RGB')
            pix = np.array(mod)
            return pix
    elif pil == True:
        rgb_scaled = (rgb*255).astype('uint8')  
        mod = Image.fromarray(rgb_scaled)
        return mod
    else:
        return rgb