%matplotlib inline
import numpy as np
import cv2
import cv
import matplotlib.pyplot as plt
import os
from IPython.html.widgets import interact, interactive, fixed
from IPython.display import clear_output, display, HTML
from IPython.html import widgets
from collections import Counter
from matplotlib import pyplot as plt

from io import BytesIO
import PIL
from IPython.display import display, Image

def display_img_array(ima, cvt=cv2.COLOR_BGR2RGB, **kwargs):
    if cvt:
        ima = cv2.cvtColor(ima, cvt)
    im = PIL.Image.fromarray(ima)
    bio = BytesIO()
    im.save(bio, format='png')
    display(Image(bio.getvalue(), format='png', **kwargs))

%qtconsole --colors=linux

from random import randint
def sim(c1, c2, D=80):
    return sum(abs(int(c1[k])-int(c2[k])) for k in range(3))<=D
imgs = [cv2.imread('img/rail{}.jpg'.format(i)) for i in range(10)]
for img in imgs:
    display_img_array(img, width=600)
    z = iter(img.flatten())
    mcolor = np.array(Counter(zip(z,z,z)).most_common(1)[0][0], dtype=np.float32)
    h,w,c = img.shape
    for i in range(h):
        for j in range(w):
            if sim(mcolor, img[i,j],65):
                img[i,j]=(0,0,0)            
    fimg = img.reshape((-1,3)).astype(dtype=np.float32)
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
    K = 7
    ret,label,center=cv2.kmeans(fimg,K, criteria,10, cv2.KMEANS_RANDOM_CENTERS)
    """
    plt.scatter(fimg[:,0], fimg[:,1])
    plt.scatter(center[:,0], center[:,1],  s=80, c='r')
    plt.axis([0,256]*2)
    plt.show()
    """
    nimg = center.astype(np.uint8)[label.flatten()].reshape(img.shape)
    random_color = np.array([[randint(0,255) for j in range(3)] for i in range(K+1)], dtype=np.uint8)
    nimg2 = random_color[label.flatten()].reshape(img.shape)        
    #display_img_array(nimg, width=600)
    #display_img_array(nimg2, width=600)

    def remove_hline(img, img2, i,j):
        if img[i,j] in (img[i-1,j], img[i+1,j]):
            return        
        if img[i, j-1] == img[i, j+1] == img[i-1, j] == img[i+1,j]:
            img2[i,j]= img[i, j-1]
            return
        j1=j
        #while j1>=0 and img[i,j1] not in (0,K,img[i-1,j1],img[i+1,j1]):
        while j1>=0 and img[i,j1]!=img[i-1,j1] and img[i-1,j1] ==  img[i+1,j1]:
            j1-=1
        j1+=1
        j2 = j
        while j2<img.shape[1] and  img[i,j2]!=img[i-1,j2] and img[i-1,j2] == img[i+1,j2]:
            j2+=1
        if j2-j1 >=4:
            for j0 in range(j1,j2):
                img2[i,j0]=img[i-1,j0]
    icenter = map(tuple, center.astype(np.uint8))
    zero = icenter.index((0,0,0))
    print "zero", zero
    
    def remove_vline(img, img2,  i,j):
        if img[i,j] in (img[i,j-1], img[i,j+1]):
            return
        if img[i, j-1] == img[i, j+1] == img[i-1, j] == img[i+1,j]:
            img2[i,j]= img[i, j-1]
            return
        i1=i
        while i1>=0 and img[i1,j]!=img[i1,j-1] and img[i1,j-1] == img[i1,j+1]:
            i1-=1
        i1+=1
        i2 = i
        #while i2<img.shape[0] and  img[i2,j] not in (zero,K, img[i2,j-1], img[i2,j+1]):
        while i2<img.shape[0] and  img[i2,j]!=img[i2,j-1] and img[i2,j-1]==img[i2,j+1]:
            i2+=1
        if i2-i1 >=4:
            for i0 in range(i1,i2):
                img2[i0,j]=img[i0,j-1]
    limg = label.reshape(img.shape[:2])
    for iterate in range(2):
        limg2 = limg.copy()
        for i in range(1,h-1):
            for j in range(1,w-1):
                remove_vline(limg, limg2, i,j)
                remove_hline(limg, limg2, i, j)                
        limg = limg2
    bw=np.array([255 if i !=zero else 0 for i in range(K)], dtype=np.uint8)
    nimg3 = bw[limg2.flatten()].reshape(img.shape[:2])    
    display_img_array(nimg3, None, width=600)

