#include "SourceStamp.h"

#include <EGL/egl.h>
#include <EGL/eglext.h>
#include <GLES3/gl3.h>
#include <png.h>

#include <algorithm>

#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <limits>
#include <map>
#include <set>
#include <string>
#include <string_view>
#include <utility>
#include <unordered_set>
#include <vector>

namespace
{
    struct FVertex
    {
        float Position[3];
        float Normal[3];
        float Height;
        float WorldPosition[3];
        float TileData[4];
        float Material[4]{};
    };

    struct FScene
    {
        std::uint32_t VoxelSizeMicrons=0;
        std::array<double,3> CameraMeters{};
        double YawDegrees=0.0,PitchDegrees=0.0;
        float MinimumHeight=std::numeric_limits<float>::infinity();
        float MaximumHeight=-std::numeric_limits<float>::infinity();
        std::vector<FVertex> Vertices;
        std::vector<std::uint32_t> Indices;
        std::vector<std::uint32_t> SurfaceIndices;
        struct FTileRange
        {
            std::array<std::int32_t,4> Key{};
            std::size_t IndexBegin=0;
            std::uint32_t IndexCount=0;
            std::vector<std::uint8_t> FaceMasks;
        };
        std::vector<FTileRange> Tiles;
        std::uint64_t Triangles=0,SkirtTriangles=0;
        std::uint64_t VisibleSkirtTriangles=0,HiddenSameLodSkirtTriangles=0;
        std::uint64_t DrawnTriangles=0;
        std::uint32_t TileCount=0;
    };

    struct FTileKeyHash
    {
        std::size_t operator()(const std::array<std::int32_t,4>& Key) const noexcept
        {
            std::size_t Hash=0xcbf29ce484222325ull;
            for(const std::int32_t Value:Key)
            {
                Hash^=static_cast<std::uint32_t>(Value);
                Hash*=0x100000001b3ull;
            }
            return Hash;
        }
    };

    std::int32_t FloorDiv(std::int64_t Value,std::int64_t Divisor)
    {
        const std::int64_t Quotient=Value/Divisor, Remainder=Value%Divisor;
        return static_cast<std::int32_t>(Quotient-((Remainder<0)?1:0));
    }

    std::uint8_t ComputeVisibleSkirtFaceMask(
        const std::array<std::int32_t,4>& Tile,
        const std::unordered_set<std::array<std::int32_t,4>,FTileKeyHash>& Desired)
    {
        const std::int64_t Extent=static_cast<std::int64_t>(32)<<Tile[3];
        const std::int64_t Origin[3]={static_cast<std::int64_t>(Tile[0])*Extent,
            static_cast<std::int64_t>(Tile[1])*Extent,static_cast<std::int64_t>(Tile[2])*Extent};
        std::uint8_t Mask=0;
        for(std::int32_t Axis=0;Axis<3;++Axis)for(std::int32_t Side=0;Side<2;++Side)
        {
            std::int64_t Probe[3]={Origin[0]+Extent/2,Origin[1]+Extent/2,
                Origin[2]+Extent/2};
            Probe[Axis]=Side?Origin[Axis]+Extent:Origin[Axis]-1;
            std::int32_t NeighborLevel=-1;
            for(std::int32_t Level=0;Level<=8;++Level)
            {
                const std::int64_t NeighborExtent=static_cast<std::int64_t>(32)<<Level;
                const std::array<std::int32_t,4> Neighbor={
                    FloorDiv(Probe[0],NeighborExtent),FloorDiv(Probe[1],NeighborExtent),
                    FloorDiv(Probe[2],NeighborExtent),Level};
                if(Desired.contains(Neighbor))
                {NeighborLevel=Level;break;}
            }
            if(NeighborLevel<0||NeighborLevel!=Tile[3])
                Mask|=static_cast<std::uint8_t>(1u<<(Axis*2+Side));
        }
        return Mask;
    }

    template<typename T> bool Read(std::ifstream& File,T& Value)
    {
        File.read(reinterpret_cast<char*>(&Value),sizeof(Value));
        return static_cast<bool>(File);
    }

