AgParameterisationStreamline.py

from pylab import *
import struct
import time
#import wavetools as wv
from scipy import interpolate,signal,ndimage
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import norm
import pylab
from IPython.display import display, Math, Latex
from scipy.ndimage import gaussian_filter1d, gaussian_filter, median_filter
import random
import scipy.ndimage
from mpl_toolkits import mplot3d
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from matplotlib import ticker
import math
import plotly.graph_objects as go
import pandas as pd
import plotly.io as pio
#from lmfit import Parameters,minimize, fit_report
from scipy.spatial.transform import Rotation as R
import math as m
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.path import Path
import matplotlib.patches as patches
from skimage import measure
from sklearn.neighbors import KernelDensity
from scipy.interpolate import RegularGridInterpolator
#from mayavi import mlab
import plotly.graph_objects as go
from scipy.spatial import Voronoi
#import pyvista as pv
import matplotlib.colors as mcolors
import copy
import glfw
from glfw.GLFW import *
import OpenGL
import open3d as o3d
import pdb
import seaborn as sns
from PIL import Image
import cv2 
import skimage
import skimage.exposure as exposure
import torch
import torch.nn as nn
from skimage.io import imread
from skimage.filters import threshold_otsu
from skimage import color
from skimage.restoration import denoise_bilateral
from skimage.segmentation import slic
from skimage.color import label2rgb
import re
from scipy.spatial.transform import Rotation as R
from scipy.spatial import KDTree
from scipy.interpolate import griddata
import scipy.fftpack as fft
from sklearn.neighbors import KDTree
from scipy.spatial import distance
from scipy.spatial import ConvexHull
from vtk.util.numpy_support import numpy_to_vtk, vtk_to_numpy
import vtk
#%%

def FS(k_x,k_y,k_z,a,C_0,C_200,C_211,C_220,C_310,C_222,C_321):
    Calc=(3
          -(np.cos(0.5*a*k_y)*np.cos(0.5*a*k_z)+np.cos(0.5*a*k_z)*np.cos(0.5*a*k_x)+np.cos(0.5*a*k_x)*np.cos(0.5*a*k_y))
          +C_200*(3-(np.cos(a*k_y)+np.cos(a*k_z)+np.cos(a*k_x)))
          +C_211*(3-((np.cos(a*k_x)*np.cos(0.5*a*k_y)*np.cos(0.5*a*k_z))+(np.cos(a*k_y)*np.cos(0.5*a*k_z)*np.cos(0.5*a*k_x))+(np.cos(a*k_z)*np.cos(0.5*a*k_x)*np.cos(0.5*a*k_y))))
          +C_220*(3-((np.cos(a*k_y)*np.cos(a*k_z))+(np.cos(a*k_z)*np.cos(a*k_x))+(np.cos(a*k_x)*np.cos(a*k_y))))
          +C_310*(6-((np.cos(1.5*a*k_x)*np.cos(0.5*a*k_y))+(np.cos(1.5*a*k_y)*np.cos(0.5*a*k_z))+(np.cos(1.5*a*k_z)*np.cos(0.5*a*k_x)))
                  -((np.cos(1.5*a*k_y)*np.cos(0.5*a*k_x))+(np.cos(1.5*a*k_z)*np.cos(0.5*a*k_y))+(np.cos(1.5*a*k_x)*np.cos(0.5*a*k_z))))
          +C_222*(1-(np.cos(a*k_x)*np.cos(a*k_y)*np.cos(a*k_z)))
          +C_321*(6-((np.cos(1.5*a*k_x)*np.cos(a*k_y)*np.cos(0.5*a*k_z))+(np.cos(1.5*a*k_y)*np.cos(a*k_z)*np.cos(0.5*a*k_x))+(np.cos(1.5*a*k_z)*np.cos(a*k_x)*np.cos(0.5*a*k_y)))
                  -((np.cos(1.5*a*k_z)*np.cos(a*k_y)*np.cos(0.5*a*k_x))+(np.cos(1.5*a*k_x)*np.cos(a*k_z)*np.cos(0.5*a*k_y))+(np.cos(1.5*a*k_y)*np.cos(a*k_x)*np.cos(0.5*a*k_z))))
          )-C_0
    return Calc

def cart2lat(k_x,k_y,k_z,revec):
    npoints=size(k_x)
    cartk=np.zeros([npoints,3])
    cartk[:,0]=k_x
    cartk[:,1]=k_y
    cartk[:,2]=k_z
    latk=np.zeros([npoints,3])
    if np.linalg.det(revec)==0:
        print('inverse revec not exist')
    else:
        #ndim=len(cartk[0])
        inv_revec=np.linalg.inv(revec)
        latk[:,0:3]=(np.matmul(revec,cartk[:,0:3].T)).T
        latk[:,0:3]=np.round(latk[:,0:3],8)
        #latk[:,4:ndim]=cartk[:,4:ndim]
    return latk

def lat2cart(latkp,revec): #%latkpt should be n *3 matrix
    eps=10**-10
    v11=len(latkp)
    v22=size(latkp[0])
    latk=np.zeros([v22,v11])
    latk[:,0]=latkp[0]
    latk[:,1]=latkp[1]
    latk[:,2]=latkp[2]
    cartkp=np.zeros([v22,v11])
    cartkp[:,0:3]=np.matmul(revec,latk[:,0:3].T).T
    if v22>=4:
        cartkp[:,4:v22]=latkp[:,4:v22]
    cartkp=np.round(cartkp,12)
    return cartkp

       
# def Plane(plane,npoints,dim):
#     k_x=
#     k_y=
#     k_z=

def PlaneR(a,npoints,plane,order):
    RefPlane=[1,1,0]
    Rm=np.zeros([3,3])
    phi=[]
    theta=[]
    psi=[]
    
    k_x = k_y = k_z =np.linspace(0,2*np.pi*a,npoints)
    
    #Calculate Rotation Matrix going from RefPlane (n) -> plane (b)
    #Normalize Miller Indices to get unit vector
    n=RefPlane/np.sqrt(sum(np.square(RefPlane)))
    b=plane/np.sqrt(sum(np.square(plane)))
    t=np.cross(n,b)/np.sqrt(sum(np.square(np.cross(n,b))))
    #Create Rotation Matrix
    Rm[:,0]=b
    Rm[:,1]=t
    Rm[:,2]=n
    #Calculate Euler Angles 
    if order == 'zxz':
        # phi = (np.arctan(-Rm[0,2] / Rm[1,2])) #*180/np.pi
        # theta = (np.arctan(-Rm[1,2] / (Rm[2,2] *np.cos(phi)))) #*180/np.pi
        # psi = (np.arctan(Rm[2,0] / Rm[2,1])) #*180/np.pi
       phi = (np.arccos(Rm[2,2]))*180/np.pi
       theta = (np.arctan(Rm[0,2] /-Rm[1,2]))*180/np.pi
       psi = (np.arctan(Rm[2,0]/Rm[2,1]))*180/np.pi
    return Rm


def Rx(theta):
  return np.matrix([[ 1, 0           , 0           ],
                   [ 0, m.cos(theta),-m.sin(theta)],
                   [ 0, m.sin(theta), m.cos(theta)]])
  
def Ry(theta):
  return np.matrix([[ m.cos(theta), 0, m.sin(theta)],
                   [ 0           , 1, 0           ],
                   [-m.sin(theta), 0, m.cos(theta)]])
  
def Rz(theta):
  return np.matrix([[ m.cos(theta), -m.sin(theta), 0 ],
                   [ m.sin(theta), m.cos(theta) , 0 ],
                   [ 0           , 0            , 1 ]])


def lod_mesh_export(mesh, lods, extension):
    mesh_lods={}
    for i in lods:
        mesh_lod = mesh.simplify_quadric_decimation(i)
        o3d.io.write_triangle_mesh("lod_"+str(i)+extension, mesh_lod)
        mesh_lods[i]=mesh_lod
    print("generation of "+str(i)+" LoD successful")
    return mesh_lods

def ChangeSampleDensity(array,whichtype,new_length):    #1 is defined length, 2 is defined density change
    if shape(array)[1] > 3:
        array=transpose(np.array(array))
    
    nx=np.arange(0, array.shape[0])   
    fit = interpolate.interp1d(nx, array, axis=0)
    new=[]
    #array[:,0]=np.linspace(np.min(array[:,0]),np.max(array[:,0]), len(array))
    if whichtype == 1:
        new = fit(np.linspace(0, array.shape[0]-1, new_length))
    if whichtype == 2:
        new = fit(np.linspace(0, array.shape[0]-1, new_length*len(array)))
    return new

def zdata(xydata,Rmatrix):
    zcol=[]
    xyzdata=[]
    for i in range(len(xydata)):
        zcol.append(np.matmul(Rmatrix[2,0:2],transpose(np.linalg.inv(Rmatrix[0:2,0:2]).dot(xydata[i,0:2])))[0,0])
    xyzdata=np.insert(xydata, 2, np.array([zcol]), axis=1)
    return xyzdata

# def symnorm(contour):
#     contour=np.round(countour,4)
#     for i in range(len(contour)):
#         indexnumber=np.where((np.round(p[:,0],4)==-np.round(p[0,0],4)) & (np.round(p[:,1],4)==-np.round(p[0,1],4)))[0][0]
        
def BZCalc(cell,lc):
    """
    Generate the Brillouin Zone of a given cell. The BZ is the Wigner-Seitz cell
    of the reciprocal lattice, which can be constructed by Voronoi decomposition
    to the reciprocal lattice.  A Voronoi diagram is a subdivision of the space
    into the nearest neighborhoods of a given set of points. 

    """
    cell= np.linalg.inv(cell).T
    cell=cell*lc  #0.81293684
    cell = np.asarray(cell, dtype=float)
    
    assert cell.shape == (3, 3)

    px, py, pz = np.tensordot(cell, np.mgrid[-1:2, -1:2, -1:2], axes=[0, 0])
    points = np.c_[px.ravel(), py.ravel(), pz.ravel()]

    vor = Voronoi(points)

    bz_facets = []
    bz_ridges = []
    bz_vertices = []

    # for rid in vor.ridge_vertices:
    #     if( np.all(np.array(rid) >= 0) ):
    #         bz_ridges.append(vor.vertices[np.r_[rid, [rid[0]]]])
    #         bz_facets.append(vor.vertices[rid])

    for pid, rid in zip(vor.ridge_points, vor.ridge_vertices):
        # WHY 13 ????
        # The Voronoi ridges/facets are perpendicular to the lines drawn between the
        # input points. The 14th input point is [0, 0, 0].
        if(pid[0] == 13 or pid[1] == 13):
            bz_ridges.append(vor.vertices[np.r_[rid, [rid[0]]]])
            bz_facets.append(vor.vertices[rid])
            bz_vertices += rid

    bz_vertices = list(set(bz_vertices))

    return vor.vertices[bz_vertices], bz_ridges, bz_facets

def BZPlanes(rlv, version): #version=0 you didnt give the reciprocal lattice vector
    if version == 0:
        rlv=reclatvec(rlv)
        print('You gave the unit cell, reciprocal lattice vector calculated')
    else: print('You gave the reciprocal lattice vector')
    
    #G=np.array([np.sum(rlv,axis=1)])
    G=np.zeros([1,3])
    b1=rlv[0,:]
    b2=rlv[1,:]
    b3=rlv[2,:]
    check=np.zeros([((3**3)-1),3])
    counter=0
    planes=[]
    
    for i in range(-1,2): 
        for j in range(-1,2): 
            for k in range(-1,2):
                if i == j == k == 0:
                    continue
                G[0,0]=b1[0]*i+b1[1]*j+b1[2]*k
                G[0,1]=b2[0]*i+b2[1]*j+b2[2]*k
                G[0,2]=b3[0]*i+b3[1]*j+b3[2]*k
                G2=np.array((G[0,0],G[0,1],G[0,2]))
                
                #print(G2)
                
                if (m.dist([0,0,0],G2)/2) <= abs(b1[2]+b2[2]+b3[2]):
                    planes.append(G2/2)
                
                check[counter,:]=G2
                counter=counter+1
    
    planes=np.array(planes)
            
        
    return planes

def BZCheck(rlv,x,y,z,version): #version=0 you didnt give the reciprocal lattice vector
    if version == 0:
        rlv=reclatvec(rlv)
        print('You gave the unit cell, reciprocal lattice vector calculated')
    else: print('You gave the reciprocal lattice vector')
    
    #G=np.array([np.sum(rlv,axis=1)])
    G=np.zeros([1,3])
    b1=rlv[0,:]
    b2=rlv[1,:]
    b3=rlv[2,:]
    check=np.zeros([((3**3)-1),3])
    counter=0
    
    for i in range(-1,2): 
        for j in range(-1,2): 
            for k in range(-1,2):
                if i == j == k == 0:
                    continue
                G[0,0]=b1[0]*i+b1[1]*j+b1[2]*k
                G[0,1]=b2[0]*i+b2[1]*j+b2[2]*k
                G[0,2]=b3[0]*i+b3[1]*j+b3[2]*k
                G2=np.array([(G[0,0],G[0,1],G[0,2])])
                check[counter,:]=G2
                counter=counter+1

    rtc=[]    
    
    for i in range(len(x)):
        test=np.array([x[i],y[i],z[i]])
        for j in range(len(check)):
            rangetest=np.matmul(test,(check[j,:]))
            rtc=np.dot(-check[j,:],check[j,:])/2            
            if abs(rangetest) > abs(rtc):
                x[i]=np.nan
                y[i]=np.nan
                z[i]=np.nan
                #FSurf[i]=np.nan
                break
            
        
    return x,y,z
        

def reclatvec(cell):
    assert cell.shape == (3, 3)
    
    rlv=np.zeros([3,3])
    a1=np.array([cell[0,:]])
    a2=np.array([cell[1,:]])
    a3=np.array([cell[2,:]])
    nm=(np.dot(a1,(np.cross(a2,a3)).T))
    
    rlv[0,:]=2*np.pi*((np.cross(a2,a3))/nm)
    rlv[1,:]=2*np.pi*((np.cross(a3,a1))/nm)
    rlv[2,:]=2*np.pi*((np.cross(a1,a2))/nm)
    
    return rlv 
    
def SF(FerSur,Ef,o,IMS):
    #IMS=0.001
    SF=np.zeros(FerSur.shape)
    FerSur=FerSur+Ef
    SF=(-1/np.pi)*(IMS)/((o-FerSur)**2+IMS**2)
    #noise2 = np.random.poisson(SF,[SF.shape[0],SF.shape[1]])
    # noise = np.random.normal(0,np.nanstd(SF),[SF.shape[0],SF.shape[1]])
    
    # pois=SF

    # for i in range(len(pois[:,0])):
    #     for j in range(len(pois[0,:])):
    #         if SF[i,j]>0:
    #             pois[i,j]=np.random.poisson(SF[i,j])
    #         if noise[i,j]>0:
    #             noise[i,j]=np.random.poisson(noise[i,j])
                

    
    # SF=(SF*pois)+noise
    return SF #pois+noise

def SFFSCont(x,y,SpecFunc,skipcon):
    
    cs=plt.contour(x,y,SpecFunc,50,cmap='viridis')
    plt.close()
    store=[]
    
    print(len(cs.allsegs))
    
    for i in range(len(cs.allsegs)):
        if len(cs.allsegs[i]) == 0:  #First few contours are interfering with BZ
            continue
        if i <= skipcon:
            continue
        o=np.concatenate(cs.allsegs[i])
        store.append(o)
        # plt.plot(o[:, 0], o[:, 1],'.', linewidth=2)
        # plt.show()
        
    store=np.concatenate(store)
    
    df = pd.DataFrame(store, columns = ["xval", "yval"])
 
    
    df2 = pd.DataFrame(df[(df['xval']>0) & (df['yval']>0)].groupby("yval").mean())
    df2 = df2.reset_index()
    df2 = pd.DataFrame(df2.groupby("xval").mean())
    df2 = df2.reset_index()      
    

    df3 = pd.DataFrame(df[(df['xval']>0) & (df['yval']<0)].groupby("yval").mean())
    df3 = df3.reset_index()
    df3 = pd.DataFrame(df3.groupby("xval").mean())
    df3 = df3.reset_index()
    df3 = df3.iloc[::-1]
            

    df4 = pd.DataFrame(df[(df['xval']<0) & (df['yval']>0)].groupby("yval").mean())
    df4 = df4.reset_index()
    df4 = pd.DataFrame(df4.groupby("xval").mean())
    df4 = df4.reset_index()


    df5 = pd.DataFrame(df[(df['xval']<0) & (df['yval']<0)].groupby("yval").mean())
    df5 = df5.reset_index()
    df5 = pd.DataFrame(df5.groupby("xval").mean())
    df5 = df5.reset_index()
    df5 = df5.iloc[::-1]
    
    df2=df2.append(df3, ignore_index = True)
    df2=df2.append(df5, ignore_index = True)
    df2=df2.append(df4, ignore_index = True)
    avg=df2.to_numpy()
    
    #avg=avg[avg[:,0].argsort()]
    
    store=np.array(store)

    return store #avg

