/* NOTE: THIS FILE CONTAINS A REDUCED PORTION OF THE PORTAL-RENDERING 2 ENGINE,
 * MOST IMPORTANTLY WITHOUT ANY PORTAL RENDERING (USES Z-BUFFERING).
 * I TRIED TO MAKE IT AS MINIMAL AS POSSIBLE...
 * BUT UNFORTUNATELY DON'T HAVE THE TIME TO MAKE IT INTO A VIDEO RIGHT NOW,
 * AS EVEN THE REDUCED SOURCE ENDED UP BEING ABOUT 2X THE LENGTH I ESTIMATED.
 *
 * THIS IS NOT THE FULL ENGINE. IT IS JUST A SNEAK PEEK.
 * Compile with -std=c++17, and link with SDL 1.2. Does not require external files.
 * Controls are Descent-style. No collision checking is implemented.
 *
 * -Bisqwit, 2016-09-22
 *   Updated 2016-09-25: Renders now whole polygons instead of triangles & tesselation
 *   Updated 2016-10-13: Removed rotation matrices; using only quaternions for simpler code
 *
 * If c++17 is not available in your compiler, watch https://www.youtube.com/watch?v=wrwwa68JXNk
 * for instructions how to make it c++14 compatible.
 */

#include <utility>     // std::swap, std::move, std::forward
#include <algorithm>   // std::swap, std::min
#include <cstdlib>     // std::abs
#include <type_traits> // std::decay
#include <array>
#include <cmath>
#include <vector>
#include <chrono>
#include <string>
#include <map>
#include <SDL.h>

