#----------------------------------------------------------------------
#  OpenGFDM - utility functions
#  Author: Andrea Pavan
#  Date: 08/12/2022
#  License: GPLv3-or-later
#----------------------------------------------------------------------



"""
PARSESU2MESH reads a su2 mesh file and extracts the points (2d and 3d)
"""
function parseSU2mesh(filename::String)
    lines = readlines(filename);
    NDIME = parse(Int, lines[1][end]);
    pointcloud = ElasticArray{Float64}(undef,NDIME,0);
    parseline = false;
    for line in lines
        if startswith(line, "NMARK=") || startswith(line, "NELEM=")
            parseline = false;
        end
        if parseline
            lineNums = parse.(Float64,split(line,' '));
            append!(pointcloud, lineNums[1:NDIME]);
        end
        if startswith(line, "NPOIN=")
            parseline = true;
        end
    end
    return pointcloud;
end


"""
MEMORYUSAGE return the memory usage of one or two variables in a convenient measurement unit
"""
function memoryUsage(var1,var2=0,digits=2)
    mem = Base.summarysize(var1);       #in bytes
    if var2 != 0
        mem += Base.summarysize(var2);
    end
    units = ["KB","MB","GB","TB","PB","EB","ZB","YB","RB","QB"];
    idx = Int(floor(floor(log10(mem))/3));
    return string(round(mem/(1000^idx),digits=digits))*" "*units[idx];       #example: "1.23 MB"
end