def outliers(cont,pointskips,tol):   #At the moment the data needs to be ordered so it follows the contour in one direction 
    distances=[]  
    cir=[]
    shiftedx=[]
    orx=[]
    ory=[]
    orx2=[]
    ory2=[]
    ol=[]
    neck=[]
    
    for i in range(len(cont)):
        if i == (len(cont)-1):
            distances.append(m.dist(cont[i,:], cont[0,:]))
            w=cont[i,0]+((cont[i,0]-cont[0,0])/2)
            w1=w
            u=cont[i,1]+((cont[i,1]-cont[0,1])/2)
            u1=u
            d=m.dist(cont[i,:], cont[0,:])
            d=abs(d)
        else:
            distances.append(m.dist(cont[i,:], cont[i+1,:]))
            w=cont[i,0]+((cont[i+1,0]-cont[i,0])/2)
            w1=w
            u=cont[i,1]+((cont[i+1,1]-cont[i,1])/2)
            u1=u
            d=m.dist(cont[i,:], cont[i+1,:])
            d=abs(d)
        ory2.append(u1)
        orx2.append(w1)
        if i == 0:
            dold=d
        if d > 50*dold:
            #print(i,d)
            if i == (len(cont)-2):
                d1=m.dist(cont[(len(cont)-1),:], cont[0,:])
                d1=abs(d1)
                d2=m.dist(cont[i,:], cont[i-1,:])
                d2=abs(d2) 
            if i == (len(cont)-1):
                d1=m.dist(cont[0,:], cont[1,:])
                d1=abs(d1)
                d2=m.dist(cont[i,:], cont[i-1,:])
                d2=abs(d2) 
            else:
                d1=m.dist(cont[i+1,:], cont[i+2,:])
                d1=abs(d1)
                d2=m.dist(cont[i,:], cont[i-1,:])
                d2=abs(d2) 
            #print((d/d2)**2, (d/d1)**2, (((((((d/d2)**2+(d/d1)**2)**2)*d2**2)**2)*d1**2)**2),i,'\n',d,d1,d2)
            if (((d/d2)**2)*((d/d1)**2)) > 1000**2 and (((((((d/d2)**2+(d/d1)**2)**2)*d2**2)**2)*d1**2)**2) > 0.9 :
                neck.append(i)
                continue
            else:
                ol.append(i)
                #print(d,d1,d2,i)
                #print((d/d2)**2, (d/d1)**2, (((((((d/d2)**2+(d/d1)**2)**2)*d2**2)**2)*d1**2)**2),i)
            
        s=w**2+u**2
        w=np.sqrt(w**2/s)*d*(w1/abs(w1))
        u=np.sqrt(u**2/s)*d*(u1/abs(u1))
            
        cir.append(u)
        shiftedx.append(w)
        ory.append(u1)
        orx.append(w1)
        dold=d

    distances=np.array(distances)
    ory=np.array(ory)
    orx=np.array(orx)
    ory2=np.array(ory2)
    orx2=np.array(orx2)
    shiftedx=np.array(shiftedx)
    cir=np.array(cir)
    neck=np.array(neck)
    
    mx=np.where(np.sqrt(cir**2+shiftedx**2) == np.amax(np.sqrt(cir**2+shiftedx**2)))[0][0]
    
    scale=((cont[mx,0]**2+cont[mx,1]**2)*1.2)/max((np.sqrt(cir**2+shiftedx**2)))
    
    s=orx**2+ory**2
    cirx=np.sqrt(orx**2/4)*(orx/abs(orx))
    ciry=np.sqrt(ory**2/4)*(ory/abs(ory))
    
    shiftedx=(np.array(shiftedx)*scale)+cirx
    cir=(np.array(cir)*scale)+ciry  
    
    s=orx2**2+ory2**2
    cirx2=np.sqrt(orx2**2/4)*(orx2/abs(orx2))
    ciry2=np.sqrt(ory2**2/4)*(ory2/abs(ory2))
         
    average=sum(distances)/len(distances)
    
    s=cont[:,0]**2+cont[:,1]**2
    w=(np.sqrt(cont[:,0]**2/s)*average*(cont[:,0]/abs(cont[:,0]))*scale)+cirx2
    u=(np.sqrt(cont[:,1]**2/s)*average*(cont[:,1]/abs(cont[:,1]))*scale)+ciry2
    
    for j in range(pointskips):
        pointsskipped=j
        for i in range(len(cont)):
            if i >= (len(cont)-pointsskipped-1):
               p=abs(m.dist(cont[i-(1+pointsskipped),:], cont[i,:]))+abs(m.dist(cont[i,:], cont[0+(pointsskipped+(1-len(cont))),:]))
            else:
               p=abs(m.dist(cont[i-(1+pointsskipped),:], cont[i,:]))+abs(m.dist(cont[i,:], cont[i+(1+pointsskipped),:]))
               
            if p > (average*2*(1+pointsskipped)):
                ol.append(i)
                ol.append(i-(1+pointsskipped))
                ol.append(i+(1+pointsskipped))
            
    ol=np.array(ol)
    
    ol = [i for i in ol if i not in neck] 
    
    for i in range(len(ol)):
        if ol[i] > (len(cont)-1):
            ol[i]=ol[i]-len(cont)

    df = pd.DataFrame(ol, columns = ["position"])
    df=df.groupby(["position"])["position"].count().reset_index(name="count")
    ol=df.to_numpy()
    olp=cont[ol[:,0],:]
    
    cmap = mcolors.LinearSegmentedColormap.from_list("mymap", ["yellow","orange","red"])
        
    plt.figure(figsize=(5, 3), dpi=200)
    plt.xlim(np.min(cont[:,0])*1.5,np.max(cont[:,0])*1.5)
    plt.ylim(np.min(cont[:,1])*1.5,np.max(cont[:,1])*1.5)
    plt.title('Extracted Contour with distance between points')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.plot(cont[:,0],cont[:,1],'.', linewidth=2,label='SF Contour')
    plt.plot([cirx[0], shiftedx[0]], [ciry[0], cir[0]],color='purple',label='Distance between points')
    plt.plot([cirx, shiftedx], [ciry, cir],color='purple')
    #plt.plot([np.zeros(len(shiftedx)), shiftedx], [np.zeros(len(cir)), cir])
    plt.plot(w,u,'.', linewidth=2, color='orange', label='Average')
    plt.scatter(olp[:,0],olp[:,1],c=ol[:,1],linewidth=2,vmin=np.min(ol[:,1]), vmax=np.max(ol[:,1]),cmap=cmap, label='SF Contour Outlier Chance')
    plt.legend(loc='lower right',prop={'size': 5})
    plt.show()
    
    threshold=np.round(tol*np.max(ol[:,1]))
    ol=ol[ol[:,1] >= threshold]
    cont=np.delete(cont, ol[:,0], axis=0)
    
    return cont,neck

def orthogonal_points(coords, s):
    """Given a tuple coords containing two 2-dimensional points and also
    given a positive number s, return two other distinct points such
    that the line segment between each output point and coords[0] is
    orthogonal (perpendicular) to the line segment between coords[0] and
    coords[1] and the distance from each output point to coords[0] is s.
    """
    (point1x, point1y), (point2x, point2y) = coords
    points_vectorx, points_vectory = point2x - point1x, point2y - point1y
    points_vector_length = math.hypot(points_vectorx, points_vectory)
    normalized_x, normalized_y = (points_vectorx * s / points_vector_length,
                                  points_vectory * s / points_vector_length)
    newpoint1x, newpoint1y = point1x + normalized_y, point1y - normalized_x
    newpoint2x, newpoint2y = point1x - normalized_y, point1y + normalized_x
    return ([newpoint1x, newpoint1y], [newpoint2x, newpoint2y])


def display_inlier_outlier(cloud, ind):
    inlier_cloud = cloud.select_by_index(ind)
    outlier_cloud = cloud.select_by_index(ind, invert=True)

    print("Showing outliers (red) and inliers (gray): ")
    outlier_cloud.paint_uniform_color([1, 0, 0])
    inlier_cloud.paint_uniform_color([0.8, 0.8, 0.8])
    o3d.visualization.draw([inlier_cloud, outlier_cloud],show_skybox=False, bg_color=(0.0, 0.0, 0.0, 1.0))

def mirrorfold(sf,x,y):
    
    #x and y are the lines at which the folding happens 
    
    dd=[]
    # Get all valid indices
    idx1 = np.flatnonzero(sf[:,0]>=x)
    idx2 = np.flatnonzero(sf[:,0]<=x)
    idx3 = np.flatnonzero(sf[:,1]>=y)
    idx4 = np.flatnonzero(sf[:,1]<=y)

    q1=sf[np.intersect1d(idx1,idx3)]
    q2=sf[np.intersect1d(idx2,idx4)]
    q3=sf[np.intersect1d(idx3,idx2)]
    q4=sf[np.intersect1d(idx4,idx1)]

    q=q1
    q=np.append(q,(q2*-1),axis=0)
    q=np.append(q,(q3*[-1,1]),axis=0)
    q=np.append(q,(q4*[1,-1]),axis=0)

    dd=q
    dd=np.append(dd,(q*-1),axis=0)
    dd=np.append(dd,(q*[-1,1]),axis=0)
    dd=np.append(dd,(q*[1,-1]),axis=0)
    
    return dd


def kdefilter(d,bndwdth):
    
    kde=KernelDensity(kernel='gaussian', bandwidth=bndwdth).fit(d) #
    score=kde.score_samples(d)

    plt.axes().set_aspect('equal')
    sns.kdeplot(d,x=d[:, 0], y=d[:, 1], cmap="Blues", fill=True, thresh=0.05)
    plt.scatter(d[:, 0], d[:, 1],c=score)
    plt.colorbar()
    plt.show()

    lim=input('Select bottom cut off value:') # Get the input
    while lim != "":        # Loop until it is a blank line
        idx = np.flatnonzero(score>=float(lim))

        plt.axes().set_aspect('equal')
        plt.scatter(d[idx, 0], d[idx, 1],c=score[idx])
        plt.colorbar()
        plt.show()
        lim = input('OPTIONS: \nSelect NEW cut off value,\ntype "redo" to recalculate the KDE for the current cut off or\nhit enter to collect data!:')   # Get the input again
        if lim == "redo":
            d=d[idx,:]
            kde2=KernelDensity(kernel='gaussian', bandwidth=bndwdth).fit(d) #
            score=kde2.score_samples(d)

            plt.axes().set_aspect('equal')
            sns.kdeplot(d,x=d[:, 0], y=d[:, 1], cmap="Blues", fill=True, thresh=0.05)
            plt.scatter(d[:, 0], d[:, 1],c=score)
            plt.colorbar()
            plt.show() 
            
            idx = np.flatnonzero(score>=np.min(score))
            
            lim=input('Select bottom cut off value or\nhit enter to collect data!:') # Get the input

    d=d[idx,:]    
    
    
    lim=input('Do you want your data downsampled? \n"y" or "n":') # Get the input
    if lim == "y":        # Loop until it is a blank line
        lim=input('How do you want to downsample? \n1: Groups points by corresponding mathcing x then y coordinates. \n2: Uses a radius about a point to average the data within that radius about a point. \n"1" or "2"?:') # Get the input
        if lim == "1":
                db4=d    
            
                df = pd.DataFrame(d, columns = ["xval", "yval"])
        
                df2 = pd.DataFrame(df[(df['xval']>0) & (df['yval']>0)].groupby("yval").mean())
                df2 = df2.reset_index()
                df2 = pd.DataFrame(df2.groupby("xval").mean())
                df2 = df2.reset_index()      
                
            
                df3 = pd.DataFrame(df[(df['xval']>0) & (df['yval']<0)].groupby("yval").mean())
                df3 = df3.reset_index()
                df3 = pd.DataFrame(df3.groupby("xval").mean())
                df3 = df3.reset_index()
                df3 = df3.iloc[::-1]
                        
            
                df4 = pd.DataFrame(df[(df['xval']<0) & (df['yval']>0)].groupby("yval").mean())
                df4 = df4.reset_index()
                df4 = pd.DataFrame(df4.groupby("xval").mean())
                df4 = df4.reset_index()
            
            
                df5 = pd.DataFrame(df[(df['xval']<0) & (df['yval']<0)].groupby("yval").mean())
                df5 = df5.reset_index()
                df5 = pd.DataFrame(df5.groupby("xval").mean())
                df5 = df5.reset_index()
                df5 = df5.iloc[::-1]
                
                df2=df2.append(df3, ignore_index = True)
                df2=df2.append(df5, ignore_index = True)
                df2=df2.append(df4, ignore_index = True)
                d=df2.to_numpy()
                
                plt.rcParams['figure.figsize'] = [10,8]
                plt.plot(db4[:, 0], db4[:, 1],'.', linewidth=1,label='Before.',)
                plt.plot(d[:, 0], d[:, 1],'.', linewidth=2,label='After',)
                plt.legend()
                plt.show()
        
        if lim == "2":
                cpoint=np.flatnonzero(score>=(max(score)*0.80))[0]
                radest=m.dist(d[cpoint,:], d[cpoint+1,:])*1


                fig1 = plt.figure()
                ax1 = fig1.add_subplot(111, aspect='equal')
                plt.scatter(d[:, 0], d[:, 1],c='black')
                plt.scatter(d[cpoint, 0], d[cpoint, 1],c='red')
                ax1.add_patch(patches.Circle((d[cpoint, 0], d[cpoint, 1]), radest, color='g', linewidth=1, fill=False))
                plt.show()
                
                fig1 = plt.figure()
                ax1 = fig1.add_subplot(111, aspect='equal')
                plt.xlim(d[cpoint, 0]-(radest*2),d[cpoint, 0]+(radest*2))
                plt.ylim(d[cpoint, 1]-(radest*2),d[cpoint, 1]+(radest*2))
                plt.scatter(d[:, 0], d[:, 1],c='black')
                plt.scatter(d[cpoint, 0], d[cpoint, 1],c='red')
                ax1.add_patch(patches.Circle((d[cpoint, 0], d[cpoint, 1]), radest, color='g', linewidth=1, fill=False))
                plt.show()
                print(radest)
                lim=input('Input Radius? \n"y" or "n":') # Get the input
                if lim=="y":
                    radestnew=input('Radius:')
                    #radestnew=radest                        
                    while radestnew !=" ":    
                            
                        fig1 = plt.figure()
                        ax1 = fig1.add_subplot(111, aspect='equal')
                        plt.xlim(d[cpoint, 0]-(radestnew*2),d[cpoint, 0]+(radestnew*2))
                        plt.ylim(d[cpoint, 1]-(radestnew*2),d[cpoint, 1]+(radestnew*2))
                        plt.scatter(d[:, 0], d[:, 1],c='black')
                        plt.scatter(d[cpoint, 0], d[cpoint, 1],c='red')
                        ax1.add_patch(patches.Circle((d[cpoint, 0], d[cpoint, 1]), radestnew, color='g', linewidth=1, fill=False))
                        ax1.add_patch(patches.Circle((d[cpoint, 0], d[cpoint, 1]), radest, color='r', linewidth=1, fill=False))
                        plt.show()
                            
                        radest=radestnew
                        
                        radestnew=input('Radius (enter to take last value):')
                            
                            
                    
                meanx=[]
                meany=[]
                dhold=np.array(d)
                deletemask=np.array(d)
                y=0

                for y in range(len(dhold)):
                    if np.isnan(dhold[y,0]) == True:
                        continue
                    x=np.where(d==dhold[y,:])[0][0]
                    maskmap = d-d[x,:]
                    mask = maskmap[:,0]**2 + maskmap[:,1]**2 <= 0.05**2
                    deletemask[mask,:]=np.nan
                    dhold=deletemask #np.delete(dhold,mask,axis=0)
                    #un-mask center and values below 0
                    #mask[centerMask]=False
                    #mask[inputMat<0]=False

                    #get the mean
                    meanx.append(np.mean(d[mask,0]))
                    meany.append(np.mean(d[mask,1]))
                    
                    deletemask[x,0],deletemask[x,1]=np.mean(d[mask,0]),np.mean(d[mask,1])
                    #y=y+1

                meanx=np.array(meanx)
                meany=np.array(meany)
                meanAll=np.zeros([len(meanx),2])
                meanAll[:,0]=meanx
                meanAll[:,1]=meany
                
                fig1 = plt.figure()
                ax1 = fig1.add_subplot(111, aspect='equal')
                plt.scatter(d[:, 0], d[:, 1],c='black',label='Old')
                plt.scatter(deletemask[:, 0], deletemask[:, 1],c='red',label='Removed')
                plt.scatter(meanAll[:, 0], meanAll[:, 1],c='blue',label='New')
                plt.legend()
                plt.show()
                        
                #lim=input('Hit enter to keep! Anything else to retry!:') # Get the input
                    
                d=meanAll

            
    lim=input('More KDE? \n"y" or "n":') # Get the input
    if lim == "y":        # Loop until it is a blank line

            kde2=KernelDensity(kernel='gaussian', bandwidth=bndwdth).fit(d) #
            score=kde2.score_samples(d)

            plt.axes().set_aspect('equal')
            sns.kdeplot(d,x=d[:, 0], y=d[:, 1], cmap="Blues", fill=True, thresh=0.05)
            plt.scatter(d[:, 0], d[:, 1],c=score)
            plt.colorbar()
            plt.show() 

            lim=input('Select bottom cut off value:') # Get the input
            while lim != "":        # Loop until it is a blank line
                idx = np.flatnonzero(score>=float(lim))
        
                plt.axes().set_aspect('equal')
                plt.scatter(d[idx, 0], d[idx, 1],c=score[idx])
                plt.colorbar()
                plt.show()
                lim = input('OPTIONS: \nSelect NEW cut off value,\ntype "redo" to recalculate the KDE for the current cut off or\nhit enter to collect data!:')   # Get the input again
                if lim == "redo":
                    d=d[idx,:]
                    kde2=KernelDensity(kernel='gaussian', bandwidth=bndwdth).fit(d) #
                    score=kde2.score_samples(d)
        
                    plt.axes().set_aspect('equal')
                    sns.kdeplot(d,x=d[:, 0], y=d[:, 1], cmap="Blues", fill=True, thresh=0.05)
                    plt.scatter(d[:, 0], d[:, 1],c=score)
                    plt.colorbar()
                    plt.show() 
                    
                    idx = np.flatnonzero(score>=np.min(score))

                    lim=input('Select bottom cut off value or\nhit enter to collect data!:') # Get the input

            d=d[idx,:]    
    
    #Neck checklist:
        #Distance between previous 3 nearest neighbours are further than the next neighbour.
        #Distance between the point and the edge of the BZ is smaller than the next point.
        #The gradient of the 3 previous neighbours  
        #It can only be a neck if the sandwhiching distance and distance ratios are small and not in the same magnitude 
            #Use   
        
        #1st check to see if a neck may be present if so then we go into more details
    distances=[]
    for i in range(len(d)):
        if i == (len(d)-1):
            distances.append(m.dist(d[i,:], d[0,:]))
        else:
            distances.append(m.dist(d[i,:], d[i+1,:]))

    distances=np.array(distances)
    distances=abs(distances)    
    sm=0
    neckdetector=False
    necklist=[]
    ratio=[]

    # for i in range(len(distances)):
    #     sm=sm+distances[i]
        
    #     if i == (len(distances)-1):
    #         if distances[0] > 20*(sm/(i+1)):
    #             neckdetector=True
    #             necklist.append(i+1)
    #             continue 
    #     else:
    #         ratio.append(distances[i+1]/(sm/(i+1)))
    #         if distances[i+1] > 20*(sm/(i+1)):
    #             neckdetector=True
    #             necklist.append(i+1)
    #             sm=sm-distances[i+1]
    #             continue
            
    #ratio=np.array(ratio)
    distanceindex=distances.argsort()
    for i in range(10):
        if distanceindex[-i]==(len(distances)-1):
            if (distances[distanceindex[-(i+1)]]/distances[distanceindex[-i]])>=0.80 or (distances[distanceindex[-i]]/distances[distanceindex[-(i-1)]])>=0.80:
                if distances[distanceindex[-i]]/(distances[distanceindex[-i]-1]+distances[0])>1.5:
                    if np.floor(np.log10(np.abs(distances[distanceindex[-i]])))>np.floor(np.log10(np.abs(distances[0]))) or np.floor(np.log10(np.abs(distances[distanceindex[-i]])))>np.floor(np.log10(np.abs(distances[distanceindex[-i]-1]))):
                                       necklist.append(distanceindex[-i])
                                       neckdetector=True

        elif i==len(distanceindex):
            if (distances[distanceindex[-1]]/distances[distanceindex[-i]])>=0.80 or (distances[distanceindex[-i]]/distances[distanceindex[-(i-1)]])>=0.80:
                if distances[distanceindex[-i]]/(distances[distanceindex[-i]-1]+distances[distanceindex[-i]+1])>1:
                    if np.floor(np.log10(np.abs(distances[distanceindex[-i]])))>np.floor(np.log10(np.abs(distances[distanceindex[-i]+1]))) or np.floor(np.log10(np.abs(distances[distanceindex[-i]])))>np.floor(np.log10(np.abs(distances[distanceindex[-i]-1]))):
                                  necklist.append(distanceindex[-i])
                                  neckdetector=True

        else:
            if (distances[distanceindex[-(i+1)]]/distances[distanceindex[-i]])>=0.80 or (distances[distanceindex[-i]]/distances[distanceindex[-(i-1)]])>=0.80:
                if distances[distanceindex[-i]]/(distances[distanceindex[-i]-1]+distances[distanceindex[-i]+1])>0.9:
                    if np.floor(np.log10(np.abs(distances[distanceindex[-i]])))>np.floor(np.log10(np.abs(distances[distanceindex[-i]+1]))) or np.floor(np.log10(np.abs(distances[distanceindex[-i]])))>np.floor(np.log10(np.abs(distances[distanceindex[-i]-1]))):
                        necklist.append(distanceindex[-i])
                        neckdetector=True
        
        
        
    #print(necklist)
    necklist=np.array(necklist)
    necklist=necklist[necklist.argsort()]
    notnecks=[]
    
    if neckdetector==True:
        
        plt.axes().set_aspect('equal')
        plt.scatter(d[:, 0], d[:, 1],c='black')
        plt.scatter(d[necklist, 0], d[necklist, 1],c='green')
        plt.show()
        
        print(len(necklist),'Potential necks found')
        lim=input('Are necks present? "y" or "n":') # Get the input
       
        if lim=="y":
            #bzp=BZPlanes(rcell*0.92, 1)
            # kde2=KernelDensity(kernel='gaussian', bandwidth=bndwdth).fit(d) #
            # score=kde2.score_samples(d2[necklist])
            # print(score)
            for i in range(len(necklist)):
                if necklist[i] == 0:            #NEED TO ADJUST THIS SO IT CHECKS THE DISTANCE OF THE PREVIOUS POINT AND THE NEXT AND SELECTS THE LARGER DISTANCE
                    p=1
                elif necklist[i] % 2 == 0:
                    p=1 # Even 
                else:
                    p=-1
                
                plt.axes().set_aspect('equal')
                #sns.kdeplot(d2,x=d2[:, 0], y=d2[:, 1], cmap="Blues", fill=True, thresh=0.05)
                plt.scatter(d[:, 0], d[:, 1],c='black')
                plt.scatter(d[necklist[i]:(necklist[i]+p), 0], d[necklist[i]:(necklist[i]+p), 1],c='green')
                plt.plot(d[necklist[i]:(necklist[i]+p), 0], d[necklist[i]:(necklist[i]+p), 1],c='red')
                plt.colorbar()
                plt.show()
                
                lim=input('Is this a neck? "y" or "n"')
                if lim == "y":
                    print('saved')
                if lim == "n":
                    notnecks.append(i)
        
        necklist=np.delete(necklist,notnecks,axis=0)
        
        print(necklist)
        
    return d #,distances,distanceindex


