quantized_mesh.h

// Copyright (c) 2018 Coronis Computing S.L. (Spain)
// All rights reserved.
//
// This file is part of EMODnet Quantized Mesh Generator for Cesium.
//
// This program 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.
//
// This program 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 for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// Author: Ricard Campos (ricardcd@gmail.com)

#ifndef EMODNET_QMGC_QUANTIZED_MESH_H
#define EMODNET_QMGC_QUANTIZED_MESH_H


#include <string>
#include <vector>
#include "cgal/cgal_utils.h"
#include <ctb.hpp>
#include "misc_utils.h"

class QuantizedMesh {

public:
    // --- Various structs as described in the documentation ---
    // (https://cesiumjs.org/data-and-assets/terrain/formats/quantized-mesh-1.0.html)
    struct Header
    {
        // The center of the tile in Earth-centered Earth-Fixed coordinates
        double CenterX; 
        double CenterY; 
        double CenterZ; 

        // The minimum and maximum heights in the area covered by this tile.
        // The minimum may be lower and the maximum may be higher than
        // the height of any vertex in this tile in the case that the min/max vertex
        // was removed during mesh simplification, but these are the appropriate
        // values to use for analysis or visualization.
        float MinimumHeight; 
        float MaximumHeight; 

        // The tile’s bounding sphere.  The X,Y,Z coordinates are again expressed
        // in Earth-centered Fixed coordinates, and the radius is in meters.
        double BoundingSphereCenterX; 
        double BoundingSphereCenterY; 
        double BoundingSphereCenterZ; 
        double BoundingSphereRadius;  

        // The horizon occlusion point, expressed in the ellipsoid-scaled Earth-centered Fixed frame.
        // If this point is below the horizon, the entire tile is below the horizon.
        // See http://cesiumjs.org/2013/04/25/Horizon-culling/ for more information.
        double HorizonOcclusionPointX; 
        double HorizonOcclusionPointY; 
        double HorizonOcclusionPointZ; 
    };

    struct VertexData
    {
        unsigned int vertexCount = 0; 
        std::vector<unsigned short> u; 
        std::vector<unsigned short> v; 
        std::vector<unsigned short> height; 
    };

    struct IndexData
    {
        unsigned int triangleCount ; 
        std::vector< unsigned int > indices ; 
    };

    struct EdgeIndices
    {
        unsigned int westVertexCount; 
        std::vector<unsigned int> westIndices; 

        unsigned int southVertexCount; 
        std::vector<unsigned int> southIndices; 

        unsigned int eastVertexCount; 
        std::vector<unsigned int> eastIndices; 

        unsigned int northVertexCount; 
        std::vector<unsigned int> northIndices; 
    };

    struct VertexNormals
    {
        std::vector<float> nx; 
        std::vector<float> ny; 
        std::vector<float> nz; 
    };

    struct WaterMask {
        std::vector<unsigned char> mask ; 
    };

    enum ExtensionIds {OCT_VERTEX_NORMALS = 1, WATER_MASK = 2};

    // --- Functions ---
    QuantizedMesh() {
        m_header = Header();
        m_vertexData = VertexData();
        m_indexData = IndexData();
        m_edgeIndices = EdgeIndices();
        m_vertexNormals = VertexNormals();
        m_waterMask = WaterMask();
        m_bytesPerIndex = 0;
    }

    QuantizedMesh(const std::string &filePath) { readFile(filePath); }

    bool readFile(const std::string &filePath);

    bool writeFile(const std::string &filePath);

    void print();

    void printHeader();

    void setHeader(const Header& header) { m_header = header; }

    void setVertexData(const VertexData& vd) { m_vertexData = vd; }

    void setIndexData(const IndexData& id) { m_indexData = id; }

    void setEdgeIndices(const EdgeIndices& ei) { m_edgeIndices = ei; }

    void setVertexNormals(const VertexNormals& vn) { m_vertexNormals = vn; }

    Header getHeader() const { return m_header; }

    VertexData getVertexData() const { return m_vertexData; }

    IndexData getIndexData() const { return m_indexData; }

    VertexNormals getVertexNormals() const { return m_vertexNormals; }

    EdgeIndices getEdgeIndices() const { return m_edgeIndices; }

    // --- Constants ---
    const unsigned short int TILE_SIZE = 65;
    static const unsigned short MAX_VERTEX_DATA = 32767;

    // --- Convenience static functions ---
    static unsigned short remapToVertexDataValue(const double& value, const double& minOr, const double& maxOr) {
        double valueRemap = remap(value, minOr, maxOr, 0, MAX_VERTEX_DATA) ;
        return static_cast<unsigned short>(valueRemap);
    }

    static double remapFromVertexDataValue(const double& value, const double& minOr, const double& maxOr) {
        return remap(value, 0, MAX_VERTEX_DATA, minOr, maxOr ) ;
    }



private:

    // --- Attributes ---
    Header m_header ;
    VertexData m_vertexData ;
    IndexData m_indexData ;
    EdgeIndices m_edgeIndices ;
    int m_bytesPerIndex ; // Number of bytes used per index
    VertexNormals m_vertexNormals ;
    WaterMask m_waterMask ;

    // --- Functions ---

    // Decode a zig-zag encoded value
    unsigned short zigZagDecode( const unsigned short &value ) {
        return (value >> 1) ^ (-(value & 1)) ;
    }

    unsigned short zigZagEncode( const short &value ) {
        return (value << 1) ^ (value >> 31) ;
    }

    void octEncode( const float xyz[3], unsigned char octxy[2] ) {
        const float invL1Norm = (1.0f) / (fabs(xyz[0]) + fabs(xyz[1]) + fabs(xyz[2]));

        if (xyz[2] < 0.0f) {
            octxy[0] = static_cast<unsigned char>((1.0f - float(fabs(xyz[1] * invL1Norm))) * sign(xyz[0]));
            octxy[1] = static_cast<unsigned char>((1.0f - float(fabs(xyz[0] * invL1Norm))) * sign(xyz[1]));
        } else {
            octxy[0] = static_cast<unsigned char>(xyz[0] * invL1Norm);
            octxy[1] = static_cast<unsigned char>(xyz[1] * invL1Norm);
        }
    }

    void octDecode( const unsigned char octxy[2], float xyz[3] ) {

        xyz[0] = uchar2float(octxy[0]);
        xyz[1] = uchar2float(octxy[1]);
        xyz[2] = 1.0f - (fabs(xyz[0]) + fabs(xyz[1]));

        if (xyz[2] < 0.0f) {
            float oldX = xyz[0];
            xyz[0] = ((1.0f) - fabs(xyz[1])) * sign(oldX);
            xyz[1] = ((1.0f) - fabs(oldX))   * sign(xyz[1]);
        }
    }

    inline float sign(float v) {
        return (v < 0.0f) ? -1.0f : 1.0f;
    }

    float uchar2float( unsigned char u ) {
        return float(clamp(int(8) * (1.0f / float((uint64_t(1) << 7) - 1)), -1.0f, 1.0f));
    }

    unsigned char float2uchar( float f ) {
        return (unsigned char)round(clamp(f, -1.0f, 1.0f) * ((uint64_t(1) << 7) - 1));
    }

    inline float clamp(float val, float low, float hi) {
        if (val <= low) {
            return low;
        } else if (val >= hi) {
            return hi;
        } else {
            return val;
        }
    }

    // Round a value
    inline float round(float f) {
        return floor(f + 0.5f);
    }

};


#endif //EMODNET_QMGC_QUANTIZED_MESH_H