"""
CARTESIANNEIGHBORSEARCH find the nearest nodes in a pointcloud using a cartesian subgrid
INPUT:
    pointcloud: 3xN Array containing the nodes coordinates [X;Y;Z]
    meshSize: Float64 target pointcloud distances
    minNeighbors: Int minimum number of neighbors required [default: 15]
OUTPUT:
    neighbors: Vector of N Vectors containing the indices of each node neighbors
    Nneighbors: Vector of size N containing the number of neighbors of each node
    cell: 3D tensor of Vectors containing the indices of each cell nodes
"""
function cartesianNeighborSearch(pointcloud,meshSize::Float64,minNeighbors::Int=15)
    #generate cartesian subgrid
    Pmax = maximum(pointcloud,dims=2);
    Pmin = minimum(pointcloud,dims=2);
    boundingBoxSize = Pmax-Pmin;
    Ncellx = max(1,Int(round(boundingBoxSize[1]/(3.5*meshSize))));
    Ncelly = max(1,Int(round(boundingBoxSize[2]/(3.5*meshSize))));
    Ncellz = max(1,Int(round(boundingBoxSize[3]/(3.5*meshSize))));
    cellsizex = boundingBoxSize[1]/Ncellx;
    cellsizey = boundingBoxSize[2]/Ncelly;
    cellsizez = boundingBoxSize[3]/Ncellz;
    cell = Array{Vector{Int},3}(undef,Ncellx,Ncelly,Ncellz);
    for i=1:Ncellx
        for j=1:Ncelly
            for k=1:Ncellz
                cell[i,j,k] = Vector{Int}(undef,0);
            end
        end
    end
    N = size(pointcloud,2);     #number of points
    #cellIdx = Vector{Vector{Int},N}(undef);     #vector containing the index of the cell
    cellIdx = Matrix{Int}(undef,3,N);
    for i=1:N
        idx1 = min(Ncellx,1+Int(floor((pointcloud[1,i]-Pmin[1])/cellsizex)));
        idx2 = min(Ncelly,1+Int(floor((pointcloud[2,i]-Pmin[2])/cellsizey)));
        idx3 = min(Ncellz,1+Int(floor((pointcloud[3,i]-Pmin[3])/cellsizez)));
        push!(cell[idx1,idx2,idx3], i);
        cellIdx[:,i] = [idx1,idx2,idx3];
    end

    #find neighbors using brute-force on a subset of cells
    neighbors = Vector{Vector{Int}}(undef,N);
    Nneighbors = Vector{Int}(undef,N);
    for i=1:N
        searchRadius = 0.6*meshSize;
        neighbors[i] = Int[];
        Nneighbors[i] = 0;
        while Nneighbors[i]<minNeighbors
            searchRadius += meshSize/2;
            #determine on which cells to search
            #=cellstosearch = [cellIdx[:,i],
                            [min(Ncellx,1+Int(floor(abs(pointcloud[1,i]+searchRadius-Pmin[1])/cellsizex))), min(Ncelly,1+Int(floor(abs(pointcloud[2,i]+searchRadius-Pmin[2])/cellsizey))), min(Ncellz,1+Int(floor(abs(pointcloud[3,i]+searchRadius-Pmin[3])/cellsizez)))],
                            [min(Ncellx,1+Int(floor(abs(pointcloud[1,i]+searchRadius-Pmin[1])/cellsizex))), min(Ncelly,1+Int(floor(abs(pointcloud[2,i]+searchRadius-Pmin[2])/cellsizey))), min(Ncellz,1+Int(floor(abs(pointcloud[3,i]-searchRadius-Pmin[3])/cellsizez)))],
                            [min(Ncellx,1+Int(floor(abs(pointcloud[1,i]+searchRadius-Pmin[1])/cellsizex))), min(Ncelly,1+Int(floor(abs(pointcloud[2,i]-searchRadius-Pmin[2])/cellsizey))), min(Ncellz,1+Int(floor(abs(pointcloud[3,i]+searchRadius-Pmin[3])/cellsizez)))],
                            [min(Ncellx,1+Int(floor(abs(pointcloud[1,i]+searchRadius-Pmin[1])/cellsizex))), min(Ncelly,1+Int(floor(abs(pointcloud[2,i]-searchRadius-Pmin[2])/cellsizey))), min(Ncellz,1+Int(floor(abs(pointcloud[3,i]-searchRadius-Pmin[3])/cellsizez)))],
                            [min(Ncellx,1+Int(floor(abs(pointcloud[1,i]-searchRadius-Pmin[1])/cellsizex))), min(Ncelly,1+Int(floor(abs(pointcloud[2,i]+searchRadius-Pmin[2])/cellsizey))), min(Ncellz,1+Int(floor(abs(pointcloud[3,i]+searchRadius-Pmin[3])/cellsizez)))],
                            [min(Ncellx,1+Int(floor(abs(pointcloud[1,i]-searchRadius-Pmin[1])/cellsizex))), min(Ncelly,1+Int(floor(abs(pointcloud[2,i]+searchRadius-Pmin[2])/cellsizey))), min(Ncellz,1+Int(floor(abs(pointcloud[3,i]-searchRadius-Pmin[3])/cellsizez)))],
                            [min(Ncellx,1+Int(floor(abs(pointcloud[1,i]-searchRadius-Pmin[1])/cellsizex))), min(Ncelly,1+Int(floor(abs(pointcloud[2,i]-searchRadius-Pmin[2])/cellsizey))), min(Ncellz,1+Int(floor(abs(pointcloud[3,i]+searchRadius-Pmin[3])/cellsizez)))],
                            [min(Ncellx,1+Int(floor(abs(pointcloud[1,i]-searchRadius-Pmin[1])/cellsizex))), min(Ncelly,1+Int(floor(abs(pointcloud[2,i]-searchRadius-Pmin[2])/cellsizey))), min(Ncellz,1+Int(floor(abs(pointcloud[3,i]-searchRadius-Pmin[3])/cellsizez)))]];=#
            cellstosearch = [cellIdx[:,i],
                            [min(Ncellx,cellIdx[1,i]+1), min(Ncelly,cellIdx[2,i]+1), min(Ncellz,cellIdx[3,i]+1)],
                            [min(Ncellx,cellIdx[1,i]+1), min(Ncelly,cellIdx[2,i]+1), cellIdx[3,i]],
                            [min(Ncellx,cellIdx[1,i]+1), min(Ncelly,cellIdx[2,i]+1), max(1,cellIdx[3,i]-1)],
                            [min(Ncellx,cellIdx[1,i]+1), cellIdx[2,i], min(Ncellz,cellIdx[3,i]+1)],
                            [min(Ncellx,cellIdx[1,i]+1), cellIdx[2,i], cellIdx[3,i]],
                            [min(Ncellx,cellIdx[1,i]+1), cellIdx[2,i], max(1,cellIdx[3,i]-1)],
                            [min(Ncellx,cellIdx[1,i]+1), max(1,cellIdx[2,i]-1), min(Ncellz,cellIdx[3,i]+1)],
                            [min(Ncellx,cellIdx[1,i]+1), max(1,cellIdx[2,i]-1), cellIdx[3,i]],
                            [min(Ncellx,cellIdx[1,i]+1), max(1,cellIdx[2,i]-1), max(1,cellIdx[3,i]-1)],
                            
                            [cellIdx[1,i], min(Ncelly,cellIdx[2,i]+1), min(Ncellz,cellIdx[3,i]+1)],
                            [cellIdx[1,i], min(Ncelly,cellIdx[2,i]+1), cellIdx[3,i]],
                            [cellIdx[1,i], min(Ncelly,cellIdx[2,i]+1), max(1,cellIdx[3,i]-1)],
                            [cellIdx[1,i], cellIdx[2,i], min(Ncellz,cellIdx[3,i]+1)],
                            [cellIdx[1,i], cellIdx[2,i], max(1,cellIdx[3,i]-1)],
                            [cellIdx[1,i], max(1,cellIdx[2,i]-1), min(Ncellz,cellIdx[3,i]+1)],
                            [cellIdx[1,i], max(1,cellIdx[2,i]-1), cellIdx[3,i]],
                            [cellIdx[1,i], max(1,cellIdx[2,i]-1), max(1,cellIdx[3,i]-1)],

                            [max(1,cellIdx[1,i]-1), min(Ncelly,cellIdx[2,i]+1), min(Ncellz,cellIdx[3,i]+1)],
                            [max(1,cellIdx[1,i]-1), min(Ncelly,cellIdx[2,i]+1), cellIdx[3,i]],
                            [max(1,cellIdx[1,i]-1), min(Ncelly,cellIdx[2,i]+1), max(1,cellIdx[3,i]-1)],
                            [max(1,cellIdx[1,i]-1), cellIdx[2,i], min(Ncellz,cellIdx[3,i]+1)],
                            [max(1,cellIdx[1,i]-1), cellIdx[2,i], cellIdx[3,i]],
                            [max(1,cellIdx[1,i]-1), cellIdx[2,i], max(1,cellIdx[3,i]-1)],
                            [max(1,cellIdx[1,i]-1), max(1,cellIdx[2,i]-1), min(Ncellz,cellIdx[3,i]+1)],
                            [max(1,cellIdx[1,i]-1), max(1,cellIdx[2,i]-1), cellIdx[3,i]],
                            [max(1,cellIdx[1,i]-1), max(1,cellIdx[2,i]-1), max(1,cellIdx[3,i]-1)]];
            unique!(cellstosearch);
            for cindex in cellstosearch
                #find the neighbors on cell[cindex]
                for j in cell[cindex[1],cindex[2],cindex[3]]
                    #if i!=j && all((abs.(pointcloud[:,i]-pointcloud[:,j]).<searchRadius))
                    if i!=j && (pointcloud[:,i]-pointcloud[:,j])'*(pointcloud[:,i]-pointcloud[:,j]).<searchRadius^2
                        push!(neighbors[i], j);
                    end
                end
            end
            unique!(neighbors[i]);
            Nneighbors[i] = length(neighbors[i]);
        end
    end
    return neighbors,Nneighbors,cell;