// Very bare-bones vector mathematics module. See https://en.wikipedia.org/wiki/Euclidean_vector
template<typename T, unsigned dim>
struct vec
{
    std::array<T,dim> d;
    template<typename...R> vec(R&&... r) : d{} { Assign<0u>(std::forward<R>(r)...); }
    template<typename F> static vec From(F&& f) { vec res; for(unsigned a=0; a<dim; ++a) res.d[a] = f(a); return res; }
    // Basic vector arithmetic operators (u+b, u-b, u*b, u/b) with a vector or scalar b, in-place or const
    #define d(o) \
        vec& operator o##= (const vec& b) { for(unsigned a=0; a<dim; ++a) d[a] o##= b.d[a]; return *this; } \
        vec& operator o##= (T b)          { for(unsigned a=0; a<dim; ++a) d[a] o##= b;      return *this; } \
        template<typename U> vec operator o(U&& b) const { return vec(*this) o##= std::forward<U>(b); }
    d(+) d(-) d(*) d(/)
    #undef d
    // Individual element access (const only)
    float operator[] (unsigned n) const { return d[n]; }
    // Normalization (û), horizontal sum, length (‖u‖), dot product (u·b), and cross product (u×b)
    vec Normalized() const                { auto l = Length(); return std::abs(l) > T(1e-10) ? *this/l : vec{}; }
    T HorizSum() const                    { T res{}; for(auto f: d) res += f; return res; }
    T Length() const                      { return std::sqrt(DotProduct(*this)); }
    T DotProduct(const vec& b) const      { return (*this * b).HorizSum(); }
    auto CrossProduct(const vec& b) const
    {
        // Three-dimensional cross-product, written in a manner that encourages SIMD optimizations:
        if constexpr (dim==3)
            return [m = vec<T,8>{d[1],d[2],d[0],0, d[2],d[0],d[1],0} * vec<T,8>{b[2],b[0],b[1],0, b[1],b[2],b[0],0}]
                { return (m.template Slice<4>(0) - m.template Slice<4>(4)).template Slice<3>(); }();
        // Two-dimensional "cross-product" is the Z component of the cross-product between two 3D vectors where Z=0.
        // Alternatively it can be thought as the dotproduct of a vector and another vector but rotated 90 degrees.
        // It is very useful in computer graphics, even though math pedantics insist it doesn't exist.
        if constexpr (dim==2) return DotProduct({b[1], -b[0]});
    }
    // Slice(): Retrieve a sub-vector (such as vec<2> containing first two coordinates of vec<3>)
    template<unsigned n> auto Slice(unsigned first=0) const { vec<T,n> res; std::copy_n(&d[first], n, &res.d[0]); return res; }
    // PartialAssign(): Assign values into a subset of coordinates within the vector
    template<unsigned first, typename...R> void PartialAssign(R&&... r) { Assign<first>(std::forward<R>(r)..., vec<T,0>{}); }
private:
    template<unsigned index, typename... R>
        void Assign(T v, R&&... r)                      { d[index]=v;                       Assign<index+1>(std::forward<R>(r)...); }
    template<unsigned index, typename U, unsigned n, typename... R>
        void Assign(const vec<U,n>& v, R&&... r)        {                                   Assign<index>(v.d, std::forward<R>(r)...); }
    template<unsigned index, typename U, std::size_t n, typename... R>
        void Assign(const std::array<U,n>& v, R&&... r) { std::copy_n(&v[0], n, &d[index]); Assign<index+n>(std::forward<R>(r)...); }
    template<unsigned index>
        void Assign() { static_assert(index==0 || index==dim, "Assign: Too few parameters"); }
    template<unsigned>
        void Assign(vec<T,0>) { } // Warningless sentinel used in PartialAssign()
};

using xy  = vec<float,2>; // 2D coordinates: A vector with two dimensions
using xyz = vec<float,3>; // 3D coordinates: A vector with three dimensions.

// Quaternion stuff. See https://en.wikipedia.org/wiki/Quaternion
struct quat: public vec<float,4> // Quaternion: real,I,J,K.
{
    template<typename... R> quat(R&&... r) : vec{std::forward<R>(r)...} {}
    static quat FromAngle(float theta, const xyz& v) { return {std::cos(theta * 0.5f), v * std::sin(theta * 0.5f)}; }
    // Inversion
    quat operator-() const { return {d[0],-d[1],-d[2],-d[3]}; }
    // Non-commutative quaternion-by-quaternion multiplication
    quat operator*(const quat& o) const
    {
        vec<float,16> m{d,d,d,d}; m *= vec<float,16>{-o, o[1],o[0],o[3],-o[2], o[2],-o[3],o[0],o[1], o[3],o[2],-o[1],o[0]};
        return { m.Slice<4>(0).HorizSum(), m.Slice<4>(4).HorizSum(), m.Slice<4>(8).HorizSum(), m.Slice<4>(12).HorizSum() };
    }
    // Rotating a vector using a quaternion. Slower, but shorter, than making an intermediate rotation matrix and using that.
    xyz Rotate(const xyz& v) const { return (*this * quat{0,v} * -*this).Slice<3>(1); }
};

// 3D Plane stuff. See https://en.wikipedia.org/wiki/Plane_%28geometry%29
struct Plane: private vec<float,4> // 3D plane: normal-vector, distance from origo
{
    template<typename... R> Plane(R&&... r) : vec{std::forward<R>(r)...} {}
    static Plane FromPoints(const xyz& a, const xyz& b, const xyz& c)
    {
        xyz normal = (b-a).CrossProduct(c-a).Normalized();
        return { normal, normal.DotProduct(a) };
    }
    float DistanceTo(const xyz& p) const { return Slice<3>().DotProduct(p) - (*this)[3]; }
};


/* Convex polygon rasterization primitive (without clipping). */
// Limitations: Polygon must be convex, and the maximum number of vertices on the same scanline is two.
// Other than that, the number of vertexes is unlimited (i.e. not only triangles).
template<bool TwosidedSupport,      // Whether we support both clockwise & ccw polygons
         bool PerspectiveCorrected, // Whether we do perspective correction. If not, GetZ is ignored.
         unsigned Nprops, // Number of properties to interpolate inside the triangle (such as uv coordinates etc.)
         typename PointType, typename XYfunc, typename Zfunc, typename PropFunc, typename Plotter>
inline void RasterizePolygon
    (const PointType* points, unsigned npoints, // List of corner vertices
     XYfunc&& GetXY,      // A functor for extracting xy<sometype> from a PointType
     Zfunc&& GetZ,        // A functor for extracting depth from a PointType. Ignored if PerspectiveCorrected==false
     PropFunc&& GetProp,  // A functor for extracting a single property from PointType for interpolation
     Plotter&& Plot)      // A functor for plotting a single pixel
{
    // Find the top-leftmost vertex.
    const PointType* topleft = points;
    const auto Less = [](auto a, auto b) { return a[1] != b[1] ? a[1] < b[1] : a[0] < b[0]; };
    for(unsigned n=1; n<npoints; ++n)
        if(Less(GetXY(points[n]), GetXY(*topleft)))
            topleft = &points[n];

    // Do the rasterization from top to bottom, keeping track of the left and right sides.
    const PointType* left[2] {topleft, topleft}, *right[2] {topleft, topleft};
    auto Prev = [b=points, e=points+npoints-1](const PointType* p) { return p != b ? p - 1 : e; };
    auto Next = [b=points, e=points+npoints-1](const PointType* p) { return p != e ? p + 1 : b; };

    // Determine whether this is a clockwise or counter-clockwise polygon, using a crossproduct.
    // In clockwise polygons, increasing vertex numbers go on the right side and decreasing
    // vertex numbers go on the left side. In counter-clockwise polygons, it's vice versa.
    const bool clockwise = TwosidedSupport
        ? (GetXY(*Prev(topleft)) - GetXY(*topleft)).CrossProduct(GetXY(*Next(topleft)) - GetXY(*topleft)) < 0
        : true; // If the template parameter is false, support only clockwise polygons for slightly faster rendering.

    constexpr unsigned YSlopeSize = (Nprops+1+PerspectiveCorrected), Z=Nprops+0, X=Nprops+PerspectiveCorrected;
    constexpr unsigned XSlopeSize = (Nprops+0+PerspectiveCorrected);

    // Create slopes for the left and right side. A slope is a pair of vectors: Current value and step-increment (second).
    struct Slope { vec<float,YSlopeSize> current, increment; } sleft, sright;
    auto Update = [&,Prev,Next](auto& slope, const PointType* (&p)[2], unsigned &next_y, const PointType* opposite, bool dir)
    {
        auto& from = p[0], &to = p[1];
        from  = to;                        // The "endpoint" for this side becomes the new startpoint
        to    = dir ? Next(to) : Prev(to); // The next vertex in chain becomes the endpoint for this side
        if(to == opposite) return false;   // Terminate, if we reach the opposite edge.
        next_y = GetXY(*to)[1] + 1;        // Remember the next Y coordinate where we must recalculate slope.

        // Create the slope.
        float z1 = 1.f;
        if constexpr(PerspectiveCorrected) z1 /= GetZ(*from);
        slope.current   = vec<float,YSlopeSize>::From([&](unsigned n) { return n<Nprops ? GetProp(*from,n)*z1 : n==X ? GetXY(*from)[0] : z1; });
        float z2 = 1.f;
        if constexpr(PerspectiveCorrected) z2 /= GetZ(*to);
        slope.increment = (vec<float,YSlopeSize>::From([&](unsigned n) { return n<Nprops ? GetProp(*to,n)*z2 : n==X ? GetXY(*to)[0] : z2; })
                             - slope.current) * (1.f / (GetXY(*to)[1] - GetXY(*from)[1]));
        return true;
    };

    // The main rasterization loop. Note that this is intentionally designed
    // such that there's only one place where Plot() is invoked. This will
    // minimize the chances that the compiler fails to inline the functor.
    for(unsigned y = GetXY(*topleft)[1], next_left = y, next_right = y;
        (y < next_left  || Update(sleft,  left,  next_left,  right[0], !clockwise))
    &&  (y < next_right || Update(sright, right, next_right, left[0],   clockwise)); ++y)
    {
        // On a single scanline, we go from the left X coordinate to the right X coordinate.
        // Take the current X coordinates from both sides
        int x = sleft.current[X], endx = sright.current[X];
        // Interpolate between the current values of the two endpoints
        auto x_current   = sleft.current.template Slice<XSlopeSize>();
        auto x_increment = (sright.current - sleft.current).template Slice<XSlopeSize>() * (1.f / (endx-x));
        for(; x < endx; ++x, x_current += x_increment)
        {
            float z = 1.f;
            if constexpr(PerspectiveCorrected) z /= x_current[Z];
            Plot(x,y, z, [&x_current,z](unsigned n) { return x_current[n] * z; });
        }
        // Progress both sides to the next scanline
        sleft.current  += sleft.increment;
        sright.current += sright.increment;
    }
}

// Clip a convex 3D polygon against a plane
template<typename Container, typename VertFun, typename GenFun1, typename GenFun2>
void ClipPolygon(const Plane& p, Container& container, VertFun&& GetVertex,
                 GenFun1&& GenerateIntermediateVertex, GenFun2&& GenerateOtherVertex)
{
    // Against each edge of the polygon
    for(auto pi = container.begin(); pi != container.end(); )
    {
        auto ni = pi;
        if(++ni == container.end()) { ni = container.begin(); }
        float outside = p.DistanceTo(GetVertex(pi)), outsidenext = p.DistanceTo(GetVertex(ni));

        // If this corner is outside the plane, keep it
        bool in  = outside <= 0, out = outside >= 0;

        // If this edge of the polygon _crosses_ the plane, generate an intersection point
        if((outside < 0 && outsidenext > 0) || (outside > 0 && outsidenext < 0))
            GenerateIntermediateVertex(pi, ni, outside / (outside - outsidenext), in,out);
        else
            GenerateOtherVertex(pi, in,out);
    }
}

template<unsigned MaxVertices = 32,
         typename PointMaker, typename Get3D,
         typename ConvertTo2D, typename Get2D, typename ForwardIterator,
         typename PlaneVectorType, // such as const std::vector<Plane>&
         typename Receiver>
void ProjectAndClipPolygon(ForwardIterator&& begin, ForwardIterator&& end,
                           PlaneVectorType&& frustum, PointMaker&& make_point,
                           Get3D&& get3d, ConvertTo2D&& convert2d, Get2D&& get2d,
                           Receiver&& receive, bool Twosided = false)
{
    if(std::distance(begin,end) < 3) return;

    // The corner vertices of this polygon; translated relative
    // to the player's position and rotated to the player's view.
    typedef std::result_of_t<PointMaker(ForwardIterator&&)> PointType;
    std::vector<PointType> points;
    points.reserve(std::distance(begin,end));

    // Collect the corner vertics of the polygon.
    // As soon as we have enough vertices, test whether the polygon is facing the camera.
    // If we render both sides of the polygon, skip the check and keep all points.
    if(Twosided)
        while(begin != end) points.emplace_back( make_point(begin++) );
    else
        for(bool whichface = false; begin != end; )
        {
            points.emplace_back( make_point(begin++) );
            // Discard polygons that are not facing the camera (back-face culling).
            if(points.size() >= 3 && !whichface)
            {
                auto p=points.end(), cur=--p, prev=--p, prev2=--p;
                float side = get3d(prev).CrossProduct(get3d(cur)).DotProduct(get3d(prev2));
                if(side < 0) { return; }         // it's facing the wrong way, so skip rendering.
                if(side > 0) { whichface=true; } // we know which face it is now
            }
        }

    // Clip the polygon against the frustum while it's still 3D.
    for(const Plane& p: frustum) // Each edge of the frustum
    {
        // Against each edge of the polygon. Transform the array of points in place.
        bool keepfirst = true;
        ClipPolygon(p, points, get3d,
                    [&](auto& pi, auto ni, float scale, bool, bool keep)
                    {
                        if(pi==points.begin()) { keepfirst=keep; keep=true; }
                        auto b = *pi + (*ni - *pi) * scale;
                        if(keep) { ++pi; pi = points.insert(pi, std::move(b)); ++pi; }
                        else     { *pi = std::move(b); ++pi; }
                    },
                    [&](auto& pi, bool,bool keep)
                    {
                        if(pi==points.begin()) { keepfirst=keep; keep=true; }
                        if(keep) ++pi; else pi = points.erase(pi);
                    });
        if(!keepfirst) points.erase(points.begin());

        // Can only render surfaces. If we don't have at least three points, it's not a surface.
        if(points.size() < 3) return;
    }

    // Perspective-project (transform into 2D).
    // After this function, get2d(point) can be used to retrieve 2D coordinates.
    for(PointType& p: points) convert2d(p);

    // Eliminate points that are too close to each others.
    // In other words, keep only points that differ from the next point by at least 1 pixel.
    // Eliminate also points that are too close to the line between their adjacent points.
recheck_points:
    if(points.size() < 3) return;
    auto prev = points.end(), now = points.begin(), next = now;--prev;++next;
    for(; now != points.end(); prev = now++)
    {
        auto d = get2d(now) - get2d(prev);
        if(d.DotProduct(d) < 1
        || std::abs(d.CrossProduct(get2d(next) - get2d(prev))) < 1 )
        {
            points.erase( points.begin() + (now-points.begin()) );
            goto recheck_points;
        }
        if(++next == points.end()) next = points.begin();
    }
    receive(points);
}


/* A view module for dealing with 2D-3D geometry conversions */
template<unsigned Width, unsigned Height>
struct View
{
    static constexpr unsigned W = Width, H = Height;
    float zbuffer[W*H];
    xy    fov, center, scale;
    View(float f) : fov(f, f*H/W), center{W*.5f,H*.5f}, scale{center[0] / std::tan((fov[0]/2) * std::atan2(0.f,-1.f)/180),
                                                              center[1] / std::tan((fov[1]/2) * std::atan2(0.f,-1.f)/180)} {}
    // Perspective projection (3D to 2D, and 2D to 3D)
    xy PerspectiveProject(const xyz& p) const            { return center + p.Slice<2>() * scale / p[2];   }
    xyz PerspectiveUnproject(const xy& p, float z) const { return xyz{(p - center) / scale, 1} * z; }

