/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2017     Statoil ASA
// 
//  ResInsight is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
// 
//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
//  FITNESS FOR A PARTICULAR PURPOSE.
// 
//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html> 
//  for more details.
//
/////////////////////////////////////////////////////////////////////////////////


#include "RigCellGeometryTools.h"
#include "cvfStructGrid.h"
#include "cvfGeometryTools.h"

#include "cafHexGridIntersectionTools/cafHexGridIntersectionTools.h"
#include "cvfBoundingBox.h"

#include "clipper/clipper.hpp"

#include <vector>
#include <array>


//--------------------------------------------------------------------------------------------------
/// Internal functions
//--------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------
/// Splits a line in a number of equal parts
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> splitLine(cvf::Vec3d ptStart, cvf::Vec3d ptEnd, size_t partCount);

//--------------------------------------------------------------------------------------------------
/// Calculates all points on a face described by edge points from all four edges.
/// The result is a grid of points including the given edge points
///
///                    edgeXPtsHigh
///                  |-------------|
///                  |             |
///   edgeYPtsLow    |             |   edgeYPtsHigh
///                  |             |
///                  |-------------|
///                    edgeXPtsLow
///
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<cvf::Vec3d>> calcFacePoints(const std::vector<cvf::Vec3d> edgeXPtsLow,
    const std::vector<cvf::Vec3d> edgeXPtsHigh,
    const std::vector<cvf::Vec3d> edgeYPtsLow,
    const std::vector<cvf::Vec3d> edgeYPtsHigh);

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigCellGeometryTools::createPolygonFromLineSegments(std::list<std::pair<cvf::Vec3d, cvf::Vec3d>> &intersectionLineSegments, std::vector<std::vector<cvf::Vec3d>> &polygons)
{
    bool startNewPolygon = true;
    while (!intersectionLineSegments.empty())
    {
        if (startNewPolygon)
        {
            std::vector<cvf::Vec3d> polygon;
            //Add first line segments to polygon and remove from list
            std::pair<cvf::Vec3d, cvf::Vec3d > linesegment = intersectionLineSegments.front();
            polygon.push_back(linesegment.first);
            polygon.push_back(linesegment.second);
            intersectionLineSegments.remove(linesegment);
            polygons.push_back(polygon);
            startNewPolygon = false;
        }

        std::vector<cvf::Vec3d>& polygon = polygons.back();

        //Search remaining list for next point...

        bool isFound = false;
        float tolerance = 0.0001f;

        for (std::list<std::pair<cvf::Vec3d, cvf::Vec3d > >::iterator lIt = intersectionLineSegments.begin(); lIt != intersectionLineSegments.end(); lIt++)
        {
            cvf::Vec3d lineSegmentStart = lIt->first;
            cvf::Vec3d lineSegmentEnd = lIt->second;
            cvf::Vec3d polygonEnd = polygon.back();
            
            double lineSegmentLength = (lineSegmentStart - lineSegmentEnd).lengthSquared();
            if (lineSegmentLength < tolerance*tolerance)
            {
                intersectionLineSegments.erase(lIt);
                isFound = true; 
                break;
            }


            double lineSegmentStartDiff = (lineSegmentStart - polygonEnd).lengthSquared();
            if (lineSegmentStartDiff < tolerance*tolerance)
            {
                polygon.push_back(lIt->second);
                intersectionLineSegments.erase(lIt);
                isFound = true;
                break;
            }

            double lineSegmentEndDiff = (lineSegmentEnd - polygonEnd).lengthSquared();
            if (lineSegmentEndDiff < tolerance*tolerance)
            {
                polygon.push_back(lIt->first);
                intersectionLineSegments.erase(lIt);
                isFound = true;
                break;
            }
        }

        if (isFound)
        {
            continue;
        }
        else
        {
            startNewPolygon = true;
        }
    }



}