def radMask(index,radius,array):
  a,b = index
  nx,ny = array.shape
  y,x = np.ogrid[-a:nx-a,-b:ny-b]
  mask = x*x + y*y <= radius*radius

  return mask

def rotate(origin, point, angle):
    """
    Rotate a point counterclockwise by a given angle around a given origin.

    The angle should be given in degrees.
    """
    angle=angle*(m.pi/180)
    ox, oy = origin
    
    dta=point.copy()
    
    for i in range(len(point)):
        px, py = point[i,:]
        qx = ox + math.cos(angle) * (px - ox) - math.sin(angle) * (py - oy)
        qy = oy + math.sin(angle) * (px - ox) + math.cos(angle) * (py - oy)
        
        dta[i,0]=qx
        dta[i,1]=qy
        
    return dta

def Downsample(d,datatype):
    
    if datatype=='2D':
        
        #meanAll is going to store ~18000 points
        meanAll=np.zeros((600,600))
        
        for x in range(600):
            for y in range(600):
                centerMask=(x,y)
                mask=radMask(centerMask,50,d)
                #un-mask center and values below 0
                mask[centerMask]=False
                #mask[d<0]=False
        
                #get the mean
                meanAll[x,y]=np.mean(d[mask])
    else:
            
        
        ymin=np.min(d[:,1])
        xmin=np.min(d[:,0])
        
        d[:,0]=d[:,0]+abs(xmin)
        d[:,1]=d[:,1]+abs(ymin)
    
        cpoint=int(np.round(len(d)/2))
        
        radest=m.dist(d[cpoint,:], d[cpoint+1,:])*1
        
        fig1 = plt.figure()
        ax1 = fig1.add_subplot(111, aspect='equal')
        plt.scatter(d[:, 0], d[:, 1],c='black')
        plt.scatter(d[cpoint, 0], d[cpoint, 1],c='red')
        ax1.add_patch(patches.Circle((d[cpoint, 0], d[cpoint, 1]), radest, color='g', linewidth=1, fill=False))
        plt.show()
        
        fig1 = plt.figure()
        ax1 = fig1.add_subplot(111, aspect='equal')
        plt.xlim(d[cpoint, 0]-(radest*2),d[cpoint, 0]+(radest*2))
        plt.ylim(d[cpoint, 1]-(radest*2),d[cpoint, 1]+(radest*2))
        plt.scatter(d[:, 0], d[:, 1],c='black')
        plt.scatter(d[cpoint, 0], d[cpoint, 1],c='red')
        ax1.add_patch(patches.Circle((d[cpoint, 0], d[cpoint, 1]), radest, color='g', linewidth=1, fill=False))
        plt.show()
        print(radest)
        lim=input('Input Radius? \n"y" or "n":') # Get the input
        if lim=="y":
            radestnew=input('Radius:')
            #radestnew=radest                        
            while radestnew !="":    
                
                radestnew=float(radestnew)
                
                fig1 = plt.figure()
                ax1 = fig1.add_subplot(111, aspect='equal')
                plt.xlim(d[cpoint, 0]-(radestnew*2),d[cpoint, 0]+(radestnew*2))
                plt.ylim(d[cpoint, 1]-(radestnew*2),d[cpoint, 1]+(radestnew*2))
                plt.scatter(d[:, 0], d[:, 1],c='black')
                plt.scatter(d[cpoint, 0], d[cpoint, 1],c='red')
                ax1.add_patch(patches.Circle((d[cpoint, 0], d[cpoint, 1]), radestnew, color='g', linewidth=1, fill=False))
                ax1.add_patch(patches.Circle((d[cpoint, 0], d[cpoint, 1]), radest, color='r', linewidth=1, fill=False))
                plt.show()
                    
                radest=radestnew
                
                radestnew=input('Radius (enter to take last value):')
                    
                    
            
        meanx=[]
        meany=[]
        dhold=np.array(d)
        deletemask=np.array(d)
        y=0
    
        for y in range(len(dhold)):
            if np.isnan(dhold[y,0]) == True:
                continue
            x=np.where(d==dhold[y,:])[0][0]
            maskmap = d-d[x,:]
            mask = maskmap[:,0]**2 + maskmap[:,1]**2 <= 0.05**2
            deletemask[mask,:]=np.nan
            dhold=deletemask #np.delete(dhold,mask,axis=0)
            #un-mask center and values below 0
            #mask[centerMask]=False
            #mask[inputMat<0]=False
    
            #get the mean
            meanx.append(np.mean(d[mask,0]))
            meany.append(np.mean(d[mask,1]))
            
            deletemask[x,0],deletemask[x,1]=np.mean(d[mask,0]),np.mean(d[mask,1])
            #y=y+1
    
        meanx=np.array(meanx)
        meany=np.array(meany)
        meanAll=np.zeros([len(meanx),2])
        meanAll[:,0]=meanx
        meanAll[:,1]=meany
        
        d[:,0]=d[:,0]-abs(xmin)
        d[:,1]=d[:,1]-abs(ymin)
        
        fig1 = plt.figure()
        ax1 = fig1.add_subplot(111, aspect='equal')
        plt.scatter(d[:, 0], d[:, 1],c='black',label='Old')
        #plt.scatter(deletemask[:, 0], deletemask[:, 1],c='red',label='Removed')
        plt.scatter(meanAll[:, 0], meanAll[:, 1],c='blue',label='New')
        plt.legend()
        plt.show()
                
        #lim=input('Hit enter to keep! Anything else to retry!:') # Get the input
        
    #d=meanAll
    
    return meanAll

def remove_rows(dblearray, trimmingarray):
    while trimmingarray.shape[0] > 2 * dblearray.shape[0]:
        trimmingarray = np.delete(trimmingarray, -1, axis=0)
    return trimmingarray
#%%

b=0.5

cell = np.array([[0.0, b, b],
                  [b, 0.0, b],
                  [b, b, 0.0]])



rcell=reclatvec(cell*7.7289957)

# rcell=np.array([[0.8129368357,0.8129368357,-0.8129368357],
# [0.8129368357,-0.8129368357,-0.8129368357],
# [-0.000000000,-0.000000000,1.625873671]])

a=7.7289957 #4.0692 Angstrom

bz=BZCalc(cell,0.81293684*1) #np.max(p2)) 0.8129


#Ag7
C_0  = 0.62443
C_200=-0.03918
C_211=-0.60245
C_220=-0.08698
C_310=-0.08429
C_222=-0.00083
C_321=0.01406

#Ag5
# C_0  =-0.89789
# C_200=-0.12030
# C_211=-0.90187
# C_220=-0.14086
# C_310=-0.09483
# C_222=0
# C_321=0

#Define Plane
npoints=64 #for 2*pi*a to 0 
plane=[1,1,0] #plane you want to observe
k_x = k_y = k_z =np.linspace(-np.pi/2,np.pi/2,npoints)
k_z=k_z*0

FermiSurface=np.zeros([npoints**2,1])


plane2=np.zeros([1*npoints**2,3])
#plane[:,2]=0



ik=0
plane=np.zeros([1*npoints**2,3])
for j in range(npoints):
    for k in range(npoints):
        #for p in range(npoints):
        plane[ik,0]=k_x[j]
        plane[ik,1]=k_y[k]
        ik+=1


plane2[:,0],plane2[:,1],plane2[:,2]=BZCheck(rcell,plane[:,0],plane[:,1],plane[:,2],1)


for i in range(len(plane)):
        FermiSurface[i]=FS(plane2[i,0],plane2[i,1],plane2[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)


Efzero1=np.nanmin(np.abs(FermiSurface))


FermiSurface=FermiSurface.reshape([npoints,npoints])

x=plane2[:,0].reshape([npoints,npoints])
y=plane2[:,1].reshape([npoints,npoints])
z=plane2[:,2].reshape([npoints,npoints])


#sfp=SF(FermiSurface,-0.82,0,-0.06)
sfp=SF(FermiSurface,0,0,-0.06)
s=SFFSCont(x,y,sfp,skipcon=10)#,skipcon=20) #skipcon=4
#d,nk=outliers(s,pointskips=5,tol=0.80) 
d=s
d=np.insert(d, 2, np.zeros(len(d)), axis=1)

# plt.axes().set_aspect('equal')
# plt.title(' [001] Norm. Theory - Spec. Func. of Parameterisation')
# plt.xlim(-np.pi/2, np.pi/2)
# plt.ylim(-np.pi/2, np.pi/2)
# plt.xlabel('100')
# plt.ylabel('010')
# cs=plt.contour(x,y,sfp,50,cmap='viridis')#,levels=0)
# #cs=plt.contour(x,y,FermiSurface,50,cmap='viridis')#,levels=0)
# #plt.xlim([lq,l3q])
# #plt.ylim([lq,l3q])
# plt.colorbar()
# plt.show()  

# plt.plot(d[:, 0], d[:, 1],'.', linewidth=2)
# plt.show()    

plt.figure(figsize=(5, 3), dpi=200)
plt.axes().set_aspect('equal')
ax = plt.subplot(111)
#ax.spines.right.set_visible(False)
#ax.spines.top.set_visible(False)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.xaxis.label.set_color('black')
ax.yaxis.label.set_color('black')
ax.tick_params(axis='x', colors='black')
ax.tick_params(axis='y', colors='black')



#title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Noisy Spectral Function Transformation\n',c='k')
#plt.title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Isosurface',c='k')
plt.xlim(-np.pi/4, np.pi/4)
plt.ylim(-np.pi/4, np.pi/4)
plt.locator_params(axis='x', nbins=3)
plt.locator_params(axis='y', nbins=3)  
plt.xlabel(r'[100] $\left(\frac{2\pi}{a}\right)$')
plt.ylabel(r'[010] $\left(\frac{2\pi}{a}\right)$')
#cs=plt.contour(x,y,sfp,50,cmap='viridis')#,levels=0)
cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)
#savefig('Huesler111.png', transparent=True,bbox_inches = "tight")
#plt.savefig('Images/Ag2DParam[001].svg',transparent=True,bbox_inches = "tight", format='svg')
show()


