"""
Created on Tue Oct 23 08:06:37 2012

@author: Jean-Patrick
"""
import os
import numpy as np
import matplotlib.pyplot as plt
import skimage.io as sio
import skimage as sk
import skimage.viewer as viewer
from skimage.draw import ellipse
from skimage.measure import find_contours, approximate_polygon, \
    subdivide_polygon
import cmath
from math import pi
import itertools
    
def addpixelborder(im):
    x,y = im.shape
    dy = np.zeros((y,))
    
    im = np.vstack((im,dy))
    im = np.vstack((dy,im))
    
    im = np.transpose(im)
    x,y = im.shape
    
    dy = np.zeros((y,))
    im = np.vstack((im,dy))
    im = np.vstack((dy,im))
    im = np.transpose(im)
    return im

def angle(p1,p2,p3):
    '''get the angle between the vectors p2p1,p2p3  with p a point
    '''
    z1=complex(p1[0],p1[1])
    z2=complex(p2[0],p2[1])
    z3=complex(p3[0],p3[1])
    v12=z1-z2
    v32=z3-z2
    #print v12,v32
    angle=(180/pi)*cmath.phase(v32/v12)
    return angle
    
def angleBetweenSegments(contour):
    angles = []
    points = list(contour)
    points.pop()
    #print 'contour in def',len(contour)
    l = len(points)
    for i in range(len(points)+1):
        #print 'i',i,' i-1',i-1,' (i+1)%l',(i+1)%l
        prv_point=points[(i-1)%l]
        mid_point=points[i % l]
        nex_point=points[(i+1)%l]
        
        angles.append(angle(prv_point,mid_point,nex_point))
    return angles

def quadArea(A,B,C,D):
    x1 = A[0]
    y1 = A[1]
    x2 = B[0]
    y2 = B[1]
    x3 = C[0]
    y3 = C[1]
    x4 = D[0]
    y4 = D[1]
    #sign +/- if direct/indirect quadrilateral
    return 0.5*(x1*y2+x2*y3+x3*y4+x4*y1-x2*y1-x3*y2-x4*y3-x1*y4)

def maximalQuad(A,B,C,D):
   maxarea = 0
   for quad in itertools.permutations((A,B,C,D)):
        area = quadArea(quad[0],quad[1],quad[2],quad[3])
        if area > maxarea:
            maxarea = area
            maxquad = quad
   return maxquad
        


user=os.path.expanduser("~")
#workdir=os.path.join(user,"QFISH","JPPAnimal","JPP52","11","DAPI","particles")
#image='part25.png'
#QFISH/CytoProject/Jpp48/1/DAPI/particles/part9.png
workdir=os.path.join(user,"QFISH","CytoProject","Jpp48","1","DAPI","particles")
image='part9.png'
complete_path=os.path.join(workdir,image)
#print complete_path
im = sio.imread(complete_path)
im = addpixelborder(im)
original_im = np.copy(im)
fig, ax1= plt.subplots(ncols=1, figsize=(10, 10))
#subplot(111,xticks=[],yticks=[])
contours = find_contours(im, 0.01)
p = subdivide_polygon(contours[0], degree=5, preserve_ends=True)
#coord = approximate_polygon(p, tolerance=1.5)
coord2 = approximate_polygon(p, tolerance=2)


ax1.imshow(im, interpolation='nearest', cmap=plt.cm.gray)

contour0=contours[0]
pf=contours[-1]
angles = angleBetweenSegments(coord2)
table = np.array([coord2[0][0],coord2[0][1],angles[0]])
for i in range(len(angles[1:])):
    current =  np.array([coord2[i][0],coord2[i][1],angles[i]])
    table = np.vstack((table,current))


pos = table[np.where(np.logical_and(table[:,2]>80, table[:,2]<145))]
neg = table[np.where(np.logical_and(table[:,2]>-130, table[:,2]<0))]
flat = table[np.where(np.logical_and(table[:,2]>=-180, table[:,2]<=-140))]

for n, contour in enumerate(contours):
    ax1.plot(contour[:, 1], contour[:, 0], linewidth=4)
    ax1.plot(coord2[:, 1], coord2[:, 0], linewidth=1, color='r')
    #ax1.plot(coord2[:, 1], coord2[:, 0], linewidth=1, color='y')