end


"""
DEFAULTCROSSSECTION draws a 2D rectangular cross-section with semicircular ends of width w and height h.
INPUT:
    w: Float64 extrusion width
    h: Float64 layer height
    resolution: Float64 contour mesh size
OUTPUT:
    section: Nx2 Float64 matrix with [X,Z] local coordinates as columns
"""
function defaultCrossSection(w::Float64, h::Float64, resolution::Float64)
    circularNsectors = max(1,round(Int,0.25*pi*h/resolution));      #number of segments in a quarter of circumference
    linearNsectors = max(1,round(Int,0.5*(w-h)/resolution));        #number of segments in a half width of the rectangle

    #section = zeros(1+4*circularNsectors+4*linearNsectors,2);
    #section = zeros(4*circularNsectors+4*linearNsectors,2);
    section = Matrix{Float64}(undef,2,4*circularNsectors+4*linearNsectors);
    sectionnormals = Matrix{Float64}(undef,2,4*circularNsectors+4*linearNsectors);
    section[:,1] = [-w/2, -h/2];
    sectionnormals[:,1] = [-1.0,0.0];
    lastIdx = 1;        #index of the last written row (updated every block, useful to repeat previous blocks)
    for i=1:circularNsectors
        section[:,lastIdx+i] = [-(w-h)/2+0.5*h*cos(pi-0.5*pi*i/circularNsectors), -h/2+0.5*h*sin(pi-0.5*pi*i/circularNsectors)];
        sectionnormals[:,lastIdx+i] = [cos(pi-0.5*pi*i/circularNsectors), sin(pi-0.5*pi*i/circularNsectors)];
    end
    lastIdx += circularNsectors;
    for i=1:linearNsectors
        section[:,lastIdx+i] = [-(w-h)/2+0.5*(w-h)*i/linearNsectors, 0.0];
        sectionnormals[:,lastIdx+i] = [0.0,1.0];
    end
    lastIdx += linearNsectors;
    j = 0;
    for i=lastIdx-1:-1:1
        j += 1;
        section[:,lastIdx+j] = [-section[1,i], section[2,i]];
        sectionnormals[:,lastIdx+j] = [-sectionnormals[1,i], sectionnormals[2,i]];
    end
    lastIdx += lastIdx-1;
    j = 0;
    for i=lastIdx-1:-1:2
        j += 1;
        section[:,lastIdx+j] = [section[1,i], -h-section[2,i]];
        sectionnormals[:,lastIdx+j] = [sectionnormals[1,i], -sectionnormals[2,i]];
    end
    lastIdx += lastIdx-1;
    return section,sectionnormals;
