/Alex/retrieve_features.py

import cv2
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
import seaborn as sns
import sys
HOME_DIR = '/Users/agalicina/Term10/genehack/'  ##change for your machine
sys.path.append(HOME_DIR+'genehack2016')
from process import *
from sklearn.metrics import roc_curve, auc
import math

def read_img(path):
  return cv2.imread(path, 0)

##########################
#   GAUSSIAN ASSESSMENT
##########################

def get_spectrum(img):
    f = np.fft.fft2(img)
    fshift = np.fft.fftshift(f)
    magnitude_spectrum = fshift
    return magnitude_spectrum
def plot_doub(num1, num2, img, spec):
    plt.subplot(num1),plt.imshow(img, cmap = 'gray')
    plt.title('Input Image'), plt.xticks([]), plt.yticks([])
    plt.subplot(num2),plt.imshow(spec, cmap = 'gray')
    plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
def plot(num1,spec):
    plt.subplot(num1),plt.imshow(spec, cmap = 'gray')
    plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])

def compute_sum_gaussian(img_mask, img_pic, sigma=10, size=150, mode='abs'):

    #img_mask = cv2.resize(cv2.imread(img_path_mask,0), (150,150))
    #img_pic  = cv2.resize(cv2.imread(img_path_pic, 0), (150,150))
    img_mask = cv2.resize(img_mask, (150,150))
    img_pic  = cv2.resize(img_pic, (150,150))

    spec_mask = get_spectrum(img_mask)
    spec_pic  = get_spectrum(img_pic)

    x = cv2.getGaussianKernel(size,sigma)
    gaussian = x*x.T

    if mode=='abs':
      res = gaussian*(np.abs(spec_pic-spec_mask))
    else:
      res = gaussian*(spec_pic)

    return np.sum(res)

def plot_gaussian_steps():
#### TO DO: check and rewrite

    img_mask = cv2.resize(cv2.imread(img_path_mask,0), (150,150))
    img_pic  = cv2.resize(cv2.imread(img_path_pic, 0), (150,150))

    spec_mask = get_spectrum(img_mask)
    spec_pic  = get_spectrum(img_pic)

    pl1 = 20*np.log(np.abs(spec_mask))
    pl2 = 20*np.log(np.abs(spec_pic))
    pl3 = 20*np.log(np.abs(spec_pic-spec_mask))
    plt.figure(figsize=(10,10))
    plot_doub(321,322,img_mask, pl1)
    plot_doub(323,324,img_pic, pl2)
    plot(326, pl3)
    plt.show()

##########################
#   CONTOUR ASSESSMENT
##########################

def find_contour(img):
    ret, thresh = cv2.threshold(img,127,255,0)
    contours, hierarchy = cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)
    return contours[0]

def get_perimeter(cnt):
    perimeter = cv2.arcLength(cnt,True)
    return perimeter


def plot_contour(img, cnt, path):
    img1 = img.copy()
    cv2.drawContours(img1,[cnt],0,(0,0,255),2)
    cv2.imwrite(path, img1)

    '''cv2.namedWindow("imc")
    cv2.imshow("imc", img1)
    cv2.waitKey(0)'''

def get_defects(cnt):
    hull = cv2.convexHull(cnt,returnPoints = False)
    defects = cv2.convexityDefects(cnt,hull)
    return defects

def plot_and_count_defects(img, defects, cnt, depth_tresh):
    img1 = img.copy()
    true_def_count = 0
    max_def = 0

    n=0
    for i in range(defects.shape[0]):
        n+=1
        s,e,f,d = defects[i,0]
        start = tuple(cnt[s][0])
        end = tuple(cnt[e][0])
        far = tuple(cnt[f][0])
        med = (min(start[0],end[0])+(abs(start[0]-end[0])/2), min(start[1],end[1])+(abs(start[1]-end[1])/2))
        depth = math.sqrt(((med[0]-far[0])**2)+((med[1]-far[1])**2))
        if n != 1:
            max_def = max(max_def, depth)
            if depth > depth_tresh:
                true_def_count += 1
                cv2.circle(img1,far,2,[0,0,255],-1) #red!
            else:
                cv2.circle(img1,far,2,[255,0,0],-1) #blue!
            cv2.line(img1,start,end,[0,255,0],1)

    '''cv2.imshow('img',img1)
    cv2.waitKey(0)
    cv2.destroyAllWindows()'''

    #cv2.imwrite(path, img1)
    return true_def_count, max_def

def solidity(cnt, area):
    hull = cv2.convexHull(cnt)
    hull_area = cv2.contourArea(hull)
    return float(area)/hull_area

def ellipse_jaccard(cnt, imgray):
    #doesn't work
    ellipse = cv2.fitEllipse(cnt)
    blank = imgray.copy()
    blank[:] = 0
    e_img = cv2.ellipse(blank, ellipse, 255, -1)
    union_img = np.logical_or(imgray, e_img)
    inter_img = np.logical_and(imgray, e_img)
    union = sum([sum(x) for x in union_img])
    inter = sum([sum(x) for x in inter_img])
    return inter/float(union)