    bool LoadScene(const std::filesystem::path& Path,FScene& Scene,std::string& Error)
    {
        std::ifstream File(Path,std::ios::binary);
        if(!File){Error="cannot open mesh scene";return false;}
        char Magic[8]{};
        File.read(Magic,sizeof(Magic));
        const bool HasMaterials=File&&std::memcmp(Magic,"VFMSH002",8)==0;
        if(!File||(!HasMaterials&&std::memcmp(Magic,"VFMSH001",8)!=0))
        {Error="mesh scene magic is invalid";return false;}
        std::uint32_t Version=0,ExpectedTiles=0;
        if(!Read(File,Version)||!Read(File,Scene.VoxelSizeMicrons)||!Read(File,ExpectedTiles)
            ||Version!=(HasMaterials?2u:1u)||Scene.VoxelSizeMicrons==0||ExpectedTiles>1000000)
        {Error="mesh scene header is malformed";return false;}
        for(double& Value:Scene.CameraMeters)if(!Read(File,Value))
        {Error="mesh scene camera is truncated";return false;}
        if(!Read(File,Scene.YawDegrees)||!Read(File,Scene.PitchDegrees))
        {Error="mesh scene camera orientation is truncated";return false;}
        Scene.TileCount=ExpectedTiles;
        std::uint64_t TotalVertices=0,TotalIndices=0;
        for(std::uint32_t Tile=0;Tile<ExpectedTiles;++Tile)
        {
            std::int32_t Key[4]{};
            std::uint32_t SurfaceVertices=0,SkirtCount=0,VertexCount=0,IndexCount=0,MaskCount=0;
            for(auto& Value:Key)if(!Read(File,Value)){Error="tile key is truncated";return false;}
            if(!Read(File,SurfaceVertices)||!Read(File,SkirtCount)||!Read(File,VertexCount)
                ||!Read(File,IndexCount)||!Read(File,MaskCount)
                ||VertexCount>10000000||IndexCount>60000000||MaskCount>20000000
                ||IndexCount%3!=0||(MaskCount!=0&&MaskCount!=IndexCount/3)
                ||Key[3]<0||Key[3]>8)
            {Error="tile mesh header is malformed";return false;}
            if(TotalVertices+VertexCount>std::numeric_limits<std::uint32_t>::max()
                ||TotalIndices+IndexCount>std::numeric_limits<std::size_t>::max())
            {Error="mesh scene exceeds supported index space";return false;}
            const std::size_t Base=Scene.Vertices.size();
            Scene.Vertices.resize(Base+VertexCount);
            Scene.Indices.reserve(Scene.Indices.size()+IndexCount);
            for(std::size_t I=Base;I<Base+VertexCount;++I)
                for(float& Value:Scene.Vertices[I].Position)if(!Read(File,Value))
                {Error="tile positions are truncated";return false;}
            for(std::size_t I=Base;I<Base+VertexCount;++I)
            {
                for(float& Value:Scene.Vertices[I].Normal)if(!Read(File,Value))
                {Error="tile normals are truncated";return false;}
                Scene.Vertices[I].Height=Scene.Vertices[I].Position[2];
                Scene.MinimumHeight=std::min(Scene.MinimumHeight,Scene.Vertices[I].Height);
                Scene.MaximumHeight=std::max(Scene.MaximumHeight,Scene.Vertices[I].Height);
            }
            if(HasMaterials)
                for(std::size_t I=Base;I<Base+VertexCount;++I)
                {
                    std::uint16_t Host=0,Ids[2]{},Weights[2]{};
                    if(!Read(File,Host)||!Read(File,Ids[0])||!Read(File,Ids[1])
                        ||!Read(File,Weights[0])||!Read(File,Weights[1]))
                    {Error="tile material attributes are truncated";return false;}
                    Scene.Vertices[I].Material[0]=static_cast<float>(Ids[0]);
                    Scene.Vertices[I].Material[1]=static_cast<float>(Weights[0])/65535.0f;
                    Scene.Vertices[I].Material[2]=static_cast<float>(Ids[1]);
                    Scene.Vertices[I].Material[3]=static_cast<float>(Weights[1])/65535.0f;
                    (void)Host;
                }
            const double TileVoxels=static_cast<double>(32u<<Key[3]);
            const float TileExtentMeters=static_cast<float>(TileVoxels
                *static_cast<double>(Scene.VoxelSizeMicrons)/1000000.0);
            const float TileOriginX=static_cast<float>(static_cast<double>(Key[0])
                *TileExtentMeters);
            const float TileOriginY=static_cast<float>(static_cast<double>(Key[1])
                *TileExtentMeters);
            for(std::size_t I=Base;I<Base+VertexCount;++I)
            {
                FVertex& Vertex=Scene.Vertices[I];
                std::copy(std::begin(Vertex.Position),std::end(Vertex.Position),
                    std::begin(Vertex.WorldPosition));
                Vertex.TileData[0]=TileOriginX;
                Vertex.TileData[1]=TileOriginY;
                Vertex.TileData[2]=TileExtentMeters;
                Vertex.TileData[3]=static_cast<float>(Key[3]);
            }
            const std::size_t IndexOffset=Scene.Indices.size();
            Scene.Indices.resize(IndexOffset+IndexCount);
            if(IndexCount)File.read(reinterpret_cast<char*>(Scene.Indices.data()+IndexOffset),
                static_cast<std::streamsize>(sizeof(std::uint32_t)*IndexCount));
            std::vector<std::uint8_t> Masks(MaskCount);
            if(MaskCount)File.read(reinterpret_cast<char*>(Masks.data()),MaskCount);
            if(!File){Error="tile index or skirt data is truncated";return false;}
            for(std::size_t Index=IndexOffset;Index<Scene.Indices.size();++Index)
            {
                if(Scene.Indices[Index]>=VertexCount)
                {Error="tile index points outside its vertex array";return false;}
                Scene.Indices[Index]+=static_cast<std::uint32_t>(Base);
            }
            Scene.Triangles+=IndexCount/3;
            Scene.SkirtTriangles+=SkirtCount;
            Scene.Tiles.push_back({{Key[0],Key[1],Key[2],Key[3]},IndexOffset,
                IndexCount,std::move(Masks)});
            TotalVertices+=VertexCount;TotalIndices+=IndexCount;
        }
        if(File.peek()!=std::char_traits<char>::eof())
        {Error="mesh scene has trailing bytes";return false;}
        if(Scene.Indices.empty()){Error="mesh scene contains no triangles";return false;}
        std::unordered_set<std::array<std::int32_t,4>,FTileKeyHash> Desired;
        Desired.reserve(Scene.Tiles.size());
        for(const auto& Tile:Scene.Tiles)Desired.insert(Tile.Key);
        std::vector<std::uint32_t> VisibleIndices;
        VisibleIndices.reserve(Scene.Indices.size());
        Scene.SurfaceIndices.reserve(Scene.Indices.size()-
            static_cast<std::size_t>(Scene.SkirtTriangles)*3);
        for(const auto& Tile:Scene.Tiles)
        {
            const std::uint8_t VisibleMask=ComputeVisibleSkirtFaceMask(Tile.Key,Desired);
            for(std::uint32_t Triangle=0;Triangle<Tile.IndexCount/3;++Triangle)
            {
                const std::size_t MaskIndex=Triangle;
                const std::uint8_t FaceMask=MaskIndex<Tile.FaceMasks.size()
                    ?Tile.FaceMasks[MaskIndex]:0;
                if(FaceMask!=0&&(FaceMask&VisibleMask)==0)
                {++Scene.HiddenSameLodSkirtTriangles;continue;}
                if(FaceMask!=0)++Scene.VisibleSkirtTriangles;
                const std::size_t Index=Tile.IndexBegin+static_cast<std::size_t>(Triangle)*3;
                VisibleIndices.insert(VisibleIndices.end(),Scene.Indices.begin()+Index,
                    Scene.Indices.begin()+Index+3);
                if(FaceMask==0)
                    Scene.SurfaceIndices.insert(Scene.SurfaceIndices.end(),
                        Scene.Indices.begin()+Index,Scene.Indices.begin()+Index+3);
            }
        }
        Scene.Indices=std::move(VisibleIndices);
        Scene.DrawnTriangles=Scene.Indices.size()/3;
        return true;
    }

    struct FClipCapResult
    {
        std::uint64_t Intersections=0,Loops=0,OpenContours=0,Triangles=0;
        std::uint64_t DegreeOneNodes=0,BranchNodes=0;
        std::uint64_t ClosedOpenContours=0;
        double OpenGapMetres=0.0,NearestOpenGapMetres=0.0;
    };