def sim(c1, c2):
    return sum(abs(int(c1[k])-int(c2[k])) for k in range(3))<=80
        
def find_line(img,i,j):    
    L=5
    for j2 in range(j, j+4):
        if not sim(img[i,j2], img[i,j2+1]):
            return False
        if sim(img[i,j2], img[i-1,j2]) or sim(img[i,j2], img[i+1,j2]):
            return False
    return True

def find_vline(img,i,j):    
    L=5
    for i2 in range(i, i+4):
        if not sim(img[i2,j], img[i2+1,j]):
            return False
        if sim(img[i2,j], img[i2,j-1]) or sim(img[i2,j], img[i2,j+1]):
            return False
    return True
                    
imgs = [cv2.imread('img/rail{}.jpg'.format(i)) for i in range(10)]        
    
for img in imgs:
    display_img_array(img)
    z = iter(img.flatten())
    mcolor = np.array(Counter(zip(z,z,z)).most_common(1)[0][0], dtype=np.float32)
    h,w,c = img.shape
    nimg = np.zeros_like(img)
    for i in range(1,h-5):
        for j in range(1,w-5):
            nimg[i,j]=img[i,j]
            if (abs(img[i,j]-mcolor)).max() >30:
                if  find_line(img, i, j):
                    nimg[i,j]=(0,255,0)
                elif  find_vline(img, i, j):
                    nimg[i,j]=(0,255,0)
                else:
                    nimg[i,j]=img[i,j]
                    #nimg[i,j]=img[i-1,j]
                    #(img[i+1,j].astype(np.uint32)+img[i-1,j].astype(np.uint32))/2
       
    display_img_array(nimg)


def good_point(img,i,j):
    for k in range(3):
        v=[img[i,j,k], img[i+1,j+1,k], img[i+1,j-1,k], img[i-1,j-1,k], img[i-1,j+1,k]]
        v = map(int, v)
        if max(v)-min(v) > 30:
            return False
    return True

def good_nbhd(img, i0, j0):
    color = img[i0,j0]
    pts= {(i0,j0)}
    newpts = {(i0,j0)}
    while newpts:
        newpts2 = set()
        for pt in newpts:
            for i in range(pt[0]-2, pt[0]+3):
                if i<0 or i >= img.shape[0]:
                    continue
                for j in range(pt[1]-2, pt[1]+3):
                    if j < 0 or j >=img.shape[1]:
                        continue
                    if (i,j) in pts:
                        continue
                    color2 = img[i,j]
                    if max(abs(int(color[k])-int(color2[k])) for k in range(3)) <=15:
                        newpts2.add((i,j))
                        pts.add((i,j))
        newpts = newpts2
    return pts
                    
imgs = [cv2.imread('img/rail{}.jpg'.format(i)) for i in range(10)]        
    
for img in imgs[:1]:
    display_img_array(img)
    z = iter(img.flatten())
    mcolor = np.array(Counter(zip(z,z,z)).most_common(1)[0][0], dtype=np.float32)
    h,w,c = img.shape
    nimg = np.zeros_like(img)
    for pt in good_nbhd(img,28,12):
            nimg[pt[0],pt[1]]=(128,0,128)
    nimg[30,10]=(255,255,255)
    display_img_array(nimg)
    """
    for i in range(1,h-5):
        for j in range(1,w-5):
            if (abs(img[i,j]-mcolor)).max() >30:
                if True or good_rect(img, i, j):
                    #nimg[i,j]=(255,255,255) 
                    #print i,j, len(good_nbhd(img,i,j))
                    for pt in good_nbhd(img,i,j):
                        #if tuple(nimg[pt[0],pt[1]]) in [(0,0,0), (0,128,0)]:                        
                            nimg[pt[0],pt[1]]=(128,0,128)
                    break
                else:
                    if tuple(nimg[i,j]) == (0,0,0):
                        nimg[i,j]=(0,128,0)
    display_img_array(nimg)
"""    

imgs = [cv2.cvtColor(cv2.imread('img/rail{}.jpg'.format(i)), cv2.COLOR_BGR2GRAY ) for i in range(10)]
for img in imgs:
    display_img_array(img, None)

# Find key points
#orb = cv2.ORB()
sift = cv2.SIFT()
img1 = imgs[3]
img2 = imgs[5]
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)
k1 = cv2.drawKeypoints(img1, kp1)
display_img_array(k1)
k2 = cv2.drawKeypoints(img2, kp2)
display_img_array(k2)