fig = plt.figure(figsize=(12, 12))#,layout='constrained')
ax = fig.add_subplot(projection='3d')
fig.tight_layout()
#ax.set_box_aspect(aspect=None, zoom=1.4)
ax.set_aspect('equal')
#fig.subplots_adjust(top=10)
#plt.tight_layout()

ax.set_xlim3d(-0.6,0.6)
ax.set_ylim3d(-0.6,0.6)
ax.set_zlim3d(-0.6,0.6)

# make the panes transparent
ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
# make the grid lines transparent
ax.xaxis._axinfo["grid"]['color'] =  (1,1,1,0)
ax.yaxis._axinfo["grid"]['color'] =  (1,1,1,0)
ax.zaxis._axinfo["grid"]['color'] =  (1,1,1,0)
# Get rid of the spines
ax.xaxis.line.set_color((1.0, 1.0, 1.0, 1.0))
ax.yaxis.line.set_color((1.0, 1.0, 1.0, 1.0))
ax.zaxis.line.set_color((1.0, 1.0, 1.0, 1.0))
# Get rid of the ticks
ax.set_xticks([]) 
ax.set_yticks([]) 
ax.set_zticks([])


ax.scatter(plane2[:,0],plane2[:,1],plane2[:,2],color='purple',s=200)
# ax.scatter(v2[:,0],v2[:,1],v2[:,2])
#ax.plot(box[:,0],box[:,1],box[:,2])
#ax.plot(linr2[:,0],linr2[:,1],linr2[:,2])

for i in range(len(bz[1])):
      ax.plot(bz[1][i][:,0],bz[1][i][:,1],bz[1][i][:,2],c='b',linewidth=3)
      ax.plot(bz[2][i][:,0],bz[2][i][:,1],bz[2][i][:,2],c='b')

plt.savefig('Images/BZSlice2.svg',transparent=True,bbox_inches = "tight", format='svg')


plt.show()

#%%

plt.axes().set_aspect('equal')
plt.title(' [110] Theory')
plt.xlim(-np.pi/2, np.pi/2)
plt.ylim(-np.pi/2, np.pi/2)
plt.xlabel('100')
plt.ylabel('010')
cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)
#plt.xlim([lq,l3q])
#plt.ylim([lq,l3q])
plt.colorbar()
plt.show()

p=cs.allsegs[1][-1]
p=p[p[:,0].argsort()]
plt.plot(p[:, 0], p[:, 1],'.', linewidth=2)
plt.show()    

p=np.insert(p, 2, np.zeros(len(p)), axis=1)



plt.figure(figsize=(5, 3), dpi=200)
plt.axes().set_aspect('equal')
ax = plt.subplot(111)
#ax.spines.right.set_visible(False)
#ax.spines.top.set_visible(False)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.xaxis.label.set_color('black')
ax.yaxis.label.set_color('black')
ax.tick_params(axis='x', colors='black')
ax.tick_params(axis='y', colors='black')


#title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Noisy Spectral Function Transformation\n',c='k')
#plt.title('2D Parameterisation of the Fermi Surface \n- Ag [110] Plane Normal -\n Isosurface (Points)',c='k')
plt.xlim(-np.pi/4, np.pi/4)
plt.ylim(-np.pi/4, np.pi/4)
plt.yticks(ticks=[-0.5,0,0.5])
plt.xticks(ticks=[-0.5,0,0.5])
#plt.locator_params(axis='x', nbins=3)
#plt.locator_params(axis='y', nbins=3)  
plt.xlabel(r'[100] $\left(\frac{2\pi}{a}\right)$')
plt.ylabel(r'[010] $\left(\frac{2\pi}{a}\right)$')
#cs=plt.contour(x,y,sfp,50,cmap='viridis')#,levels=0)
#cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)
plt.plot(p[:, 0], p[:, 1],'.', markersize=0.5)
#savefig('Huesler111.png', transparent=True,bbox_inches = "tight")
#plt.savefig('Images/Ag2DParam[001]Points.svg',transparent=True,bbox_inches = "tight", format='svg')
show()


# t=p.copy()

# t=mirrorfold(p[:,0:2],0,0)

# plt.rcParams['figure.figsize'] = [10,8]
# plt.plot(t[:, 0], t[:, 1],'.', linewidth=1,label='Comb.',)
# plt.plot(p[:, 0], p[:, 1],'.', linewidth=2,label='Full Orig Spec Func',)
# plt.legend()
# plt.show()   

# ==============================SPECTRAL FUNCTION==============================
# d=mirrorfold(s,0,0)
# 
# plt.rcParams['figure.figsize'] = [10,8]
# plt.plot(d[:, 0], d[:, 1],'.', linewidth=1,label='Comb.',)
# plt.plot(s[:, 0], s[:, 1],'.', linewidth=2,label='Full Orig Spec Func',)
# plt.legend()
# plt.show()   
# 
# d=kdefilter(d,0.1)
# 
# d=np.insert(d, 2, np.zeros(len(d)), axis=1)
# =============================================================================


#%%



#=============================Probably a placeholder==========================
# FermiSurface001=np.zeros([len(p)])
# for i in range(len(p)): 
#     FermiSurface001[i]=FS(p[i,0],p[i,1],p[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)
#=============================================================================

# =============================================================================
# plt.axes().set_aspect('equal')
# plt.title('Comparison of the Eigenvalues and extracted Spectral function')    
# plt.plot(p[:, 0], p[:, 1],'.', linewidth=4,c='black',label='Eigenvalues')
# plt.scatter(d[:, 0], d[:, 1],c='orange',label='Spectral Function')
# plt.legend()
# plt.show()
# 
# =============================================================================

#%%
x,y,z,A=loadtxt('FERMISURF-2D.OUT',skiprows=1,unpack=True)
x=x.reshape([64,64])[0]
y=y.reshape([64,64])#[0]
z=z.reshape([64,64])[:,0]
A=A.reshape([64,64])

z2=[]
x2=[]
for i in range(len(x)-1):
    x2.append(-x[-(i+1)])
    z2.append(-z[-(i+1)])


for i in range(len(x)-1):
    x2.append(x[i])
    z2.append(z[i])

    
x=np.array(x2)
z=np.array(z2)

#A=np.matrix(A)
A2=A.tolist()
eA=np.lib.pad(A2,((0,62),(62,0)),'reflect')

plt.axes().set_aspect('equal')
plt.title(' [110] Theory')
plt.xlim(-(x[-1]/2), (x[-1]/2))
plt.ylim(-(z[-1]/2), (z[-1]/2))
#plt.xlim(-x[-1], x[-1])
#plt.ylim(-z[-1], z[-1])
plt.xlabel('100')
plt.ylabel('010')
#plt.contourf((x-(x[-1]/2)),(z-(z[-1]/2)),A,50,cmap='viridis')
cs=plt.contour(x,z,eA,50,cmap='viridis',levels=0)
#plt.xlim([lq,l3q])
#plt.ylim([lq,l3q])
plt.colorbar()
plt.show()

fig, axes = plt.subplots(1,2, figsize=(8, 4), sharex=False, sharey=False)
#plt.axes(aspect='equal')
ax = axes.ravel()

#fig.suptitle('[001] Theoretical 2D Fermi Surface Conotour',fontsize=15,y=0.98)

ax[0].contour(x+0.025404275999999948,z,eA,50,cmap='viridis') #+0.025404275999999948 for centering purposes
ax[0].set_xlim(-(x[-1]/2), (x[-1]/2))
ax[0].set_ylim(-(z[-1]/2), (z[-1]/2))
#ax[0].set_title('No Isovalue level')
# Hide X and Y axes label marks
ax[0].xaxis.set_tick_params(labelbottom=False)
ax[0].yaxis.set_tick_params(labelleft=False)

# Hide X and Y axes tick marks
ax[0].set_xticks([])
ax[0].set_yticks([])
#ax[0].axis('off')

ax[1].contour(x+0.025404275999999948,z,eA,50,cmap='viridis',levels=0)
ax[1].set_aspect('equal')
ax[1].set_xlim(-(x[-1]/2), (x[-1]/2))
ax[1].set_ylim(-(z[-1]/2), (z[-1]/2))
#ax[1].set_title("Isovalue level = 0")
ax[1].xaxis.set_tick_params(labelbottom=False)
ax[1].yaxis.set_tick_params(labelleft=False)

# Hide X and Y axes tick marks
ax[1].set_xticks([])
ax[1].set_yticks([])
#ax[0].axis('off')

fig.tight_layout()
plt.savefig('Images/AgContouringComparison.svg',transparent=True,bbox_inches = "tight", format='svg')
plt.show()

#%%
pelk=cs.allsegs[1][-1]
pelk=pelk[pelk[:,0].argsort()]
plt.plot(pelk[:, 0], pelk[:, 1],'.', linewidth=2)
plt.show()    



plt.axes().set_aspect('equal')
plt.xlim(-np.pi/2, np.pi/2)
plt.ylim(-np.pi/2, np.pi/2)
plt.title('Comparison of the Parameterised Eigenvalues and Elk evaluated 2D Fermi surface isosurface')
plt.plot(pelk[:, 0]+abs(pelk[np.where(pelk==np.nanmin(pelk[:,1]))[0][0],0]), pelk[:, 1],'o', linewidth=1,c='blue',label='Elk')    
plt.plot(p[:, 0], p[:, 1],'.', linewidth=4,c='black',label='Eigenvalues')
#plt.scatter(d[:, 0], d[:, 1],c='orange',label='Spectral Function')
#plt.contour(x,z,eA,50,cmap='viridis',levels=0,label='Elk')
plt.legend()
plt.show()



plt.figure(figsize=(7, 4), dpi=200)
plt.axes().set_aspect('equal')
ax = plt.subplot(111)
#ax.spines.right.set_visible(False)
#ax.spines.top.set_visible(False)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.xaxis.label.set_color('black')
ax.yaxis.label.set_color('black')
ax.tick_params(axis='x', colors='black')
ax.tick_params(axis='y', colors='black')


#title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Noisy Spectral Function Transformation\n',c='k')
#plt.title('2D Comparison of the Fermi Surface Isosurface \n- Ag [001] Plane Normal -\nParameterised Eigenvalues vs Elk Fermi Surface',c='k')
plt.xlim(-np.pi/4, np.pi/4)
plt.ylim(-np.pi/4, np.pi/4)
plt.yticks(ticks=[-0.5,0,0.5])
plt.xticks(ticks=[-0.5,0,0.5])
#plt.locator_params(axis='x', nbins=3)
#plt.locator_params(axis='y', nbins=3)  
plt.xlabel(r'[100] $\left(\frac{2\pi}{a}\right)$')
plt.ylabel(r'[010] $\left(\frac{2\pi}{a}\right)$')
#cs=plt.contour(x,y,sfp,50,cmap='viridis')#,levels=0)
plt.scatter(x=(pelk[:, 0]+abs(pelk[np.where(pelk==np.nanmin(pelk[:,1]))[0][0],0])), y=pelk[:, 1],s=1,c='blue',label='Elk')    
plt.scatter(x=p[:, 0], y=p[:, 1],s=0.5,c='black',label='Parameterisation')
plt.legend(loc='lower right',prop={'size': 6.5})

#savefig('Huesler111.png', transparent=True,bbox_inches = "tight")
#plt.savefig('Images/Ag2DParamVsElk.svg',transparent=True,bbox_inches = "tight", format='svg')
show()


# stop
#%%

phi = 0*(m.pi/180)
theta = 45*(m.pi/180)
psi = 0*(m.pi/180)
print("phi =", phi)
print("theta  =", theta)
print("psi =", psi)
  
ik=0
plane=np.zeros([1*npoints**2,3])
for j in range(npoints):
    for k in range(npoints):
        #for p in range(npoints):
        plane[ik,0]=k_x[j]
        plane[ik,1]=k_y[k]
        ik+=1

R = Rz(psi) * Rx(theta) * Rz(phi)
print(np.round(R, decimals=2))

axis=np.zeros([npoints**2,3])


FermiSurface=np.zeros([npoints**2])

for ik in range(npoints**2):
    axis[ik,:]=np.matmul(R,plane[ik,:])

x=axis[:,0].reshape([npoints,npoints])
y=axis[:,1].reshape([npoints,npoints])
z=axis[:,2].reshape([npoints,npoints])

x2=x
y2=y
z2=z

plane2[:,0],plane2[:,1],plane2[:,2]=BZCheck(rcell,axis[:,0],axis[:,1],axis[:,2],1)


for i in range(len(plane2)):
        FermiSurface[i]=FS(plane2[i,0],plane2[i,1],plane2[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)

Efzero2=np.nanmin(np.abs(FermiSurface))

# p2=transpose(np.matmul(R,transpose(p)))
# FermiSurface011=np.zeros([len(p)])
# for i in range(len(p)): 
#     FermiSurface011[i]=FS(p2[i,0],p2[i,1],p2[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)

FermiSurface=FermiSurface.reshape([npoints,npoints])

sfp=SF(FermiSurface,-0,0,-0.06)
s=SFFSCont(x,y,sfp,skipcon=6)  #skipcon=6
#d2,nk=outliers(s,pointskips=0,tol=0.95)
d2=s

plt.axes().set_aspect('equal')
plt.title(' [110] Theory Converted into Spectral Function')
plt.xlim(-np.pi/2, np.pi/2)
plt.ylim(-np.pi/2, np.pi/2)
plt.xlabel('100')
plt.ylabel('010')
cs=plt.contour(x,y,sfp,50,cmap='viridis')
#plt.xlim([lq,l3q])
#plt.ylim([lq,l3q])
plt.colorbar()
plt.show()


# plt.plot(s[:, 0], s[:, 1],'.', linewidth=2)
# plt.show()    

plt.plot(d2[:, 0], d2[:, 1],'.', linewidth=2)
plt.show()   

d2=zdata(d2,R)

#%%
# =============================================================================
# 
# d2=mirrorfold(s,0,0)
# 
# plt.rcParams['figure.figsize'] = [10,8]
# plt.plot(d2[:, 0], d2[:, 1],'.', linewidth=1,label='Comb.',)
# plt.plot(s[:, 0], s[:, 1],'.', linewidth=2,label='Full Orig',)
# plt.legend()
# plt.show()   
# 
# d2=kdefilter(d2,0.1)
# 
# d2=zdata(d2,R)
# =============================================================================
#%%


plt.axes().set_aspect('equal')
plt.title(' [110] Theory')
plt.xlim(-np.pi/2, np.pi/2)
plt.ylim(-np.pi/2, np.pi/2)
plt.xlabel('100')
plt.ylabel('010')
cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)
#plt.xlim([lq,l3q])
#plt.ylim([lq,l3q])
plt.colorbar()
plt.show()
    
p2=np.concatenate(cs.allsegs[1])
p2=p2[p2[:,0].argsort()]
plt.plot(p2[:, 0], p2[:, 1],'.', linewidth=2)
plt.show()    

p2=zdata(p2,R)

# FermiSurface011=np.zeros([len(p2)])
# for i in range(len(p2)): 
#     FermiSurface011[i]=FS(p2[i,0],p2[i,1],p2[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)

plt.figure(figsize=(5, 3), dpi=200)
plt.axes().set_aspect('equal')
ax = plt.subplot(111)
#ax.spines.right.set_visible(False)
#ax.spines.top.set_visible(False)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.xaxis.label.set_color('black')
ax.yaxis.label.set_color('black')
ax.tick_params(axis='x', colors='black')
ax.tick_params(axis='y', colors='black')


#title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Noisy Spectral Function Transformation\n',c='k')
#plt.title('2D Parameterisation of the Fermi Surface \n- Ag [110] Plane Normal -\n Isosurface',c='k')
plt.xlim(-np.pi/4, np.pi/4)
plt.ylim(-np.pi/4, np.pi/4)
plt.yticks(ticks=[-0.5,0,0.5])
plt.xticks(ticks=[-0.5,0,0.5])
#plt.locator_params(axis='x', nbins=3)
#plt.locator_params(axis='y', nbins=3)  
plt.xlabel(r'[0$\overline{1}$1] $\left(\frac{2\pi}{a}\right)$')
plt.ylabel(r'[100] $\left(\frac{2\pi}{a}\right)$')
#cs=plt.contour(x,y,sfp,50,cmap='viridis')#,levels=0)
cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)

#savefig('Huesler111.png', transparent=True,bbox_inches = "tight")
#plt.savefig('Images/Ag2DParam[110].svg',transparent=True,bbox_inches = "tight", format='svg')
show()


plt.figure(figsize=(5, 3), dpi=200)
plt.axes().set_aspect('equal')
ax = plt.subplot(111)
#ax.spines.right.set_visible(False)
#ax.spines.top.set_visible(False)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.xaxis.label.set_color('black')
ax.yaxis.label.set_color('black')
ax.tick_params(axis='x', colors='black')
ax.tick_params(axis='y', colors='black')