ax1.scatter(contour0[0][1],contour0[0][0],s=150, c='m',marker='o')
ax1.scatter(pf[0][1],pf[0][0],s=150, c='y',marker='*')

#print table[pos]
ax1.scatter(pos[:,1], pos[:,0], s=200, c='r')
ax1.scatter(neg[:,1], neg[:,0], s=200, c='y')
ax1.scatter(flat[:,1], flat[:,0], s=200, c='g')

plt.show()

neg[:,0:3]

neg[0,0:2]

coord = []
for r in range(0, neg[:,0:2].shape[0]):
    coord.append(neg[r,0:2])
quad_it = itertools.combinations(coord,4)

q=next(quad_it)
q

maxQuad = maximalQuad(*q)
print maxQuad
print quadArea(*maxQuad)


from matplotlib.patches import Polygon
p = Polygon(maxQuad, alpha=0.5, color='b' )
fig, ax1= plt.subplots(ncols=1, figsize=(8, 8))
ax1.imshow(np.transpose(im), interpolation='nearest', cmap=plt.cm.gray)
ax1.scatter(neg[:,0], neg[:,1], s=100, c='y')
ax1.add_artist(p)

from skimage.draw import polygon
print im.shape

overlapp = np.zeros(im.shape, dtype=int)
print overlapp.shape
x = np.asarray([int(r[0]) for r in maxQuad])
y = np.asarray([int(r[1]) for r in maxQuad])
print x
print y
rr,cc = polygon(x, y)
overlapp[rr,cc] = True
cut = np.logical_and(im>0,np.logical_not(overlapp))
subplot(121, xticks=[],yticks=[])
figsize(12,12)
imshow(cut, interpolation ='nearest')
subplot(122, xticks=[],yticks=[])
title('overlapp domain')
imshow(overlapp, interpolation ='nearest')

import mahotas as mh
labcut = sk.morphology.label(cut, neighbors=4)
labsize = mh.labeled.labeled_size(labcut)
#print labsize
#print labsize < 2
#print np.where(labsize < 2)
smallimage = labcut==np.where(labsize < 2)
labcut = mh.labeled.remove_regions(labcut,4)
mh.labeled.relabel(labcut, inplace=True)
overlapp = (labcut.max()+1)* (smallimage+overlapp > 0)
figsize(4,4)
subplot(121,xticks=[],yticks=[])
title(str(mh.labeled.labeled_size(labcut)[1:]))
#print np.histogram(labcut, bins=1)
imshow(labcut, interpolation = 'nearest')
subplot(122,xticks=[],yticks=[])
title('overlapp lab:'+str(overlapp.max()))
imshow(smallimage+overlapp, interpolation = 'nearest')

from matplotlib import colors
testim = np.array([[0,1,1,0],[0,2,0,0],[0,0,3,3],[4,0,0,0]])
elements = [labcut==l for l in range(1,labcut.max()+1)]
#print len (elements)
figsize(6,6)
    
elements2=[]
mycmap = colors.ListedColormap(['black','red','cyan','blue','green','yellow'])
bounds = [0,1,2,3,4,5]
mynorm = colors.BoundaryNorm(bounds, mycmap.N+1)

for imlab in range(1,labcut.max()+1):
    cur = labcut == imlab
    cur = cur * imlab
    #print imlab, cur.min(), cur.max()
    elements2.append(cur)
    subplot(2,2,imlab, xticks=[],yticks=[])
    title(str(cur.max())+' label='+str(imlab))
    #print cur.max()
    imshow(cur, interpolation='nearest',cmap=mycmap, norm=mynorm)

    
    

candidates_it = itertools.combinations(elements2,2)
candidates = [c for c in candidates_it]
assemble = []
for c in candidates:
    image = np.zeros(im.shape, dtype = int)
    for i in c:
        image = image + i
    image = image+overlapp
    #image = image > 0
    assemble.append(image)