end



"""
BUILDQUADTREE builds a quadtree from a pointcloud
INPUT:
    pointcloud: 2xN matrix containing the coordinates of the points
OUTPUT:
    quadtree: vector containing the parent of each cell
    quadtreeSize: vector containing the size of each cell
    quadtreeCenter: matrix containing the coordinates of each cell center
    quadtreePoints: vector containing the indices of the points inside each cell
    quadtreeNpoints: vector containing the number of points inside each cell
"""
function buildQuadtree(pointcloud)
    Pmax = maximum(pointcloud,dims=2);
    Pmin = minimum(pointcloud,dims=2);
    boundingBoxSize = Pmax-Pmin;
    quadtree = Vector{Int}(undef,0);      #vector containing the parent of each cell
    quadtreeSize = Vector{Float64}(undef,0);      #vector containing the size of each cell
    quadtreeCenter = ElasticArray{Float64}(undef,2,0);        #matrix containing the coordinates of each cell center
    quadtreePoints = Vector{Vector{Int}}(undef,0);        #vector containing the indices of the boundary points inside each cell
    quadtreeNpoints = Vector{Int}(undef,0);       #vector containing the number of points inside each cell
    push!(quadtree,0);
    push!(quadtreeSize,maximum(boundingBoxSize));
    append!(quadtreeCenter, Pmin.+(maximum(boundingBoxSize)/2));
    push!(quadtreePoints,collect(range(1,size(pointcloud,2))));
    push!(quadtreeNpoints,size(pointcloud,2));

    maximumNpoints = size(pointcloud,2);
    while maximumNpoints>1
        #println("Maximum Npoints: ",maximumNpoints)
        n = length(quadtree);
        for i=1:n
            if quadtreeNpoints[i]>1 || quadtreeSize[i]>meshSize
                #cell division
                divideQuadtreeCell!(i,quadtree,quadtreeSize,quadtreeCenter,quadtreePoints,quadtreeNpoints);
            end
        end
        maximumNpoints = maximum(quadtreeNpoints);
    end
    return quadtree,quadtreeSize,quadtreeCenter,quadtreePoints,quadtreeNpoints;
end


function divideQuadtreeCell!(i,quadtree,quadtreeSize,quadtreeCenter,quadtreePoints,quadtreeNpoints)
    push!(quadtree, i);
    push!(quadtree, i);
    push!(quadtree, i);
    push!(quadtree, i);
    push!(quadtreeSize, quadtreeSize[i]/2);
    push!(quadtreeSize, quadtreeSize[i]/2);
    push!(quadtreeSize, quadtreeSize[i]/2);
    push!(quadtreeSize, quadtreeSize[i]/2);
    append!(quadtreeCenter, quadtreeCenter[:,i]+[quadtreeSize[i]/4,quadtreeSize[i]/4]);
    append!(quadtreeCenter, quadtreeCenter[:,i]+[quadtreeSize[i]/4,-quadtreeSize[i]/4]);
    append!(quadtreeCenter, quadtreeCenter[:,i]+[-quadtreeSize[i]/4,quadtreeSize[i]/4]);
    append!(quadtreeCenter, quadtreeCenter[:,i]+[-quadtreeSize[i]/4,-quadtreeSize[i]/4]);
    PPpoints = Vector{Int}(undef,0);
    PMpoints = Vector{Int}(undef,0);
    MPpoints = Vector{Int}(undef,0);
    MMpoints = Vector{Int}(undef,0);
    for j in quadtreePoints[i]
        relpos = pointcloud[:,j]-quadtreeCenter[:,i];
        if relpos[1]>=0 && relpos[2]>=0
            push!(PPpoints,j);
        elseif relpos[1]>=0 && relpos[2]<=0
            push!(PMpoints,j);
        elseif relpos[1]<=0 && relpos[2]>=0
            push!(MPpoints,j);
        elseif relpos[1]<=0 && relpos[2]<=0
            push!(MMpoints,j);
        end
    end
    push!(quadtreePoints, PPpoints);
    push!(quadtreePoints, PMpoints);
    push!(quadtreePoints, MPpoints);
    push!(quadtreePoints, MMpoints);
    push!(quadtreeNpoints, length(PPpoints));
    push!(quadtreeNpoints, length(PMpoints));
    push!(quadtreeNpoints, length(MPpoints));
    push!(quadtreeNpoints, length(MMpoints));
    quadtreePoints[i] = [];
    quadtreeNpoints[i] = -1;