# Match Key points
if 0: # Brute Force Match
    #bf = cv2.BFMatcher(cv2.NORM_HAMMING) #, crossCheck=True)
    bf = cv2.BFMatcher(cv2.NORM_L2) #, crossCheck=True)
    matches = bf.knnMatch(des1, des2, k=2)
else: #Flann Match
    index_params= dict(algorithm = 0, tree=5)
    search_params = dict(checks=50) 
    flann = cv2.FlannBasedMatcher(index_params, search_params)
    matches = flann.knnMatch(des1, des2, k=2)

# Draw Matches
H1, W1 = img1.shape[:2]
slide_h, slide_w = H2,W2 = img2.shape[:2]
#print img.shape, frame.shape
from random import randint
_img2 = cv2.resize(img2, ( 800, H2*800/W2))
NH1 = H1*800/W1
_img1 = cv2.resize(img1, (800, NH1))
outimg = np.concatenate((_img1, _img2), axis=0)
def draw_match_line(pt1, pt2):
    draw1(pt1)
    draw2(pt2)
    pt1 = pt1[0]*800/W1, pt1[1]*800/W1
    pt2 = pt2[0]*800/W2, pt2[1]*800/W2+NH1
    poly =  [pt1, pt2]
    cv2.polylines(outimg, np.int32([poly]), True, (randint(0,255), randint(0,255), randint(0,255)), 1)
def draw1(pt):
    pt = np.float32(pt)*800/W1
    pt = tuple(np.int32(pt))
    cv2.circle(outimg, pt, 3, (0,0,255))
def draw2(pt2):
    pt2 = int(pt2[0]*800/W2), int(pt2[1]*800/W2+NH1)
    cv2.circle(outimg, pt2, 3, (0,0,255))
    #cv2.polylines(outimg, np.int32([poly]), True, (0,255,255), 3)
m2 = []
for x in matches:
    #if len(x)==2 and x[0].distance> 0.6*x[1].distance:
    #    continue
    if len(x)==0:
        continue
    m2.append(x[0])
m2.sort(key=lambda x:x.distance)
for x in m2[:20]:
    #print x.queryIdx, x.trainIdx
    pt1 = kp1[x.queryIdx].pt
    pt2 = kp2[x.trainIdx].pt
    # print kp1[x.queryIdx].angle,  kp2[x.trainIdx].angle
    print pt1, pt2    
    draw_match_line(pt1,pt2)
display_img_array(outimg, cvt=None)

imgs = [cv2.cvtColor(cv2.imread('img/rail{}.jpg'.format(i)), cv2.COLOR_BGR2GRAY ) for i in range(10)]
for img in imgs:
    display_img_array(img, None)
    thresh = cv2.adaptiveThreshold(img,  255, 
                                                  cv2.ADAPTIVE_THRESH_MEAN_C, 
                                                  cv2.THRESH_BINARY, 
                                                  11, 
                                                  -5.0)
    display_img_array(thresh, None)
    ret,thresh = cv2.threshold(img, 110,255, cv2.THRESH_BINARY_INV)
    #display_img_array(thresh, None)
    #ret,thresh = cv2.threshold(thresh, 30,255, cv2.THRESH_TOZERO)
    contours, hierarchy = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
    im = np.zeros_like(img)
    cv2.drawContours(im,contours,-1,255,1)    
    #display_img_array(im, None)

%qtconsole --colors=linux

imgs = [cv2.imread('img/rail{}.jpg'.format(i)) for i in range(10)]
from collections import Counter
for img in imgs:
    z = iter(img.flatten())
    mcolor = np.array(Counter(zip(z,z,z)).most_common(1)[0][0], dtype=np.float32)-128.0
    mv = mcolor / np.linalg.norm(mcolor)
    print mv
    #nimg = abs((img - np.kron(np.dot(img, mv), mv).reshape(60,200,3))).astype(np.uint8)
    #nimg = cv2.cvtColor(nimg, cv2.COLOR_BGR2GRAY)
    #ret,thresh = cv2.threshold(nimg, 10,255, cv2.THRESH_BINARY_INV)
    h,w,c = img.shape
    nimg = np.zeros((h,w), dtype=np.uint8)
    print h,w,c,nimg.shape
    display_img_array(img)
    for i in range(h):
        for j in range(w):
            v = img[i,j].astype(np.float32)-128.0
            n = np.linalg.norm(v)    
            r = np.dot(v,mv)/n
            if all(v[i]>=mcolor[i]+5 for i in range(3)):
                nimg[i,j] = 255
            else:
                nimg[i,j]= 0
            #2000*(r-0.9)
            #0 if n ==0 else 200*r
    #nimg = np.array([[img[i,j]-np.dot(img[i,j],mv)*mv for j in range(w)] for i in range(h)], dtype=np.uint8)
    print nimg.min(), nimg.max()
    display_img_array(nimg, None)
    

    