    FClipCapResult AddClipPlaneCaps(FScene& Scene,const std::array<double,4>& InputPlane)
    {
        FClipCapResult Result;
        const double NormalLength=std::hypot(InputPlane[0],InputPlane[1],InputPlane[2]);
        if(!(NormalLength>0.0))return Result;
        const std::array<double,3> N={InputPlane[0]/NormalLength,
            InputPlane[1]/NormalLength,InputPlane[2]/NormalLength};
        const double D=InputPlane[3]/NormalLength;
        const std::array<double,3> Reference=std::abs(N[2])<0.85
            ?std::array<double,3>{0.0,0.0,1.0}:std::array<double,3>{0.0,1.0,0.0};
        const auto Cross=[](const std::array<double,3>& A,const std::array<double,3>& B)
        {
            return std::array<double,3>{A[1]*B[2]-A[2]*B[1],
                A[2]*B[0]-A[0]*B[2],A[0]*B[1]-A[1]*B[0]};
        };
        const auto Normalize=[](const std::array<double,3>& A)
        {
            const double Length=std::hypot(A[0],A[1],A[2]);
            return Length>0.0?std::array<double,3>{A[0]/Length,A[1]/Length,A[2]/Length}
                :std::array<double,3>{1.0,0.0,0.0};
        };
        const auto Dot=[](const std::array<double,3>& A,const std::array<double,3>& B)
        {return A[0]*B[0]+A[1]*B[1]+A[2]*B[2];};
        const std::array<double,3> U=Normalize(Cross(N,Reference));
        const std::array<double,3> V=Normalize(Cross(N,U));
        const std::array<double,3> Origin={-D*N[0],-D*N[1],-D*N[2]};
        using FPointKey=std::pair<std::int64_t,std::int64_t>;
        std::map<FPointKey,std::size_t> NodeForPoint;
        std::vector<std::array<double,2>> PlanePoints;
        std::vector<std::array<double,3>> WorldPoints;
        std::set<std::pair<std::size_t,std::size_t>> Segments;
        const auto AddPoint=[&](const std::array<double,3>& Point)
        {
            const std::array<double,3> Relative={Point[0]-Origin[0],
                Point[1]-Origin[1],Point[2]-Origin[2]};
            const std::array<double,2> PlanePoint={Dot(Relative,U),Dot(Relative,V)};
            constexpr double Quantization=1000.0;
            const FPointKey Key{static_cast<std::int64_t>(std::llround(PlanePoint[0]*Quantization)),
                static_cast<std::int64_t>(std::llround(PlanePoint[1]*Quantization))};
            const auto Found=NodeForPoint.find(Key);
            if(Found!=NodeForPoint.end())return Found->second;
            const std::size_t Index=PlanePoints.size();
            NodeForPoint.emplace(Key,Index);PlanePoints.push_back(PlanePoint);
            WorldPoints.push_back(Point);return Index;
        };
        for(std::size_t Triangle=0;Triangle+2<Scene.Indices.size();Triangle+=3)
        {
            std::array<std::array<double,3>,3> Points{};
            std::array<double,3> Distances{};
            for(std::size_t Vertex=0;Vertex<3;++Vertex)
            {
                const auto& Source=Scene.Vertices[Scene.Indices[Triangle+Vertex]];
                Points[Vertex]={Source.WorldPosition[0],Source.WorldPosition[1],Source.WorldPosition[2]};
                Distances[Vertex]=Dot(N,Points[Vertex])+D;
            }
            std::array<std::array<double,3>,3> Hits{};std::size_t HitCount=0;
            for(std::size_t Edge=0;Edge<3;++Edge)
            {
                const std::size_t A=Edge,B=(Edge+1)%3;
                const double DA=Distances[A],DB=Distances[B];
                if(std::abs(DA)<=1.0e-9)
                {
                    bool Duplicate=false;
                    for(std::size_t I=0;I<HitCount;++I)
                        Duplicate=Duplicate||std::hypot(Hits[I][0]-Points[A][0],
                            Hits[I][1]-Points[A][1],Hits[I][2]-Points[A][2])<1.0e-7;
                    if(!Duplicate&&HitCount<Hits.size())Hits[HitCount++]=Points[A];
                }
                if((DA< -1.0e-9&&DB>1.0e-9)||(DA>1.0e-9&&DB< -1.0e-9))
                {
                    const double T=DA/(DA-DB);
                    const std::array<double,3> Hit={
                        std::lerp(Points[A][0],Points[B][0],T),
                        std::lerp(Points[A][1],Points[B][1],T),
                        std::lerp(Points[A][2],Points[B][2],T)};
                    bool Duplicate=false;
                    for(std::size_t I=0;I<HitCount;++I)
                        Duplicate=Duplicate||std::hypot(Hits[I][0]-Hit[0],
                            Hits[I][1]-Hit[1],Hits[I][2]-Hit[2])<1.0e-7;
                    if(!Duplicate&&HitCount<Hits.size())Hits[HitCount++]=Hit;
                }
            }
            if(HitCount!=2)continue;
            const std::size_t A=AddPoint(Hits[0]),B=AddPoint(Hits[1]);
            if(A==B)continue;
            Segments.emplace(std::min(A,B),std::max(A,B));
        }
        Result.Intersections=Segments.size();
        if(Segments.empty())return Result;
        std::vector<std::vector<std::size_t>> Adjacency(PlanePoints.size());
        for(const auto& [A,B]:Segments)
        {Adjacency[A].push_back(B);Adjacency[B].push_back(A);}
        std::vector<std::size_t> DegreeOne;
        for(std::size_t Node=0;Node<Adjacency.size();++Node)
        {
            const auto& Neighbours=Adjacency[Node];
            if(Neighbours.size()==1){++Result.DegreeOneNodes;DegreeOne.push_back(Node);}
            else if(Neighbours.size()>2)++Result.BranchNodes;
        }
        struct FEndpointPair {double Distance=0.0;std::size_t A=0,B=0;};
        std::vector<FEndpointPair> EndpointPairs;
        for(std::size_t I=0;I<DegreeOne.size();++I)
            for(std::size_t J=I+1;J<DegreeOne.size();++J)
            {
                const std::size_t A=DegreeOne[I],B=DegreeOne[J];
                EndpointPairs.push_back({std::hypot(PlanePoints[A][0]-PlanePoints[B][0],
                    PlanePoints[A][1]-PlanePoints[B][1]),A,B});
            }
        std::sort(EndpointPairs.begin(),EndpointPairs.end(),[](const auto& A,const auto& B)
            {return A.Distance<B.Distance;});
        if(!EndpointPairs.empty())Result.NearestOpenGapMetres=EndpointPairs.front().Distance;
        std::set<std::size_t> PairedEndpoints;
        for(const FEndpointPair& Pair:EndpointPairs)
        {
            if(Pair.Distance>1.0)break;
            if(PairedEndpoints.contains(Pair.A)||PairedEndpoints.contains(Pair.B))continue;
            const auto ClosingEdge=std::minmax(Pair.A,Pair.B);
            if(!Segments.emplace(ClosingEdge.first,ClosingEdge.second).second)continue;
            Adjacency[Pair.A].push_back(Pair.B);Adjacency[Pair.B].push_back(Pair.A);
            PairedEndpoints.insert(Pair.A);PairedEndpoints.insert(Pair.B);
            Result.OpenGapMetres=std::max(Result.OpenGapMetres,Pair.Distance);
            ++Result.ClosedOpenContours;
        }
        std::vector<std::map<std::size_t,std::size_t>> Continuation(Adjacency.size());
        for(std::size_t Node=0;Node<Adjacency.size();++Node)
        {
            const auto& Neighbours=Adjacency[Node];
            if(Neighbours.size()==2)
            {
                Continuation[Node][Neighbours[0]]=Neighbours[1];
                Continuation[Node][Neighbours[1]]=Neighbours[0];
                continue;
            }
            if(Neighbours.size()<3)continue;
            struct FTurnPair {double Dot=1.0;std::size_t A=0,B=0;};
            std::vector<FTurnPair> Pairs;
            for(std::size_t I=0;I<Neighbours.size();++I)
                for(std::size_t J=I+1;J<Neighbours.size();++J)
                {
                    const auto& A=PlanePoints[Neighbours[I]];
                    const auto& B=PlanePoints[Neighbours[J]];
                    const double AX=A[0]-PlanePoints[Node][0];
                    const double AY=A[1]-PlanePoints[Node][1];
                    const double BX=B[0]-PlanePoints[Node][0];
                    const double BY=B[1]-PlanePoints[Node][1];
                    const double LengthA=std::hypot(AX,AY),LengthB=std::hypot(BX,BY);
                    if(LengthA<=1.0e-9||LengthB<=1.0e-9)continue;
                    Pairs.push_back({(AX*BX+AY*BY)/(LengthA*LengthB),
                        Neighbours[I],Neighbours[J]});
                }
            std::sort(Pairs.begin(),Pairs.end(),[](const auto& A,const auto& B)
                {return A.Dot<B.Dot;});
            std::set<std::size_t> Paired;
            for(const FTurnPair& Pair:Pairs)
            {
                if(Pair.Dot>0.25||Paired.contains(Pair.A)||Paired.contains(Pair.B))continue;
                Continuation[Node][Pair.A]=Pair.B;Continuation[Node][Pair.B]=Pair.A;
                Paired.insert(Pair.A);Paired.insert(Pair.B);
            }
        }
        std::set<std::pair<std::size_t,std::size_t>> UsedSegments;
        std::vector<std::vector<std::size_t>> Loops;
        for(const auto& Segment:Segments)
        {
            if(UsedSegments.contains(Segment))continue;
            std::vector<std::size_t> Loop{Segment.first};
            std::size_t Previous=Segment.first,Current=Segment.second;
            UsedSegments.insert(Segment);bool Closed=false,Branched=false;
            for(std::size_t Guard=0;Guard<=Segments.size();++Guard)
            {
                if(Current==Loop.front()){Closed=true;break;}
                Loop.push_back(Current);
                const auto ContinuationIt=Continuation[Current].find(Previous);
                if(ContinuationIt==Continuation[Current].end())
                {Branched=true;break;}
                const std::size_t Next=ContinuationIt->second;
                const auto NextEdge=std::minmax(Current,Next);
                if(UsedSegments.contains(NextEdge))
                {Branched=true;break;}
                UsedSegments.insert(NextEdge);Previous=Current;Current=Next;
            }
            if(Closed&&Loop.size()>=3)Loops.push_back(std::move(Loop));
            else if(Branched||!Closed)++Result.OpenContours;
        }
        Result.Loops=Loops.size();
        if(Loops.empty())return Result;

        struct FPlaneEdge {std::array<double,2> A{},B{};};
        std::vector<FPlaneEdge> Edges;
        std::vector<double> Ys;
        for(const auto& Loop:Loops)
            for(std::size_t I=0;I<Loop.size();++I)
            {
                const auto& A=PlanePoints[Loop[I]];
                const auto& B=PlanePoints[Loop[(I+1)%Loop.size()]];
                if(std::hypot(A[0]-B[0],A[1]-B[1])<1.0e-6)continue;
                Edges.push_back({A,B});Ys.push_back(A[1]);Ys.push_back(B[1]);
            }
        std::sort(Ys.begin(),Ys.end());
        Ys.erase(std::unique(Ys.begin(),Ys.end(),[](double A,double B)
            {return std::abs(A-B)<1.0e-6;}),Ys.end());
        const auto AppendTriangle=[&](const std::array<double,2>& A,
            const std::array<double,2>& B,const std::array<double,2>& C)
        {
            const std::array<std::array<double,2>,3> PlaneVertices={A,B,C};
            const std::uint32_t Base=static_cast<std::uint32_t>(Scene.Vertices.size());
            for(const auto& PlaneVertex:PlaneVertices)
            {
                FVertex Vertex{};
                for(std::size_t Axis=0;Axis<3;++Axis)
                {
                    const double World=Origin[Axis]+U[Axis]*PlaneVertex[0]
                        +V[Axis]*PlaneVertex[1]+N[Axis]*0.001;
                    Vertex.Position[Axis]=static_cast<float>(World);
                    Vertex.WorldPosition[Axis]=static_cast<float>(World);
                    Vertex.Normal[Axis]=static_cast<float>(-N[Axis]);
                }
                Vertex.Height=Vertex.Position[2];
                Scene.Vertices.push_back(Vertex);
            }
            Scene.Indices.insert(Scene.Indices.end(),{Base,Base+1,Base+2});
            ++Result.Triangles;
        };
        for(std::size_t Slab=0;Slab+1<Ys.size();++Slab)
        {
            const double Y0=Ys[Slab],Y1=Ys[Slab+1];
            if(Y1-Y0<1.0e-7)continue;
            const double Mid=(Y0+Y1)*0.5;
            struct FCrossing {double X=0.0;const FPlaneEdge* Edge=nullptr;};
            std::vector<FCrossing> Crossings;
            for(const FPlaneEdge& Edge:Edges)
            {
                const double Low=std::min(Edge.A[1],Edge.B[1]);
                const double High=std::max(Edge.A[1],Edge.B[1]);
                if(Mid<=Low||Mid>=High)continue;
                const double T=(Mid-Edge.A[1])/(Edge.B[1]-Edge.A[1]);
                Crossings.push_back({std::lerp(Edge.A[0],Edge.B[0],T),&Edge});
            }
            std::sort(Crossings.begin(),Crossings.end(),[](const auto& A,const auto& B)
                {return A.X<B.X;});
            const auto XAt=[](const FPlaneEdge& Edge,double Y)
            {return std::lerp(Edge.A[0],Edge.B[0],(Y-Edge.A[1])/(Edge.B[1]-Edge.A[1]));};
            for(std::size_t I=0;I+1<Crossings.size();I+=2)
            {
                const FPlaneEdge& Left=*Crossings[I].Edge;
                const FPlaneEdge& Right=*Crossings[I+1].Edge;
                const std::array<double,2> A={XAt(Left,Y0),Y0};
                const std::array<double,2> B={XAt(Right,Y0),Y0};
                const std::array<double,2> C={XAt(Right,Y1),Y1};
                const std::array<double,2> E={XAt(Left,Y1),Y1};
                if(std::abs(A[0]-B[0])+std::abs(E[0]-C[0])<1.0e-7)continue;
                AppendTriangle(A,B,C);AppendTriangle(A,C,E);
            }
        }
        Scene.DrawnTriangles=Scene.Indices.size()/3;
        return Result;
    }

