
import numpy
import math
import cmath
from  matplotlib.pyplot import *
from PIL import Image
img = Image.open("../../../../simul/image/bords-kapla.jpg")
array=numpy.array(img)*1.0
figure(figsize=(8,6))
imshow(array,cmap=cm.gray)
            

(Ny,Nx) = array.shape
rho= 1.0
theta=1.0
Ntheta = int(180.0/theta)
Nrho = int(math.floor(math.sqrt(Nx*Nx+Ny*Ny))/rho)
dtheta = math.pi/Ntheta
drho = math.floor(math.sqrt(Nx*Nx+Ny*Ny))/Nrho
accum = numpy.zeros((Ntheta,Nrho))
            

for j in range(Ny):
    for i in range(Nx):
        if array[j][i]!=0:
            for i_theta in range(Ntheta):
                theta = i_theta*dtheta
                rho = i*math.cos(theta)+(Ny-j)*math.sin(theta)
                i_rho = int(rho/drho)
                if (i_rho>0) and (i_rho<Nrho):
                    accum[i_theta][i_rho] += 1
            

figure(figsize=(12,6))
imshow(accum,cmap=cm.gray)
            

seuil=130
accum_seuil = accum.copy()
for i_theta in range(Ntheta):
    for i_rho in range(Nrho):
        if accum[i_theta][i_rho]<seuil:
            accum_seuil[i_theta][i_rho] = 0
            
figure(figsize=(12,6))
imshow(accum_seuil,cmap=cm.gray)
            

lignes = []
for i_theta in range(Ntheta):
    for i_rho in range(Nrho):
        if accum_seuil[i_theta][i_rho]!=0:
            lignes.append((i_rho*drho,i_theta*dtheta))
            

figure(figsize=(8,6))
axis([0,Nx,0,Ny])
for rho,theta in lignes:
    a = math.cos(theta)
    b = math.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))
    plot([x1,x2],[y1,y2],color="b")

            

def coloriage_pixel_pile(image,result,shape,seuil,valeur,pile,i,j):
    image[j][i] = valeur
    result[j][i] = 255
    voisins = [(i+1,j),(i-1,j),(i,j-1),(i,j+1)]
    for pixel in voisins:
        (k,l) = pixel
        if k>=0 and k<shape[1] and l>=0 and l<shape[0]:
            if image[l][k]>seuil:
                image[l][k] = valeur
                pile.append(pixel)
                  

def analyse(image,seuil):
    shape = image.shape
    Nx = shape[1]
    Ny = shape[0]
    compteur = 0
    positions = []
    result = numpy.zeros((Ny,Nx),dtype=numpy.uint8)
    for y in range(Ny):
        for x in range(Nx):
            if image[y][x] > seuil:
                compteur += 1
                pile = [(x,y)]
                si = 0
                sj = 0
                npix = 0
                while len(pile)>0:
                    (i,j) = pile.pop()
                    si += i
                    sj += j
                    npix += 1
                    coloriage_pixel_pile(image,result,shape,seuil,0,pile,i,j)
                xb = si*1.0/npix
                yb = sj*1.0/npix
                positions.append((xb,yb,npix))
    print("Nombre de taches : %d"%compteur)
    positions = sorted(positions, key=lambda positions:-positions[2])
    for p in positions:
        print("x=%f, y=%f, npix=%d"%(p[0],p[1],p[2]))
    return (positions,result)
            

(positions,result) = analyse(accum,130)
shape = accum.shape
figure(figsize=(6,6))
for p in positions:
    plot([p[0]],[shape[0]-p[1]],'ob')
axis([0,shape[1],0,shape[0]])
            

figure(figsize=(8,6))
axis([0,Nx,0,Ny])
for point in positions:
    i_rho = point[0]
    i_theta = point[1]
    rho = i_rho*drho
    theta = i_theta*dtheta
    a = math.cos(theta)
    b = math.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))
    plot([x1,x2],[y1,y2],color="b")

            

figure(figsize=(8,6))
axis([0,Nx,0,Ny])
for point in positions[0:8]:
    i_rho = point[0]
    i_theta = point[1]
    rho = i_rho*drho
    theta = i_theta*dtheta
    a = math.cos(theta)
    b = math.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))
    plot([x1,x2],[y1,y2],color="b")

            
