import scipy.io as sio
from scipy import interpolate
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from copy import deepcopy

mapa = sio.loadmat("data_map.mat")

mapa1 = mapa['data_map']

x = mapa1[:, 0]
y = mapa1[:, 1]
z = mapa1[:, 2]
z1 = [0. for i in range(len(x))]


# [X, Y] = np.meshgrid(x, y)
# [Z, Z] = np.meshgrid(z1, z1)
#
# for i in range(len(x)):
#     Z[i][i] = float(z[i])
#
# for j in range(0, len(x)-1):
#     for i in range(0, len(x)-j-1):
#         Z[i][j+i+1] = (Z[i][j+i] + Z[i+1][j+i+1])/2
#         Z[j + i + 1][i] = (Z[i][j + i] + Z[i + 1][j + i + 1]) / 2

x1 = np.arange(np.round(min(x), 3), max(x)+0.005, 0.001)
y1 = np.arange(min(y)-0.001, max(y)+0.001, 0.001)
z2 = [0. for i in range(len(x1))]

startx = np.round(min(x), 3)
starty = min(y)-0.001

[X1, Y1] = np.meshgrid(x1, y1)
[Z1, Z2] = np.meshgrid(z2, z2)

for i in range(len(x)):
    iks = (np.round(x[i],3) - startx)*1000
    iks1 = int(round(iks))
    igrek = (np.round(y[i], 3) - starty) * 1000
    igrek1 = int(round(igrek))
    if Z1[igrek1][iks1] == 0:
        Z1[igrek1][iks1] = z[i]
        Z2[igrek1][iks1] += 1
    else:
        Z1[igrek1][iks1] = (Z1[igrek1][iks1]*Z2[igrek1][iks1] + z[i])/(Z2[igrek1][iks1]+1)
        Z2[igrek1][iks1] += 1

iter = 0
przerwa = 6
zrobiono = False

def znajdz_k(lista):
    index = -1
    for i in range(len(lista)):
        if lista[i] != 0:
            index = i
            break
    return index

def znajdz_r(lista):
    index = -1
    for i in range(len(lista)):
        if lista[len(lista)-i-1] != 0:
            index = len(lista)-i-1
            break
    return index


while przerwa < 300:
    operacje = 0
    if iter % 2 == 0:
        for j in range(len(x1)):
            kon = False
            idx1 = znajdz_k(Z1[0:len(x1), j])
            if idx1 != -1:
                idx2 = znajdz_k(Z1[idx1+1:len(x1), j]) + idx1+1
            else:
                kon = True
            while not(kon):
                if idx1 == -1 or idx2 == -1 or idx1 == idx2:
                    kon = True
                elif np.abs(idx2-idx1) > przerwa+1 or np.abs(idx2-idx1) == 1:
                    idx1 = deepcopy(idx2)
                    idx2 = znajdz_k(Z1[idx1 + 1:len(x1), j]) + idx1 + 1
                else:
                    roznica = (Z1[idx2][j] - Z1[idx1][j])/(idx2 - idx1)
                    for k in range(1, idx2 - idx1):
                        Z1[k+idx1][j] = Z1[idx1+k-1][j] + roznica
                    operacje += 1
                    idx1 = deepcopy(idx2)
                    idx2 = znajdz_k(Z1[idx1+1:len(x1), j]) + idx1+1
    else:
        for j in range(len(x1)):
            kon = False
            idx1 = znajdz_k(Z1[j, 0:len(x1)])
            if idx1 != -1:
                idx2 = znajdz_k(Z1[j, idx1+1:len(x1)]) + idx1+1
            else:
                kon = True
            while not(kon):
                if idx1 == -1 or idx2 == -1 or idx1 == idx2:
                    kon = True
                elif np.abs(idx2-idx1) > przerwa+1 or np.abs(idx2-idx1) == 1:
                    idx1 = deepcopy(idx2)
                    idx2 = znajdz_k(Z1[j, idx1 + 1:len(x1)]) + idx1 + 1
                else:
                    roznica = (Z1[j][idx2] - Z1[j][idx1])/(idx2 - idx1)
                    for k in range(1, idx2 - idx1):
                        Z1[j][k+idx1] = Z1[j][idx1+k-1] + roznica
                    operacje += 1
                    idx1 = deepcopy(idx2)
                    idx2 = znajdz_k(Z1[j, idx1+1:len(x1)]) + idx1+1
    iter += 1
    if operacje == 0:
        if przerwa < 48:
            przerwa += 6
        else:
            przerwa += 48

interx = []
intery = []


operacje = 1


while operacje > 0:
    operacje = 0
    for j in range(1, len(x1)-1):
        for i in range(1, len(x1)-1):
            if Z1[i][j] == 0:
                if Z1[i][j-1] != 0 and Z1[i-1][j] != 0 :
                    Z1[i][j] = (Z1[i][j-1] + Z1[i-1][j])/2
                    operacje +=1
                elif Z1[i][j-1] != 0 and Z1[i+1][j] != 0 :
                    Z1[i][j] = (Z1[i][j - 1] + Z1[i + 1][j]) / 2
                    operacje += 1
                elif Z1[i][j + 1] != 0 and Z1[i + 1][j] != 0:
                    Z1[i][j] = (Z1[i][j + 1] + Z1[i + 1][j]) / 2
                    operacje += 1
                elif Z1[i][j + 1] != 0 and Z1[i - 1][j] != 0:
                    Z1[i][j] = (Z1[i][j + 1] + Z1[i - 1][j]) / 2
                    operacje += 1

for j in range(1):
    for i in range(1, len(x1)-1):
        if Z1[i][j] == 0:
            Z1[i][j] = Z1[i][j+1]

for j in range(417,len(x1)):
    for i in range(1, len(x1)-1):
        if Z1[i][j] == 0:
            Z1[i][j] = Z1[i][j-1]

for i in range(len(x1)):
    if Z1[0][i] == 0:
        Z1[0][i] = Z1[1][i]
    if Z1[len(x1)-1][i] == 0:
        Z1[len(x1) - 1][i] = Z1[len(x1)-2][i]

print(min(z)," ", max(z))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z)

ax.plot_surface(X1, Y1, Z1,cmap=cm.coolwarm,
                       linewidth=0, antialiased=True)

plt.show()