#title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Noisy Spectral Function Transformation\n',c='k')
#plt.title('2D Parameterisation of the Fermi Surface \n- Ag [110] Plane Normal -\n Isosurface (Points)',c='k')
plt.xlim(-np.pi/4, np.pi/4)
plt.ylim(-np.pi/4, np.pi/4)
plt.yticks(ticks=[-0.5,0,0.5])
plt.xticks(ticks=[-0.5,0,0.5])
#plt.locator_params(axis='x', nbins=3)
#plt.locator_params(axis='y', nbins=3)  
plt.xlabel(r'[0$\overline{1}$1] $\left(\frac{2\pi}{a}\right)$')
plt.ylabel(r'[100] $\left(\frac{2\pi}{a}\right)$')
#cs=plt.contour(x,y,sfp,50,cmap='viridis')#,levels=0)
#cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)
plt.plot(p2[:, 0], p2[:, 1],'.', markersize=0.5)
#savefig('Huesler111.png', transparent=True,bbox_inches = "tight")
#plt.savefig('Images/Ag2DParam[110]Points.svg',transparent=True,bbox_inches = "tight", format='svg')
show()

#%%

phi = 0*(m.pi/180)
theta = 54.74*(m.pi/180)
psi = 45*(m.pi/180)
print("phi =", phi)
print("theta  =", theta)
print("psi =", psi)
  
ik=0
plane=np.zeros([1*npoints**2,3])
for j in range(npoints):
    for k in range(npoints):
        #for p in range(npoints):
        plane[ik,0]=k_x[j]
        plane[ik,1]=k_y[k]
        ik+=1
  
R = Rz(psi) * Rx(theta) * Rz(phi)
print(np.round(R, decimals=2))

axis=np.zeros([npoints**2,3])


FermiSurface=np.zeros([npoints**2])

for ik in range(npoints**2):
    axis[ik,:]=np.matmul(R,plane[ik,:])
    
x=axis[:,0].reshape([npoints,npoints])
y=axis[:,1].reshape([npoints,npoints])
z=axis[:,2].reshape([npoints,npoints])    

plane2[:,0],plane2[:,1],plane2[:,2]=BZCheck(rcell,axis[:,0],axis[:,1],axis[:,2],1)


for i in range(len(plane2)):
        FermiSurface[i]=FS(plane2[i,0],plane2[i,1],plane2[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)

Efzero3=np.nanmin(np.abs(FermiSurface))


# p2=transpose(np.matmul(R,transpose(p)))
# FermiSurface011=np.zeros([len(p)])
# for i in range(len(p)): 
#     FermiSurface011[i]=FS(p2[i,0],p2[i,1],p2[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)

FermiSurface=FermiSurface.reshape([npoints,npoints])

sfp=SF(FermiSurface,-0,0,-0.06)
s=SFFSCont(x,y,sfp,skipcon=6)  #skipcon=6
#d3,nk=outliers(s,pointskips=0,tol=0.95)
d3=s

# plt.axes().set_aspect('equal')
# plt.title(' [110] Theory Converted into Spectral Function')
# plt.xlim(-np.pi/2, np.pi/2)
# plt.ylim(-np.pi/2, np.pi/2)
# plt.xlabel('100')
# plt.ylabel('010')
# cs=plt.contour(x,y,sfp,50,cmap='viridis')
# #plt.xlim([lq,l3q])
# #plt.ylim([lq,l3q])
# plt.colorbar()
# plt.show()


# plt.plot(s[:, 0], s[:, 1],'.', linewidth=2)
# plt.show()    

plt.plot(d3[:, 0], d3[:, 1],'.', linewidth=2)
plt.show() 

d3=zdata(d3,R)
#%%

# =============================================================================
# d3=rotate([0,0], s, 45)
# d3_2=d3
# d3=mirrorfold(d3,0,0)
# 
# #plt.rcParams['figure.figsize'] = [10,8]
# plt.plot(d3[:, 0], d3[:, 1],'.', linewidth=1,label='Comb.',)
# plt.plot(s[:, 0], s[:, 1],'.', linewidth=2,label='Full Orig',)
# #plt.plot(d3[:, 0], d3[:, 1],'.', linewidth=2,label='t',)
# plt.legend()
# plt.show()   
# 
# d3=rotate([0,0], d3, -45)
# 
# 
# d3=kdefilter(d3,0.09)
#
# d3=zdata(d3,R)
# =============================================================================


#%%
plt.axes().set_aspect('equal')
plt.title(' [111] Theory')
plt.xlim(-np.pi/2, np.pi/2)
plt.ylim(-np.pi/2, np.pi/2)
plt.xlabel('110')
plt.ylabel('011')
cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)
#plt.xlim([lq,l3q])
#plt.ylim([lq,l3q])
plt.colorbar()
plt.show()
    

p3=np.concatenate(cs.allsegs[1]) #max(cs.allsegs[1], key=len)
p3=p3[p3[:,0].argsort()]
#p3=rotate([0,0], p3[:,0:2], 45)
plt.plot(p3[:, 0], p3[:, 1],'.', linewidth=2)
plt.show()    


p3=zdata(p3,R)
  
    
plt.figure(figsize=(5, 3), dpi=200)
plt.axes().set_aspect('equal')
ax = plt.subplot(111)
#ax.spines.right.set_visible(False)
#ax.spines.top.set_visible(False)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.xaxis.label.set_color('black')
ax.yaxis.label.set_color('black')
ax.tick_params(axis='x', colors='black')
ax.tick_params(axis='y', colors='black')


#title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Noisy Spectral Function Transformation\n',c='k')
#plt.title('2D Parameterisation of the Fermi Surface \n- Ag [111] Plane Normal -\n Isosurface',c='k')
plt.xlim(-np.pi/4, np.pi/4)
plt.ylim(-np.pi/4, np.pi/4)
plt.yticks(ticks=[-0.5,0,0.5])
plt.xticks(ticks=[-0.5,0,0.5])
#plt.locator_params(axis='x', nbins=3)
#plt.locator_params(axis='y', nbins=3)  
plt.xlabel(r'[0$\overline{1}$1] $\left(\frac{2\pi}{a}\right)$')
plt.ylabel(r'[1$\overline{1}$0] $\left(\frac{2\pi}{a}\right)$')
#cs=plt.contour(x,y,sfp,50,cmap='viridis')#,levels=0)
cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)

#savefig('Huesler111.png', transparent=True,bbox_inches = "tight")
#plt.savefig('Images/Ag2DParam[111].svg',transparent=True,bbox_inches = "tight", format='svg')
show()


plt.figure(figsize=(5, 3), dpi=200)
plt.axes().set_aspect('equal')
ax = plt.subplot(111)
#ax.spines.right.set_visible(False)
#ax.spines.top.set_visible(False)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.xaxis.label.set_color('black')
ax.yaxis.label.set_color('black')
ax.tick_params(axis='x', colors='black')
ax.tick_params(axis='y', colors='black')


#title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Noisy Spectral Function Transformation\n',c='k')
#plt.title('2D Parameterisation of the Fermi Surface \n- Ag [111] Plane Normal -\n Isosurface (Points)',c='k')
plt.xlim(-np.pi/4, np.pi/4)
plt.ylim(-np.pi/4, np.pi/4)
plt.yticks(ticks=[-0.5,0,0.5])
plt.xticks(ticks=[-0.5,0,0.5])
#plt.locator_params(axis='x', nbins=3)
#plt.locator_params(axis='y', nbins=3)  
plt.xlabel(r'[0$\overline{1}$1] $\left(\frac{2\pi}{a}\right)$')
plt.ylabel(r'[1$\overline{1}$0] $\left(\frac{2\pi}{a}\right)$')
#cs=plt.contour(x,y,sfp,50,cmap='viridis')#,levels=0)
#cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)
plt.plot(p3[:, 0], p3[:, 1],'.', markersize=0.5)
#savefig('Huesler111.png', transparent=True,bbox_inches = "tight")
#plt.savefig('Images/Ag2DParam[111]Points.svg',transparent=True,bbox_inches = "tight", format='svg')
show()

# stop

#%%

#ax.imshow(FermiSurface, cmap=plt.cm.gray)
# cont=measure.find_contours(FermiSurface, 0.0)
# cont=max(cont,key=len)
# cont=np.round(cont)

# FermSing=np.zeros([npoints,npoints])
# FermSing[:] = np.nan
# x3d=np.zeros([len(cont),1])
# y3d=np.zeros([len(cont),1])
# z3d=np.zeros([len(cont),1])

# for i in range(len(cont)):
#     FermSing[int(cont[i,0]),int(cont[i,1])]=FermiSurface[int(cont[i,0]),int(cont[i,1])]
#     x3d[i]=x[int(cont[i,0]),0]
#     y3d[i]=y[0,int(cont[i,1])]
#     z3d[i]=FermiSurface[int(cont[i,0]),int(cont[i,1])]
    
    
# b=0.5

# # cell = np.array([[0.0, b, b],
# #                  [b, 0.0, b],
# #                  [b, b, 0.0]])

# # rcell=reclatvec(cell*7.7289957)

# bz=BZCalc(cell,0.81293684)

# fig = plt.figure(figsize=(12, 12))
# ax = fig.add_subplot(projection='3d')
# ax.set_xlim(-np.pi/2,np.pi/2) 
# ax.set_ylim(-np.pi/2,np.pi/2) 
# ax.set_zlim(-np.pi/2,np.pi/2)

# for i in range(len(bz[1])):
#     ax.plot(bz[1][i][:,0],bz[1][i][:,1],bz[1][i][:,2])


# plt.show()

#%%
# p=d.copy()
# p2=d2.copy()
# p3=d3.copy()

psym=[]
p2sym=[]
p3sym=[]

# p_hld=p.copy()
# p2_hld=p2.copy()
# p3_hld=p3.copy()

# p[:, [0, 2]] = p[:, [2, 0]]
# p2[:, [0, 2]] = p2[:, [2, 0]]
# p3[:, [0, 2]] = p3[:, [2, 0]]

o=0.5


cell = np.array([[0.0, o, o],
                  [o, 0.0, o],
                  [o, o, 0.0]])



rcell=reclatvec(cell*7.7289957)
rcell=rcell#/abs(np.min(rcell))

myLine = ""
for line in (open("SYMLAT.OUT", "r").read()).splitlines():
    if ": nsymlat" in line:
        myLine = line
        break

symnum=int(re.findall("-?\d+", myLine)[0])
symlist=list((np.array(range(0,symnum,1))*5)+3)

#U=rcell2/abs(rcell2[0,0])
#O=np.linalg.inv(U)



for k in range(len(symlist)):
    j=symlist[k]
    with open("SYMLAT.OUT") as fp:
        for i, line in enumerate(fp):
            if i == j:
                a=np.array([list(map(int,list(re.findall("-?\d+",line)))),list(map(int,list(re.findall("-?\d+",fp.readline())))),list(map(int,list(re.findall("-?\d+",fp.readline()))))])
                #b=np.matrix(lat2cart(a,revec))
                
                #a[:, [0, 2]] = a[:, [2, 0]]
                b=np.matrix(a) #np.cross(np.cross(O.T,a/0.5),O)
                print(b)
            
                
                for l in range(len(p)-1):
                    n=np.matrix(p[l,:])
                    n=np.matrix(rcell)*np.matrix(n).T  #cart2lat(p[l,0],p[l,1],p[l,2],rcell)[0]
                    n=(b)*n
                    n=np.linalg.inv(rcell)*n #n=np.linalg.inv(cell)*n #lat2cart(n.T,rcell)
                    psym.append((np.array(n)).T)

                    
                    
                for l in range(len(p2)-1):
                    n=np.matrix(p2[l,:])
                    n=np.matrix(rcell)*np.matrix(n).T  #cart2lat(p[l,0],p[l,1],p[l,2],rcell)[0]
                    n=(b)*n
                    n=np.linalg.inv(rcell)*n #n=np.linalg.inv(cell)*n #lat2cart(n.T,rcell)
                    p2sym.append((np.array(n)).T)
                    
                    
                    
                for l in range(len(p3)-1):
                    n=np.matrix(p3[l,:])
                    n=np.matrix(rcell)*np.matrix(n).T  #cart2lat(p[l,0],p[l,1],p[l,2],rcell)[0]
                    n=(b)*n
                    n=np.linalg.inv(rcell)*n #n=np.linalg.inv(cell)*n #lat2cart(n.T,rcell)
                    p3sym.append((np.array(n)).T)

                
psym=np.vstack(psym)
p2sym=np.vstack(p2sym)
p3sym=np.vstack(p3sym)


print(p[0],'\n',psym[105])

# p=p_hld.copy()
# p2=p2_hld.copy()
# p3=p3_hld.copy()




#%%
# %matplotlib qt

# o=0.5

# # cell2 = np.array([[b, b, 0.0],     #Avec used in Elk
# #                   [b, -b, 0.0],
# #                   [b, 0.0, b]])

# # rcell2=reclatvec(cell2*7.7289957)


# cell = np.array([[0.0, o, o],
#                   [o, 0.0, o],
#                   [o, o, 0.0]])



# rcell=reclatvec(cell*7.7289957)
# rcell=rcell#/abs(np.min(rcell))

# linr=np.array([[0,0,0],[1,0,0],[0,0,0],[0,1,0],[0,0,0],[0,0,1]])
# linr2=b*np.matrix([[0,0,0],[1,0,0],[0,0,0],[0,1,0],[0,0,0],[0,0,1]]).T
# linr2=np.array(linr2.T)

# normfac=np.matrix([np.sqrt(sum(np.array(b)[:, 0]**2)),np.sqrt(sum(np.array(b)[:, 1]**2)),np.sqrt(sum(np.array(b)[:, 2]**2))])
# #normfacrot=np.matrix([np.sqrt(sum(np.array(rot)[:, 0]**2)),np.sqrt(sum(np.array(rot)[:, 1]**2)),np.sqrt(sum(np.array(rot)[:, 2]**2))])
# #normfac=np.sqrt(np.array(b)[:, 0]**2 + np.array(b)[:, 1]**2 + np.array(b)[:, 2]**2)

# b2=np.matrix([ [-1, 0, 0], [0, 0, -1], [1, 1, 1]]) #np.matrix([ [2/3, 1/3, 2/3], [1/3, 2/3, -2/3], [-2/3, 2/3, 1/3]]) #np.matrix([ [0, -1, 0], [-1, 0, 0], [1, 1, 1]])

# rot=np.matrix(p2.copy())
# #tst=np.sqrt(np.array(rot)[:, 0]**2 + np.array(rot)[:, 1]**2 + np.array(rot)[:, 2]**2)
# rot=np.matrix(rcell)*np.matrix(rot).T
# rot=b2*rot
# rot=np.matrix(np.linalg.inv(rcell))*rot
# rot=np.array(rot.T)
# #tst2=np.sqrt(np.array(rot)[:, 0]**2 + np.array(rot)[:, 1]**2 + np.array(rot)[:, 2]**2)
# #rot=rot*np.max(tst)/np.max(tst2)
# # rot=np.matrix(rcell)*np.matrix(rot).T  #cart2lat(p[l,0],p[l,1],p[l,2],rcell)[0]
# # rot=(b2)*rot
# # rot=np.linalg.inv(rcell)*rot #n=np.linalg.inv(cell)*n #lat2cart(n.T,rcell)
# # rot=np.array(rot.T)


fig = plt.figure(figsize=(10, 10))#,layout='constrained')
ax = fig.add_subplot(projection='3d')
fig.tight_layout()
ax.set_box_aspect(aspect=None, zoom=0.97)
ax.set_aspect('equal')
#fig.subplots_adjust(top=10)
#plt.tight_layout()

# make the panes transparent
ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
# make the grid lines transparent
ax.xaxis._axinfo["grid"]['color'] =  (1,1,1,0)
ax.yaxis._axinfo["grid"]['color'] =  (1,1,1,0)
ax.zaxis._axinfo["grid"]['color'] =  (1,1,1,0)
# Get rid of the spines
#ax.xaxis.line.set_color((1.0, 1.0, 1.0, 1.0))
#ax.yaxis.line.set_color((1.0, 1.0, 1.0, 1.0))
#ax.zaxis.line.set_color((1.0, 1.0, 1.0, 1.0))
# Get rid of the ticks
ax.set_xticks([]) 
ax.set_yticks([]) 
ax.set_zticks([])

# ax.set_xlim(-np.pi/4,np.pi/4) 
# ax.set_ylim(-np.pi/4,np.pi/4) 
# ax.set_zlim(-np.pi/4,np.pi/4)
ax.set_xlabel(r'(100)', fontsize = 35)#,labelpad=20) #$\left(\frac{2\pi}{a}\right)$
ax.set_ylabel(r'(010)', fontsize = 35)#,labelpad=20)
ax.set_zlabel(r'(001)', fontsize = 35)#,labelpad=20)

ax.set_xlim3d(-0.5, 0.5)
ax.set_ylim3d(-0.5, 0.5)
ax.set_zlim3d(-0.5, 0.5)

#ax.set_title('3D Fermi Surface Parameterisation Contours \n- Ag Originally [001] -\n Symmetry Rotations Applied',c='k',fontsize = 40)
ax.scatter(psym[:,0],psym[:,1],psym[:,2],label='(001)')#='Rotated Contours')
#ax.scatter(p3[:,2],p3[:,1],p3[:,0],c='red',s=100,label='Original Contour')
# #ax.scatter(rot[:,2],rot[:,1],rot[:,0])
# # ax.plot(linr[:,0],linr[:,1],linr[:,2])
# # ax.plot(linr2[:,0],linr2[:,1],linr2[:,2])



