#! /usr/bin/env python
""" 
This script aligns a set of Spatial Transcriptomics datasets using the first
image as reference. It provides different alignment methods (manual
and automatic) and it outputs the aligned images, the aligned counts
matrices and the alignment spot coordinates.

@Author Jose Fernandez Navarro <jc.fernandez.navarro@gmail.com>
"""
import tissue_recognition as tr
import numpy as np
from matplotlib import pyplot as plt
import matplotlib
from matplotlib.image import AxesImage
import cv2
import os
import pandas as pd
import sys
import argparse
from sklearn.neighbors import NearestNeighbors
from skimage.measure import regionprops
from skimage import filters

def best_fit_transform(A, B):
    """
    Calculates the least-squares best-fit transform that maps 
    corresponding points A to B in m spatial dimensions
    Input:
      A: Nxm numpy array of corresponding points
      B: Nxm numpy array of corresponding points
    Returns:
      T: (m+1)x(m+1) homogeneous transformation matrix that maps A on to B
      R: mxm rotation matrix
      t: mx1 translation vector
    """
    assert A.shape == B.shape
    # get number of dimensions
    m = A.shape[1]
    # translate points to their centroids
    centroid_A = np.mean(A, axis=0)
    centroid_B = np.mean(B, axis=0)
    AA = A - centroid_A
    BB = B - centroid_B
    # rotation matrix
    H = np.dot(AA.T, BB)
    U, S, Vt = np.linalg.svd(H)
    R = np.dot(Vt.T, U.T)
    # special reflection case
    if np.linalg.det(R) < 0:
       Vt[m-1,:] *= -1
       R = np.dot(Vt.T, U.T)
    # translation
    t = centroid_B.T - np.dot(R,centroid_A.T)
    # homogeneous transformation
    T = np.identity(m+1)
    T[:m, :m] = R
    T[:m, m] = t
    return T, R, t

def nearest_neighbor(src, dst):
    """
    Find the nearest (Euclidean) neighbor in dst for each point in src
    Input:
        src: Nxm array of points
        dst: Nxm array of points
    Returns:
        distances: Euclidean distances of the nearest neighbor
        indices: dst indices of the nearest neighbor
    """
    assert src.shape == dst.shape
    neigh = NearestNeighbors(n_neighbors=1)
    neigh.fit(dst)
    distances, indices = neigh.kneighbors(src, return_distance=True)
    return distances.ravel(), indices.ravel()

def icp(A, B, init_pose=None, max_iterations=50, tolerance=0.001):
    """
    The Iterative Closest Point method: finds best-fit transform
    that maps points A on to points B
    Input:
        A: Nxm numpy array of source mD points
        B: Nxm numpy array of destination mD point
        init_pose: (m+1)x(m+1) homogeneous transformation
        max_iterations: exit algorithm after max_iterations
        tolerance: convergence criteria
    Returns:
        T: final homogeneous transformation that maps A on to B
        distances: Euclidean distances (errors) of the nearest neighbor
        i: number of iterations to converge
    """
    assert A.shape == B.shape
    # get number of dimensions
    m = A.shape[1]
    # make points homogeneous, copy them to maintain the originals
    src = np.ones((m+1,A.shape[0]))
    dst = np.ones((m+1,B.shape[0]))
    src[:m,:] = np.copy(A.T)
    dst[:m,:] = np.copy(B.T)
    # apply the initial pose estimation
    if init_pose is not None:
        src = np.dot(init_pose, src)
    prev_error = 0
    for i in range(max_iterations):
        # find the nearest neighbors between the current source and destination points
        distances, indices = nearest_neighbor(src[:m,:].T, dst[:m,:].T)
        # compute the transformation between the current source and nearest destination points
        T,_,_ = best_fit_transform(src[:m,:].T, dst[:m,indices].T)
        # update the current source
        src = np.dot(T, src)
        # check error
        mean_error = np.mean(distances)
        if np.abs(prev_error - mean_error) < tolerance:
            break
        prev_error = mean_error
    # calculate final transformation
    T,_,_ = best_fit_transform(A, src[:m,:].T)
    return T, distances, i