end



"""
BUILDOCTREE builds an octree from a pointcloud
INPUT:
    pointcloud: 3xN matrix containing the coordinates of the points
OUTPUT:
    octree: vector containing the parent of each cell
    octreeSize: vector containing the size of each cell
    octreeCenter: matrix containing the coordinates of each cell center
    octreePoints: vector containing the indices of the points inside each cell
    octreeNpoints: vector containing the number of points inside each cell
"""
function buildOctree(pointcloud)
    Pmax = maximum(pointcloud,dims=2);
    Pmin = minimum(pointcloud,dims=2);
    boundingBoxSize = Pmax-Pmin;
    octree = Vector{Int}(undef,0);      #vector containing the parent of each cell
    octreeSize = Vector{Float64}(undef,0);      #vector containing the size of each cell
    octreeCenter = ElasticArray{Float64}(undef,3,0);        #matrix containing the coordinates of each cell center
    octreePoints = Vector{Vector{Int}}(undef,0);        #vector containing the indices of the boundary points inside each cell
    octreeNpoints = Vector{Int}(undef,0);       #vector containing the number of points inside each cell
    push!(octree,0);
    push!(octreeSize,maximum(boundingBoxSize));
    append!(octreeCenter, Pmin.+(maximum(boundingBoxSize)/2));
    push!(octreePoints,collect(range(1,size(pointcloud,2))));
    push!(octreeNpoints,size(pointcloud,2));

    maximumNpoints = size(pointcloud,2);
    while maximumNpoints>1
        #println("Maximum Npoints: ",maximumNpoints)
        n = length(octree);
        for i=1:n
            if octreeNpoints[i]>1 || octreeSize[i]>meshSize
                #cell division
                divideOctreeCell!(i,octree,octreeSize,octreeCenter,octreePoints,octreeNpoints);
            end
        end
        maximumNpoints = maximum(octreeNpoints);
    end
    return octree,octreeSize,octreeCenter,octreePoints,octreeNpoints;
end