# #ax.scatter(rot[:,0],rot[:,1],rot[:,2])
# #ax.plot(linr2[:,0],linr2[:,1],linr2[:,2])
# #ax.plot(linr[:,0],linr[:,1],linr[:,2])
# # #ax.scatter(axis[:,0],axis[:,1],FermiSurface)
ax.scatter(p2sym[:,0],p2sym[:,1],p2sym[:,2],label='(110)')
# #ax.scatter(p2[:,0],p2[:,1],p2[:,2])
# # ax.scatter(p2[:,1],p2[:,0],p2[:,2])
# # ax.scatter(p2[:,1],-p2[:,0],-p2[:,2])
#ax.scatter(p2[:,2],p2[:,1],p2[:,0])
#ax.scatter(p2[:,2],-p2[:,1],-p2[:,0])

ax.scatter(p3sym[:,0],p3sym[:,1],p3sym[:,2],label='(111)')
# # ax.scatter(p3[:,1],p3[:,0],p3[:,2])


# #ax.scatter(ptest[:,0],ptest[:,1],ptest[:,2])

# for i in range(len(bz[1])):
#       ax.plot(bz[1][i][:,0],bz[1][i][:,1],bz[1][i][:,2])
#       ax.plot(bz[2][i][:,0],bz[2][i][:,1],bz[2][i][:,2])
#ax.view_init(elev=25,azim=45)
#pylab.legend(loc='upper right',prop={'size': 25},markerscale=3,bbox_to_anchor=(1.2,1))
plt.savefig('Images/Ag3DSymAll.svg',transparent=True,bbox_inches = "tight", format='svg')
plt.show()



# plt.figure(figsize=(5, 3), dpi=200)
# plt.axes().set_aspect('equal')
# ax = plt.subplot(111)
# #ax.spines.right.set_visible(False)
# #ax.spines.top.set_visible(False)
# ax.spines['bottom'].set_color('black')
# ax.spines['left'].set_color('black')
# ax.xaxis.label.set_color('black')
# ax.yaxis.label.set_color('black')
# ax.tick_params(axis='x', colors='black')
# ax.tick_params(axis='y', colors='black')


# #title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Noisy Spectral Function Transformation\n',c='k')
# plt.title('2D Parameterisation of the Fermi Surface \n- Ag [110] Plane Normal -\n Isosurface',c='k')
# plt.xlim(-np.pi/2, np.pi/2)
# plt.ylim(-np.pi/2, np.pi/2)
# plt.xlabel(r'[100] $\left(\frac{2\pi}{a}\right)$')
# plt.ylabel(r'[011] $\left(\frac{2\pi}{a}\right)$')
# #cs=plt.contour(x,y,sfp,50,cmap='viridis')#,levels=0)
# cs=plt.contour(x,y,FermiSurface,50,cmap='viridis',levels=0)

# #savefig('Huesler111.png', transparent=True,bbox_inches = "tight")
# ##plt.savefig('Images/Ag2DParam[110].svg',transparent=True,bbox_inches = "tight", format='svg')
# show()

stop
#%matplotlib inline
#%%

# =============================================================================
# fig = plt.figure(figsize=(12, 12))
# ax = fig.add_subplot(projection='3d')
# ax.set_xlim(-np.pi/2,np.pi/2) 
# ax.set_ylim(-np.pi/2,np.pi/2) 
# ax.set_zlim(-np.pi/2,np.pi/2)
# 
# ax.scatter(d[:,0],d[:,1],d[:,2])
# ax.scatter(d[:,2],d[:,1],d[:,0])
# 
# #ax.scatter(axis[:,0],axis[:,1],FermiSurface)
# ax.scatter(d2[:,0],d2[:,1],d2[:,2])
# ax.scatter(d2[:,0],-d2[:,1],d2[:,2])
# 
# 
# ax.scatter(d3[:,0],d3[:,1],d3[:,2])
# ax.scatter(d3[:,1],d3[:,0],d3[:,2])
# 
# #ax.scatter(ptest[:,0],ptest[:,1],ptest[:,2])
# 
# # for i in range(len(bz[1])):
# #     ax.plot(bz[1][i][:,0],bz[1][i][:,1],bz[1][i][:,2])
# #     ax.plot(bz[2][i][:,0],bz[2][i][:,1],bz[2][i][:,2])
# 
# plt.show()
# =============================================================================

#%matplotlib inline
#%%
# p=transpose(np.array([list(p[:,0]),list(p[:,1]),list(p[:,2])]))
# p_2=transpose(np.array([list(p[:,2]),list(p[:,1]),list(p[:,0])]))

# p2_1=transpose(np.array([list(p2[:,0]),list(p2[:,1]),list(p2[:,2])]))
# p2_2=transpose(np.array([list(p2[:,0]),list(-p2[:,1]),list(p2[:,2])]))
# p2_3=transpose(np.array([list(p2[:,1]),list(p2[:,0]),list(p2[:,2])]))
# p2_4=transpose(np.array([list(p2[:,1]),list(-p2[:,0]),list(-p2[:,2])]))
# p2_5=transpose(np.array([list(p2[:,2]),list(p2[:,1]),list(p2[:,0])]))
# p2_6=transpose(np.array([list(p2[:,2]),list(-p2[:,1]),list(-p2[:,0])]))




# p3_1=transpose(np.array([list(p3[:,0]),list(p3[:,1]),list(p3[:,2])]))
# p3_2=transpose(np.array([list(p3[:,1]),list(p3[:,0]),list(p3[:,2])]))

#p3_2=transpose(np.array([list(p3[:,1]),list(p3[:,0]),list(p3[:,2])]))

p_t=transpose(np.array([list(psym[:,0]),list(psym[:,1]),list(psym[:,2])]))
p2_t=transpose(np.array([list(p2sym[:,0]),list(p2sym[:,1]),list(p2sym[:,2])]))
p3_t=transpose(np.array([list(p3sym[:,0]),list(p3sym[:,1]),list(p3sym[:,2])]))

#point_cloud=np.concatenate((p, p_2,p2_1,p2_2,p2_3,p2_4,p2_5,p2_6,p3_1,p3_2), axis=0)

point_cloud=np.concatenate((p_t, p2_t,p3_t), axis=0)
# point_cloud=np.concatenate((p_t, p2_t), axis=0) 
# point_cloud=p2_t.copy()


pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(point_cloud[:,:3])

# point_cloud2=point_cloud


N=np.zeros([len(point_cloud),3])
pcd2 = pcd #.voxel_down_sample(voxel_size=0.05)

b=0.5

cell = np.array([[0.0, b, b],
                  [b, 0.0, b],
                  [b, b, 0.0]])

rcell=reclatvec(cell*7.7289957)

b1=rcell[0,:]
b2=rcell[1,:]
b3=rcell[2,:]

#minfrac=calc_fraction_to_edge(point_cloud2, rcell)

rcellrg=reclatvec(cell*7.7289957*0.95)
point_cloud=np.asarray(pcd2.points)
point_cloudshft=point_cloud.copy()

for i in range(-1,2): 
    for j in range(-1,2): 
        for k in range(-1,2):
            if i == j == k == 0:
                  continue
            shifted = np.empty_like(point_cloud)            
            #shifted=np.array([point_cloud[:,0]+(i*b1[0]),point_cloud[:,1]+(j*b2[0]),point_cloud[:,2]+(k*b3[0])])
            for iv in range(3):
                shifted[:,iv]=point_cloud[:,iv]+i*b1[iv]+j*b2[iv]+k*b3[iv] 
            shifted[:, 0], shifted[:, 1], shifted[:, 2]=BZCheck(rcellrg,shifted[:,0],shifted[:,1],shifted[:,2],1)
            #shifted=np.reshape(shifted, (len(point_cloud),3))
            point_cloudshft=np.append(point_cloudshft,shifted,axis=0)
            # point_cloud2shft=BZCheck(rcellrg,point_cloud2shft[:,0],point_cloud2shft[:,1],point_cloud2shft[:,2],1)
            print(i,j,k)
            
# # create a boolean mask for rows containing only NaN values
mask = np.isnan(point_cloudshft).all(axis=1)
# # select the rows you want to keep
point_cloudshft = point_cloudshft[~mask]

pcd2 = o3d.geometry.PointCloud()
pcd2.points = o3d.utility.Vector3dVector(point_cloudshft[:,:3])


#o3d.visualization.draw([pcd2])


o3d.io.write_point_cloud('Ag3con.xyz', pcd2)

stop
#%%
# import os
# avgs=[]
# recon_dir = 'ParametricGaussRecon/results/AgSFtoit2/recon/'

# for file_name in os.listdir(recon_dir):
#     if not file_name.endswith('.ply'):
#         continue  # skip non-PLY files
    
#     file_path = os.path.join(recon_dir, file_name)
    
#     try:
#         mesh = o3d.io.read_triangle_mesh(file_path)
#         mesh.compute_vertex_normals()
    

#mesh = o3d.io.read_triangle_mesh('ParametricGaussRecon/results/AgSFtoit2/recon/AgSFtoit2_k_7_min_0.04_max_0.015_alpha_2.0_depth_min_1_depth_max_1_recon.ply')  #BEST USE WITH SMOOTHNING
mesh = o3d.io.read_triangle_mesh('ParametricGaussRecon/results/Ag1con/recon/Ag1con_k_7_min_0.03_max_0.015_alpha_2.0_depth_min_1_depth_max_1_recon.ply') 


#mesh = o3d.io.read_triangle_mesh('ParametricGaussRecon/results/AgSF/recon/AgSF_k_7_min_0.06_max_0.015_alpha_4.0_depth_min_1_depth_max_1_recon.ply')
mesh.compute_vertex_normals()
mesh = mesh.subdivide_midpoint(number_of_iterations=2)


mesh = mesh.filter_smooth_simple(number_of_iterations=10)
mesh = mesh.filter_smooth_laplacian(number_of_iterations=12) #THIS IS FOR MIN 0.04 WITH ALPHA 2
mesh.compute_vertex_normals()

#mesh = mesh.filter_smooth_taubin(number_of_iterations=2000)#filter_smooth_laplacian(number_of_iterations=10) #
#mesh = mesh.filter_smooth_simple(number_of_iterations=20)
#mesh.compute_vertex_normals()

#o3d.visualization.draw_geometries([mesh])

#%%
##################################################################################

# b=0.5*7.7289957

# cell = np.array([[0.0, b, b],
#                   [b, 0.0, b],
#                   [b, b, 0.0]])



# rcell=reclatvec(cell)

bz=BZCalc(cell,0.81293684*1) #np.max(p2)) 0.8129



count = 0
for listElem in bz[1]:
    count += len(listElem) 


points=np.zeros([count,3])
lines=np.array([(0,0)])

count = 0
for i in range(len(bz[1])):
    for j in range(len(bz[1][i])):
        points[count,:]=bz[1][i][j,:]
        if j == (len(bz[1][i])-1):
            lines=np.append(lines,[(count,count-j)],axis=0)
        else:
            lines=np.append(lines,[(count,count+1)],axis=0)
        count += 1

colors = [[0, 0, 1] for i in range(len(lines))]

line_set = o3d.geometry.LineSet(
    points=o3d.utility.Vector3dVector(points),
    lines=o3d.utility.Vector2iVector(lines),
)
line_set.colors = o3d.utility.Vector3dVector(colors)


#o3d.visualization.draw([{'name': 'difference', 'geometry': mesh}]+[pcd]+[line_set])

bzplane=BZPlanes(rcell*1, 1) #0.96

#o3d.visualization.draw([pcd2]+[line_set])


mesh2 = o3d.t.geometry.TriangleMesh.from_legacy(mesh)
#mesh2.scale(1.225, center=mesh.get_center())


bznorms = np.divide(bzplane, abs(bzplane), out=np.zeros_like(bzplane), where=bzplane!=0)
bznorms = bznorms*-1
for i in range(len(bzplane)):
    mesh2 = mesh2.clip_plane(point=bzplane[i,:], normal=bznorms[i,:])


# mesh2.scale(1/1.1, center=mesh.get_center())

# color1 = [1.0, 0.0, 0.0]  # red
# color2 = [0.0, 1.0, 0.0]  # green
# color3 = [0.0, 0.0, 1.0]  # blue

o3d.visualization.draw([line_set])

#o3d.visualization.draw([mesh2]+[line_set]+[pcd2])

#%%
################################################################################      

from scipy.spatial.transform import Rotation
from scipy.spatial.distance import cdist

def distance_stats(arr):
    distances = np.linalg.norm(arr, axis=1)
    avg_distance = np.mean(distances)
    max_distance = np.max(distances)
    min_distance = np.min(distances)
    return avg_distance, max_distance, min_distance

# Define your plane equation (point and normal)
plane_point = [0, 0, 0]
plane_normal = [0, 0, 1]

# necks=[]
# for i in range(len(bzplane)):
#     holder = mesh2.slice_plane(bzplane[i,:], bznorms[i,:])   
#     holder = o3d.t.geometry.LineSet.to_legacy(holder)
#     necks.append(np.asarray(holder.points))

# for i in range(2):
#     for arr in necks:
#         if len(arr) == 0:
#             necks.remove(arr)

mesho=o3d.t.geometry.TriangleMesh.from_legacy(mesh)

necks = []
ddis=[]
rdis=[]

for i in range(len(bzplane)):
    holder = mesho.slice_plane(bzplane[i, :], bznorms[i, :])
    #o3d.visualization.draw([mesho]+[holder])
    holder = o3d.t.geometry.LineSet.to_legacy(holder)
    if len(holder.points) == 0:
        continue
    holder=np.asarray(holder.points)
    
# =============================================================================
    center = np.mean(holder, axis=0)  # calculate the center of the circle
    distances = np.linalg.norm(holder - center, axis=1)  # calculate the distances to the center
    rdis.append(np.mean(distances)) 
    
    # Calculate the longest distance between any two points in holder
    dist_matrix = cdist(holder, holder)
    max_distance = np.max(dist_matrix)
    ddis.append(max_distance)
# =============================================================================
    
    holder = holder - np.mean(holder, axis=0)  # center the points around (0, 0, 0)
    normal = bznorms[i, :]
    z_axis = np.array([0, 0, 1])  # the target z-axis
    v = np.cross(z_axis, normal)  # the rotation axis
    c = np.dot(z_axis, normal)  # the cosine of the angle between z_axis and normal
    s = np.linalg.norm(v)  # the sine of the angle between z_axis and normal
    if s > 0:
        v /= s
    kmat = np.matrix([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])  # the cross-product matrix
    rmat = np.eye(3) + kmat + kmat @ kmat * (1 - c) / (s ** 2)  # the rotation matrix
    holder = np.asarray(np.matrix(holder) @ rmat)  # apply the rotation matrix
    phi = np.arctan2(rmat[1, 0], rmat[0, 0])
    theta = np.arccos(rmat[2, 2])
    psi = np.arctan2(-rmat[2, 1], rmat[2, 0])
    R = Rz(psi) * Rx(theta) * Rz(phi)  # calculate the rotation matrix to lie flat on the x-y plane
    holder = np.asarray(np.matrix(holder) @ R)  # apply the rotation matrix to lie flat on the x-y plane
    #holder[:, 2] = np.mean(holder[:, 2])  # set all z components to be the same
    
# =============================================================================
#     center = np.mean(holder, axis=0)  # calculate the center of the circle
#     distances = np.linalg.norm(holder - center, axis=1)  # calculate the distances to the center
#     rdis.append(np.max(distances)) 
#     
#     # Calculate the longest distance between any two points in holder
#     dist_matrix = cdist(holder, holder)
#     max_distance = np.max(dist_matrix)
#     ddis.append(max_distance)
# =============================================================================
    
    necks.append(holder)

dist_stats = [np.array(distance_stats(arr)) for arr in necks] 
overall_max = np.max([stats[1] for stats in dist_stats]) #we take the maximum because the exp will take the extremal
# overall_max = np.partition([stats[1] for stats in dist_stats], -4)[-4:]
# overall_max = np.mean(overall_max)

contours = mesh2.slice_plane([0,0,0], [0,0,1])   #, np.linspace(0,0.001))
#o3d.visualization.draw([mesh2]+[contours])
contours = o3d.t.geometry.LineSet.to_legacy(contours)
contours=np.asarray(contours.points)


# fig, ax = plt.subplots()
# ax.scatter(contours[:,0],contours[:,1])
# plt.show()

# fig = plt.figure(figsize=(12, 12))
# ax = fig.add_subplot(projection='3d')

# for cr in necks:
#     ax.scatter(cr[:,0],cr[:,1],cr[:,2])

# plt.show()

print(overall_max)
print(np.mean(rdis))


################################################################