def plot_images(image_list, is_gray=False, filename=None):
    """
    Helper function to debug images
    Input:
        image_list: list of images (OpenCV)
        is_gray: whether the images are in grayscale or not
        filename: name of the output file

    Returns: None

    """
    columns = int(round(np.sqrt(len(image_list))))
    rows = int(np.ceil(np.sqrt(len(image_list))))
    fig = plt.figure(figsize = (16, 16))
    for i in range(1, columns*rows + 1):
        try:
            img = image_list[i - 1]
            ax = fig.add_subplot(rows, columns, i)
            if is_gray:
                plt.gray()
            ax.imshow(img)
        except:
            continue
    if filename is not None:
        plt.savefig(filename)
    else:
        plt.show()
    plt.clf()
    plt.cla()
    
def get_binary_masks(image_list):
    """
    Obtains a binary mask of the image using
    the tissue recognition library (automatic).
    Returns the mask (1 inside, 0 outside).
    Input:
        - A list of images
    Returns:
        - A list of binary masks
    """
    assert(len(image_list) > 0)
    binary_masks = list()
    for img in image_list:
        msk = np.zeros(img.shape[0:2], dtype=np.uint8)
        tr.recognize_tissue(img, msk)
        binary_masks.append(tr.get_binary_mask(msk))
    return binary_masks
    
def mask_images(image_list, binary_masks):
    """
    Applies a binary masks to a list of images to keep only
    the pixels inside tissue.
    Input:
        - A list of images
        - A list of binary masks
    Returns:
        - A list of masked images
    """
    assert(len(image_list) == len(binary_masks))
    masked_images = list()
    for img, bin_msk in zip(image_list, binary_masks):
        h, w, z = img.shape
        bin_msk = bin_msk > 0.5
        comb_img = 255 - np.zeros((h, w, 3), dtype=np.uint8)
        comb_img[bin_msk] = img[bin_msk]
        masked_images.append(comb_img)
    return masked_images

def center_images(image_list, binary_masks):
    """
    Center images and binary masks applying the center of mass.
    Input:
        - A list of images
        - A list of binary masks
    Returns:
        - A list of centered images
        - A list of centered binary masks
        - A list of affine transformations (applied to center)
    """
    assert(len(image_list) == len(binary_masks))
    centered_images = list()
    centered_maks = list()
    transf_matrices = list()
    for image, bin_msk in zip(image_list, binary_masks):
        h, w, _ = image.shape
        threshold_value = filters.threshold_otsu(bin_msk)
        labeled_foreground = (bin_msk > threshold_value).astype(int)
        properties = regionprops(labeled_foreground, bin_msk)
        center_of_mass = properties[0].centroid
        translation = np.float32([[1, 0, (w / 2) - center_of_mass[1]],
                                  [0, 1, (h / 2) - center_of_mass[0]]])
        img_centered = cv2.warpAffine(image, translation, (w, h), 
                                      borderValue=(255, 255, 255))
        mask_centered = cv2.warpAffine(bin_msk, translation, (w, h), 
                                      borderValue=(0, 0, 0))
        centered_images.append(img_centered)
        transf_matrices.append(translation)
        centered_maks.append(mask_centered)
    return centered_images, centered_maks, transf_matrices

def add_borders(image_list, size=1000):
    """
    Adds a white border to images
    Input:
        - A list of images
        - Amount of border to add (in pixels)
    Returns:
        - A list of images with borders added (white)
    """
    images_border = list()
    for image in image_list:
        border_image = cv2.copyMakeBorder(image, 
                                          size, size, size, size,
                                          cv2.BORDER_CONSTANT,
                                          value=(255, 255, 255))
        images_border.append(border_image)
    return images_border

def remove_borders(image_list, size=1000):
    """
    Removes a white border of images
    Input:
        - A list of images
        - Amount of border to remove (in pixels)
    Returns:
        - A list of images with borders removed (white)
    """
    images_no_border = list()
    for image in image_list:
        h, w, _ = image.shape
        no_border_image = image[size:h-size, size:w-size]
        images_no_border.append(no_border_image)
    return images_no_border
    
def mse(imageA, imageB):
    """
    The 'Mean Squared Error' between the two images is the
    sum of the squared difference between the two images;
    NOTE: the two images must have the same dimension
    returns the MSE, the lower the error, the more "similar"
    the two images are
    """
    assert(imageA.shape == imageB.shape)
    err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2)
    err /= float(imageA.shape[0] * imageA.shape[1])
    return err