from skimage import morphology
import skimage
imgs = [cv2.imread('img/rail{}.jpg'.format(i)) for i in range(20)]
for img in imgs:
    display_img_array(img, width=600)
    z = iter(img.flatten())
    mcolor = np.array(Counter(zip(z,z,z)).most_common(1)[0][0], dtype=np.float32)
    h,w,c = img.shape
    nimg=np.zeros((h,w,c), dtype=np.uint8)
    for i in range(h):
        for j in range(w):
            #nimg[i,j]=abs(img[i,j]-mcolor).astype(dtype=np.uint8)
            if sim(mcolor, img[i,j], 120):
                nimg[i,j]=(0,0,0)
            else:
                nimg[i,j]=img[i,j]
    display_img_array(nimg, width=600)
    z = iter(nimg[40:,:,:].flatten())
    result = Counter(zip(z,z,z)).most_common(2)
    print result
    mcolor2 = np.array(result[1][0], dtype=np.float32)
    nimg2=np.zeros((h,w,c), dtype=np.uint8)
    print mcolor2
    for i in range(h):
        for j in range(w):
            #r=abs(nimg[i,j]-mcolor2).astype(dtype=np.uint8)            
            if tuple(nimg[i,j])==(0,0,0) or sim(mcolor2, nimg[i,j], 40): # or sim(mcolor3, nimg[i,j], 40):
                nimg2[i,j]=(0,0,0)
            else:
                nimg2[i,j]=(255,255,255)
    display_img_array(nimg2, width=600)
    nimg2 = cv2.cvtColor(nimg2, cv2.COLOR_BGR2GRAY)
    nimg3=np.zeros((h,w), dtype=np.uint8)
       
    for i in range(5,h-5):
        for j in range(5,w-5):
            r = nimg2[i-2:i+3, j-2:j+3].sum()
            if r <= 255*2:
                nimg3[i,j]=0
            elif  r >= 255*17:
                nimg3[i,j]=255
            else:
                nimg3[i,j]=nimg2[i,j]
    def remove_hline(img, i,j):
        if img[i-1,j]  >=127 or img[i+1,j]>=127:
            return
        if img[i, j-1]<127 and img[i, j+1]<127:
            img[i,j]=0
            return
        j1=j
        while j1>=0 and img[i,j1]>=127 and img[i-1,j1]  <127 and img[i+1,j1] <127:
            j1-=1
        j1+=1
        j2 = j
        while j2<img.shape[1] and  img[i,j2]>=127 and img[i-1,j2]  <127 and img[i+1,j2] <127:
            j2+=1
        if j2-j1 >=4:
            for j0 in range(j1,j2):
                img[i,j0]=0

    def remove_vline(img, i,j):
        if img[i,j-1]  >=127 or img[i,j+1]>=127:
            return
        if img[i-1, j]<127 and img[i+1, j]<127:
            img[i,j]=0
            return
        i1=i
        while i1>=0 and img[i1,j]>=127 and img[i1,j-1]  <127 and img[i1,j+1] <127:
            i1-=1
        i1+=1
        i2 = i
        while i2<img.shape[0] and  img[i2,j]>=127 and img[i2,j-1]  <127 and img[i2,j+1] <127:
            i2+=1
        if i2-i1 >=4:
            for i0 in range(i1,i2):
                img[i0,j]=0
    for i in range(2, h-2):
        for j in range(2,w-2):
            remove_hline(nimg3, i,j)
            remove_vline(nimg3, i,j)
    ximg3= img.copy() #cv2.cvtColor(nimg2, cv2.COLOR_GRAY2BGR)
    itv = None
    result =[]
    for j in range(2,w-2):
        if nimg3[:,j].max() ==0:
            cv2.polylines(ximg3, np.int32([[[j,0],[j,h]]]), False, (0,0,128))
            if itv:
                result.append(itv)
                itv = None
        else:
            if itv is None:
                itv=(j,j)
            else:
                itv = (itv[0], j)
    if itv is not None:
        result.append(itv)
    last = None
    result2=[]
    for itv in result:        
        if itv[1]-itv[0]>=16:
            if last and last[1]-last[0]>=9:
                result2.append(last)
            if itv[1]-itv[0]>=28:
                mid = int((itv[0]+itv[1])/2)
                result2.append( (itv[0],mid) )
                result2.append( (mid,itv[1]) )
            else:
                result2.append(itv)
            last = None
        elif last and last[1]+2 >= itv[0] and 30>itv[1]-last[0]>=10:
            result2.append((last[0], itv[1]))
            last = None
        elif itv[1]-itv[0]>=1:
            if last and last[1]-last[0] >=9:
                result2.append(last)                
            last = itv
        else:
            if last and last[1]-last[0] >=9:
                result2.append(last)     
            last = None   
    if last and last[1]-last[0] >=9:
        result2.append(last)   
    print result
    print result2
    for itv in result2:
        for i in range(h):
            if nimg3[i,itv[0]:itv[1]+1].max() ==0:
                cv2.polylines(ximg3, np.int32([[[itv[0],i],[itv[1]+1,i]]]), False, (0,0,128))
    
    display_img_array(ximg3, width=600)
    nimg4 = np.zeros_like(nimg3)
    for i in range(3, h-3):
        for j in range(3,w-3):
            s = nimg3[i-2:i+3, j-1:j+2].astype(int).sum()
            if s > 255*2:
                nimg4[i,j]=128
                if nimg3[i,j]:
                    nimg4[i,j]=255
    display_img_array(nimg4, None, width=600)
    #L = morphology.label(nimg4)
    L=skimage.measure.label(nimg4)
    print "Number of components:", np.max(L)
    print Counter(L.flatten()).most_common()
    contours, hierarchy = cv2.findContours(nimg4,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
    im = np.zeros_like(nimg4)
    cv2.drawContours(im,contours,-1,255,1)    
    display_img_array(im, None, width=600)
    


def sim(c1, c2, D=80):
    return sum(abs(int(c1[k])-int(c2[k])) for k in range(3))<=D
imgs = [cv2.imread('img/rail{}.jpg'.format(i)) for i in range(10)]
for img in imgs:
    display_img_array(img,  width=800)
    z = iter(img.flatten())
    mcolor = np.array(Counter(zip(z,z,z)).most_common(1)[0][0], dtype=np.float32)
    h,w,c = img.shape
    nimg=np.zeros((h,w,c), dtype=np.uint8)
    for i in range(h-5):
        for j in range(w-5):
            if sim(mcolor, img[i,j],65):
                nimg[i,j]=(0,0,0)
            else:
                nimg[i,j]=(255,255,255)
    #display_img_array(nimg)
    img = nimg
    nimg=cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    
    def remove_hline(img, i,j):
        if img[i-1,j]  >=127 or img[i+1,j]>=127:
            return
        if img[i, j-1]<127 and img[i, j+1]<127:
            img[i,j]=0
            return
        j1=j
        while j1>=0 and img[i,j1]>=127 and img[i-1,j1]  <127 and img[i+1,j1] <127:
            j1-=1
        j1+=1
        j2 = j
        while j2<img.shape[1] and  img[i,j2]>=127 and img[i-1,j2]  <127 and img[i+1,j2] <127:
            j2+=1
        if j2-j1 >=4:
            for j0 in range(j1,j2):
                img[i,j0]=0

    def remove_vline(img, i,j):
        if img[i,j-1]  >=127 or img[i,j+1]>=127:
            return
        if img[i-1, j]<127 and img[i+1, j]<127:
            img[i,j]=0
            return
        i1=i
        while i1>=0 and img[i1,j]>=127 and img[i1,j-1]  <127 and img[i1,j+1] <127:
            i1-=1
        i1+=1
        i2 = i
        while i2<img.shape[0] and  img[i2,j]>=127 and img[i2,j-1]  <127 and img[i2,j+1] <127:
            i2+=1
        if i2-i1 >=4:
            for i0 in range(i1,i2):
                img[i0,j]=0
                
    for i in range(2,h-2):
        for j in range(2,w-2):
            if nimg[i,j]==0:
                continue
            remove_hline(nimg, i,j)          
            remove_vline(nimg,i,j)
    display_img_array(nimg, None, width=800)

import math
imgs = [cv2.imread('img/rail{}.jpg'.format(i)) for i in range(10)]
gimgs= [cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) for img in imgs]
for img, img0 in zip(gimgs, imgs):
    cimg = cv2.Canny(img, 50, 200,1)
    display_img_array(img, None)
    display_img_array(cimg, None)
    lines = cv2.HoughLinesP(cimg, 1, math.pi/180, 20,3,3)
    for l in lines[0]: 
        cv2.polylines(img0, np.int32([[l[:2],l[2:]]]), False, (0,0,255), 1)
    display_img_array(img0)