function divideOctreeCell!(i,octree,octreeSize,octreeCenter,octreePoints,octreeNpoints)
    push!(octree, i);
    push!(octree, i);
    push!(octree, i);
    push!(octree, i);
    push!(octree, i);
    push!(octree, i);
    push!(octree, i);
    push!(octree, i);
    push!(octreeSize, octreeSize[i]/2);
    push!(octreeSize, octreeSize[i]/2);
    push!(octreeSize, octreeSize[i]/2);
    push!(octreeSize, octreeSize[i]/2);
    push!(octreeSize, octreeSize[i]/2);
    push!(octreeSize, octreeSize[i]/2);
    push!(octreeSize, octreeSize[i]/2);
    push!(octreeSize, octreeSize[i]/2);
    append!(octreeCenter, octreeCenter[:,i]+[octreeSize[i]/4,octreeSize[i]/4,octreeSize[i]/4]);
    append!(octreeCenter, octreeCenter[:,i]+[octreeSize[i]/4,octreeSize[i]/4,-octreeSize[i]/4]);
    append!(octreeCenter, octreeCenter[:,i]+[octreeSize[i]/4,-octreeSize[i]/4,octreeSize[i]/4]);
    append!(octreeCenter, octreeCenter[:,i]+[octreeSize[i]/4,-octreeSize[i]/4,-octreeSize[i]/4]);
    append!(octreeCenter, octreeCenter[:,i]+[-octreeSize[i]/4,octreeSize[i]/4,octreeSize[i]/4]);
    append!(octreeCenter, octreeCenter[:,i]+[-octreeSize[i]/4,octreeSize[i]/4,-octreeSize[i]/4]);
    append!(octreeCenter, octreeCenter[:,i]+[-octreeSize[i]/4,-octreeSize[i]/4,octreeSize[i]/4]);
    append!(octreeCenter, octreeCenter[:,i]+[-octreeSize[i]/4,-octreeSize[i]/4,-octreeSize[i]/4]);
    PPPpoints = Vector{Int}(undef,0);
    PPMpoints = Vector{Int}(undef,0);
    PMPpoints = Vector{Int}(undef,0);
    PMMpoints = Vector{Int}(undef,0);
    MPPpoints = Vector{Int}(undef,0);
    MPMpoints = Vector{Int}(undef,0);
    MMPpoints = Vector{Int}(undef,0);
    MMMpoints = Vector{Int}(undef,0);
    for j in octreePoints[i]
        relpos = pointcloud[:,j]-octreeCenter[:,i];
        if relpos[1]>=0 && relpos[2]>=0 && relpos[3]>=0
            push!(PPPpoints,j);
        elseif relpos[1]>=0 && relpos[2]>=0 && relpos[3]<=0
            push!(PPMpoints,j);
        elseif relpos[1]>=0 && relpos[2]<=0 && relpos[3]>=0
            push!(PMPpoints,j);
        elseif relpos[1]>=0 && relpos[2]<=0 && relpos[3]<=0
            push!(PMMpoints,j);
        elseif relpos[1]<=0 && relpos[2]>=0 && relpos[3]>=0
            push!(MPPpoints,j);
        elseif relpos[1]<=0 && relpos[2]>=0 && relpos[3]<=0
            push!(MPMpoints,j);
        elseif relpos[1]<=0 && relpos[2]<=0 && relpos[3]>=0
            push!(MMPpoints,j);
        elseif relpos[1]<=0 && relpos[2]<=0 && relpos[3]<=0
            push!(MMMpoints,j);
        end
    end
    push!(octreePoints, PPPpoints);
    push!(octreePoints, PPMpoints);
    push!(octreePoints, PMPpoints);
    push!(octreePoints, PMMpoints);
    push!(octreePoints, MPPpoints);
    push!(octreePoints, MPMpoints);
    push!(octreePoints, MMPpoints);
    push!(octreePoints, MMMpoints);
    push!(octreeNpoints, length(PPPpoints));
    push!(octreeNpoints, length(PPMpoints));
    push!(octreeNpoints, length(PMPpoints));
    push!(octreeNpoints, length(PMMpoints));
    push!(octreeNpoints, length(MPPpoints));
    push!(octreeNpoints, length(MPMpoints));
    push!(octreeNpoints, length(MMPpoints));
    push!(octreeNpoints, length(MMMpoints));
    octreePoints[i] = [];
    octreeNpoints[i] = -1;
end