//==================================================================================================
/// 
//==================================================================================================
void RigCellGeometryTools::findCellLocalXYZ(const std::array<cvf::Vec3d, 8>& hexCorners, cvf::Vec3d& localXdirection, cvf::Vec3d& localYdirection, cvf::Vec3d& localZdirection)
{
    cvf::ubyte faceVertexIndices[4];
    cvf::StructGridInterface::FaceEnum face;

    face = cvf::StructGridInterface::NEG_I;
    cvf::StructGridInterface::cellFaceVertexIndices(face, faceVertexIndices);
    cvf::Vec3d faceCenterNegI = cvf::GeometryTools::computeFaceCenter(hexCorners[faceVertexIndices[0]], hexCorners[faceVertexIndices[1]], hexCorners[faceVertexIndices[2]], hexCorners[faceVertexIndices[3]]);
    //TODO: Should we use face centroids instead of face centers?

    face = cvf::StructGridInterface::POS_I;
    cvf::StructGridInterface::cellFaceVertexIndices(face, faceVertexIndices);
    cvf::Vec3d faceCenterPosI = cvf::GeometryTools::computeFaceCenter(hexCorners[faceVertexIndices[0]], hexCorners[faceVertexIndices[1]], hexCorners[faceVertexIndices[2]], hexCorners[faceVertexIndices[3]]);

    face = cvf::StructGridInterface::NEG_J;
    cvf::StructGridInterface::cellFaceVertexIndices(face, faceVertexIndices);
    cvf::Vec3d faceCenterNegJ = cvf::GeometryTools::computeFaceCenter(hexCorners[faceVertexIndices[0]], hexCorners[faceVertexIndices[1]], hexCorners[faceVertexIndices[2]], hexCorners[faceVertexIndices[3]]);

    face = cvf::StructGridInterface::POS_J;
    cvf::StructGridInterface::cellFaceVertexIndices(face, faceVertexIndices);
    cvf::Vec3d faceCenterPosJ = cvf::GeometryTools::computeFaceCenter(hexCorners[faceVertexIndices[0]], hexCorners[faceVertexIndices[1]], hexCorners[faceVertexIndices[2]], hexCorners[faceVertexIndices[3]]);

    cvf::Vec3d faceCenterCenterVectorI = faceCenterPosI - faceCenterNegI;
    cvf::Vec3d faceCenterCenterVectorJ = faceCenterPosJ - faceCenterNegJ;

    localZdirection.cross(faceCenterCenterVectorI, faceCenterCenterVectorJ);
    localZdirection.normalize();

    cvf::Vec3d crossPoductJZ;
    crossPoductJZ.cross(faceCenterCenterVectorJ, localZdirection);
    localXdirection = faceCenterCenterVectorI + crossPoductJZ;
    localXdirection.normalize();

    cvf::Vec3d crossPoductIZ;
    crossPoductIZ.cross(faceCenterCenterVectorI, localZdirection);
    localYdirection = faceCenterCenterVectorJ - crossPoductIZ;
    localYdirection.normalize();

    //TODO: Check if we end up with 0-vectors, and handle this case...
}