    EGLDisplay MakeDisplay()
    {
        using TGetPlatformDisplay=EGLDisplay(*)(EGLenum,void*,const EGLint*);
        auto GetPlatformDisplay=reinterpret_cast<TGetPlatformDisplay>(
            eglGetProcAddress("eglGetPlatformDisplayEXT"));
#ifdef EGL_PLATFORM_SURFACELESS_MESA
        if(GetPlatformDisplay)
        {
            const EGLDisplay Display=GetPlatformDisplay(EGL_PLATFORM_SURFACELESS_MESA,
                EGL_DEFAULT_DISPLAY,nullptr);
            if(Display!=EGL_NO_DISPLAY)return Display;
        }
#endif
        return eglGetDisplay(EGL_DEFAULT_DISPLAY);
    }

    GLuint CompileShader(GLenum Type,const char* Source)
    {
        const GLuint Shader=glCreateShader(Type);
        glShaderSource(Shader,1,&Source,nullptr);
        glCompileShader(Shader);
        GLint Success=GL_FALSE;
        glGetShaderiv(Shader,GL_COMPILE_STATUS,&Success);
        if(Success==GL_FALSE)
        {
            GLint Length=0;glGetShaderiv(Shader,GL_INFO_LOG_LENGTH,&Length);
            std::string Log(static_cast<std::size_t>(std::max(Length,1)),0);
            glGetShaderInfoLog(Shader,Length,nullptr,Log.data());
            std::fprintf(stderr,"[OffscreenGL] shader compile failed: %s\n",Log.c_str());
            glDeleteShader(Shader);return 0;
        }
        return Shader;
    }

