from scipy import ndimage
import math
import matplotlib.pyplot as plt
from skimage import filters
import numpy as np
from skimage import feature

img = ndimage.imread("/home/tpolgrabia/Pobrane/CSFID/tracks_cropped/00003.jpg")
otsu_lvl = filters.threshold_otsu(img)
img_mask = img >= otsu_lvl
img_mask_fil = True ^ ndimage.minimum_filter(img_mask, 3)
img_label, nb_labels = ndimage.label(img_mask_fil)
sizes = ndimage.sum(img_mask_fil, img_label, range(0, nb_labels+1))
remove_labels = sizes < 100
lbls = True & remove_labels
lbls[0] = False

print("Remove labels: {}".format(ndimage.sum(remove_labels)))

remove_mask = remove_labels[img_label]
img_label[remove_mask] = 0

idx = 0
n2 = nb_labels + 1 - np.sum(remove_labels)
centers = np.zeros((nb_labels+1,2))
(n,m) = img.shape
corners = np.zeros((n,m), dtype=np.bool)

for i in range(0, nb_labels+1):
    if remove_labels[i]:
        continue

    seg = img_label == i
    slice_y, slice_x = ndimage.find_objects(seg)[0]
    seg_window = seg[slice_y, slice_x]
    edge = feature.canny(seg_window)
    plt.imshow(edge, cmap="Greys_r")
    plt.show()
    (en,em) = edge.shape
    distances = np.zeros((en,em))

    centers[idx,:] = ndimage.center_of_mass(seg)
    for j in range(0, en):
        for k in range(0, em):
            if edge[j,k]:
                diff = centers[i,:] - np.array([slice_y.start, slice_x.start])
                dx = k - diff[1]
                dy = j - diff[0]
                d = math.sqrt(dx*dx+dy*dy)
                distances[j,k] = d
 #   max_distances = ndimage.maximum_filter(distances, 20)
 #   t_corners = np.abs(max_distances - distances) < 1e-6
 #   print(t_corners)
    # plt.imshow(t_corners, cmap="Greys_r")
    # plt.show()
#    corners[slice_y,slice_x] = corners[slice_y,slice_x] + edge

    idx += 1

(n,m) = img.shape

mx = np.tile(np.arange(0,m),n).reshape((n,m))
my = np.tile(np.arange(0,n),m).reshape((m,n)).T
xs2 = mx[corners]
ys2 = my[corners]

xs = centers[:,1]
ys = centers[:,0]

plt.imshow(img_label, cmap="Greys_r")
# plt.plot(xs, ys, "r+")
plt.plot(xs2,ys2,"b+")
plt.show()