//--------------------------------------------------------------------------------------------------
///  
//--------------------------------------------------------------------------------------------------
double RigCellGeometryTools::polygonLengthInLocalXdirWeightedByArea(const std::vector<cvf::Vec3d>& polygonToCalcLengthOf)
{
    //Find bounding box
    cvf::BoundingBox polygonBBox;
    for (cvf::Vec3d nodeCoord : polygonToCalcLengthOf) polygonBBox.add(nodeCoord);
    cvf::Vec3d bboxCorners[8];
    polygonBBox.cornerVertices(bboxCorners);


    //Split bounding box in multiple polygons (2D)
    int resolutionOfLengthCalc = 20;
    double widthOfPolygon = polygonBBox.extent().y() / resolutionOfLengthCalc;

    std::vector<double> areasOfPolygonContributions;
    std::vector<double> lengthOfPolygonContributions;

    cvf::Vec3d directionOfLength(1, 0, 0);

    for (int i = 0; i < resolutionOfLengthCalc; i++)
    {
        cvf::Vec3d pointOnLine1(bboxCorners[0].x(), bboxCorners[0].y() + i*widthOfPolygon, 0);
        cvf::Vec3d pointOnLine2(bboxCorners[0].x(), bboxCorners[0].y() + (i + 1)*widthOfPolygon, 0);

        std::pair<cvf::Vec3d, cvf::Vec3d> line1 = getLineThroughBoundingBox(directionOfLength,
                                                                            polygonBBox,
                                                                            pointOnLine1);
        std::pair<cvf::Vec3d, cvf::Vec3d> line2 = getLineThroughBoundingBox(directionOfLength,
                                                                            polygonBBox,
                                                                            pointOnLine2);

        std::vector<cvf::Vec3d> polygon;
        polygon.push_back(line1.first);
        polygon.push_back(line1.second);
        polygon.push_back(line2.second);
        polygon.push_back(line2.first);

        //Use clipper to find overlap between bbpolygon and fracture
        std::vector<std::vector<cvf::Vec3d> > clippedPolygons = intersectPolygons(polygonToCalcLengthOf, polygon);

        double area = 0;
        double length = 0;
        cvf::Vec3d areaVector = cvf::Vec3d::ZERO;

        //Calculate length (max-min) and area    
        for (std::vector<cvf::Vec3d> clippedPolygon : clippedPolygons)
        {
            areaVector = cvf::GeometryTools::polygonAreaNormal3D(clippedPolygon);
            area += areaVector.length();
            length += (getLengthOfPolygonAlongLine(line1, clippedPolygon)
                     + getLengthOfPolygonAlongLine(line2, clippedPolygon)) / 2;
        }
        areasOfPolygonContributions.push_back(area);
        lengthOfPolygonContributions.push_back(length);
    }

    //Calculate area-weighted length average.  
    double totalArea = 0.0;
    double totalAreaXlength = 0.0;

    for (size_t i = 0; i < areasOfPolygonContributions.size(); i++)
    {
        totalArea += areasOfPolygonContributions[i];
        totalAreaXlength += (areasOfPolygonContributions[i] * lengthOfPolygonContributions[i]);
    }

    double areaWeightedLength = totalAreaXlength / totalArea;
    return areaWeightedLength;
}


double clipperConversionFactor = 10000; //For transform to clipper int 

ClipperLib::IntPoint toClipperPoint(const cvf::Vec3d& cvfPoint)
{
    int xInt = cvfPoint.x()*clipperConversionFactor;
    int yInt = cvfPoint.y()*clipperConversionFactor;
    int zInt = cvfPoint.z()*clipperConversionFactor;
    return  ClipperLib::IntPoint(xInt, yInt, zInt); 
}

cvf::Vec3d fromClipperPoint(const ClipperLib::IntPoint& clipPoint)
{
    double zDValue; 
    
    if (clipPoint.Z == std::numeric_limits<int>::max()) 
    {
        zDValue = HUGE_VAL;
    }
    else 
    {
        zDValue = clipPoint.Z;
    }

    return cvf::Vec3d (clipPoint.X, clipPoint.Y, zDValue ) /clipperConversionFactor;
}

//--------------------------------------------------------------------------------------------------
///  
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::intersectPolygons(const std::vector<cvf::Vec3d>& polygon1,
                                                                              const std::vector<cvf::Vec3d>& polygon2)
{
    std::vector<std::vector<cvf::Vec3d> > clippedPolygons;
    
    // Convert to int for clipper library and store as clipper "path"
    ClipperLib::Path polygon1path;
    for (const cvf::Vec3d& v : polygon1)
    {
        polygon1path.push_back(toClipperPoint(v));
    }

    ClipperLib::Path polygon2path;
    for (const cvf::Vec3d& v : polygon2)
    {
        polygon2path.push_back(toClipperPoint(v));
    }

    ClipperLib::Clipper clpr;
    clpr.AddPath(polygon1path, ClipperLib::ptSubject, true);
    clpr.AddPath(polygon2path, ClipperLib::ptClip, true);

    ClipperLib::Paths solution;
    clpr.Execute(ClipperLib::ctIntersection, solution, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);

    // Convert back to std::vector<std::vector<cvf::Vec3d> >
    for (ClipperLib::Path pathInSol : solution)
    {
        std::vector<cvf::Vec3d> clippedPolygon;
        for (ClipperLib::IntPoint IntPosition : pathInSol)
        {
            clippedPolygon.push_back(fromClipperPoint(IntPosition));
        }
        clippedPolygons.push_back(clippedPolygon);
    }

    return clippedPolygons;
}