    GLuint MakeProgram()
    {
        static constexpr char VertexSource[]=R"GLSL(#version 300 es
precision highp float;
layout(location=0) in vec3 inCameraPosition;
layout(location=1) in vec3 inNormal;
layout(location=2) in float inHeight;
layout(location=3) in vec3 inWorldPosition;
layout(location=4) in vec4 inTileData;
layout(location=5) in vec4 inMaterial;
out vec3 vNormal;
out float vHeight;
out vec3 vWorldPosition;
out vec4 vTileData;
out vec4 vMaterial;
uniform float uTanHalfHorizontal;
uniform float uTanHalfVertical;
uniform float uNear;
uniform float uFar;
void main()
{
    float z=max(inCameraPosition.z,0.0001);
    float clipZ=((uFar+uNear)/(uFar-uNear))*z-(2.0*uFar*uNear)/(uFar-uNear);
    gl_Position=vec4(inCameraPosition.x/uTanHalfHorizontal,
        inCameraPosition.y/uTanHalfVertical,clipZ,z);
    vNormal=normalize(inNormal);
    vHeight=inHeight;
    vWorldPosition=inWorldPosition;
    vTileData=inTileData;
    vMaterial=inMaterial;
})GLSL";
        static constexpr char FragmentSource[]=R"GLSL(#version 300 es
precision highp float;
in vec3 vNormal;
in float vHeight;
in vec3 vWorldPosition;
in vec4 vTileData;
in vec4 vMaterial;
out vec4 outColor;
uniform int uReadableMode;
uniform int uTileEdges;
uniform int uHeightTint;
uniform float uHeightCenter;
uniform float uHeightRange;
uniform float uTileEdgeWidth;
uniform vec3 uLightDirection;
uniform int uHeadlamp;
uniform int uGrazingLight;
uniform int uClipEnabled;
uniform vec4 uClipPlane;
void main()
{
    if(uClipEnabled!=0&&dot(vec4(vWorldPosition,1.0),uClipPlane)<0.0)discard;
    vec3 n=normalize(vNormal);
    vec3 color;
    float illumination=1.0;
    if(uReadableMode==0)
    {
        vec3 normalColor=abs(n)*0.78+vec3(0.12);
        float band=mod(floor((vHeight+10000.0)/80.0),2.0);
        vec3 bandColor=mix(vec3(0.24,0.31,0.37),vec3(0.39,0.42,0.35),band);
        color=mix(bandColor,normalColor,0.68);
    }
    else
    {
        float heightT=clamp(0.5+(vHeight-uHeightCenter)/max(uHeightRange,0.001),0.0,1.0);
        float smoothHeightT=heightT*heightT*(3.0-2.0*heightT);
        vec3 heightTint=uHeightTint!=0
            ?mix(vec3(0.62,0.78,1.0),vec3(1.0,0.82,0.62),smoothHeightT)
            :vec3(1.0);
        float diffuse=max(dot(n,normalize(uLightDirection)),0.0);
        illumination=uGrazingLight!=0?0.32+1.20*diffuse:
            uHeadlamp!=0?0.78+0.34*diffuse:0.62+0.38*diffuse;
        color=vec3(0.64,0.62,0.56)*heightTint*illumination;
    }
    if(uReadableMode!=0)
    {
        vec3 resourceColor=vec3(0.0);
        float resourceWeight=0.0;
        if(vMaterial.x>0.5)
        {
            float w=clamp(vMaterial.y,0.0,1.0);
            vec3 oreColor=vMaterial.x<1.5?vec3(0.94,0.36,0.12):
                vMaterial.x<2.5?vec3(0.20,0.48,0.92):vec3(0.25,0.88,0.84);
            resourceColor+=oreColor*w;resourceWeight+=w;
        }
        if(vMaterial.z>0.5)
        {
            float w=clamp(vMaterial.w,0.0,1.0);
            vec3 oreColor=vMaterial.z<1.5?vec3(0.94,0.36,0.12):
                vMaterial.z<2.5?vec3(0.20,0.48,0.92):vec3(0.25,0.88,0.84);
            resourceColor+=oreColor*w;resourceWeight+=w;
        }
        if(resourceWeight>0.0)
            color=mix(color,resourceColor/max(resourceWeight,0.001),
                clamp(resourceWeight,0.0,0.95))*illumination;
    }
    if(uTileEdges!=0&&vTileData.z>0.0)
    {
        vec2 local=mod(vWorldPosition.xy-vTileData.xy,vTileData.z);
        vec2 edge=min(local,vec2(vTileData.z)-local);
        float distanceToTileEdge=min(edge.x,edge.y);
        float antialias=max(fwidth(distanceToTileEdge),0.005);
        float line=1.0-smoothstep(uTileEdgeWidth-antialias,
            uTileEdgeWidth+antialias,distanceToTileEdge);
        vec3 edgeColor=mix(vec3(0.16,0.22,0.25),vec3(0.30,0.23,0.14),
            clamp(vTileData.w/8.0,0.0,1.0));
        color=mix(color,edgeColor,0.32*line);
    }
    outColor=vec4(color,1.0);
})GLSL";
        const GLuint VS=CompileShader(GL_VERTEX_SHADER,VertexSource);
        const GLuint FS=CompileShader(GL_FRAGMENT_SHADER,FragmentSource);
        if(!VS||!FS)return 0;
        const GLuint Program=glCreateProgram();
        glAttachShader(Program,VS);glAttachShader(Program,FS);glLinkProgram(Program);
        glDeleteShader(VS);glDeleteShader(FS);
        GLint Success=GL_FALSE;glGetProgramiv(Program,GL_LINK_STATUS,&Success);
        if(Success==GL_FALSE)
        {
            GLint Length=0;glGetProgramiv(Program,GL_INFO_LOG_LENGTH,&Length);
            std::string Log(static_cast<std::size_t>(std::max(Length,1)),0);
            glGetProgramInfoLog(Program,Length,nullptr,Log.data());
            std::fprintf(stderr,"[OffscreenGL] program link failed: %s\n",Log.c_str());
            glDeleteProgram(Program);return 0;
        }
        return Program;
    }

    bool WritePng(const std::filesystem::path& Path,int Width,int Height,
                  const std::vector<std::uint8_t>& Pixels)
    {
        std::error_code Error;
        if(!Path.parent_path().empty())std::filesystem::create_directories(Path.parent_path(),Error);
        if(Error)return false;
        FILE* File=std::fopen(Path.c_str(),"wb");
        if(!File)return false;
        png_structp Png=png_create_write_struct(PNG_LIBPNG_VER_STRING,nullptr,nullptr,nullptr);
        png_infop Info=Png?png_create_info_struct(Png):nullptr;
        if(!Png||!Info){if(Png)png_destroy_write_struct(&Png,nullptr);std::fclose(File);return false;}
        if(setjmp(png_jmpbuf(Png)))
        {png_destroy_write_struct(&Png,&Info);std::fclose(File);return false;}
        png_init_io(Png,File);
        png_set_IHDR(Png,Info,static_cast<png_uint_32>(Width),static_cast<png_uint_32>(Height),
            8,PNG_COLOR_TYPE_RGBA,PNG_INTERLACE_NONE,PNG_COMPRESSION_TYPE_DEFAULT,PNG_FILTER_TYPE_DEFAULT);
        png_write_info(Png,Info);
        const std::size_t Stride=static_cast<std::size_t>(Width)*4;
        for(int Y=Height-1;Y>=0;--Y)
            png_write_row(Png,const_cast<png_bytep>(Pixels.data()+static_cast<std::size_t>(Y)*Stride));
        png_write_end(Png,nullptr);
        png_destroy_write_struct(&Png,&Info);
        const bool Good=std::fclose(File)==0;
        return Good;
    }

    bool Render(FScene& Scene,const std::filesystem::path& Output,int Width,int Height,
                double HorizontalFovDegrees,double HeightRangeMeters,
                double HeightCenterMeters,bool ReadableMode,
                bool HeightTint,bool Headlamp,bool TileEdges,
                double TileEdgeWidthMeters,bool HasLightDirection,
                const std::array<double,3>& LightDirection,bool ClipPlaneEnabled,
                const std::array<double,4>& ClipPlane)
    {
        EGLDisplay Display=MakeDisplay();
        EGLint Major=0,Minor=0;
        if(Display==EGL_NO_DISPLAY||eglInitialize(Display,&Major,&Minor)!=EGL_TRUE)
        {std::fprintf(stderr,"[OffscreenEGL] eglInitialize failed (0x%x)\n",eglGetError());return false;}
        if(eglBindAPI(EGL_OPENGL_ES_API)!=EGL_TRUE)
        {std::fprintf(stderr,"[OffscreenEGL] eglBindAPI failed (0x%x)\n",eglGetError());eglTerminate(Display);return false;}
        const EGLint ConfigAttributes[]={EGL_SURFACE_TYPE,EGL_PBUFFER_BIT,
            EGL_RENDERABLE_TYPE,EGL_OPENGL_ES3_BIT,EGL_RED_SIZE,8,EGL_GREEN_SIZE,8,
            EGL_BLUE_SIZE,8,EGL_ALPHA_SIZE,8,EGL_DEPTH_SIZE,24,EGL_NONE};
        EGLConfig Config{};EGLint ConfigCount=0;
        if(eglChooseConfig(Display,ConfigAttributes,&Config,1,&ConfigCount)!=EGL_TRUE||ConfigCount==0)
        {std::fprintf(stderr,"[OffscreenEGL] no GLES3 pbuffer config (0x%x)\n",eglGetError());eglTerminate(Display);return false;}
        const EGLint SurfaceAttributes[]={EGL_WIDTH,Width,EGL_HEIGHT,Height,EGL_NONE};
        EGLSurface Surface=eglCreatePbufferSurface(Display,Config,SurfaceAttributes);
        const EGLint ContextAttributes[]={EGL_CONTEXT_CLIENT_VERSION,3,EGL_NONE};
        EGLContext Context=eglCreateContext(Display,Config,EGL_NO_CONTEXT,ContextAttributes);
        if(Surface==EGL_NO_SURFACE||Context==EGL_NO_CONTEXT
            ||eglMakeCurrent(Display,Surface,Surface,Context)!=EGL_TRUE)
        {std::fprintf(stderr,"[OffscreenEGL] context/surface creation failed (0x%x)\n",eglGetError());
            if(Context!=EGL_NO_CONTEXT)eglDestroyContext(Display,Context);
            if(Surface!=EGL_NO_SURFACE)eglDestroySurface(Display,Surface);
            eglTerminate(Display);return false;}

        const GLubyte* Renderer=glGetString(GL_RENDERER);
        const GLubyte* GlVersion=glGetString(GL_VERSION);
        std::printf("[OffscreenEGL] renderer=%s version=%s pbuffer=%dx%d\n",
            Renderer?reinterpret_cast<const char*>(Renderer):"unknown",
            GlVersion?reinterpret_cast<const char*>(GlVersion):"unknown",Width,Height);
        GLuint ColorTexture=0,DepthBuffer=0,Framebuffer=0;
        glGenTextures(1,&ColorTexture);glBindTexture(GL_TEXTURE_2D,ColorTexture);
        glTexImage2D(GL_TEXTURE_2D,0,GL_RGBA8,Width,Height,0,GL_RGBA,GL_UNSIGNED_BYTE,nullptr);
        glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER,GL_NEAREST);
        glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_NEAREST);
        glGenRenderbuffers(1,&DepthBuffer);glBindRenderbuffer(GL_RENDERBUFFER,DepthBuffer);
        glRenderbufferStorage(GL_RENDERBUFFER,GL_DEPTH_COMPONENT24,Width,Height);
        glGenFramebuffers(1,&Framebuffer);glBindFramebuffer(GL_FRAMEBUFFER,Framebuffer);
        glFramebufferTexture2D(GL_FRAMEBUFFER,GL_COLOR_ATTACHMENT0,GL_TEXTURE_2D,ColorTexture,0);
        glFramebufferRenderbuffer(GL_FRAMEBUFFER,GL_DEPTH_ATTACHMENT,GL_RENDERBUFFER,DepthBuffer);
        if(glCheckFramebufferStatus(GL_FRAMEBUFFER)!=GL_FRAMEBUFFER_COMPLETE)
        {std::fprintf(stderr,"[OffscreenGL] framebuffer is incomplete\n");return false;}

        constexpr double Pi=3.14159265358979323846;
        const double Yaw=Scene.YawDegrees*Pi/180.0,Pitch=Scene.PitchDegrees*Pi/180.0;
        const std::array<double,3> Forward={std::sin(Yaw)*std::cos(Pitch),
            std::cos(Yaw)*std::cos(Pitch),std::sin(Pitch)};
        const std::array<double,3> Right={std::cos(Yaw),-std::sin(Yaw),0.0};
        const std::array<double,3> Up={Right[1]*Forward[2]-Right[2]*Forward[1],
            Right[2]*Forward[0]-Right[0]*Forward[2],Right[0]*Forward[1]-Right[1]*Forward[0]};
        const double TanHalfHorizontal=std::tan(HorizontalFovDegrees*Pi/360.0);
        const double Aspect=static_cast<double>(Width)/Height;
        const double TanHalfVertical=TanHalfHorizontal/Aspect;
        const double Near=0.1,Far=20000.0;
        for(FVertex& Vertex:Scene.Vertices)
        {
            const double DX=Vertex.Position[0]-Scene.CameraMeters[0];
            const double DY=Vertex.Position[1]-Scene.CameraMeters[1];
            const double DZ=Vertex.Position[2]-Scene.CameraMeters[2];
            const float CameraX=static_cast<float>(DX*Right[0]+DY*Right[1]+DZ*Right[2]);
            const float CameraY=static_cast<float>(DX*Up[0]+DY*Up[1]+DZ*Up[2]);
            const float CameraZ=static_cast<float>(DX*Forward[0]+DY*Forward[1]+DZ*Forward[2]);
            Vertex.Position[0]=CameraX;Vertex.Position[1]=CameraY;Vertex.Position[2]=CameraZ;
        }
        const GLuint Program=MakeProgram();
        if(!Program)return false;
        GLuint Vao=0,Vbo=0,Ebo=0;
        glGenVertexArrays(1,&Vao);glBindVertexArray(Vao);
        glGenBuffers(1,&Vbo);glBindBuffer(GL_ARRAY_BUFFER,Vbo);
        glBufferData(GL_ARRAY_BUFFER,static_cast<GLsizeiptr>(Scene.Vertices.size()*sizeof(FVertex)),
            Scene.Vertices.data(),GL_STATIC_DRAW);
        glGenBuffers(1,&Ebo);glBindBuffer(GL_ELEMENT_ARRAY_BUFFER,Ebo);
        glBufferData(GL_ELEMENT_ARRAY_BUFFER,static_cast<GLsizeiptr>(Scene.Indices.size()*sizeof(std::uint32_t)),
            Scene.Indices.data(),GL_STATIC_DRAW);
        glEnableVertexAttribArray(0);glVertexAttribPointer(0,3,GL_FLOAT,GL_FALSE,sizeof(FVertex),
            reinterpret_cast<void*>(offsetof(FVertex,Position)));
        glEnableVertexAttribArray(1);glVertexAttribPointer(1,3,GL_FLOAT,GL_FALSE,sizeof(FVertex),
            reinterpret_cast<void*>(offsetof(FVertex,Normal)));
        glEnableVertexAttribArray(2);glVertexAttribPointer(2,1,GL_FLOAT,GL_FALSE,sizeof(FVertex),
            reinterpret_cast<void*>(offsetof(FVertex,Height)));
        glEnableVertexAttribArray(3);glVertexAttribPointer(3,3,GL_FLOAT,GL_FALSE,sizeof(FVertex),
            reinterpret_cast<void*>(offsetof(FVertex,WorldPosition)));
        glEnableVertexAttribArray(4);glVertexAttribPointer(4,4,GL_FLOAT,GL_FALSE,sizeof(FVertex),
            reinterpret_cast<void*>(offsetof(FVertex,TileData)));
        glEnableVertexAttribArray(5);glVertexAttribPointer(5,4,GL_FLOAT,GL_FALSE,sizeof(FVertex),
            reinterpret_cast<void*>(offsetof(FVertex,Material)));
        glViewport(0,0,Width,Height);glDisable(GL_CULL_FACE);glEnable(GL_DEPTH_TEST);
        // A covered empty ray represents open air. Keep it visibly blue so it cannot
        // be mistaken for the old near-black result from a scene that lacked tiles.
        glDepthFunc(GL_LEQUAL);glClearColor(0.30f,0.48f,0.62f,1.0f);
        glClearDepthf(1.0f);glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
        glUseProgram(Program);
        glUniform1f(glGetUniformLocation(Program,"uTanHalfHorizontal"),static_cast<float>(TanHalfHorizontal));
        glUniform1f(glGetUniformLocation(Program,"uTanHalfVertical"),static_cast<float>(TanHalfVertical));
        glUniform1f(glGetUniformLocation(Program,"uNear"),static_cast<float>(Near));
        glUniform1f(glGetUniformLocation(Program,"uFar"),static_cast<float>(Far));
        glUniform1i(glGetUniformLocation(Program,"uReadableMode"),ReadableMode?1:0);
        glUniform1i(glGetUniformLocation(Program,"uTileEdges"),TileEdges?1:0);
        glUniform1i(glGetUniformLocation(Program,"uHeightTint"),HeightTint?1:0);
        glUniform1f(glGetUniformLocation(Program,"uHeightCenter"),
            static_cast<float>(HeightCenterMeters));
        glUniform1f(glGetUniformLocation(Program,"uHeightRange"),
            static_cast<float>(HeightRangeMeters));
        glUniform1f(glGetUniformLocation(Program,"uTileEdgeWidth"),
            static_cast<float>(TileEdgeWidthMeters));
        glUniform1i(glGetUniformLocation(Program,"uClipEnabled"),ClipPlaneEnabled?1:0);
        glUniform1i(glGetUniformLocation(Program,"uHeadlamp"),Headlamp?1:0);
        glUniform1i(glGetUniformLocation(Program,"uGrazingLight"),HasLightDirection?1:0);
        glUniform4f(glGetUniformLocation(Program,"uClipPlane"),
            static_cast<float>(ClipPlane[0]),static_cast<float>(ClipPlane[1]),
            static_cast<float>(ClipPlane[2]),static_cast<float>(ClipPlane[3]));
        const std::array<double,3> ActiveLightDirection=HasLightDirection?LightDirection:
            Headlamp?std::array<double,3>{-Forward[0],-Forward[1],-Forward[2]}:
                std::array<double,3>{-0.42,-0.58,0.70};
        glUniform3f(glGetUniformLocation(Program,"uLightDirection"),
            static_cast<float>(ActiveLightDirection[0]),
            static_cast<float>(ActiveLightDirection[1]),
            static_cast<float>(ActiveLightDirection[2]));
        glDrawElements(GL_TRIANGLES,static_cast<GLsizei>(Scene.Indices.size()),GL_UNSIGNED_INT,nullptr);
        glFinish();
        std::vector<std::uint8_t> Pixels(static_cast<std::size_t>(Width)*Height*4);
        glReadPixels(0,0,Width,Height,GL_RGBA,GL_UNSIGNED_BYTE,Pixels.data());
        const GLenum GlError=glGetError();
        const bool PngGood=GlError==GL_NO_ERROR&&WritePng(Output,Width,Height,Pixels);
        if(!PngGood)std::fprintf(stderr,"[OffscreenGL] render or PNG write failed (GL=0x%x)\n",GlError);
        else
        {
            std::error_code SizeError;
            const auto Bytes=std::filesystem::file_size(Output,SizeError);
            std::printf("[OffscreenCapture] source_triangles=%llu drawn_triangles=%llu tiles=%u "
                "skirts=%llu visible_skirt_triangles=%llu hidden_same_lod_skirt_triangles=%llu "
                "shading=%s height_tint=%s tile_edges=%s edge_width_m=%.3f png_bytes=%llu status=PASS path=%s\n",
                static_cast<unsigned long long>(Scene.Triangles),
                static_cast<unsigned long long>(Scene.DrawnTriangles),Scene.TileCount,
                static_cast<unsigned long long>(Scene.SkirtTriangles),
                static_cast<unsigned long long>(Scene.VisibleSkirtTriangles),
                static_cast<unsigned long long>(Scene.HiddenSameLodSkirtTriangles),
                ReadableMode?"readable":"normal",HeightTint?"on":"off",
                TileEdges?"on":"off",TileEdgeWidthMeters,
                static_cast<unsigned long long>(SizeError?0:Bytes),Output.string().c_str());
        }
        glDeleteProgram(Program);glDeleteBuffers(1,&Ebo);glDeleteBuffers(1,&Vbo);
        glDeleteVertexArrays(1,&Vao);glDeleteFramebuffers(1,&Framebuffer);
        glDeleteRenderbuffers(1,&DepthBuffer);glDeleteTextures(1,&ColorTexture);
        eglMakeCurrent(Display,EGL_NO_SURFACE,EGL_NO_SURFACE,EGL_NO_CONTEXT);
        eglDestroyContext(Display,Context);eglDestroySurface(Display,Surface);eglTerminate(Display);
        return PngGood;
    }
}

