segmentacion.py

# -*- coding: utf-8 -*-
"""
Created on Fri Aug 22 11:28:11 2014

@author: german
"""


from PyQt4.QtOpenGL import *
import numpy as np
import scipy.signal as signal
from drlse import *
import cv2
from skimage import morphology
from alineamiento import *
import copy
import shutil
from scipy.signal.signaltools import correlate2d as c2d

###################################################################################
###################################################################################
def gauss_kern():
    """ Returns a normalized 2D gauss kernel array for convolutions """
    h1 = 15
    h2 = 15
    x, y = np.mgrid[0:h2, 0:h1] 
    x = x-h2/2
    y = y-h1/2
    sigma = 1.5
    g = np.exp( -( x**2 + y**2 ) / (2*sigma**2) );
    return g / g.sum()
###################################################################################
def levelset(img,point):
    try:
        Img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    except:
        Img=img
    g = gauss_kern()
    Img_smooth = signal.convolve(Img,g,mode='same')
    Iy,Ix=np.gradient(Img_smooth)
    
    f=Ix**2+Iy**2
    g=1. / (1.+f)  # edge indicator function.
    epsilon=1.5 # the papramater in the definition of smoothed Dirac function
    timestep=5 # time step
    mu=0.2/timestep  # coefficient of the internal (penalizing) 
                      # energy term P(\phi)
                      # Note: the product timestep*mu must be less 
                      # than 0.25 for stability!
    
    lam=5 # coefficient of the weighted length term Lg(\phi)
    alf=-3 # coefficient of the weighted area term Ag(\phi);
          # Note: Choose a positive(negative) alf if the 
          # initial contour is outside(inside) the object.
    nrow, ncol=Img.shape
   
    c0=2
    a=6
    initialLSF=c0*np.ones((nrow,ncol))    
    initialLSF[point[1]-a:point[1]+a, point[0]-a:point[0]+a]=-c0    
    u=initialLSF 
    
    for n in range(500):    
        u=evolution(u, g ,lam, mu, alf, epsilon, timestep, 1)
    mask=np.uint8(u)
    ret,mask= cv2.threshold(mask,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
    contours,hierarchy = cv2.findContours(copy.deepcopy(mask),cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE) 
    largos=[]
    for i in contours:
        largos.append(len(i))
    cnt=largos.index(max(largos))
    cnt=contours[cnt]
    
    row, col=mask.shape
    for i in range(col):
        for j in range(row):
            test=cv2.pointPolygonTest(cnt,(i,j),False)
            if test<0:
                mask[j][i]=0
#            else:
#                mask[j][i]=255
#                
#    cv2.imshow('mask',mask)
#    cv2.waitKey(0)
    
    return mask
###################################################################################
def binary(img,umbral1=120,umbral2=0):
#    img=cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

    ret,img1=cv2.threshold(img,umbral1,255,cv2.THRESH_BINARY_INV)
    ret,img2=cv2.threshold(img,umbral2,255,cv2.THRESH_BINARY)
    
    img=cv2.bitwise_and(img1,img2)
    
    return img
###################################################################################
###################################################################################    
def Point(img, imgNext0,cnt,ptP,ptR,):
   
    if len(imgNext0.shape)>2:
        imgNext=cv2.cvtColor(imgNext0, cv2.COLOR_RGB2GRAY)
    else:
        imgNext=copy.deepcopy(imgNext0)
    
    
    difx=ptP[0]
    dify=ptP[1]
    
    imgcrop=crop(img,ptP[0],ptP[1],ptR[0],ptR[1])
#    skel, dist = morphology.medial_axis(imgcrop, mask=None, return_distance=True)
    
    
    r,c = imgcrop.shape[:2]
    dist=np.zeros(imgcrop.shape[:2],np.float32)   
    
    for i in range(r):
        for j in range(c):
            dist.itemset((i,j),cv2.pointPolygonTest(cnt,(j+difx,i+dify),True))
        
    maximo=np.amax(dist)
    positions=np.where(dist==maximo)
    pos=(int(positions[1][0]+difx),int(positions[0][0]+dify))
    posANT=copy.deepcopy(pos)
    print pos, imgNext[pos[1]][pos[0]],dist[positions[0][0]][positions[1][0]]
    print '---'
    
    continuar = False
    if (imgNext[pos[1]][pos[0]]>0)  or dist[positions[0][0]][positions[1][0]]<=0:
        print pos, imgNext[pos[1]][pos[0]],dist[positions[0][0]][positions[1][0]]
        continuar = True
    
        
    while (continuar==True):
        print '...'
        dist[positions[0][0]][positions[1][0]]=0
        maximo=np.amax(dist)
        positions=np.where(dist==maximo)
        pos=(int(positions[1][0]+difx),int(positions[0][0]+dify))
        print pos, imgNext[pos[1]][pos[0]],dist[positions[0][0]][positions[1][0]]
        if (imgNext[pos[1]][pos[0]]<=0 and dist[positions[0][0]][positions[1][0]]>0) :
            continuar = False
            print 'encontro pto'
            break
        if dist[positions[0][0]][positions[1][0]]<=0:
            pos=copy.deepcopy(posANT)
            continuar = False
            print 'NO encontro pto'
            break 
        
    return pos
###################################################################################
###################################################################################        
def copyMaskFiles(src,dst):    
    try:
        shutil.rmtree(dst)
    except:
        pass    
    try:            
        shutil.copytree(src,dst)                
    except:
        pass    
###################################################################################
###################################################################################    

def CorrIndex(img1,img2, centroid):
   if len(img1.shape)>2:
        im1=cv2.cvtColor(img1, cv2.COLOR_RGB2GRAY)
        im2=cv2.cvtColor(img2, cv2.COLOR_RGB2GRAY)
   else:
       im1=copy.deepcopy(img1)
       im2=copy.deepcopy(img2)
   a=50    
   im1=crop(im1,centroid[0]-a,centroid[1]-a,centroid[0]+a,centroid[1]+a)
   im2=crop(im2,centroid[0]-a,centroid[1]-a,centroid[0]+a,centroid[1]+a)
       
   im1=(im1 - im1.mean()) / im1.std()
   im2=(im2 - im2.mean()) / im2.std()
   
   corrIndex1=c2d(im1, im1, mode='same')
   corrIndex2=c2d(im1, im2, mode='same')
   
   return corrIndex1.max() , corrIndex2.max()