from skimage import morphology
import skimage
import skimage.feature
TP = [cv2.cvtColor(cv2.imread("data/{}a.png".format(fn)), cv2.COLOR_BGR2GRAY) for fn in [6,8,9]]
imgs = [cv2.imread('img/rail{}.jpg'.format(i)) for i in range(20)]
for img in imgs:
    #display_img_array(img, width=600)
    z = iter(img.flatten())
    mcolor = np.array(Counter(zip(z,z,z)).most_common(1)[0][0], dtype=np.float32)
    h,w,c = img.shape
    nimg=np.zeros((h,w,c), dtype=np.uint8)
    for i in range(h):
        for j in range(w):
            #nimg[i,j]=abs(img[i,j]-mcolor).astype(dtype=np.uint8)
            if sim(mcolor, img[i,j], 120):
                nimg[i,j]=(0,0,0)
            else:
                nimg[i,j]=img[i,j]
    # display_img_array(nimg, width=600)
    z = iter(nimg[40:,:,:].flatten())
    result = Counter(zip(z,z,z)).most_common(2)
    print result
    mcolor2 = np.array(result[1][0], dtype=np.float32)
    nimg2=np.zeros((h,w,c), dtype=np.uint8)    
    for i in range(h):
        for j in range(w):
            #r=abs(nimg[i,j]-mcolor2).astype(dtype=np.uint8)            
            if tuple(nimg[i,j])==(0,0,0) or sim(mcolor2, nimg[i,j], 40): # or sim(mcolor3, nimg[i,j], 40):
                nimg2[i,j]=(0,0,0)
            else:
                nimg2[i,j]=(255,255,255)
    #display_img_array(nimg2, width=600)
    nimg2 = cv2.cvtColor(nimg2, cv2.COLOR_BGR2GRAY)

    nimg3=np.zeros((h,w), dtype=np.uint8)
       
    for i in range(5,h-5):
        for j in range(5,w-5):
            r = nimg2[i-2:i+3, j-2:j+3].sum()
            if r <= 255*2:
                nimg3[i,j]=0
            elif  r >= 255*17:
                nimg3[i,j]=255
            else:
                nimg3[i,j]=nimg2[i,j]
    def remove_hline(img, i,j):
        if img[i-1,j]  >=127 or img[i+1,j]>=127:
            return
        if img[i, j-1]<127 and img[i, j+1]<127:
            img[i,j]=0
            return
        j1=j
        while j1>=0 and img[i,j1]>=127 and img[i-1,j1]  <127 and img[i+1,j1] <127:
            j1-=1
        j1+=1
        j2 = j
        while j2<img.shape[1] and  img[i,j2]>=127 and img[i-1,j2]  <127 and img[i+1,j2] <127:
            j2+=1
        if j2-j1 >=4:
            for j0 in range(j1,j2):
                img[i,j0]=0

    def remove_vline(img, i,j):
        if img[i,j-1]  >=127 or img[i,j+1]>=127:
            return
        if img[i-1, j]<127 and img[i+1, j]<127:
            img[i,j]=0
            return
        i1=i
        while i1>=0 and img[i1,j]>=127 and img[i1,j-1]  <127 and img[i1,j+1] <127:
            i1-=1
        i1+=1
        i2 = i
        while i2<img.shape[0] and  img[i2,j]>=127 and img[i2,j-1]  <127 and img[i2,j+1] <127:
            i2+=1
        if i2-i1 >=4:
            for i0 in range(i1,i2):
                img[i0,j]=0
    for i in range(2, h-2):
        for j in range(2,w-2):
            remove_hline(nimg3, i,j)
            remove_vline(nimg3, i,j)
    for i, t1 in enumerate(TP):
        res = cv2.matchTemplate(nimg3, t1, cv2.TM_CCORR_NORMED)
        h,w=t1.shape
        #min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(res)    
        #tl = np.int32(max_loc)
        #print "%d: "%i, tl, max_val
        loc =  skimage.feature.peak_local_max(res, 5, 0.75)
        for tl in loc:
            tl=(tl[1], tl[0])
            cv2.rectangle(img, tuple(tl), (tl[0]+w, tl[1]+h), [(0,0,255), (0,255,0), (255,0,0)][i] )   
        
    display_img_array(img, width=600)