int main(int argc,char** argv)
{
    VoxelForgePrintSourceStamp();
    if(argc==2&&std::string_view(argv[1])=="--print-source-stamp")return 0;
    if(argc<3)
    {
        std::fprintf(stderr,"Usage: %s scene.vfmesh output.png [--size=WIDTHxHEIGHT] [--fov-deg=90] [--height-range-m=N] [--height-center-m=N] [--camera-m=X,Y,Z] [--yaw-deg=N] [--pitch-deg=N] [--shade=readable|normal] [--height-tint=on|off] [--headlamp=on|off] [--light-direction=X,Y,Z] [--tile-edges=on|off] [--tile-skirts=on|off] [--tile-edge-width-m=N] [--clip-plane=NX,NY,NZ,D] [--focus-region=X,Y,Z,R]\n",argv[0]);
        return 2;
    }
    int Width=2560,Height=1440;double Fov=90.0,HeightRangeMeters=64.0;
    double HeightCenterMeters=0.0;bool HasHeightCenter=false;
    bool HasYaw=false,HasPitch=false;double Yaw=0.0,Pitch=0.0;
    bool HasCamera=false;std::array<double,3> Camera{};
    bool HasLightDirection=false;std::array<double,3> LightDirection{};
    bool ReadableMode=true,HeightTint=true,TileEdges=false,TileSkirts=true,Headlamp=false;
    double TileEdgeWidthMeters=0.15;
    bool ClipPlaneEnabled=false;std::array<double,4> ClipPlane{};
    bool FocusRegionEnabled=false;std::array<double,3> FocusCenter{};double FocusRadius=0.0;
    for(int I=3;I<argc;++I)
    {
        const std::string_view Arg(argv[I]);
        if(Arg.starts_with("--size="))
        {
            const std::string Value(Arg.substr(7));
            if(std::sscanf(Value.c_str(),"%dx%d",&Width,&Height)!=2||Width<=0||Height<=0)
            {std::fprintf(stderr,"Invalid image size\n");return 2;}
        }
        else if(Arg.starts_with("--fov-deg="))
        {
            char* End=nullptr;const std::string Value(Arg.substr(10));Fov=std::strtod(Value.c_str(),&End);
            if(End==Value.c_str()||*End!='\0'||Fov<=1.0||Fov>=179.0)
            {std::fprintf(stderr,"Invalid horizontal FOV\n");return 2;}
        }
        else if(Arg.starts_with("--height-range-m="))
        {
            char* End=nullptr;const std::string Value(Arg.substr(17));
            HeightRangeMeters=std::strtod(Value.c_str(),&End);
            if(End==Value.c_str()||*End!='\0'||!std::isfinite(HeightRangeMeters)
                ||HeightRangeMeters<=0.0||HeightRangeMeters>1000000.0)
            {std::fprintf(stderr,"Invalid height tint range\n");return 2;}
        }
        else if(Arg.starts_with("--height-center-m="))
        {
            char* End=nullptr;const std::string Value(Arg.substr(18));
            HeightCenterMeters=std::strtod(Value.c_str(),&End);
            if(End==Value.c_str()||*End!='\0'||!std::isfinite(HeightCenterMeters))
            {std::fprintf(stderr,"Invalid height tint center\n");return 2;}
            HasHeightCenter=true;
        }
        else if(Arg.starts_with("--camera-m="))
        {
            const std::string Value(Arg.substr(11));
            if(std::sscanf(Value.c_str(),"%lf,%lf,%lf",&Camera[0],&Camera[1],&Camera[2])!=3
                ||!std::isfinite(Camera[0])||!std::isfinite(Camera[1])||!std::isfinite(Camera[2]))
            {std::fprintf(stderr,"Invalid camera position\n");return 2;}
            HasCamera=true;
        }
        else if(Arg.starts_with("--focus-region="))
        {
            const std::string Value(Arg.substr(15));
            if(std::sscanf(Value.c_str(),"%lf,%lf,%lf,%lf",&FocusCenter[0],&FocusCenter[1],
                    &FocusCenter[2],&FocusRadius)!=4
                ||!std::isfinite(FocusCenter[0])||!std::isfinite(FocusCenter[1])
                ||!std::isfinite(FocusCenter[2])||!std::isfinite(FocusRadius)
                ||FocusRadius<=0.0)
            {std::fprintf(stderr,"Invalid focus region\n");return 2;}
            FocusRegionEnabled=true;
        }
        else if(Arg.starts_with("--yaw-deg=")||Arg.starts_with("--pitch-deg="))
        {
            const bool IsYaw=Arg.starts_with("--yaw-deg=");
            const std::size_t Prefix=IsYaw?10:12;
            const std::string Value(Arg.substr(Prefix));char* End=nullptr;
            const double Parsed=std::strtod(Value.c_str(),&End);
            if(End==Value.c_str()||*End!='\0'||!std::isfinite(Parsed))
            {std::fprintf(stderr,"Invalid camera angle\n");return 2;}
            if(IsYaw){Yaw=Parsed;HasYaw=true;}else{Pitch=Parsed;HasPitch=true;}
        }
        else if(Arg.starts_with("--shade="))
        {
            const std::string_view Mode=Arg.substr(8);
            if(Mode=="readable")ReadableMode=true;
            else if(Mode=="normal")ReadableMode=false;
            else {std::fprintf(stderr,"Invalid shading mode\n");return 2;}
        }
        else if(Arg.starts_with("--height-tint="))
        {
            const std::string_view Mode=Arg.substr(14);
            if(Mode=="on")HeightTint=true;
            else if(Mode=="off")HeightTint=false;
            else {std::fprintf(stderr,"Invalid height-tint mode\n");return 2;}
        }
        else if(Arg.starts_with("--headlamp="))
        {
            const std::string_view Mode=Arg.substr(11);
            if(Mode=="on")Headlamp=true;
            else if(Mode=="off")Headlamp=false;
            else {std::fprintf(stderr,"Invalid headlamp mode\n");return 2;}
        }
        else if(Arg.starts_with("--light-direction="))
        {
            const std::string Value(Arg.substr(18));
            if(std::sscanf(Value.c_str(),"%lf,%lf,%lf",&LightDirection[0],
                &LightDirection[1],&LightDirection[2])!=3
                ||!std::isfinite(LightDirection[0])||!std::isfinite(LightDirection[1])
                ||!std::isfinite(LightDirection[2])
                ||std::hypot(LightDirection[0],LightDirection[1],LightDirection[2])<=0.0)
            {std::fprintf(stderr,"Invalid light direction; use X,Y,Z with non-zero length\n");return 2;}
            HasLightDirection=true;
        }
        else if(Arg.starts_with("--tile-edges="))
        {
            const std::string_view Mode=Arg.substr(13);
            if(Mode=="on")TileEdges=true;
            else if(Mode=="off")TileEdges=false;
            else {std::fprintf(stderr,"Invalid tile-edge mode\n");return 2;}
        }
        else if(Arg.starts_with("--tile-skirts="))
        {
            const std::string_view Mode=Arg.substr(14);
            if(Mode=="on")TileSkirts=true;
            else if(Mode=="off")TileSkirts=false;
            else {std::fprintf(stderr,"Invalid tile-skirt mode\n");return 2;}
        }
        else if(Arg.starts_with("--tile-edge-width-m="))
        {
            char* End=nullptr;const std::string Value(Arg.substr(20));
            TileEdgeWidthMeters=std::strtod(Value.c_str(),&End);
            if(End==Value.c_str()||*End!='\0'||!std::isfinite(TileEdgeWidthMeters)
                ||TileEdgeWidthMeters<=0.0||TileEdgeWidthMeters>5.0)
            {std::fprintf(stderr,"Invalid tile-edge width\n");return 2;}
        }
        else if(Arg.starts_with("--clip-plane="))
        {
            const std::string Value(Arg.substr(13));
            if(std::sscanf(Value.c_str(),"%lf,%lf,%lf,%lf",&ClipPlane[0],&ClipPlane[1],
                &ClipPlane[2],&ClipPlane[3])!=4
                ||!std::isfinite(ClipPlane[0])||!std::isfinite(ClipPlane[1])
                ||!std::isfinite(ClipPlane[2])||!std::isfinite(ClipPlane[3])
                ||std::hypot(ClipPlane[0],ClipPlane[1],ClipPlane[2])<=0.0)
            {std::fprintf(stderr,"Invalid clip plane; use NX,NY,NZ,D with a non-zero normal\n");return 2;}
            ClipPlaneEnabled=true;
        }
        else {std::fprintf(stderr,"Unknown option: %.*s\n",static_cast<int>(Arg.size()),Arg.data());return 2;}
    }
    FScene Scene;std::string Error;
    if(!LoadScene(argv[1],Scene,Error))
    {std::fprintf(stderr,"[OffscreenCapture] FAIL: %s\n",Error.c_str());return 2;}
    if(HasCamera)Scene.CameraMeters=Camera;
    if(HasYaw)Scene.YawDegrees=Yaw;
    if(HasPitch)Scene.PitchDegrees=Pitch;
    if(!HasHeightCenter)HeightCenterMeters=Scene.CameraMeters[2];
    const std::size_t CapIndexBegin=Scene.Indices.size();
    FClipCapResult ClipCap;
    if(ClipPlaneEnabled)
    {
        ClipCap=AddClipPlaneCaps(Scene,ClipPlane);
        const bool Capped=ClipCap.Loops>0&&ClipCap.Triangles>0;
        std::printf("[OffscreenClipCap] segments=%llu loops=%llu open_contours=%llu closed_open_contours=%llu degree_one=%llu branch_nodes=%llu open_gap_m=%.3f nearest_open_gap_m=%.3f cap_triangles=%llu status=%s\n",
            static_cast<unsigned long long>(ClipCap.Intersections),
            static_cast<unsigned long long>(ClipCap.Loops),
            static_cast<unsigned long long>(ClipCap.OpenContours),
            static_cast<unsigned long long>(ClipCap.ClosedOpenContours),
            static_cast<unsigned long long>(ClipCap.DegreeOneNodes),
            static_cast<unsigned long long>(ClipCap.BranchNodes),
            ClipCap.OpenGapMetres,ClipCap.NearestOpenGapMetres,
            static_cast<unsigned long long>(ClipCap.Triangles),Capped?"PASS":"FAIL");
        if(!Capped)return 3;
    }
    if(!TileSkirts)
    {
        const std::size_t SurfaceTriangleCount=Scene.SurfaceIndices.size()/3;
        std::vector<std::uint32_t> SurfaceAndCaps=std::move(Scene.SurfaceIndices);
        // Caps were appended after the original mesh indices. Keep them while dropping
        // the LOD seam skirts, which are useful for streamed play but obscure close captures.
        if(ClipPlaneEnabled&&Scene.Indices.size()>CapIndexBegin)
            SurfaceAndCaps.insert(SurfaceAndCaps.end(),Scene.Indices.begin()+CapIndexBegin,
                Scene.Indices.end());
        Scene.Indices=std::move(SurfaceAndCaps);
        Scene.DrawnTriangles=Scene.Indices.size()/3;
        std::printf("[OffscreenMeshDisplay] tile_skirts=off kept_surface_triangles=%llu "
            "cap_triangles=%llu status=PASS\n",
            static_cast<unsigned long long>(SurfaceTriangleCount),
            static_cast<unsigned long long>(Scene.Indices.size()/3-
                SurfaceTriangleCount));
    }
    if(FocusRegionEnabled)
    {
        std::vector<std::uint32_t> FocusedIndices;
        FocusedIndices.reserve(Scene.Indices.size());
        const double RadiusSquared=FocusRadius*FocusRadius;
        for(std::size_t Triangle=0;Triangle+2<Scene.Indices.size();Triangle+=3)
        {
            std::array<double,3> Center{};
            for(std::size_t Corner=0;Corner<3;++Corner)
            {
                const auto& Vertex=Scene.Vertices[Scene.Indices[Triangle+Corner]];
                for(std::size_t Axis=0;Axis<3;++Axis)
                    Center[Axis]+=Vertex.WorldPosition[Axis]/3.0;
            }
            const double DX=Center[0]-FocusCenter[0],DY=Center[1]-FocusCenter[1],
                DZ=Center[2]-FocusCenter[2];
            if(DX*DX+DY*DY+DZ*DZ>RadiusSquared)continue;
            FocusedIndices.insert(FocusedIndices.end(),Scene.Indices.begin()+Triangle,
                Scene.Indices.begin()+Triangle+3);
        }
        if(FocusedIndices.empty())
            std::printf("[OffscreenFocusRegion] center_m=(%.3f,%.3f,%.3f) radius_m=%.3f "
                "triangles=0 status=EMPTY_FALLBACK full_scene_triangles=%llu\n",
                FocusCenter[0],FocusCenter[1],FocusCenter[2],FocusRadius,
                static_cast<unsigned long long>(Scene.DrawnTriangles));
        else
        {
            Scene.Indices=std::move(FocusedIndices);
            Scene.DrawnTriangles=Scene.Indices.size()/3;
            std::printf("[OffscreenFocusRegion] center_m=(%.3f,%.3f,%.3f) radius_m=%.3f "
                "triangles=%llu status=PASS\n",FocusCenter[0],FocusCenter[1],FocusCenter[2],
                FocusRadius,static_cast<unsigned long long>(Scene.DrawnTriangles));
        }
    }
    std::printf("[OffscreenScene] voxels_per_meter=%.6f camera_m=(%.3f,%.3f,%.3f) yaw=%.2f pitch=%.2f height_m=%.3f..%.3f height_center_m=%.3f shade=%s height_tint=%s tile_edges=%s vertices=%zu triangles=%llu\n",
        1000000.0/Scene.VoxelSizeMicrons,Scene.CameraMeters[0],Scene.CameraMeters[1],Scene.CameraMeters[2],
        Scene.YawDegrees,Scene.PitchDegrees,Scene.MinimumHeight,Scene.MaximumHeight,
        HeightCenterMeters,
        ReadableMode?"readable":"normal",HeightTint?"on":"off",
        TileEdges?"on":"off",Scene.Vertices.size(),
        static_cast<unsigned long long>(Scene.Triangles));
    return Render(Scene,argv[2],Width,Height,Fov,HeightRangeMeters,HeightCenterMeters,ReadableMode,
        HeightTint,Headlamp,TileEdges,
        TileEdgeWidthMeters,HasLightDirection,LightDirection,ClipPlaneEnabled,ClipPlane)?0:3;
}