r=np.mean(rdis) #overall_max #(np.linalg.norm(pointspk[2] - pointspk[3])+np.linalg.norm(pointspk[0] - pointspk[1])+np.linalg.norm(pointspk[4] - pointspk[5]))/3    #0.08856#

circx=np.linspace(-r,r,len(contours[:,0]))
circy=np.sqrt(r**2-circx**2)#*(circx/abs(circx))
circx=np.append(circx,circx,axis=0)
circy=np.append(circy,-circy,axis=0)

circle=circx.copy()
circle=np.c_[circle,circy]

# plt.scatter(circle[:,0],circle[:,1],s=0.2)
# plt.scatter(contours[:,0],contours[:,1])
# plt.show()
###################################################################        
import shapely
from shapely.geometry import Polygon, Point

# Define the shapes as arrays of points that define the contour of the shape
shape1_points = contours[:,0:2] #Belly
shape2_points = circle.copy()   #Neck


# Create polygons from the points
# shape1_poly = Polygon(shape1_points)
# shape2_poly = Polygon(shape2_points)


center = shape1_points.mean(0)
angle = np.arctan2(*(shape1_points - center).T[::-1])
index = np.argsort(angle)
shape1_poly = Polygon(shape1_points[index])


center = shape2_points.mean(0)
angle = np.arctan2(*(shape2_points - center).T[::-1])
index = np.argsort(angle)
shape2_poly = Polygon(shape2_points[index])

# points = shapely.MultiPoint(shape1_poly.exterior.coords[1:])
# convex_hull = points.convex_hull

# points2 = shapely.MultiPoint(shape2_poly.exterior.coords[1:])
# convex_hull2 = points2.convex_hull

# Determine the bounding box that encloses both shapes
x_min, x_max = np.min(np.concatenate([shape1_points[:, 0], shape2_points[:, 0]])) * 1.2, np.max(np.concatenate([shape1_points[:, 0], shape2_points[:, 0]])) * 1.2
y_min, y_max = np.min(np.concatenate([shape1_points[:, 1], shape2_points[:, 1]])) * 1.2, np.max(np.concatenate([shape1_points[:, 1], shape2_points[:, 1]])) * 1.2

# Generate a large number of random points within the bounding box
num_points = 300000
x_vals = np.random.uniform(x_min, x_max, num_points)
y_vals = np.random.uniform(y_min, y_max, num_points)

# Determine which points fall within each shape
mask1 = np.array([shape1_poly.contains(Point(x, y)) for x, y in zip(x_vals, y_vals)])
mask2 = np.array([shape2_poly.contains(Point(x, y)) for x, y in zip(x_vals, y_vals)])

# Calculate the fraction of points in each shape
frac1 = np.sum(mask1) / num_points
frac2 = np.sum(mask2) / num_points

print('Belly to neck MC',frac1/frac2)
print('Belly to neck polygon area',shape1_poly.area/shape2_poly.area)
#%%
#Plot the random points and the shapes

plt.figure(figsize=(8, 5), dpi=200)
plt.axes().set_aspect('equal')
ax = plt.subplot(111)
#ax.spines.right.set_visible(False)
#ax.spines.top.set_visible(False)
ax.spines['bottom'].set_color('black')
ax.spines['left'].set_color('black')
ax.xaxis.label.set_color('black')
ax.yaxis.label.set_color('black')
ax.tick_params(axis='x', colors='white')
ax.tick_params(axis='y', colors='white')


#title('2D Parameterisation of the Fermi Surface \n- Ag [001] Plane Normal -\n Noisy Spectral Function Transformation\n',c='k')
plt.title('Belly-to-Neck Ratio\nMonte-Carlo Method',c='k')
ax.set_xlim(x_min, x_max)
ax.set_ylim(y_min, y_max)
ax.plot(x_vals[mask1], y_vals[mask1], '.', label='Belly')
ax.plot(x_vals[mask2], y_vals[mask2], '.', label='Neck')
ax.plot(x_vals[~mask1 & ~mask2], y_vals[~mask1 & ~mask2], '.', label='Outside')
ax.plot(*shape1_poly.exterior.xy, 'k--', label='Belly Outline')
ax.plot(*shape2_poly.exterior.xy, 'r--', label='Neck Outline')
ax.legend(loc='lower right',prop={'size': 8},bbox_to_anchor=(1.35,0))
plt.savefig('Images/AgMonCrlo.svg',transparent=True,bbox_inches = "tight", format='svg')
plt.show()
#savefig('Huesler111.png', transparent=True,bbox_inches = "tight")
#show()
#%%



volume = ConvexHull(bz[0]).volume
#meshvolume=mesh.get_volume()
#print((meshvolume/volume)*100,'%')


mesh3=o3d.t.geometry.TriangleMesh.to_legacy(mesh2)
volume2 = ConvexHull(np.asarray(mesh3.vertices)).volume
print(((volume2/volume))*100,'%')

print((0.98/2)*100,'+-',(0.003/2)*100,'% - Exp')

# avgs.append(file_name)
# avgs.append(np.mean(rdis))
# avgs.append(frac1/frac2)
# avgs.append(shape1_poly.area/shape2_poly.area)
# avgs.append(((volume2/volume))*100)

    # except Exception as e:
    #     print(f"Skipping {file_name} due to error: {e}")
    #     continue

stop


#########################################################################################
#########################################################################################
#########################################################################################
#########################################################################################
#########################################################################################
#%%
import scipy.stats as stats

#Define Plane
npoints=64 #for 2*pi*a to 0 
plane=[1,1,0] #plane you want to observe
k_x = k_y = k_z =np.linspace(-np.pi/2,np.pi/2,npoints)
k_z=k_z*0

a=7.7289957 #4.0692 Angstrom

FermiSurface100=np.zeros([npoints**2,1])


ik=0
plane=np.zeros([1*npoints**2,3])
for j in range(npoints):
    for k in range(npoints):
        #for p in range(npoints):
        plane[ik,0]=k_x[j]
        plane[ik,1]=k_y[k]
        ik+=1
        