    // Viewport geometry control
    std::array<xy,4> GetEdges() const                    { return { xy{0,0}, xy{W-1,0}, xy{W-1,H-1}, xy{0,H-1} }; }
    std::vector<Plane> GetFrustum() const         // Generates clipping planes for any arbitrary viewport
    {
        constexpr float znear = 16, zcorner = 1;  // Near clipping plane distance, and arbitrary z value.
        std::vector<Plane> result;
        auto e = GetEdges();
        for(unsigned n=0; n<4; ++n)
            result.push_back(Plane::FromPoints( xyz{0,0,0}, PerspectiveUnproject(e[n], zcorner), PerspectiveUnproject(e[(n+1)%4], zcorner) ));
        result.push_back(Plane::FromPoints( xyz{0,0,znear}, xyz{1,0,znear}, xyz{0,1,znear} ));
        return result;
    }

    // Z-buffer control
    void ResetZ()                    { for(auto& z: zbuffer) z = 1e30f; }
    float ReadZ(unsigned c) const    { return zbuffer[c]; }
    void WriteZ(unsigned c, float z) { zbuffer[c] = z; }
    unsigned Coordinate(unsigned x,unsigned y) const     { return y*W + x; }
};


/* Level data */
/* Vertices are 3 letters each (x,y,z). */
static const char vlist[] = "FHHKHHKAHFAHFHBFBBHBBHHBFABHABHEBIEBIHBIABKEBKHBKABFHEHHEHHHFHDHHDIHEIHHKHEIHDKHDKCDKC"
"HKCBKCFKAFKCEKAEKRQMRQMNQKNQKRAKNAMNAMRAKREMREKRDMRDMNGKNGMNFKNFMNEKNEKQEKQDKNDKQCKRCKQAFLHMLHMIHFIHFLBFKBIKBILBFJ"
"BIJBFIBHIBHJBIIBLIBLLBMIBMLBFLEKLEKLHKLGMLGKLFMLFMLEHIHHIEFIEHIDFIDIIHIIELIHLIEMIEMIDIIDKIDKIBLIDGXDIXDIRDGRDMXDGX"
"AGRAMXAGRCJRCJRAGSDGSAGSCFNEFNAFLAMLAFNCKNCJNAJNCILALLAMHHMCHLCDLHDLEDMEDMHDMCDLHELHHMCFLCFLCEFRCJQCFQCFRAJQAIXEMX"
"EIRELEBIEALEAIKAIJAIIAIHALIAFJAFBAHBAHJAHIAHHAHEAFXDFSDFXAFSAFSCLCBLICLHCLECLCCMAHLAFMAFFRDIQDFQDFKAMECMHCMIAIQELAE";
/* Edges are convex polygons where each corner is 4 letters: 2 hex digits (vertex number) followed by U and V coordinates. */
static const std::string edges[] {    "03HA02HF01AF00AA","06AC09BC08BA05AA","0BAB0DEB09EA0AAA","0EAC10EC0DEA0BAA","12AC15BC14BA11AA",
"17AB16DB12DA13AA","16AB0CDB07DA12AA","18AC1ABC19BA16AA","08GA10GF02AF03AA","04AG08HG03HA00AA","21CA10CD1DAD20AA","25EA24EC23AC22AA",
"29AC28EC27EA26AA","23AC2BMC2AMA22AA","2FKA2EKC24AC25AA","33BA32BC30AC31AA","2AAM33EM25EA22AA","35AB36DB33DA34AA","39AD27DD36DA35AA",
"24EA28EQ29AQ23AA","3DDA3CDH3BAH3AAA","40AD43BD42BA3FAA","4EAF4DDF4CDA3AAA","3BAC50BC4FBA4EAA","52AC53BC4DBA51AA","58BA57BC55AC56AA",
"55DA5ADB59AB54AA","5CDA5DDB3CAB5BAA","57BA5EBF5DAF55AA","45CA47CB5FAB57AA","48CA4ACB5EAB62AA","3EAG44DG3DDA3AAA","3CDA4ADG4BAG3BAA",
"66GA65GC64AC63AA","6AAG29GG69GA68AA","67AG6ADG68DA63AA","6FAC69BC6BBA70AA","2DGA29GD6AAD67AA","4CCA53CH32AH71AA","28AH74CH73CA72AA",
"33AF76CF75CA71AA","76AB27CB77CA78AA","72AE73CE4CCA71AA","53CA74CE28AE32AA","1CFA7CFC7BAC01AA","80AB82CB7DCA7FAA","83AB7EBB1ABA18AA",
"7BAB81EB7EEA84AA","1BBA7DBB87AB20AA","7DCA82CB85AB86AA","7CFA82FE81AE7BAA","75DA78DE89AE8AAA","6DAE77EE72EA8BAA","6BAB69CB8BCA88AA",
"8BAC72EC75EA88AA","78DA77DC8CAC89AA","8FGA2BGE8EAE8DAA","8EAE67BE64BA8DAA","64AB65GB8FGA8DAA","2BGA2DGB67AB8EAA","7AAD92HD91HA79AA",
"91BA92BD90AD0BAA","93AB94BB43BA40AA","95AB96BB0CBA47AA","90EA92EB97AB48AA","9BAC9AIC99IA98AA","46AC9BBC98BA42AA","99BA9ABC06AC05AA",
"98AB99IB05IA42AA","07BA9DBB9CAB45AA","06DA9ADB9EAB0AAA","A0FA6EFB63AB9FAA","68AB6FFBA2FAA1AA","63AB68DBA1DA9FAA","A2CA6FCB70ABA3AA",
"A1ADA2FDA0FA9FAA","1ABA7EBB62AB60AA","90ABA4CB1DCA0EAA","1DCAA4CB7DAB1BAA","7EBAA6BBA5AB62AA","7DCAA8CBA7AB7FAA","A8GAA4GB48ABA5AA",
"17BA84BD5BAD59AA","5CAD83BD16BA5AAA","5AAD16BD17BA59AA","84BA83BD5CAD5BAA","02CAA9CC7CAC1CAA","85ABABCBAACA86AA","1FCAABCCA9AC02AA",
"A9CAABCC85AC7CAA","00BA13BC54AC3DAA","55AC12BC11BA56AA","56AD11BD00BA3DAA","13BA12BD55AD54AA","AEBAADBD65ADACAA","8ABA37BF35AFAEAA",
"88AB8ABBAEBAACAA","14BA15BC57AC58AA","44AC04BC14BA58AA","15BA07BC45AC57AA","19BA1ABC60AC5FAA","47AC0CBC19BA5FAA","4FCA50CC2EAC2FAA",
"30AC52CC51CA31AA","31AB51CB4FCA2FAA","50CA52CB30AB2EAA","79AD93BDAFBA73AA","AFBA93BD40AD3FAA","73ABAFBB3FBA3EAA","96AB91DB9EDA9DAA",
"0CAB96BB9DBA07AA","9EBA91BB0BAB0AAA","B1ABB0DBA7DAA6AA","81ABB1BBA6BA7EAA","A7BAB0BB80AB7FAA","80DAB0DBB1AB81AA","74ABB2DB97DA7AAA",
"97BAB2BB4AAB48AA","4ADAB2DB74AB4BAA","B3BA34BC2AAC8FAA","ADBA35BC34ACB3AA","65ABADBBB3BA8FAA","26AB39BB8CBA6DAA","8CCA39CB37AB89AA",
"87ABB4CB21CA20AA","21BAB4BBAAAB1FAA","AACAB4CB87AB86AA","ACBA66BB6EABA0AA","70AB6BBB88BAA3AA","A3AB88BBACBAA0AA","94AB95BB9CBA9BAA",
"43AB94BB9BBA46AA","9CBA95BB47AB45AA" };
/* Note that this program is fully capable of rendering arbitrary 3D geometry. This is just a sample scene. */


/* And finally the main drawing module */
int main()
{
    // Initialize video
    //View<1280,720> view(90.f);
    View<212,120> view(90.f);
    SDL_Surface* surface = SDL_SetVideoMode(view.W, view.H, 32,0);

    // Give each vertex a random RGB color
    std::vector<unsigned> rndcolors(std::string(vlist).size()/3);
    for(auto& r: rndcolors) r = std::rand() & 0xFFFFFF;

    // Generate a single grayscale texture
    std::vector<unsigned char> texture(32*32);
    for(unsigned y=0; y<32; ++y)
        for(unsigned x=0; x<32; ++x)
            texture[y*32+x] = (x<1 || y<1 || (x+1)>=32 || (y+1)>=32) ? 204 : (255-38*std::pow((std::rand()%100)/100.0, 2.0));

    // Initialize the list of held keys and the pre-assigned keyboard mappings
    bool keys[16]{}, quit=false;
    const std::map<unsigned,bool*> known_keys{ {SDLK_ESCAPE,&quit},
        {SDLK_KP8,&keys[0]},{SDLK_UP,  &keys[0]},{SDLK_KP2,&keys[1]},{SDLK_DOWN, &keys[1]},{SDLK_LALT,  &keys[14]},
        {SDLK_KP4,&keys[2]},{SDLK_LEFT,&keys[2]},{SDLK_KP6,&keys[3]},{SDLK_RIGHT,&keys[3]},{SDLK_RALT,  &keys[14]},
        {SDLK_KP7,&keys[4]},{SDLK_KP9, &keys[5]},{SDLK_KP1,&keys[6]},{SDLK_KP3,  &keys[7]},{SDLK_LSHIFT,&keys[15]},
        {SDLK_KP_MINUS,&keys[8]},{SDLK_KP_PLUS,&keys[9]},{'a',&keys[10]},{'z',&keys[11]},{'q',&keys[12]},{'e',&keys[13]} };

    // Initialize camera
    xyz camera{39091.7,44185.6,34086.9}, movement{0,0,0}, rotation{0,0,0};
    quat cameradir{0.129572,0.177557,0.686448,-0.693163};

    // Main loop
    auto begin          = std::chrono::system_clock::now();
    unsigned framecount = 0;
    bool     rendered   = false;
    while(!quit)
    {
        fprintf(stderr, "    xyz camera{%g,%g,%g}, movement{0,0,0}, rotation{0,0,0};\nquat cameradir{%g,%g,%g,%g};\n",
            camera[0],camera[1],camera[2],cameradir[0],cameradir[1],cameradir[2],cameradir[3]);
        // If we are ahead the time budget (10ms per frame (100fps)), spend time doing something useful
        while(std::chrono::duration<double>(std::chrono::system_clock::now() - begin).count() < framecount*0.010)
        {
            if(rendered) { SDL_Delay(10); continue; } // If the frame has already been rendered, sleep instead
            // Clear the screen
            view.ResetZ();
            memset(surface->pixels, 255, view.W*view.H*4); // Optional
            // Draw all polygons, using this plot-pixel function:
            auto plot = [&](unsigned x,unsigned y, float z, auto&& prop)
                        {
                            unsigned coord = view.Coordinate(x,y);
                            if(coord > view.W*view.H || z > view.ReadZ(coord)) return;
                            view.WriteZ(coord, z);
                            auto dprop = [=](unsigned propno) // Round the float prop into integer with Bayer 4x4 ordered-dithering
                            {
                                // Dithering produces visually similar results as bilinear interpolation, but is much faster
                                constexpr unsigned long long bayer4x4 = 0x0819C4D53B2AF7E6; // Array of 16 4-bit values
                                return (unsigned(prop(propno)*16.f) + ((bayer4x4 >> ((y&3)*16 + (x&3)*4)) & 15)) >> 4;
                            };
                            unsigned u = dprop(0), v = dprop(1), r = dprop(2), g = dprop(3), b = dprop(4);
                            unsigned t = texture[(u&31)*32 + (v&31)]; // Sample the texture
                            unsigned pix = (unsigned(r*t/256) << 16) | (unsigned(g*t/256) << 8) | (unsigned(b*t/256) << 0);
                            ((unsigned*)surface->pixels)[coord] = pix;
                        };
            // Project, clip, tesselate and draw each polygon.
            // Technically this could be done in parallel, if Z-buffer accesses were atomic.
            constexpr bool Twosided = false;
            for(const auto& e: edges)
                ProjectAndClipPolygon(e.begin(), e.begin()+e.size()/4, view.GetFrustum(),
                                      [&](auto ei) -> vec<float,8> // Return xyz coordinates, uv coordinates and rgb color
                                      {
                                          std::size_t a = (ei-e.begin())*4;
                                          unsigned vertno = std::stoi(e.substr(a,2),nullptr,16),  c = rndcolors[vertno];
                                          return { cameradir.Rotate(xyz{ vlist[vertno*3+0]*512.f,
                                                                         vlist[vertno*3+1]*512.f,
                                                                         vlist[vertno*3+2]*512.f } - camera),
                                                   (e[a+2]-65)*32.f, (e[a+3]-65)*32.f, (c>>16)&0xFF, (c>>8)&0xFF, c&0xFF };
                                      },
                                      [](auto pi) { return pi->template Slice<3>(); },
                                      [&](auto& p) { p.template PartialAssign<0>(view.PerspectiveProject(p.template Slice<3>())); },
                                      [](auto pi) { return pi->template Slice<2>(); },
                                      [&](auto&& points)
                                      {
                                          RasterizePolygon<Twosided,true,5>(&points[0], points.size(),
                                              [](const auto& p) { return vec<int,2>{ p.template Slice<2>() } ; },
                                              [](const auto& p) { return p[2]; },
                                              [](const auto& p, unsigned prop) { return p[3+prop]; },
                                              plot);
                                      }, Twosided);
            SDL_Flip(surface);
            rendered = true;
        }

        // Wait for events
        for(SDL_Event ev; SDL_PollEvent(&ev); )
        {
            auto i = known_keys.find((ev.type == SDL_KEYDOWN || ev.type == SDL_KEYUP)
                                      ? ev.key.keysym.sym : (ev.type == SDL_QUIT) ? SDLK_ESCAPE : SDLK_UNKNOWN);
            if(i != known_keys.end()) *i->second = (ev.type == SDL_KEYDOWN);
        }

        // Convert keyboard status into relative rotation and movement vectors (un-normalized)
        xyz new_rota_vec{ keys[14] ? 0.f : (keys[0]-keys[1]), // pitch fwd
                          keys[14] ? 0.f : (keys[2]-keys[3]), // turn left
                          float(keys[4]||keys[12]) - (keys[5]||keys[13]) }; // bank left
        xyz new_move_vec{ float(keys[6] || (keys[14] && keys[3])) - (keys[7] || (keys[14] && keys[2])), // left
                          float(keys[9] || (keys[14] && keys[1])) - (keys[8] || (keys[14] && keys[0])), // up
                          float(keys[10] - keys[11]) }; // forward

        // Apply rotation with hysteresis
        rotation = rotation * 0.9f + new_rota_vec * 0.1f;
        if(rotation.Length() > 1e-3f)
            cameradir = cameradir * quat::FromAngle(rotation.Length()*0.02f, (-cameradir).Rotate(rotation.Normalized()));

        // Apply movement with hysteris
        movement = movement * 0.8f + (-cameradir).Rotate(new_move_vec.Normalized() * 32.f) * 0.2f;
        camera   += movement;

        // Update rendering statistics. Schedule frame for rendering if the camera is moving.
        ++framecount;
        if(movement.Length() > 1e-3f || rotation.Length() > 1e-3f) rendered = false;
    }
    SDL_Quit();
}