"""
QUADRANTNEIGHBORSEARCH find the nearest nodes in a pointcloud using the four quadrant neighborhood approach
INPUT:
    pointcloud: 3xN Array containing the nodes coordinates [X;Y;Z]
    meshSize: Float64 target pointcloud distances
OUTPUT:
    neighbors: Vector of N Vectors containing the indices of each node neighbors
    Nneighbors: Vector of size N containing the number of neighbors of each node
"""
function quadrantNeighborSearch(pointcloud,meshSize::Float64)
    #generate cartesian subgrid
    Pmax = maximum(pointcloud,dims=2);
    Pmin = minimum(pointcloud,dims=2);
    boundingBoxSize = Pmax-Pmin;
    Ncellx = max(1,Int(round(boundingBoxSize[1]/(3*meshSize))));
    Ncelly = max(1,Int(round(boundingBoxSize[2]/(3*meshSize))));
    Ncellz = max(1,Int(round(boundingBoxSize[3]/(3*meshSize))));
    cellsizex = boundingBoxSize[1]/Ncellx;
    cellsizey = boundingBoxSize[2]/Ncelly;
    cellsizez = boundingBoxSize[3]/Ncellz;
    cell = Array{Vector{Int},3}(undef,Ncellx,Ncelly,Ncellz);
    for i=1:Ncellx
        for j=1:Ncelly
            for k=1:Ncellz
                cell[i,j,k] = Vector{Int}(undef,0);
            end
        end
    end
    N = size(pointcloud,2);     #number of points
    #cellIdx = Vector{Vector{Int},N}(undef);     #vector containing the index of the cell
    cellIdx = Matrix{Int}(undef,3,N);
    for i=1:N
        idx1 = min(Ncellx,1+Int(floor((pointcloud[1,i]-Pmin[1])/cellsizex)));
        idx2 = min(Ncelly,1+Int(floor((pointcloud[2,i]-Pmin[2])/cellsizey)));
        idx3 = min(Ncellz,1+Int(floor((pointcloud[3,i]-Pmin[3])/cellsizez)));
        push!(cell[idx1,idx2,idx3], i);
        cellIdx[:,i] = [idx1,idx2,idx3];
    end

    #find candidates using brute-force on a subset of cells
    minNeighbors = 60;
    neighbors = Vector{Vector{Int}}(undef,N);
    candidates = Vector{Vector{Int}}(undef,N);
    Nneighbors = Vector{Int}(undef,N);
    for i=1:N
        searchRadius = 1.1*meshSize;
        neighbors[i] = Int[];
        candidates[i] = Int[];
        Nneighbors[i] = 0;
        while Nneighbors[i]<minNeighbors
            searchRadius += meshSize/2;
            #determine on which cells to search
            cellstosearch = [cellIdx[:,i],
                            [min(Ncellx,cellIdx[1,i]+1), min(Ncelly,cellIdx[2,i]+1), min(Ncellz,cellIdx[3,i]+1)],
                            [min(Ncellx,cellIdx[1,i]+1), min(Ncelly,cellIdx[2,i]+1), cellIdx[3,i]],
                            [min(Ncellx,cellIdx[1,i]+1), min(Ncelly,cellIdx[2,i]+1), max(1,cellIdx[3,i]-1)],
                            [min(Ncellx,cellIdx[1,i]+1), cellIdx[2,i], min(Ncellz,cellIdx[3,i]+1)],
                            [min(Ncellx,cellIdx[1,i]+1), cellIdx[2,i], cellIdx[3,i]],
                            [min(Ncellx,cellIdx[1,i]+1), cellIdx[2,i], max(1,cellIdx[3,i]-1)],
                            [min(Ncellx,cellIdx[1,i]+1), max(1,cellIdx[2,i]-1), min(Ncellz,cellIdx[3,i]+1)],
                            [min(Ncellx,cellIdx[1,i]+1), max(1,cellIdx[2,i]-1), cellIdx[3,i]],
                            [min(Ncellx,cellIdx[1,i]+1), max(1,cellIdx[2,i]-1), max(1,cellIdx[3,i]-1)],
                            
                            [cellIdx[1,i], min(Ncelly,cellIdx[2,i]+1), min(Ncellz,cellIdx[3,i]+1)],
                            [cellIdx[1,i], min(Ncelly,cellIdx[2,i]+1), cellIdx[3,i]],
                            [cellIdx[1,i], min(Ncelly,cellIdx[2,i]+1), max(1,cellIdx[3,i]-1)],
                            [cellIdx[1,i], cellIdx[2,i], min(Ncellz,cellIdx[3,i]+1)],
                            [cellIdx[1,i], cellIdx[2,i], max(1,cellIdx[3,i]-1)],
                            [cellIdx[1,i], max(1,cellIdx[2,i]-1), min(Ncellz,cellIdx[3,i]+1)],
                            [cellIdx[1,i], max(1,cellIdx[2,i]-1), cellIdx[3,i]],
                            [cellIdx[1,i], max(1,cellIdx[2,i]-1), max(1,cellIdx[3,i]-1)],

                            [max(1,cellIdx[1,i]-1), min(Ncelly,cellIdx[2,i]+1), min(Ncellz,cellIdx[3,i]+1)],
                            [max(1,cellIdx[1,i]-1), min(Ncelly,cellIdx[2,i]+1), cellIdx[3,i]],
                            [max(1,cellIdx[1,i]-1), min(Ncelly,cellIdx[2,i]+1), max(1,cellIdx[3,i]-1)],
                            [max(1,cellIdx[1,i]-1), cellIdx[2,i], min(Ncellz,cellIdx[3,i]+1)],
                            [max(1,cellIdx[1,i]-1), cellIdx[2,i], cellIdx[3,i]],
                            [max(1,cellIdx[1,i]-1), cellIdx[2,i], max(1,cellIdx[3,i]-1)],
                            [max(1,cellIdx[1,i]-1), max(1,cellIdx[2,i]-1), min(Ncellz,cellIdx[3,i]+1)],
                            [max(1,cellIdx[1,i]-1), max(1,cellIdx[2,i]-1), cellIdx[3,i]],
                            [max(1,cellIdx[1,i]-1), max(1,cellIdx[2,i]-1), max(1,cellIdx[3,i]-1)]];
            unique!(cellstosearch);
            for cindex in cellstosearch
                #find the neighbors on cell[cindex]
                for j in cell[cindex[1],cindex[2],cindex[3]]
                    #if i!=j && all((abs.(pointcloud[:,i]-pointcloud[:,j]).<searchRadius))
                    if i!=j && (pointcloud[:,i]-pointcloud[:,j])'*(pointcloud[:,i]-pointcloud[:,j]).<searchRadius^2
                        push!(candidates[i], j);
                    end
                end
            end
            unique!(candidates[i]);
            Nneighbors[i] = length(candidates[i]);
        end

        #extract neighbors from the candidates list
        #neighbors[i] = candidates[i];
        candidatesDistances = Vector{Float64}(undef,Nneighbors[i]);
        candidatesQuadrant = Vector{Int}(undef,Nneighbors[i]);
        for j=1:Nneighbors[i]
            relpos = pointcloud[:,j]-pointcloud[:,candidates[i][j]];
            candidatesDistances[j] = relpos'*relpos;
            if relpos[1]>=0 && relpos[2]>=0 && relpos[3]>=0
                candidatesQuadrant[j] = 1;
            elseif relpos[1]>=0 && relpos[2]>=0 && relpos[3]<=0
                candidatesQuadrant[j] = 2;
            elseif relpos[1]>=0 && relpos[2]<=0 && relpos[3]>=0
                candidatesQuadrant[j] = 3;
            elseif relpos[1]>=0 && relpos[2]<=0 && relpos[3]<=0
                candidatesQuadrant[j] = 4;
            elseif relpos[1]<=0 && relpos[2]>=0 && relpos[3]>=0
                candidatesQuadrant[j] = 5;
            elseif relpos[1]<=0 && relpos[2]>=0 && relpos[3]<=0
                candidatesQuadrant[j] = 6;
            elseif relpos[1]<=0 && relpos[2]<=0 && relpos[3]>=0
                candidatesQuadrant[j] = 7;
            elseif relpos[1]<=0 && relpos[2]<=0 && relpos[3]<=0
                candidatesQuadrant[j] = 8;
            end
        end
        sortidx = sortperm(candidatesDistances);
        candidatesQuadrant = candidatesQuadrant[sortidx];
        candidates[i] = candidates[i][sortidx];
        quadrantPoints1 = findall(candidatesQuadrant.==1);
        if length(quadrantPoints1)>=2
            push!(neighbors[i], candidates[i][quadrantPoints1[1]]);
            push!(neighbors[i], candidates[i][quadrantPoints1[2]]);
        end
        quadrantPoints2 = findall(candidatesQuadrant.==2);
        if length(quadrantPoints2)>=2
            push!(neighbors[i], candidates[i][quadrantPoints2[1]]);
            push!(neighbors[i], candidates[i][quadrantPoints2[2]]);
        end
        quadrantPoints3 = findall(candidatesQuadrant.==3);
        if length(quadrantPoints3)>=2
            push!(neighbors[i], candidates[i][quadrantPoints3[1]]);
            push!(neighbors[i], candidates[i][quadrantPoints3[2]]);
        end
        quadrantPoints4 = findall(candidatesQuadrant.==4);
        if length(quadrantPoints4)>=2
            push!(neighbors[i], candidates[i][quadrantPoints4[1]]);
            push!(neighbors[i], candidates[i][quadrantPoints4[2]]);
        end
        quadrantPoints5 = findall(candidatesQuadrant.==5);
        if length(quadrantPoints5)>=2
            push!(neighbors[i], candidates[i][quadrantPoints5[1]]);
            push!(neighbors[i], candidates[i][quadrantPoints5[2]]);
        end
        quadrantPoints6 = findall(candidatesQuadrant.==6);
        if length(quadrantPoints6)>=2
            push!(neighbors[i], candidates[i][quadrantPoints6[1]]);
            push!(neighbors[i], candidates[i][quadrantPoints6[2]]);
        end
        quadrantPoints7 = findall(candidatesQuadrant.==7);
        if length(quadrantPoints7)>=2
            push!(neighbors[i], candidates[i][quadrantPoints7[1]]);
            push!(neighbors[i], candidates[i][quadrantPoints7[2]]);
        end
        quadrantPoints8 = findall(candidatesQuadrant.==8);
        if length(quadrantPoints8)>=2
            push!(neighbors[i], candidates[i][quadrantPoints8[1]]);
            push!(neighbors[i], candidates[i][quadrantPoints8[2]]);
        end
        if length(neighbors[i])<16
            if i>length(boundaryNodes)
                println("Warning point #",i)
            end
            #neighbors[i] = candidates[i];
            while length(neighbors[i])<16
                push!(neighbors[i], candidates[i][rand(1:length(candidates[i]))]);
                unique!(neighbors[i]);
            end
        end
        Nneighbors[i] = length(neighbors[i]);

    end
    return neighbors,Nneighbors;
end