//-------------------------------------------------------------------------------------------------- 
///  
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<cvf::Vec3d>>
    RigCellGeometryTools::subtractPolygons(const std::vector<cvf::Vec3d>&              sourcePolygon,
                                           const std::vector<std::vector<cvf::Vec3d>>& polygonsToSubtract)
{
    ClipperLib::Clipper clpr;

    {
        // Convert to int for clipper library and store as clipper "path"
        ClipperLib::Path polygon1path;
        for (const auto& v : sourcePolygon)
        {
            polygon1path.push_back(toClipperPoint(v));
        }
        clpr.AddPath(polygon1path, ClipperLib::ptSubject, true);
    }

    for (const auto& path : polygonsToSubtract)
    {
        ClipperLib::Path polygon2path;
        for (const auto& v : path)
        {
            polygon2path.push_back(toClipperPoint(v));
        }

        clpr.AddPath(polygon2path, ClipperLib::ptClip, true);
    }

    ClipperLib::Paths solution;
    clpr.Execute(ClipperLib::ctDifference, solution, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);

    std::vector<std::vector<cvf::Vec3d>> clippedPolygons;

    // Convert back to std::vector<std::vector<cvf::Vec3d> >
    for (ClipperLib::Path pathInSol : solution)
    {
        std::vector<cvf::Vec3d> clippedPolygon;
        for (ClipperLib::IntPoint IntPosition : pathInSol)
        {
            clippedPolygon.push_back(fromClipperPoint(IntPosition));
        }

        clippedPolygons.push_back(clippedPolygon);
    }

    return clippedPolygons;
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void fillInterpolatedSubjectZ(ClipperLib::IntPoint& e1bot,
                              ClipperLib::IntPoint& e1top,
                              ClipperLib::IntPoint& e2bot,
                              ClipperLib::IntPoint& e2top,
                              ClipperLib::IntPoint& pt)
{
    ClipperLib::IntPoint ePLbot;
    ClipperLib::IntPoint ePLtop;

    if (e1top.Z == std::numeric_limits<int>::max()) 
    {
        ePLtop = e2top;
        ePLbot = e2bot;
    }
    else
    {
        ePLtop = e1top;
        ePLbot = e1bot;
    }

    double ePLXRange = (ePLtop.X - ePLbot.X);
    double ePLYRange = (ePLtop.Y - ePLbot.Y);

    double ePLLength =  sqrt(ePLXRange*ePLXRange + ePLYRange*ePLYRange);
    
    if (ePLLength <= 1)
    {
        pt.Z = ePLbot.Z;
        return;
    }

    double ePLBotPtXRange = pt.X - ePLbot.X;
    double ePLBotPtYRange = pt.Y - ePLbot.Y;

    double ePLBotPtLength =  sqrt(ePLBotPtXRange*ePLBotPtXRange + ePLBotPtYRange*ePLBotPtYRange);

    double fraction = ePLBotPtLength/ePLLength;

    pt.Z = std::nearbyint( ePLbot.Z + fraction*(ePLtop.Z - ePLbot.Z) );
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void fillUndefinedZ(ClipperLib::IntPoint& e1bot,
                              ClipperLib::IntPoint& e1top,
                              ClipperLib::IntPoint& e2bot,
                              ClipperLib::IntPoint& e2top,
                              ClipperLib::IntPoint& pt)
{
   pt.Z = std::numeric_limits<int>::max();
}

//--------------------------------------------------------------------------------------------------
/// Assumes x.y plane polygon. Polyline might have a Z, and the returned Z is the polyline Z, interpolated if it is clipped.
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolylineByPolygon(const std::vector<cvf::Vec3d>& polyLine, 
                                                                                  const std::vector<cvf::Vec3d>& polygon, 
                                                                                  ZInterpolationType interpolType)
{
    std::vector<std::vector<cvf::Vec3d> > clippedPolyline;

    //Adjusting polygon to avoid clipper issue with interpolating z-values when lines crosses though polygon vertecies
    std::vector<cvf::Vec3d> adjustedPolygon = ajustPolygonToAvoidIntersectionsAtVertex(polyLine, polygon);

    //Convert to int for clipper library and store as clipper "path"
    ClipperLib::Path polyLinePath;
    for (const cvf::Vec3d& v : polyLine)
    {
        polyLinePath.push_back(toClipperPoint(v));
    }

    ClipperLib::Path polygonPath;
    for (const cvf::Vec3d& v : adjustedPolygon)
    {
        ClipperLib::IntPoint intp = toClipperPoint(v);
        intp.Z = std::numeric_limits<int>::max();
        polygonPath.push_back(intp);
    }

    ClipperLib::Clipper clpr;
    clpr.AddPath(polyLinePath, ClipperLib::ptSubject, false);
    clpr.AddPath(polygonPath, ClipperLib::ptClip, true);

    if ( interpolType == INTERPOLATE_LINE_Z )
    {
        clpr.ZFillFunction(&fillInterpolatedSubjectZ);
    }
    else if ( interpolType == USE_HUGEVAL )
    {
        clpr.ZFillFunction(&fillUndefinedZ);
    }

    ClipperLib::PolyTree solution;
    clpr.Execute(ClipperLib::ctIntersection, solution, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);

    // We only expect open paths from this method (unless the polyline is self intersecting, a condition we do not handle)
    ClipperLib::Paths solutionPaths;
    ClipperLib::OpenPathsFromPolyTree(solution, solutionPaths);

    //Convert back to std::vector<std::vector<cvf::Vec3d> >
    for (ClipperLib::Path pathInSol : solutionPaths)
    {
        std::vector<cvf::Vec3d> clippedPolygon;
        for (ClipperLib::IntPoint IntPosition : pathInSol)
        {
            clippedPolygon.push_back(fromClipperPoint(IntPosition));
        }
        clippedPolyline.push_back(clippedPolygon);
    }

    return clippedPolyline;
}

//--------------------------------------------------------------------------------------------------
///  
//--------------------------------------------------------------------------------------------------
std::pair<cvf::Vec3d, cvf::Vec3d> RigCellGeometryTools::getLineThroughBoundingBox(const cvf::Vec3d& lineDirection,
                                                                                  const cvf::BoundingBox& polygonBBox,
                                                                                  const cvf::Vec3d& pointOnLine)
{
    cvf::Vec3d bboxCorners[8];
    polygonBBox.cornerVertices(bboxCorners);

    cvf::Vec3d startPoint = pointOnLine;
    cvf::Vec3d endPoint = pointOnLine;
    cvf::Vec3d lineDir = lineDirection;

    //To avoid doing many iterations in loops below linedirection should be quite large. 
    lineDir.normalize();
    lineDir = lineDir * polygonBBox.extent().length() / 5;

    //Extend line in positive direction
    while (polygonBBox.contains(startPoint))
    {
        startPoint = startPoint + lineDir;
    }
    //Extend line in negative direction
    while (polygonBBox.contains(endPoint))
    {
        endPoint = endPoint - lineDir;
    }

    std::pair<cvf::Vec3d, cvf::Vec3d> line;
    line = { startPoint, endPoint };
    return line;
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
double RigCellGeometryTools::getLengthOfPolygonAlongLine(const std::pair<cvf::Vec3d, cvf::Vec3d>& line, const std::vector<cvf::Vec3d>& polygon)
{
    cvf::BoundingBox lineBoundingBox;

    for (cvf::Vec3d polygonPoint : polygon)
    {
        cvf::Vec3d pointOnLine = cvf::GeometryTools::projectPointOnLine(line.first, line.second, polygonPoint, nullptr);
        lineBoundingBox.add(pointOnLine);
    }

    double length = lineBoundingBox.extent().length();
    
    return length;
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> RigCellGeometryTools::unionOfPolygons(const std::vector<std::vector<cvf::Vec3d>>& polygons)
{
    // Convert to int for clipper library and store as clipper "path"
    std::vector<ClipperLib::Path> polygonPaths;
    for (const std::vector<cvf::Vec3d>& polygon : polygons)
    {
        polygonPaths.emplace_back();
        auto& p = polygonPaths.back();
        for (const cvf::Vec3d& pp : polygon)
        {
            p.push_back(toClipperPoint(pp));
        }
    }

    ClipperLib::Clipper clpr;
    clpr.AddPaths(polygonPaths, ClipperLib::ptSubject, true);

    ClipperLib::Paths solution;
    clpr.Execute(ClipperLib::ctUnion, solution, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);

    // Convert back to std::vector<std::vector<cvf::Vec3d> >
    std::vector<cvf::Vec3d> unionPolygon;
    for (ClipperLib::Path pathInSol : solution)
    {
        std::vector<cvf::Vec3d> clippedPolygon;
        for (ClipperLib::IntPoint IntPosition : pathInSol)
        {
            unionPolygon.push_back(fromClipperPoint(IntPosition));
        }
    }

    return unionPolygon;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> splitLine(cvf::Vec3d ptStart, cvf::Vec3d ptEnd, size_t partCount)
{
	std::vector<cvf::Vec3d> pts = { ptStart };

	for (int i = 1; i < partCount; i++)
	{
		pts.push_back(cvf::Vec3d(ptStart.x() + (ptEnd.x() - ptStart.x()) * i / partCount,
								 ptStart.y() + (ptEnd.y() - ptStart.y()) * i / partCount,
								 ptStart.z() + (ptEnd.z() - ptStart.z()) * i / partCount));
	}
	pts.push_back(ptEnd);
	return pts;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<cvf::Vec3d>>
    calcFacePoints(const std::vector<cvf::Vec3d> edgeXPtsLow, const std::vector<cvf::Vec3d> edgeXPtsHigh,
                   const std::vector<cvf::Vec3d> edgeYPtsLow, const std::vector<cvf::Vec3d> edgeYPtsHigh)
{
    CVF_ASSERT(edgeXPtsLow.size() == edgeXPtsHigh.size() && edgeYPtsLow.size() == edgeYPtsHigh.size());

    size_t xSize = edgeXPtsLow.size();
    size_t ySize = edgeYPtsLow.size();

    std::vector<std::vector<cvf::Vec3d>> pts;

    // Add low edge points
    pts.push_back(edgeXPtsLow);

    // Interior points
    for (int y = 1; y < ySize - 1; y++)
    {
        //for (int x = 0; x < xSize; x++)
        {
            auto interiorPts = splitLine(edgeYPtsLow[y], edgeYPtsHigh[y], xSize - 1);
            pts.push_back(interiorPts);
        }
    }

    // Add low edge points
    pts.push_back(edgeXPtsHigh);
    return pts;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d>
RigCellGeometryTools::createHexCornerCoords(std::array<cvf::Vec3d, 8> mainCellCorners, size_t nx, size_t ny, size_t nz)
{
	std::array<std::pair<size_t, size_t>, 12> edgeCorners =
	{
		std::make_pair(0, 1), std::make_pair(3, 2), std::make_pair(4, 5), std::make_pair(7, 6), // X
		std::make_pair(0, 3), std::make_pair(4, 7), std::make_pair(1, 2), std::make_pair(5, 6), // Y
		std::make_pair(0, 4), std::make_pair(1, 5), std::make_pair(3, 7), std::make_pair(2, 6), // Z
	};

	std::array<size_t, 3> nxyz = { nx, ny, nz };
	std::array<std::vector<cvf::Vec3d>, 12> edgePoints;

	for (int i = 0; i < 12; i++)
	{
		int partCountsIndex = i / 4;
		edgePoints[i] = splitLine(mainCellCorners[edgeCorners[i].first], mainCellCorners[edgeCorners[i].second], nxyz[partCountsIndex]);
	}
	

    // lowIJ, highIJ, lowJK, highKJ, 
    std::vector<std::vector<std::vector<cvf::Vec3d>>> nodes;
    nodes.reserve((nx + 1)*(ny + 1)*(nz + 1));

    auto xyFacePtsLow = calcFacePoints(edgePoints[0], edgePoints[1], edgePoints[4], edgePoints[6]);
    auto xyFacePtsHigh = calcFacePoints(edgePoints[2], edgePoints[3], edgePoints[5], edgePoints[7]);
    auto yzFacePtsLow  = calcFacePoints(edgePoints[4], edgePoints[5], edgePoints[8], edgePoints[10]);
    auto yzFacePtsHigh = calcFacePoints(edgePoints[6], edgePoints[7], edgePoints[9], edgePoints[11]);
    auto xzFacePtsLow  = calcFacePoints(edgePoints[0], edgePoints[2], edgePoints[8], edgePoints[9]);
    auto xzFacePtsHigh = calcFacePoints(edgePoints[1], edgePoints[3], edgePoints[10], edgePoints[11]);

    nodes.push_back(xyFacePtsLow);

    for (int z = 1; z < nz; z++)
    {
        auto xyFacePoints = calcFacePoints(xzFacePtsLow[z], xzFacePtsHigh[z], yzFacePtsLow[z], yzFacePtsHigh[z]);
        nodes.push_back(xyFacePoints);
    }

    nodes.push_back(xyFacePtsHigh);
	
    std::vector<cvf::Vec3d> coords;
    coords.reserve(nx*ny*nz * 8);

    for (int z = 1; z < nz + 1; z++)
    {
        for (int y = 1; y < ny + 1; y++)
        {
            for (int x = 1; x < nx + 1; x++)
            {
                std::array<cvf::Vec3d, 8> cs;

                cs[0] = nodes[z - 1][y - 1][x - 1];
                cs[1] = nodes[z - 1][y - 1][x];
                cs[2] = nodes[z - 1][y][x];
                cs[3] = nodes[z - 1][y][x - 1];

                cs[4] = nodes[z][y - 1][x - 1];
                cs[5] = nodes[z][y - 1][x];
                cs[6] = nodes[z][y][x];
                cs[7] = nodes[z][y][x - 1];

                coords.insert(coords.end(), cs.begin(), cs.end());
            }
        }
    }
    return coords;
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> RigCellGeometryTools::ajustPolygonToAvoidIntersectionsAtVertex(const std::vector<cvf::Vec3d>& polyLine, 
                                                                                       const std::vector<cvf::Vec3d>& polygon)
{
    std::vector<cvf::Vec3d> adjustedPolygon; 

    double treshold =  (1.0 / 10000.0) * 5; //5 times polygonScaleFactor for converting to int for clipper

    for (cvf::Vec3d polygonPoint : polygon)
    {
        
        for (size_t i = 0; i < polyLine.size() - 1; i++)
        {
            cvf::Vec3d linePoint1(polyLine[i].x(), polyLine[i].y(), 0.0);
            cvf::Vec3d linePoint2(polyLine[i + 1].x(), polyLine[i + 1].y(), 0.0);

            double pointDistanceFromLine = cvf::GeometryTools::linePointSquareDist(linePoint1, linePoint2, polygonPoint);
            if (pointDistanceFromLine < treshold)
            {
                //calculate new polygonPoint
                cvf::Vec3d directionOfLineSegment = linePoint2 - linePoint1;
                //finding normal to the direction of the line segment in the XY plane (z=0)
                cvf::Vec3d normalToLine(-directionOfLineSegment.y(), directionOfLineSegment.x(), 0.0);
                normalToLine.normalize();
                polygonPoint = polygonPoint + normalToLine * 0.005;
            }
        }
        adjustedPolygon.push_back(polygonPoint);
    }
    
    return adjustedPolygon;
}