def calculate_all(path, depth_tresh):
    img = cv2.imread(path)
    imgray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    cnt = find_contour(imgray)
    area = cv2.contourArea(cnt)
    defects = get_defects(cnt)
    true_def_count, max_def = plot_and_count_defects(imgray, defects, cnt, depth_tresh)
    s = solidity(cnt, area)
    #jac = ellipse_jaccard(cnt, imgray)
    return true_def_count, max_def, s, 0


#########################
#   Hough transform
#########################

def get_Hough_lines_lens(img, par_len=3, par_gap=0.5):

    edges = cv2.Canny(img,0,0)
    P = get_perimeter(find_contour(img))

    minLineLength = par_len*P
    maxLineGap = par_gap*P

    lines = cv2.HoughLinesP(edges, 1, np.pi/90, 10, 100, minLineLength=minLineLength, maxLineGap=maxLineGap)
    len_sum = 0
    if lines is None:
        return 0

    for x1,y1,x2,y2 in lines[0]:
        len_sum += math.sqrt(((x1-x2)**2)+((y1-y2)**2))

    return len_sum/float(P)

def plot_Hough(img, saveto, par_len=3, par_gap=0.5):

    edges = cv2.Canny(img,0,0)
    P = get_perimeter(find_contour(img))

    minLineLength = par_len*P
    maxLineGap = par_gap*P

    img1 = img.copy()
    img1.fill(1)

    len_sum = 0
    lines = cv2.HoughLinesP(edges, 1, np.pi/90, 10, 100, minLineLength=minLineLength, maxLineGap=maxLineGap)
    if lines is None:
        return 0
    for x1,y1,x2,y2 in lines[0]:
        cv2.line(img1,(x1,y1),(x2,y2),(0,255,0),1)
        len_sum += math.sqrt(((x1-x2)**2)+((y1-y2)**2))

    plt.clf()
    plt.figure(figsize=(10,10))
    plt.subplot(121)
    plt.imshow(img1)
    plt.title(
            'Hough transformation with:\n minLen = {} ({} of P), \n maxGap = {} ({} of P), \n len_sum = {}, norm = {}'.format(
                    minLineLength, par_len, maxLineGap, par_gap, len_sum, len_sum/float(P))),\
            plt.xticks([]), plt.yticks([])
    plt.subplot(122)
    plt.imshow(img, cmap = 'gray')
    plt.title('Input Image'), plt.xticks([]), plt.yticks([])
    plt.savefig(saveto)

#########################
#   Intensity measurements
#########################

def get_intensity_max(img):
  return np.max(img)

def get_intensity_sum(img):
  return np.sum(img)

def get_intensity_mean(img):
  cnt = find_contour_my(img)
  return np.sum(img)/get_area(cnt)


#########################
#     Apply to images and compute features FINALLY!
#########################

### Load images from file, apply function and return array with values
#img_dir = HOME_DIR+'processed_images'

def my_apply(func, img_dir = 'processed_images/', all_imgs=None, y_test=None, **kwargs_dict):
  if all_imgs is None:
    all_imgs, y_test = load_data(img_dir)
  y_score = [func(all_imgs[i], **kwargs_dict) for i in range(len(y_test))]
  return [np.array(l) for l in (all_imgs, y_score, y_test)] # (imgs, ys_score, ys_test)

def my_apply_tolist(func, all_imgs=None, y_test=None, **kwargs_dict):
  y_score = [func(all_imgs[i], **kwargs_dict) for i in range(len(y_test))]
  return [np.array(l) for l in (all_imgs, y_score, y_test)] # (imgs, ys_score, ys_test)

def plot_hist(y_test, y_score, saveto, **kwargs):
  '''
  :param x: all values after you applied your function
  :param y: all categories list. must be 0 or 1
  :return: histogram saved into file

  Example usage:
   rf.plot_hist([1,2,3,4,5,6], [1,0,1,0,1,0], 'ex.png', bins=1, hist_kws={"alpha":0.5})

  '''
  #df = pd.DataFrame({'x':x, 'y':y})
  #sns_plot =
  plt.clf()
  plt.close()
  plt.figure()
  sns.set(style="white", palette="muted", color_codes=True)

  sns_plot = sns.distplot(np.array(y_score)[(np.array(y_test)==1)], color='green', **kwargs)
  sns.distplot(np.array(y_score)[np.array(y_test)==0], color='red', **kwargs)

  plt.savefig(saveto)

def plot_roc(y_test, y_score, saveto):

  # Compute ROC curve and ROC area for each class
  fpr = dict()
  tpr = dict()
  roc_auc = dict()
  fpr, tpr, _ = roc_curve(y_test, y_score)
  roc_auc = auc(fpr, tpr)

  ##############################################################################
  # Plot of a ROC curve for a specific class
  plt.figure()
  plt.plot(fpr, tpr, label='ROC curve (area = %0.2f)' % roc_auc)
  plt.plot([0, 1], [0, 1], 'k--')
  plt.xlim([0.0, 1.0])
  plt.ylim([0.0, 1.05])
  plt.xlabel('False Positive Rate')
  plt.ylabel('True Positive Rate')
  plt.title('ROC')
  plt.legend(loc="lower right")
  plt.savefig(saveto)
  return roc_auc