figsize(5,5)
convexhull00 = sk.morphology.convex_hull_image(assemble[0])
contour00 = mh.bwperim(convexhull00)
subplot(111, xticks=[],yticks=[])
imshow(2*contour00+assemble[0], interpolation = 'nearest')



def contourNegCorners(im, angle0=-130, angle1 = 0):
    contours = find_contours(im, 0.01)
    subpol = subdivide_polygon(contours[0], degree=5, preserve_ends=True)
    coord = approximate_polygon(subpol, tolerance=1.5)
    angles = angleBetweenSegments(coord)
    table = np.array([coord[0][0],coord[0][1],angles[0]])
    for i in range(len(angles[1:])):
        current =  np.array([coord2[i][0],coord2[i][1],angles[i]])
        table = np.vstack((table,current))
    neg = table[np.where(np.logical_and(table[:,2]>angle0, table[:,2]<angle1))]
    return neg

from mahotas.polygon import convexhull as mhcvxh
square = np.array([[1, 1, 1],[1, 1, 1], [1, 1, 1]])
print len(assemble)
n=1
contour = mh.bwperim(im > 0)
figsize(12,12)
chromosomes = {}
for im in assemble:
    subplot(2,3,n,xticks=[],yticks=[])
    cvxh = sk.morphology.convex_hull_image(im)
    ratio =100 * ((1.0*np.sum(im>0)) / (1.0*np.sum(cvxh > 0)))
    lab,_ = mh.label(im>0, Bc=square)
    prop = sk.measure.regionprops(lab,properties=['Eccentricity','ConvexArea','Area'])
    Ecc = prop[0]['Eccentricity']
    Area = prop[0]['Area']
    CVXArea = prop[0]['ConvexArea']
    #print Ecc, str(100.0*Area/CVXArea)[0:5]
    labim, nim = mh.label(im)
    negcorners = contourNegCorners(im>0)
    #
    cont = sk.measure.find_contours(im,0.01)
    pol = subdivide_polygon(cont[0], degree=5, preserve_ends=True)
    coord_pol = approximate_polygon(pol, tolerance=3)
    angles_pol = angleBetweenSegments(coord_pol)
    #print angles_pol

    table_im = np.array([coord_pol[0][0],coord_pol[0][1],angles_pol[0]])
    for i in range(len(angles_pol[1:])):
        current_im =  np.array([coord_pol[i][0],coord_pol[i][1],angles_pol[i]])
        table_im = np.vstack((table_im,current_im))
    neg_dom = table_im[np.where(np.logical_and(table_im[:,2]>-110, table_im[:,2]<0))]
   
    #Build a dict to filter assembled chromosomes
    chromosomes[ratio] = im    
    
    title('cc='+str(nim)+'  r='+str(ratio)[0:5]+' Ecc:'+str(100*Ecc)[0:4]+'  corner='+str(len(neg_dom)))
    imshow(im, interpolation = 'nearest',cmap=mycmap)
    convexhull = mhcvxh(im>0)
    plot(convexhull[:,1], convexhull[:,0],linewidth=1, color='w')
    plot(coord_pol[:, 1], coord_pol[:, 0], linewidth=2, color='c')
    if len(neg_dom)>0:
        scatter(neg_dom[:,1],neg_dom[:,0],s=100,c='white')
    #scatter(negcorners[:,1], negcorners[:,0], s=100, c='y')
    n =n +1

chromkeys = chromosomes.keys()
chromkeys.sort()
twoHighestRatio = chromkeys[-2:]
r1=twoHighestRatio[0]
r2=twoHighestRatio[1]
chrom1=chromosomes[r1]
chrom2=chromosomes[r2]
figsize(10,10)
subplot(221,xticks=[],yticks=[])
imshow(chrom1, interpolation='nearest')
subplot(222,xticks=[],yticks=[])
segmented1 = (chrom1>0)*original_im
imshow(segmented1, interpolation='nearest', cmap=pylab.cm.gray)
subplot(223,xticks=[],yticks=[])
imshow(chrom2, interpolation='nearest')
subplot(224,xticks=[],yticks=[])
segmented2 = (chrom2>0)*original_im
imshow(segmented2, interpolation='nearest', cmap=pylab.cm.gray)