for i in range(len(plane)):
        FermiSurface100[i]=FS(plane[i,0],plane[i,1],plane[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)
        
mesh3=o3d.t.geometry.TriangleMesh.to_legacy(mesh2)

Efpoints=np.asarray((mesh3.sample_points_uniformly(number_of_points=100000)).points)



#Efpoints=np.asarray(mesh3.vertices)

indices = np.where(np.abs(Efpoints[:,2]) <= 0.001)[0]
closest_points = Efpoints[indices]
Efzero=np.mean([Efzero1,Efzero2,Efzero3])#np.min(np.abs(FermiSurface100))

FermiSurface=np.zeros([len(Efpoints),1])
chisq=np.zeros([len(Efpoints),1])
for i in range(len(Efpoints)):
        FermiSurface[i]=FS(Efpoints[i,0],Efpoints[i,1],Efpoints[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)
        chisq[i]=(((FermiSurface[i])-Efzero)**2)/Efzero

chisqsum=np.sum(chisq)

# FermiSurface=np.zeros([len(closest_points),1])
# chisq=np.zeros([len(closest_points),1])
# for i in range(len(closest_points)):
#         FermiSurface[i]=FS(closest_points[i,0],closest_points[i,1],closest_points[i,2],a,C_0,C_200,C_211,C_220,C_310,C_222,C_321)
#         chisq[i]=(((FermiSurface[i])-Efzero)**2)/Efzero

plt.scatter(plane[:,0],FermiSurface100.flatten())
plt.scatter(Efpoints[indices,0],FermiSurface.flatten()[indices])

#%%
# # Compute the normal vector of the plane defined by the points in array1
# v1 = plane[1] - plane[0]
# v2 = Efpoints[2] - Efpoints[0]
# normal = np.cross(v1, v2)
# normal = normal / np.linalg.norm(normal)

# # Compute the distance from each point in array2 to the plane defined by array1
# distances = np.abs(np.dot(Efpoints - plane[0], normal))

# # Specify a distance threshold
# threshold = 0.01

# # Find the indices of the points in array2 that are within the distance threshold from the plane
# indices = np.where(distances <= threshold)[0]
# closest_points = Efpoints[indices]



plt.scatter(range(len(FermiSurface)),FermiSurface,s=0.3)
plt.show()

plt.scatter(range(len(FermiSurface)),chisq,s=0.3)
plt.show()

print('The reduced Chi squared value is',np.sum(chisq)/len(chisq))
print('variance per point',np.sum(FermiSurface**2)/len(FermiSurface))

#Efzero=0

t_stat, p_value = stats.ttest_1samp(FermiSurface, Efzero)
print("t-statistic:", t_stat)
print("p-value:", p_value)

alpha = 0.05
t_critical = stats.t.ppf(1 - alpha / 2, df=len(FermiSurface) - 1)
margin_of_error = t_critical * np.std(FermiSurface, ddof=1) / np.sqrt(len(FermiSurface))
confidence_interval = (np.mean(FermiSurface) - margin_of_error, np.mean(FermiSurface) + margin_of_error)

print("Confidence interval:", confidence_interval)
print('The standard deviation is:',np.sqrt(np.sum((FermiSurface-np.mean(FermiSurface))**2)/len(FermiSurface)))


#######################################################################################

q1, q3 = np.percentile(FermiSurface, [25, 75])
iqr = q3 - q1

# Define the upper and lower bounds
lower_bound = q1 - (3 * iqr)
upper_bound = q3 + (3 * iqr)

# Identify the outliers
outliers = np.logical_or(FermiSurface < lower_bound, FermiSurface > upper_bound)

# Remove the outliers
FermiSurface_clean = FermiSurface[~outliers]

plt.scatter(range(len(FermiSurface)),FermiSurface,s=0.3)
plt.scatter(range(len(FermiSurface_clean)),FermiSurface_clean,s=0.3)
plt.show()

chisq=np.zeros([len(Efpoints),1])
for i in range(len(FermiSurface_clean)):
        chisq[i]=(((FermiSurface_clean[i])-Efzero)**2)/Efzero
        
print('The reduced Chi squared value is',np.sum(chisq)/len(chisq))

sample_size=len(FermiSurface_clean) #round(len(FermiSurface_clean)*0.9)
alpha = 0.05
subsample = np.random.choice(FermiSurface_clean, size=sample_size, replace=False)
t_stat, p_value = stats.ttest_1samp(subsample, Efzero, alternative='two-sided')

t_critical = stats.t.ppf(1 - alpha / 2, df=sample_size - 1)
margin_of_error = t_critical * np.std(subsample, ddof=1) / np.sqrt(sample_size)
confidence_interval = (np.mean(subsample) - margin_of_error, np.mean(subsample) + margin_of_error)

# Print the results
print("Sample size:", sample_size)
print("Data points removed:", len(FermiSurface)-sample_size)
print("t-statistic:", t_stat)
print("p-value:", p_value)
print("Confidence interval:", confidence_interval)

print('The standard deviation is:',np.sqrt(np.sum((FermiSurface_clean-np.mean(FermiSurface_clean))**2)/len(FermiSurface_clean)))
print('variance per point',np.sum(FermiSurface_clean**2)/len(FermiSurface_clean))



#%%

x,y,z,A=loadtxt('FERMISURF.OUT',skiprows=1,unpack=True)

v1=np.empty((1,3), float)
v2=np.empty((1,3), float)
ith=[]

vlst=np.array([x,y,z]).T

# fsdim=int(np.cbrt(len(x)))
# x=x.reshape([fsdim,fsdim])[0]
# y=y.reshape([fsdim,fsdim])#[0]
# z=z.reshape([fsdim,fsdim])[:,0]
# A=A.reshape([fsdim,fsdim])

for i in range(len(A)):
    v2=np.append(v2,np.array([[x[i],y[i],z[i]]]),axis=0) #v1.append(np.array([(x[i]),(y[i]),(z[i]),(A[i])]))
    if -0.001 <= A[i] <= 0.001:
        v1=np.append(v1,np.array([[x[i],y[i],z[i]]]),axis=0) #v1.append(np.array([(x[i]),(y[i]),(z[i]),(A[i])]))
        ith.append(i)
           
    
v1=np.delete(v1,0,axis=0)
v2=np.delete(v2,0,axis=0)
#%%
rlb=0.81293684
st=[-1,-1,-1]
st2=[-rlb,-rlb,-rlb]
# linr=np.array([st,cell[0,:],st,cell[1,:],st,cell[2,:]])

linr2=np.array([st2,rcell[0,:],st2,rcell[1,:],st2,rcell[2,:]])


box=np.array([st,[1,-1,-1],st,[-1,1,-1],st,[-1,-1,1]])*rlb
points=np.array([x,y,z]).T

#box=np.array((np.matrix(box)*np.matrix([rlb,rlb,rlb]).T).T)
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(projection='3d')
# ax.set_xlim(np.min(x),np.max(x)) 
# ax.set_ylim(np.min(y),np.max(y)) 
# ax.set_zlim(np.min(z),np.max(z))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')


ax.scatter(points[:,0], points[:,1], points[:,2], c=A, s=0.1)

ax.plot(box[:,0],box[:,1],box[:,2])
ax.plot(linr2[:,0],linr2[:,1],linr2[:,2])
# ax.plot(linr[:,0],linr[:,1],linr[:,2])

for i in range(len(bz[1])):
      ax.plot(bz[1][i][:,0],bz[1][i][:,1],bz[1][i][:,2])
      ax.plot(bz[2][i][:,0],bz[2][i][:,1],bz[2][i][:,2])

plt.show()


#%%
zind=np.where(abs(v2[:,2])<=0.01)[0]
plt.scatter(plane[:,0],FermiSurface100.flatten()/np.max(FermiSurface100),s=0.1)
plt.scatter(v2[zind,0],A.flatten()[zind]/np.max(A),s=0.2)

#%%
points=np.array([x,y,z]).T

pontsv2=BZCheck(rcell,points[:,0],points[:,1],points[:,2],1)

indx=np.where(points==np.array(pontsv2).T)[0]

indx=np.unique(indx)
#%%
interp=griddata(points[indx,:],A[indx],Efpoints)

#inte=np.reshape(interp, (int(np.sqrt(len(plane2))), int(np.sqrt(len(plane2)))))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(Efpoints[:,0], Efpoints[:,1], Efpoints[:,2],s=0.01)
ax.scatter(Efpoints[:,0], Efpoints[:,1], Efpoints[:,2], c=interp, s=0.1)
plt.show()

#%%

ster=np.sqrt((np.sum((interp-np.min(abs(A)))**2))/len(interp))
ster2=np.sqrt((np.sum((interp-0)**2))/len(interp))

from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error


variance = np.var(interp)
std_deviation = np.std(interp)

plt.hist(interp, bins=1000)
plt.axvline(x=np.min(abs(A)), color='r', linestyle='dashed', linewidth=2)
plt.show()



#########################################################################################
#########################################################################################
#########################################################################################
#########################################################################################
#########################################################################################

#%%
mesh3=o3d.t.geometry.TriangleMesh.to_legacy(mesh2)

vertices = np.asarray(mesh3.vertices)
triangles = np.asarray(mesh3.triangles)

polydata = vtk.vtkPolyData()
points = vtk.vtkPoints()
cells = vtk.vtkCellArray()

for vertex in vertices:
    points.InsertNextPoint(vertex)

for triangle in triangles:
    cells.InsertNextCell(3)
    for point_idx in triangle:
        cells.InsertCellPoint(point_idx)

polydata.SetPoints(points)
polydata.SetPolys(cells)


# Fill the holes on the mesh
fill_holes = vtk.vtkFillHolesFilter()
fill_holes.SetInputData(polydata)
fill_holes.Update()
filled_mesh = fill_holes.GetOutput()


# ==========================VISUALISATION===================================
# mapper = vtk.vtkPolyDataMapper()
# mapper.SetInputData(filled_mesh)  #where you input the mesh

# actor = vtk.vtkActor()
# actor.SetMapper(mapper)

# renderer = vtk.vtkRenderer()
# renderer.AddActor(actor)

# render_window = vtk.vtkRenderWindow()
# render_window.AddRenderer(renderer)

# interactor = vtk.vtkRenderWindowInteractor()
# interactor.SetRenderWindow(render_window)

# interactor.Start()
# =============================================================================

# Clean the mesh to make it watertight
cleaner = vtk.vtkCleanPolyData()
cleaner.SetInputData(filled_mesh)
cleaner.Update()
clean_mesh = cleaner.GetOutput()

mass_props = vtk.vtkMassProperties()
mass_props.SetInputData(clean_mesh)
mass_props.Update()

# Get the volume of the mesh
trimmedvolume = mass_props.GetVolume()

print((trimmedvolume/volume)*100,'%')

################################################################################

vertices = np.asarray(mesh.vertices)
triangles = np.asarray(mesh.triangles)

polydata = vtk.vtkPolyData()
points = vtk.vtkPoints()
cells = vtk.vtkCellArray()

for vertex in vertices:
    points.InsertNextPoint(vertex)

for triangle in triangles:
    cells.InsertNextCell(3)
    for point_idx in triangle:
        cells.InsertCellPoint(point_idx)

polydata.SetPoints(points)
polydata.SetPolys(cells)


# Fill the holes on the mesh
fill_holes = vtk.vtkFillHolesFilter()
fill_holes.SetInputData(polydata)
fill_holes.Update()
filled_mesh = fill_holes.GetOutput()

# Clean the mesh to make it watertight
cleaner = vtk.vtkCleanPolyData()
cleaner.SetInputData(filled_mesh)
cleaner.Update()
clean_mesh = cleaner.GetOutput()

mass_props = vtk.vtkMassProperties()
mass_props.SetInputData(clean_mesh)
mass_props.Update()

# Get the volume of the mesh
trimmedvolume2 = mass_props.GetVolume()

print((1-(trimmedvolume/trimmedvolume2))*100,'%')
#%%



#%%
################################################################################
################################################################################
################################################################################
################################################################################
#%%



def pick_points(pcd):
    print("")
    print(
        "1) Please pick at least three correspondences using [shift + left click]"
    )
    print("   Press [shift + right click] to undo point picking")
    print("2) After picking points, press 'Q' to close the window")
    vis = o3d.visualization.VisualizerWithEditing()
    vis.create_window()
    vis.add_geometry(pcd)
    vis.run()  # user picks points
    vis.destroy_window()
    print("")
    coords=[np.asarray(pcd.points)[i] for i in vis.get_picked_points()]
    return coords


pcd22 = o3d.geometry.PointCloud()
pcd22.points = o3d.utility.Vector3dVector(mesh2.vertex.positions.numpy())


pointspk=pick_points(pcd22)

#%%


def calculate_diameters_on_plane(mesh):
    points=mesh2.vertex.positions.numpy()
    
    intersections = []
    
    sorted_points = points[np.argsort(np.abs(points[:, 2]))]
    lowest_value = sorted_points[0, 2]
    intersections = sorted_points[sorted_points[:, 2] == lowest_value]

    # Calculate the maximum and minimum distances between all pairs of vertices
    mn = mx = []
    max_diameter = 0
    min_diameter = float('inf')
    for i in range(len(intersections)):
        d = np.linalg.norm(intersections[i])
        if d > max_diameter:
            max_diameter = d
            mx=i
        if d < min_diameter:
            min_diameter = d
            mn=i

    return max_diameter, min_diameter, intersections[mx], intersections[mn], np.array(intersections)


def neckpoints(mesh):
    points=mesh2.vertex.positions.numpy()
    
    # Find the intersection points
    intersections = []
    
    def euclidian_distance(a, b):
      return np.linalg.norm(a - b)
    
    points = np.vstack(sorted(points, key=lambda point: euclidian_distance(point,[0,0,0]),reverse=True))
    

    maxdis = np.linalg.norm(points[0])
    for i in range(len(points)):
        if np.round(np.linalg.norm(points[i,:]),decimals=2) == np.round(maxdis,decimals=2):
            intersections.append(points[i,:])
            
    intersections=np.array(intersections)
    
    return intersections

dis1,dis2,mxi,mni,pnts=calculate_diameters_on_plane(mesh2)

plt.scatter(pnts[:,0],pnts[:,1])
plt.scatter(mxi[0],mxi[1])
plt.scatter(mni[0],mni[1])
plt.show()

pstv=[]
for i in range(len(pnts)):
    if np.all(pnts[i,0:2]>=0):
        pstv.append(i)
        
bellyarea=(np.trapz(pnts[pstv,1],pnts[pstv,0]))*4

plt.scatter(pnts[pstv,0],pnts[pstv,1])
plt.show()

neckarea=np.pi*(np.linalg.norm(pointspk[2] - pointspk[3])+np.linalg.norm(pointspk[0] - pointspk[1]))/4**2

nk=neckpoints(mesh2)


nkv=[]
for i in range(len(nk)):
    if np.all(nk[i,:]>=0):
        nkv.append(i)

plt.scatter(nk[nkv,0],nk[nkv,1])
plt.show()

poly=nk[nkv,:]
#%%

def project_polygon(coords):
    # Calculate the normal vector of the plane defined by the polygon
    n = np.cross(coords[1] - coords[0], coords[2] - coords[0])
    n /= np.linalg.norm(n)

    # Calculate the rotation matrix that rotates the normal to [0, 0, 1]
    v = np.array([0, 0, 1])
    c = np.dot(n, v)
    if c < 1 - 1e-6:
        s = np.sqrt(1 - c**2)
        u = np.cross(n, v)
        ux, uy, uz = u
        R = np.array([[c + ux**2*(1-c), ux*uy*(1-c) - uz*s, ux*uz*(1-c) + uy*s],
                      [uy*ux*(1-c) + uz*s, c + uy**2*(1-c), uy*uz*(1-c) - ux*s],
                      [uz*ux*(1-c) - uy*s, uz*uy*(1-c) + ux*s, c + uz**2*(1-c)]])
    else:
        R = np.eye(3)

    # Apply the rotation
    rotated_coords = np.dot(R, coords.T).T
    
    # Project onto the x-y plane
    projected_coords = rotated_coords[:, :2]
    
    return projected_coords

def area_ratio(pnts, proj):
    area_pnts = 0.5 * abs(sum(pnts[i][0]*pnts[i+1][1] - pnts[i+1][0]*pnts[i][1] for i in range(-1, len(pnts)-1)))
    area_proj = 0.5 * abs(sum(proj[i][0]*proj[i+1][1] - proj[i+1][0]*proj[i][1] for i in range(-1, len(proj)-1)))
    if area_pnts > area_proj:
        return area_proj / area_pnts
    else:
        return area_pnts / area_proj
    
proj=project_polygon(nk[nkv,:])
proj = proj - np.mean(proj, axis=0)

plt.scatter(proj[:,0],proj[:,1])
plt.scatter(pnts[:,0],pnts[:,1])
plt.scatter(pnts[:,0]/29,pnts[:,1]/29)
plt.show()


rat=area_ratio(pnts, proj)
#%%

x,y,z,nx,ny,nz = np.loadtxt('point_cloudwd.xyz',skiprows=0,usecols=(0,1,2,3,4,5),unpack=True)
normals=transpose(np.array([list(nx),list(ny),list(nz)]))

for i in range(len(point_cloud)):
    norm=np.linalg.norm(normals[i,:])   #np.sqrt(sum(point_cloud[i,:]**2))
    N[i,0],N[i,1],N[i,2]=normals[i,:]/norm#


pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(point_cloud[:,:3])
pcd.normals = o3d.utility.Vector3dVector(N[:,:3])


o3d.visualization.draw_geometries([pcd],point_show_normal=True)
#%%

# =============================================================================
# p=transpose(np.array([list(d[:,0]),list(d[:,1]),list(d[:,2])]))
# p_2=transpose(np.array([list(d[:,2]),list(d[:,1]),list(d[:,0])]))
# 
# p2_1=transpose(np.array([list(d2[:,0]),list(d2[:,1]),list(d2[:,2])]))
# 
# p2_2=transpose(np.array([list(d2[:,0]),list(-d2[:,1]),list(d2[:,2])]))
# 
# p3_1=transpose(np.array([list(d3[:,0]),list(d3[:,1]),list(d3[:,2])]))
# 
# p3_2=transpose(np.array([list(d3[:,1]),list(d3[:,0]),list(d3[:,2])]))
# 
# #p3_2=transpose(np.array([list(p3[:,1]),list(p3[:,0]),list(p3[:,2])]))
# 
# point_cloud=np.concatenate((p, p_2,p2_1,p2_2,p3_1,p3_2), axis=0)
# 
# N=np.zeros([len(point_cloud),3])
# =============================================================================

stop
#%% DEBUG



# box = o3d.geometry.TriangleMesh.create_box()
# box = o3d.t.geometry.TriangleMesh.from_legacy(box)
# sphere = o3d.geometry.TriangleMesh.create_sphere(0.8)
# sphere = o3d.t.geometry.TriangleMesh.from_legacy(sphere)

# ans = box.boolean_difference(sphere)

# o3d.visualization.draw([{'name': 'difference', 'geometry': ans}])



# pcd = o3d.geometry.PointCloud()
# pcd.points = o3d.utility.Vector3dVector(point_cloud[:,:3])
# pcd.colors = o3d.utility.Vector3dVector(point_cloud[:,:3])
# pcd.estimate_normals()#search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=10, max_nn=1),fast_normal_computation=False)
# #pcd.orient_normals_to_align_with_direction()
# pcd.orient_normals_consistent_tangent_plane(2)



# o3d.visualization.draw_geometries([pcd],point_show_normal=True)

# o3d.visualization.draw([{'name': 'difference', 'geometry': ans}])


# o3d.visualization.draw_geometries([pcd],point_show_normal=True)


#%%
b=0.5

cell = np.array([[0.0, b, b],
                 [b, 0.0, b],
                 [b, b, 0.0]])

rcell=reclatvec(cell*7.7289957)

bz=BZCalc(cell,0.81293684*1) #np.max(p2)) 0.8129

count = 0
for listElem in bz[1]:
    count += len(listElem) 


points=np.zeros([count,3])
lines=np.array([(0,0)])

count = 0
for i in range(len(bz[1])):
    for j in range(len(bz[1][i])):
        points[count,:]=bz[1][i][j,:]
        if j == (len(bz[1][i])-1):
            lines=np.append(lines,[(count,count-j)],axis=0)
        else:
            lines=np.append(lines,[(count,count+1)],axis=0)
        count += 1


print("Let's draw a box using o3d.geometry.LineSet.")

colors = [[1, 0, 0] for i in range(len(lines))]

line_set = o3d.geometry.LineSet(
    points=o3d.utility.Vector3dVector(points),
    lines=o3d.utility.Vector2iVector(lines),
)
line_set.colors = o3d.utility.Vector3dVector(colors)


#point_cloud_bu=point_cloud
#point_cloud=point_cloud2
#point_cloud=point_cloud_bu

# pcd = o3d.geometry.PointCloud()
# pcd.points = o3d.utility.Vector3dVector(point_cloud[:,:3])
# pcd.colors = o3d.utility.Vector3dVector(point_cloud[:,:3])

# =============================================================================
# for i in range(len(point_cloud)):
#     norm=np.linalg.norm(point_cloud[i,:])   #np.sqrt(sum(point_cloud[i,:]**2))
#     N[i,0],N[i,1],N[i,2]=point_cloud[i,:]/norm#
# pcd.normals = o3d.utility.Vector3dVector(N[:,:3])
# 
# pcd = pcd.voxel_down_sample(voxel_size=0.07)
# #pcd, ind = pcd.remove_radius_outlier(nb_points=15, radius=0.07)
# pcd, ind = pcd.remove_statistical_outlier(16, 2, print_progress=False)
# #pcd = pcd.voxel_down_sample(voxel_size=0.3)
# 
# 
# o3d.visualization.draw([pcd.paint_uniform_color([1, 1, 1])],show_skybox=False, bg_color=(0.0, 0.0, 0.0, 1.0))
# 
# #pcd, ind = pcd.remove_statistical_outlier(50, 2, print_progress=False)
# pcd = pcd.voxel_down_sample(voxel_size=0.05)
# #pcd, ind = pcd.remove_statistical_outlier(25, 1.5, print_progress=False)
# #pcd = pcd.voxel_down_sample(voxel_size=0.03)
# #pcd = pcd.voxel_down_sample(voxel_size=0.1)
# 
# o3d.visualization.draw([pcd.paint_uniform_color([1, 1, 1])],show_skybox=False, bg_color=(0.0, 0.0, 0.0, 1.0))
# 
# pcd, ind = pcd.remove_statistical_outlier(20, 2, print_progress=False)
# 
# 
# o3d.visualization.draw([pcd.paint_uniform_color([1, 1, 1])],show_skybox=False, bg_color=(0.0, 0.0, 0.0, 1.0))
# #display_inlier_outlier(cl, ind)
# 
# 
# #pcd = pcd.voxel_down_sample(voxel_size=0.2)
# #pcd, ind = pcd.remove_statistical_outlier(5, 2, print_progress=False)
# 
# =============================================================================
# pcd.estimate_normals()  #search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=10, max_nn=1),fast_normal_computation=False)

# #o3d.visualization.draw([pcd])

# #pcd.orient_normals_to_align_with_direction()
# pcd.orient_normals_consistent_tangent_plane(2)
# pcd.colors = o3d.utility.Vector3dVector(point_cloud[:,:3])

o3d.visualization.draw([pcd]+[line_set])
#%%
b=0.5

cell = np.array([[0.0, b, b],
                 [b, 0.0, b],
                 [b, b, 0.0]])

rcell=reclatvec(cell*7.7289957)

bz=BZCalc(cell,0.81293684*1) #np.max(p2)) 0.8129

bzplane=BZPlanes(rcell, 1)

mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(pcd, scale=2, depth=15 ,linear_fit=False) #,depth=12

count = 0
for listElem in bz[1]:
    count += len(listElem) 


points=np.zeros([count,3])
lines=np.array([(0,0)])

count = 0
for i in range(len(bz[1])):
    for j in range(len(bz[1][i])):
        points[count,:]=bz[1][i][j,:]
        if j == (len(bz[1][i])-1):
            lines=np.append(lines,[(count,count-j)],axis=0)
        else:
            lines=np.append(lines,[(count,count+1)],axis=0)
        count += 1



o3d.visualization.draw([{'name': 'difference', 'geometry': mesh}]+[pcd]+[line_set])
#%%
bzplane=BZPlanes(rcell*0.98, 1) #0.96
mesh2 = o3d.t.geometry.TriangleMesh.from_legacy(mesh)
bznorms = np.divide(bzplane, abs(bzplane), out=np.zeros_like(bzplane), where=bzplane!=0)
bznorms = bznorms*-1
for i in range(len(bzplane)):
    mesh2 = mesh2.clip_plane(point=bzplane[i,:], normal=bznorms[i,:])

o3d.visualization.draw([mesh2]+[line_set])

#%%
tetra_mesh, pt_map = o3d.geometry.TetraMesh.create_from_point_cloud(pcd)

mesh_in = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(
    pcd, 0.7, tetra_mesh, pt_map)
mesh_in.compute_vertex_normals()
o3d.visualization.draw([mesh_in])

#%%

mesh_out = mesh_in.filter_smooth_taubin(number_of_iterations=1)
mesh_out.compute_vertex_normals()
o3d.visualization.draw([mesh_out])

#%%

bzplane=BZPlanes(rcell*0.99, 1) #0.96
mesh2 = o3d.t.geometry.TriangleMesh.from_legacy(mesh_in)
bznorms = np.divide(bzplane, abs(bzplane), out=np.zeros_like(bzplane), where=bzplane!=0)
bznorms = bznorms*-1
for i in range(len(bzplane)):
    mesh2 = mesh2.clip_plane(point=bzplane[i,:], normal=bznorms[i,:])

o3d.visualization.draw([mesh2]+[line_set])

#%%

mesh_out = mesh_in.filter_smooth_laplacian(number_of_iterations=1)
mesh_out.compute_vertex_normals()
o3d.visualization.draw([mesh_out])

#%%

print(
    f'The mesh has {len(mesh_in.vertices)} vertices and {len(mesh_in.triangles)} triangles'
)
o3d.visualization.draw([mesh_in])
mesh_in = mesh_in.subdivide_loop(number_of_iterations=2)
print(
    f'After subdivision it has {len(mesh_in.vertices)} vertices and {len(mesh_in.triangles)} triangles'
)
o3d.visualization.draw([mesh_in])

#%% RECONSTRUCTION

b=0.5

cell = np.array([[0.0, b, b],
                 [b, 0.0, b],
                 [b, b, 0.0]])

rcell=reclatvec(cell*7.7289957)

bz=BZCalc(cell,0.81293684*1) #np.max(p2)) 0.8129

count = 0
for listElem in bz[1]:
    count += len(listElem) 


points=np.zeros([count,3])
lines=np.array([(0,0)])

count = 0
for i in range(len(bz[1])):
    for j in range(len(bz[1][i])):
        points[count,:]=bz[1][i][j,:]
        if j == (len(bz[1][i])-1):
            lines=np.append(lines,[(count,count-j)],axis=0)
        else:
            lines=np.append(lines,[(count,count+1)],axis=0)
        count += 1


print("Let's draw a box using o3d.geometry.LineSet.")

colors = [[1, 0, 0] for i in range(len(lines))]

line_set = o3d.geometry.LineSet(
    points=o3d.utility.Vector3dVector(points),
    lines=o3d.utility.Vector2iVector(lines),
)
line_set.colors = o3d.utility.Vector3dVector(colors)

mesh_in = o3d.io.read_triangle_mesh('point_cloud_recon.ply')

bzplane=BZPlanes(rcell*0.99, 1) #0.96
mesh2 = o3d.t.geometry.TriangleMesh.from_legacy(mesh_in)
bznorms = np.divide(bzplane, abs(bzplane), out=np.zeros_like(bzplane), where=bzplane!=0)
bznorms = bznorms*-1
for i in range(len(bzplane)):
    mesh2 = mesh2.clip_plane(point=bzplane[i,:], normal=bznorms[i,:])

o3d.visualization.draw([mesh2]+[line_set])