AdaptivePlanarGridSystem.cpp

#include "AdaptivePlanarGridSystem.h"
 
#include "Foundation/Logging/Logger.h"
 
#include "Components/Appearance/MaterialComponent.h"
#include "Components/Core/SubMeshRenderComponent.h"
#include "Components/Core/MeshComponent.h"
#include "Components/Core/TransformComponent.h"
#include "Components/Core/LodComponent.h"
#include "Components/Core/SceneComponent.h"
 
#include "Systems/Core/CameraSystem.h"
#include "Systems/Core/TransformSystem.h"
#include "Systems/Core/SubMeshRenderSystem.h"
 
#include "Resources/Mesh.h"
#include "Resources/MeshManager.h"
#include "Resources/VertexFormats.h"
#include "Resources/MaterialManager.h"
 
#include "Utilities/FrustumClassification.h"
#include "Utilities/Simplex.h"
 
#include "Context.h"
#include "EntityStore.h"
 
#include <cmath>
 
using std::atan2;
using std::round;
using std::sort;
using std::numeric_limits;
using std::vector;
using std::equal;
using std::move;
 
using glm::max;
using glm::min;
using glm::vec2;
using glm::vec3;
using glm::vec4;
using glm::bvec3;
using glm::dvec2;
using glm::dvec3;
using glm::ivec2;
using glm::ivec3;
using glm::mat4;
using glm::dmat2;
 
using glm::abs;
using glm::sin;
using glm::cos;
using glm::mix;
using glm::dot;
using glm::clamp;
using glm::normalize;
using glm::any;
using glm::lessThan;
 
// Node hierarchy
// --------------
//
// The adaptive planar grid uses a small node hierarchy to manage the tiles:
//
// Entity A [AdapativePlanarGridComponent-entity]
//    |
//    +-- Entity B [Group]
//    |      |
//    |      +-- Entity C1 [MeshPart] (Represents a tile)
//    |      |
//    |      +-- Entity C2 [MeshPart] (Represents a tile)
//    |      |
//    |     ...
//    |      |
//    |      +-- Entity Cm [MeshPart] (Represents a tile)
//    |
//    |
//    +-- Entity D [Default] (Holds debug graphics)
//
// Entity A is the entity that contains an AdaptivePlanarGridComponent, and is
// created by the user. First time the system sees it, it creates this
// hierarchy.
//
// Entity B is a group that contains all the tiles.
//
// Entities C1 through Cm are a pool of entities, where C1 through Cn <Cm are
// currently in used, and are set as visible. Each of these represents a
// visible tile in the scene. Entities Cn+1 through Cm are invisible and are
// just kept around until they are required by the scene. The number of tiles
// is limited, and thus, this pool is also limited. Further, there is no
// fixed relationship between a tile entity and a tile at a given position.
// Tile entities are reorganized each time the camera moves.
//
// Transforms
// ----------
//
//  World space -+-> Root space -+-> Cam space -+-> Clip space
//               |               |              |
//               |               |              +--- projection matrix
//               |               |
//               |               +--- Camera world-to-local
//               |
//               +--- entity A local-to-world
//
// Entity A defines its own world space. We have introduced the term root
// space, which is the world space for the node hierarchy. In the case that
// Entity A's transform is the identity, these two spaces coincide. If
// AdaptivePlanarGridComponent::positionable is true, the world space of
// the grid relative to the root may be freely specified. If it is false,
// that transform is hijacked, and is effectively a unity transform, i.e.,
// root and world is the same coordinate systems.
//
// All properties of AdaptiveGridCOmponent (e.g. extent) refers to Entity A's
// world space.
//
// AdaptivePlanarGridData::currentExtent is the intersection between the
// current view frustum and the extent. Thus, it is the subset of the
// extent that is relevant for the current frustum, and forms the basis
// for the refinement hierarchy. When the extent changes, it is made slightly
// larger, to avoid changing the position of the refinement lines all the time.
//
// AdaptivePlanarGridData::currentCoordinates are the coordinates used when
// current extent was defined.
 
namespace
{
  Cogs::Logging::Log logger = Cogs::Logging::getLogger("AdaptivePlanarGridSystem");
 
  using namespace Cogs::Core;
 
  size_t maxTiles = 1000;
  int maxLevel = 20;
  std::vector<uint8_t> simplexScratch;
  std::vector<AdaptivePlanarGridQuadTreeNode> quadTree;
  std::vector<AdaptivePlanarGridIndex> leaves;
  std::list<AdaptivePlanarGridRefinableNode> refineQueue;
 
 
  struct Triangle
  {
    float key;
    ivec3 ix;
  };
 
  struct Edge
  {
    vec3 a;
    vec3 b;
  };
 
  vec3 euclidean(const vec4& h)
  {
    return (1.f / h.w)*vec3(h);
  }
 
  struct CameraState
  {
    glm::mat4 worldToView;
    glm::mat4 worldToClip;
    glm::mat4 viewToWorld;
    glm::vec3 worldUnitXY;
    glm::vec2 viewportSize;
    float nearDistance = 0.0f;
    float farDistance = 0.0f;
  };
 
  // Find aabbox of intersection volume between the extent bounding box and the frustums.
  void analyzeView(Context* context,
                   const AdaptivePlanarGridComponent & gridComp,
                   AdaptivePlanarGridData & gridData,
                   const std::vector<CameraState>& cameras,
                   bool updateDebugGraphics)
  {
    if (!std::isfinite(gridComp.extentMin.x) || !std::isfinite(gridComp.extentMin.y) ||
        !std::isfinite(gridComp.extentMax.x) || !std::isfinite(gridComp.extentMax.y))
    {
      gridData.currentExtentMin = dvec2(0.0);
      gridData.currentExtentMax = dvec2(0.0);
      return;
    }
 
    const auto bboxMin = glm::dvec3(gridComp.extentMin, 0.0) + glm::dvec3(gridComp.displaceMin);
    const auto bboxMax = glm::dvec3(gridComp.extentMax, 0.0) + glm::dvec3(gridComp.displaceMax);
    auto bboxSize = bboxMax - bboxMin;
 
    dvec2 smallestExtentMin(std::numeric_limits<double>::max());
    dvec2 smallestExtentMax(-std::numeric_limits<double>::max());
 
    for (auto & camera : cameras) {
 
      // Find intersection shape of extent bounding box and the frustum of P.
      auto rows = glm::dmat4(glm::transpose(camera.worldToClip));
      auto l = rows[3] + rows[0];
      auto r = rows[3] - rows[0];
      auto b = rows[3] + rows[1];
      auto t = rows[3] - rows[1];
 
      double restrictions[7 * 4] =
      {
        1.0, 0.0, 0.0, bboxSize.x,
        0.0, 1.0, 0.0, bboxSize.y,
        0.0, 0.0, 1.0, bboxSize.z,
        -l.x, -l.y, -l.z, l.w + glm::dot(glm::dvec3(l),bboxMin),
        -r.x, -r.y, -r.z, r.w + glm::dot(glm::dvec3(r),bboxMin),
        -b.x, -b.y, -b.z, b.w + glm::dot(glm::dvec3(b),bboxMin),
        -t.x, -t.y, -t.z, t.w + glm::dot(glm::dvec3(t),bboxMin)
      };
 
      double objectives[4 * 3] = {
        1.0, 0.0, 0.0,
        -1.0, 0.0, 0.0,
        0.0, 1.0, 0.0,
        0.0, -1.0, 0.0
      };
 
      double solutions[4 * 3];
      auto rv = simplexMaximize(context,
                                (uintptr_t)simplexScratch.data(),
                                solutions, 3,
                                objectives, 4,
                                restrictions, 7,
                                150);
 
      if (rv == SimplexResult::NoFeasibleSolutions) {           // No intersection found
      }
      else if (rv != SimplexResult::OptimalSolutionFound) {   // Simplex method failed, just use full extent
        smallestExtentMin = dvec2(bboxMin);
        smallestExtentMax = dvec2(bboxMax);
      }
      else {                                                  // Sane solution.
 
        // Calc aabbox of the extreme points.
        glm::dvec2 minWorld = dvec2(solutions[0 * 3 + 0], solutions[0 * 3 + 1]);
        glm::dvec2 maxWorld = minWorld;
        for (int i = 1; i < 4; i++) {
          auto p = dvec2(solutions[i * 3 + 0], solutions[i * 3 + 1]);
          minWorld = glm::min(minWorld, p);
          maxWorld = glm::max(maxWorld, p);
        }
 
        // Clamp aabbox to aabox of full extent and update smallestExtent.
        smallestExtentMin = glm::min(smallestExtentMin, glm::max(minWorld + glm::dvec2(bboxMin), glm::dvec2(bboxMin)));
        smallestExtentMax = glm::max(smallestExtentMax, glm::min(maxWorld + glm::dvec2(bboxMin), glm::dvec2(bboxMax)));
      }
 
    }
 
    if ((smallestExtentMax.x - smallestExtentMin.x < std::numeric_limits<float>::epsilon())
        || (smallestExtentMax.y - smallestExtentMin.y < std::numeric_limits<float>::epsilon()))
    {
      gridData.currentExtentMin = dvec2(0.0);
      gridData.currentExtentMax = dvec2(0.0);
    }
    else {
 
      if ((gridData.currentExtentMin.x <= smallestExtentMin.x)
          && (gridData.currentExtentMin.y <= smallestExtentMin.y)
          && (smallestExtentMax.x <= gridData.currentExtentMax.x)
          && (smallestExtentMax.y <= gridData.currentExtentMax.y))
      {
        // New extent is inside current extent, check if new extent isn't too small
        dvec2 A = gridData.currentExtentMax - gridData.currentExtentMin;
        dvec2 B = smallestExtentMax - smallestExtentMin;
        if (B.x > 0.5*A.x && B.y > 0.5*A.y) {
          // It is all right, don't change anything
          return;
        }
      }
      // New extent doesn't match current extent particularly well.
      // Use new extent, but add some extra room so we don't change the extent all the time.
      // Still, make sure that we adhere to prescribed extent
      dvec2 delta = smallestExtentMax - smallestExtentMin;
      gridData.currentExtentMin = glm::max(smallestExtentMin - 0.25 * delta, glm::dvec2(bboxMin));
      gridData.currentExtentMax = glm::min(smallestExtentMax + 0.25 * delta, glm::dvec2(bboxMax));
    }
 
    if (updateDebugGraphics) {
      std::vector<glm::vec3> p;
      p.push_back(glm::vec3(smallestExtentMin.x, smallestExtentMin.y, 0.f));
      p.push_back(glm::vec3(smallestExtentMax.x, smallestExtentMin.y, 0.f));
 
      p.push_back(glm::vec3(smallestExtentMax.x, smallestExtentMin.y, 0.f));
      p.push_back(glm::vec3(smallestExtentMax.x, smallestExtentMax.y, 0.f));
 
      p.push_back(glm::vec3(smallestExtentMax.x, smallestExtentMax.y, 0.f));
      p.push_back(glm::vec3(smallestExtentMin.x, smallestExtentMax.y, 0.f));
 
      p.push_back(glm::vec3(smallestExtentMin.x, smallestExtentMax.y, 0.f));
      p.push_back(glm::vec3(smallestExtentMin.x, smallestExtentMin.y, 0.f));
 
      gridData.debugMesh->setPositions(std::move(p));
      gridData.debugMesh->primitiveType = Cogs::PrimitiveType::LineList;
    }
  }
 
 
 
  int isBoxInsideAnyFrustum(const std::vector<CameraState>& cameras, const glm::vec3& boxMin, const glm::vec3& boxMax, const glm::vec2& nearestWorld, const float tileNorm, const double tileScreenCoverageThreshold)
  {
    bool isInsideAny = false;
    for (auto & camera : cameras) {
 
      // Cull tiles that are outside the frustum
      int planes = 0;
 
      for (int i = 0; i < 8; i++) {
        vec4 p((i & 1) ? boxMin.x : boxMax.x,
               (i & 2) ? boxMin.y : boxMax.y,
               (i & 4) ? boxMin.z : boxMax.z,
               1.f);
        vec4 v = camera.worldToView * p;
        vec4 c = camera.worldToClip * p;
        planes |= (v.z < 0.f ? 1 : 0) | (-camera.farDistance <= v.z ? 2 : 0) |
                  (c.x <= c.w ? 4 : 0) | (-c.w <= c.x ? 8 : 0) |
                  (c.y <= c.w ? 16 : 0) | (-c.w <= c.y ? 32 : 0);
      }
 
      if (planes == 63) {
        // Tile is inside frustum, check refinement condition.
 
        vec4 nearestView = camera.worldToView * vec4(nearestWorld, 0.f, 1.f);
        float d = max(camera.nearDistance, -nearestView.z / nearestView.w);
        if(camera.farDistance < d) continue;
 
        const vec3 cameraOriginWorld(camera.viewToWorld[3]);
        const vec3 cameraZAxisWorld(camera.viewToWorld[2]);
 
        vec3 tileReferenceL = cameraOriginWorld - d * normalize(cameraZAxisWorld) - tileNorm * camera.worldUnitXY;
        vec3 tileReferenceR = cameraOriginWorld - d * normalize(cameraZAxisWorld) + tileNorm * camera.worldUnitXY;
        vec2 tileRefNDCL = vec2(euclidean(camera.worldToClip*vec4(tileReferenceL, 1.f)));
        vec2 tileRefNDCR = vec2(euclidean(camera.worldToClip*vec4(tileReferenceR, 1.f)));
 
        vec2 tileRefScreen = camera.viewportSize * (tileRefNDCR - tileRefNDCL);
        double tileScreenCoverage = double(tileRefScreen.x)*double(tileRefScreen.y);
 
        if (tileScreenCoverageThreshold < tileScreenCoverage) return 2;
        isInsideAny = true;
      }
    }
    return isInsideAny ? 1 : 0;
  }
 
  void refineTiles(const AdaptivePlanarGridComponent & gridComp,
                   AdaptivePlanarGridData & gridData,
                   const std::vector<CameraState>& cameras,
                   const glm::mat4 lodRefTransform,
                   const float levelOfDetailTolerance)
  {
    vec3 cameraOriginWorld(lodRefTransform[3]);
 
    int gridResolution = 1 << gridData.tileResolutionLog2;
 
    // levelOfDetailTolerance should be in pixels squared.
    // Multiply by grid resolution to get triangle quad size.
    // Divide by four since we let clip be in [-1,1].
    double tileScreenCoverageThreshold = 0.25 * double(gridResolution * gridResolution * levelOfDetailTolerance);
 
    // Breadth-first subdivision of domain
    leaves.clear();
    quadTree.clear();
    refineQueue.clear();
    refineQueue.push_back(AdaptivePlanarGridRefinableNode{ { 0, 0, 0 }, NoParent, 0 });
 
    while (!refineQueue.empty()) {
      const auto & currentNode = refineQueue.front();
 
      double delta = 1.0 / double(1 << currentNode.ix.level);
 
      vec2 tileMinWorld = mix(gridData.currentExtentMin, gridData.currentExtentMax, delta * dvec2(currentNode.ix.xIndex, currentNode.ix.yIndex));
      vec2 tileMaxWorld = mix(gridData.currentExtentMin, gridData.currentExtentMax, delta * dvec2(currentNode.ix.xIndex + 1, currentNode.ix.yIndex + 1));
 
      float tileNorm = max(tileMaxWorld.x - tileMinWorld.x, tileMaxWorld.y - tileMinWorld.y) / gridResolution;
 
      // Grow bbox to accommodate whatever displacement that the vertex shader might produce.
      vec3 displacedTileMinWorld = vec3(tileMinWorld, 0.f) + gridComp.displaceMin;
      vec3 displacedTileMaxWorld = vec3(tileMaxWorld, 0.f) + gridComp.displaceMax;
 
      // Point in tile nearest lod-reference
      vec2 nearestWorld = clamp(vec2(cameraOriginWorld), tileMinWorld, tileMaxWorld);
 
      auto s = isBoxInsideAnyFrustum(cameras, displacedTileMinWorld, displacedTileMaxWorld, nearestWorld, tileNorm, tileScreenCoverageThreshold);
 
      //if (isBoxInsideAnyFrustum(cameras, displacedTileMinWorld, displacedTileMaxWorld)) {
      if(s != 0) {
 
        // Connect the parent to this node that we will populate in this iteration.
        uint16_t currIndex = static_cast<uint16_t>(quadTree.size());
        if (currentNode.parent != NoParent) {
          quadTree[currentNode.parent].children[currentNode.isWhichChild] = currIndex;
        }
 
 
        bool tileTooLarge = s == 2;
        bool belowMaxTileLimit = leaves.size() + refineQueue.size() < size_t(maxTiles);
        bool belowMaxLevelLimit = currentNode.ix.level + 1 < maxLevel;
 
        if (belowMaxTileLimit && tileTooLarge && belowMaxLevelLimit) {
          quadTree.push_back(
            AdaptivePlanarGridQuadTreeNode{ {
                AdaptivePlanarGridQuadTreeNode::DeadEnd,
                AdaptivePlanarGridQuadTreeNode::DeadEnd,
                AdaptivePlanarGridQuadTreeNode::DeadEnd,
                AdaptivePlanarGridQuadTreeNode::DeadEnd }
            });
 
          const uint32_t x = currentNode.ix.xIndex << 1;
          const uint32_t y = currentNode.ix.yIndex << 1;
          const uint8_t level = currentNode.ix.level + 1;
          const int parent = currIndex;
 
          refineQueue.push_back(AdaptivePlanarGridRefinableNode{ { x, y, level, }, parent, 0 });
          refineQueue.push_back(AdaptivePlanarGridRefinableNode{ { x + 1, y, level, }, parent, 1 });
          refineQueue.push_back(AdaptivePlanarGridRefinableNode{ { x, y + 1, level, }, parent, 2 });
          refineQueue.push_back(AdaptivePlanarGridRefinableNode{ { x + 1, y + 1, level, }, parent, 3 });
        }
        else {
          quadTree.push_back(
            AdaptivePlanarGridQuadTreeNode{ {
                AdaptivePlanarGridQuadTreeNode::Leaf,
                AdaptivePlanarGridQuadTreeNode::Leaf,
                AdaptivePlanarGridQuadTreeNode::Leaf,
                AdaptivePlanarGridQuadTreeNode::Leaf }
            });
 
          leaves.push_back(currentNode.ix);
        }
 
      }
 
      refineQueue.pop_front();
    }
  }
 
 
  void updateTileMaterialInstances(const AdaptivePlanarGridComponent & /*gridComp*/, AdaptivePlanarGridData & gridData)
  {
    auto sceneComp = gridData.tileEntitiesGroup->getComponent<SceneComponent>();
 
    for (auto & c : sceneComp->children) {
      auto mshRndComp = c->getComponent<SubMeshRenderComponent>();
      mshRndComp->subMesh = gridData.camYaw;
    }
 
    auto & matData = gridData.materialData;
 
    for (size_t i = 0; i < gridData.tiles.size(); i++) {
      auto & tile = gridData.tiles[i];
      tile.materialInstance->setVec4Property(matData.lodLevelsKey,
                                             vec4(tile.adjacentLevels[0],
                                                  tile.adjacentLevels[1],
                                                  tile.adjacentLevels[2],
                                                  tile.adjacentLevels[3]));
      tile.materialInstance->setVec4Property(matData.texCoordWorldTransformKey, tile.texCoordWorldTransform);
      tile.materialInstance->setVec4Property(matData.texCoordUnitTransformKey, tile.texCoordUnitTransform);
      tile.materialInstance->setVec4Property(matData.texCoordLinearTransformKey, glm::make_vec4(glm::value_ptr(gridData.texCoordLinearTransform)));
    }
  }
 
  void updatePools(Context * context, Cogs::ComponentModel::Entity * container, AdaptivePlanarGridData &gridData)
  {
    const size_t oldSize = gridData.materialPool.size();
    gridData.materialPool.resize(std::max(gridData.materialPool.size(), gridData.tiles.size()));
    gridData.entityPool.resize(std::max(gridData.entityPool.size(), gridData.tiles.size()));
 
    for (size_t i = oldSize; i < gridData.materialPool.size(); ++i) {
      gridData.materialPool[i] = context->materialInstanceManager->createMaterialInstance(gridData.materialData.material);
      if (gridData.materialData.initMaterialInstanceCallback) {
        gridData.materialData.initMaterialInstanceCallback(static_cast<MaterialInstance*>(gridData.materialPool[i].get()), container, gridData.materialData.initMaterialInstanceData);
      }
    }
 
    for (size_t i = oldSize; i < gridData.entityPool.size(); ++i) {
      gridData.entityPool[i] = context->store->createChildEntity("SubMeshPart", gridData.tileEntitiesGroup.get(), "Pooled Entity " + std::to_string(i));
      auto entity = gridData.entityPool[i];
      auto sceneComp = entity->getComponent<SceneComponent>();
      sceneComp->setChanged();
      auto meshComp = entity->getComponent<MeshComponent>();
      meshComp->meshHandle = gridData.gridMesh;
      meshComp->setChanged();
    }
 
    // Keep inactive entities from showing up in the scene.
    for (size_t i = gridData.tiles.size(); i < gridData.entityPool.size(); ++i) {
      auto & entity = gridData.entityPool[i];
 
      auto renderComponent = entity->getComponent<SubMeshRenderComponent>();
      renderComponent->setVisible(false);
    }
  }
 
  void updateTileEntities(Context* context, const AdaptivePlanarGridComponent & gridComp, AdaptivePlanarGridData & gridData, const glm::mat4 & worldToRoot)
  {
    updatePools(context, gridComp.getContainer(), gridData);
 
    glm::dvec2 fullWidth = gridComp.extentMax - gridComp.extentMin;
 
    glm::vec2 unityMin = glm::vec2((gridData.currentExtentMin - gridComp.extentMin) / fullWidth);
    glm::vec2 unityMax = glm::vec2((gridData.currentExtentMax - gridComp.extentMin) / fullWidth);
    glm::dvec3 origin = context->transformSystem->getOrigin();
 
    bool originChanged = origin != gridData.currentOrigin;
    gridData.currentOrigin = origin;
 
    // Create new child items
    for (size_t i = 0; i < gridData.tiles.size(); i++) {
      auto & tile = gridData.tiles[i];
 
      const float delta = 1.f / (1 << tile.ix.level);
      const vec2 minLocal = delta * vec2(tile.ix.xIndex, tile.ix.yIndex);
      const vec2 maxLocal = minLocal + delta;
 
      auto & entity = gridData.entityPool[i];
      tile.materialInstance = gridData.materialPool[i];
 
      auto meshRendComp = entity->getComponent<SubMeshRenderComponent>();
      meshRendComp->materials.resize(gridData.tessellations);
 
      for (int k = 0; k < gridData.tessellations; k++) {
        meshRendComp->materials[k] = tile.materialInstance;
      }
 
      meshRendComp->setVisible(true);
      meshRendComp->layer = gridComp.layer;
      meshRendComp->objectId = gridData.objectId;
      meshRendComp->subMesh = -1;
      meshRendComp->setChanged();
 
      const vec2 minWorld = mix(vec2(gridData.currentExtentMin), vec2(gridData.currentExtentMax), minLocal);
      const vec2 maxWorld = mix(vec2(gridData.currentExtentMin), vec2(gridData.currentExtentMax), maxLocal);
      const vec2 minUnit = mix(unityMin, unityMax, minLocal);
      const vec2 maxUnit = mix(unityMin, unityMax, maxLocal);
 
      auto transformComponent = entity->getComponent<TransformComponent>();
      vec3 newPosition = vec3(minWorld, 0.0);
      vec3 newScale = vec3(maxWorld - minWorld, 1.f);
 
      if (transformComponent->position != newPosition ||
          transformComponent->scale != newScale ||
          originChanged) {
 
        transformComponent->position = newPosition;
        transformComponent->scale = newScale;
 
        context->transformSystem->updateTransformData(*transformComponent);
      }
 
      vec2 texCoordWorldShiftFromTilePos = gridData.texCoordLinearTransform*minWorld;
      vec2 texCoordUnitShiftFromTilePos = gridData.texCoordLinearTransform*minUnit;
 
      if (std::isfinite(gridData.texCoordPeriod.x)) {
        texCoordWorldShiftFromTilePos.x = glm::mod(texCoordWorldShiftFromTilePos.x, gridData.texCoordPeriod.x);
        texCoordUnitShiftFromTilePos.x = glm::mod(texCoordUnitShiftFromTilePos.x, gridData.texCoordPeriod.x);
      }
      if (std::isfinite(gridData.texCoordPeriod.y)) {
        texCoordWorldShiftFromTilePos.y = glm::mod(texCoordWorldShiftFromTilePos.y, gridData.texCoordPeriod.y);
        texCoordUnitShiftFromTilePos.y = glm::mod(texCoordUnitShiftFromTilePos.y, gridData.texCoordPeriod.y);
      }
 
      tile.texCoordWorldTransform = vec4(maxWorld - minWorld, texCoordWorldShiftFromTilePos);
      tile.texCoordUnitTransform = vec4(maxUnit - minUnit, texCoordUnitShiftFromTilePos);
 
      const auto worldBox0 = glm::vec3(minWorld, 0) + gridComp.displaceMin;
      const auto worldBox1 = glm::vec3(maxWorld, 0) + gridComp.displaceMax;
      vec3 rootBox0(std::numeric_limits<float>::max());
      vec3 rootBox1(-std::numeric_limits<float>::max());
      for (int ii = 0; ii < 8; ii++) {
        glm::bvec3 end((ii & 1) != 0,
                       -(ii & 2) != 0,
                       -(ii & 4) != 0);
        const auto ch = worldToRoot * vec4(mix(worldBox0, worldBox1, end), 1.f);
        const auto c = (1.f / ch.w)*vec3(ch);
        rootBox0 = min(rootBox0, c);
        rootBox1 = max(rootBox1, c);
      }
      auto renderComponent = entity->getComponent<SubMeshRenderComponent>();
      Cogs::Geometry::BoundingBox box{ rootBox0, rootBox1 };
      context->subMeshRenderSystem->setBounds(renderComponent, box, context->transformSystem->getLocalToWorld(transformComponent));
    }
  }
 
  void setupGridMesh(Context* context, const int tileResolutionLog2, AdaptivePlanarGridData & gridData)
  {
    if (gridData.tileResolutionLog2 == tileResolutionLog2) return;
 
    gridData.gridMesh = context->meshManager->create();
    gridData.tileResolutionLog2 = tileResolutionLog2;
    const int gridResolution = 1 << tileResolutionLog2;
 
    vector<Triangle> T(2 * gridResolution * gridResolution);
    ivec2 indices[2][2][3] = {
      { { ivec2(0, 0), ivec2(1, 0), ivec2(0, 1) }, { ivec2(1, 0), ivec2(1, 1), ivec2(0, 1) } },
      { { ivec2(0, 0), ivec2(1, 0), ivec2(1, 1) }, { ivec2(0, 0), ivec2(1, 1), ivec2(0, 1) } }
    };
    gridData.gridMesh->clearIndexes();
 
    vector<uint32_t> I(2 * 3 * gridResolution * gridResolution);
    for (int d = 0; d < gridData.tessellations; d++) {
      const float theta = (2.0f * glm::pi<float>() * d) / float(gridData.tessellations);
      vec2 dir(cos(theta), sin(theta));
      int e = abs(dot(vec2(1, 1), dir)) < abs(dot(vec2(-1, 1), dir)) ? 1 : 0;
      for (int j = 0; j < gridResolution; j++) {
        for (int i = 0; i < gridResolution; i++) {
          for (int t = 0; t < 2; t++) {
            vec2 key(0.f, 0.f);
            for (int k = 0; k < 3; k++) {
              ivec2 p = ivec2(i, j) + indices[e][t][k];
              key += vec2(p);
              T[2 * (j*gridResolution + i) + t].ix[k] = p.y*(gridResolution + 1) + p.x;
            }
            T[2 * (j*gridResolution + i) + t].key = dot(dir, key);
          }
        }
      }
      sort(T.begin(), T.end(), [](const Triangle& a, const Triangle& b) ->bool { return a.key > b.key; });
      for (size_t i = 0; i < T.size(); i++) {
        I[3 * i + 0] = T[i].ix[0];
        I[3 * i + 1] = T[i].ix[1];
        I[3 * i + 2] = T[i].ix[2];
      }
      gridData.gridMesh->addSubMesh(std::span(I), Cogs::PrimitiveType::TriangleList);
    }
 
    vec2 s(1.f / gridResolution, 1.f / gridResolution);
    auto vertices = gridData.gridMesh->mapPositions(0, (gridResolution + 1)*(gridResolution + 1));
    for (int j = 0; j <= gridResolution; j++) {
      for (int i = 0; i <= gridResolution; i++) {
        float edge = 0;
        if (j == 0) { edge = 3 + (i < gridResolution / 2 ? 0.f : 0.5f); }
        else if (j == gridResolution) { edge = 4 + (i < gridResolution / 2 ? 0.f : 0.5f); }
        else if (i == 0) { edge = 1 + (j < gridResolution / 2 ? 0.f : 0.5f); }
        else if (i == gridResolution) { edge = 2 + (j < gridResolution / 2 ? 0.f : 0.5f); }
        vertices[(gridResolution + 1)*j + i].x = s.x*i;
        vertices[(gridResolution + 1)*j + i].y = s.y*j;
        vertices[(gridResolution + 1)*j + i].z = edge;
      }
    }
 
    gridData.gridMesh->setBounds(Cogs::Geometry::BoundingBox{ { 0.f, 0.f, 0.f },{ 1.0f, 1.0f, 1.0f } });
  }
 
  void setupProxyMesh(Context* context, const AdaptivePlanarGridComponent & gridComp, AdaptivePlanarGridData& gridData)
  {
    // Create dummy geometry with insufficient vertices for a single triangle,
    // so it will produce no rendering output. But since we set the bounding
    // box, it will push the near and far planes so the child tiles won't be
    // culled unnecessary.
    gridData.proxyMesh = context->meshManager->create();
    vector<uint32_t> I(1);
    gridData.proxyMesh->clearIndexes();
    gridData.proxyMesh->addSubMesh(std::span(I), Cogs::PrimitiveType::TriangleList);
    gridData.proxyMesh->setBounds(Cogs::Geometry::BoundingBox{ { 0.f, 0.f, 0.f },{ 1.f, 1.f, 1.f } });
 
    auto meshComponent = gridComp.getComponent<MeshComponent>();
    if (meshComponent) {
      meshComponent->meshHandle = gridData.proxyMesh;
      meshComponent->setChanged();
    }
  }
 
  void setupDebugMesh(Context * context, const AdaptivePlanarGridComponent & gridComp, AdaptivePlanarGridData & gridData)
  {
    gridData.debugGraphics = context->store->createChildEntity("Default", gridComp.getContainer(), "Debug Graphics");
    gridData.debugMesh = context->meshManager->create();
    auto matComp = gridData.debugGraphics->getComponent<MaterialComponent>();
    matComp->diffuseColor = vec4(1.f, 0.f, 0.f, 1.f);
    matComp->enableLighting = false;
  }
 
  void updateTileTopology(AdaptivePlanarGridData & gridData)
  {
    gridData.tiles.clear();
    gridData.tiles.reserve(leaves.size());
 
    for (size_t i = 0; i < leaves.size(); i++) {
      const AdaptivePlanarGridIndex & leaf = leaves[i];
 
      // Determine refinement level of neighbors
      int levelSize = 1 << leaf.level;
      int edgeLevels[4];
 
      int shifts[4][2] = {
        { -1, 0 },
        {  1, 0 },
        {  0, -1 },
        {  0,  1 }
      };
 
      for (int e = 0; e < 4; e++) {
        int xIndex = int(leaf.xIndex) + shifts[e][0];
        int yIndex = int(leaf.yIndex) + shifts[e][1];
        int edgeLevel = leaf.level;
 
        if ((0 <= xIndex) && (xIndex < levelSize) && (0 <= yIndex) && (yIndex < levelSize)) {
 
          int currLevel = 0;
          int currQuadNode = 0;
          while (currLevel < maxLevel) {
 
            if (currQuadNode == AdaptivePlanarGridQuadTreeNode::DeadEnd) {
              // child has been culled, we can use whatever ref level we want
              break;
            }
            else if (currQuadNode == AdaptivePlanarGridQuadTreeNode::Leaf) {
              // we prematurely hit a leaf, restrict refinement level to this node
              edgeLevel = currLevel - 1;
              break;
            }
 
            // We have a valid node, we may descend
            if (currLevel == leaf.level) {
              // we have found the node, and if it is a leaf, we have its
              // refinement level. Otherwise, if it is an internal node,
              // it has more refined children that will match this level
              // at the boundary.
              //assert(quadTree[currQuadNode].debugIndex.xIndex == xIndex);
              //assert(quadTree[currQuadNode].debugIndex.yIndex == yIndex);
              break;
            }
 
            // which bit to inspect to determine which child to descend into
            int currBit = leaf.level - currLevel - 1;
            int whichChild = 2 * ((yIndex >> currBit) & 1) + ((xIndex >> currBit) & 1);
            currQuadNode = quadTree[currQuadNode].children[whichChild];
            currLevel++;
          }
        }
 
        edgeLevels[e] = max(0, 5 - max(0, leaf.level - edgeLevel));
      }
 
      gridData.tiles.push_back(AdaptivePlanarGridTile{ leaf,{
        uint8_t(edgeLevels[0]),
        uint8_t(edgeLevels[1]),
        uint8_t(edgeLevels[2]),
        uint8_t(edgeLevels[3]) },
        vec4(0.f),
        vec4(0.f),
        MaterialInstanceHandle::NoHandle });
    }
  }
 
 
}
 
Cogs::Core::AdaptivePlanarGridSystem::AdaptivePlanarGridSystem(Memory::Allocator * allocator, SizeType /*capacity*/) :
  ComponentSystemWithDataPool(allocator, 16)
{
  simplexScratch.resize(simplexScratchSize(3, 4, 8));
}
 
 
 
 
 
 
 
 
void Cogs::Core::AdaptivePlanarGridSystem::registerMaterial(const AdaptivePlanarGridComponent * gridComp, MaterialHandle material, AdaptivePlanarGridMaterialData::InitMaterialInstanceCallback * initMaterialInstanceCallback, void * initMaterialInstanceData)
{
  auto & gridData = getData(gridComp);
 
  auto & materialData = gridData.materialData;
 
  materialData.material = material;
  materialData.lodLevelsKey = material->getVec4Key("lodLevels");
  materialData.texCoordWorldTransformKey = material->getVec4Key("texCoordWorldTransform");
  materialData.texCoordUnitTransformKey = material->getVec4Key("texCoordUnitTransform");
  materialData.texCoordLinearTransformKey = material->getVec4Key("texCoordLinearTransform");
  materialData.initMaterialInstanceCallback = initMaterialInstanceCallback;
  materialData.initMaterialInstanceData = initMaterialInstanceData;
 
  auto * container = gridComp->getContainer();
  for (auto & item : gridData.materialPool) {
    item = context->materialInstanceManager->createMaterialInstance(gridData.materialData.material);
    if (initMaterialInstanceCallback) {
      initMaterialInstanceCallback(static_cast<MaterialInstance*>(item.get()), container, initMaterialInstanceData);
    }
  }
}
 
 
void Cogs::Core::AdaptivePlanarGridSystem::setTexCoordTransform(const AdaptivePlanarGridComponent * gridComp, const glm::mat2 transform, const glm::vec2 period)
{
  auto & gridData = getData(gridComp);
  gridData.texCoordLinearTransform = transform;
  gridData.texCoordPeriod = period;
}
 
 
void Cogs::Core::AdaptivePlanarGridSystem::update(Context * context)
{
  for (const auto & gridComp : pool) {
    const TransformComponent * lodRefTrComp = nullptr;
    if (auto l = gridComp.lodReference.lock(); l) {
      lodRefTrComp = l->getComponent<TransformComponent>();
    }
    if (!lodRefTrComp) {
      lodRefTrComp = context->cameraSystem->getMainCamera()->getComponent<TransformComponent>();
    }
    if(!lodRefTrComp) continue;
 
    std::vector<CameraComponent*> cameraComps;
    for (auto & weak : gridComp.cameras) {
      if (auto entity = weak.lock(); entity) {
        if (auto * comp = entity->getComponent<CameraComponent>(); comp) {
          if (((comp->flags & CameraFlags::EnableRender) != 0) && ((comp->layerMask & gridComp.layer) != 0)) {
            cameraComps.push_back(comp);
          }
        }
      }
    }
 
    if (cameraComps.empty()) {
      if (auto * comp = lodRefTrComp->getComponent<CameraComponent>(); comp) {
        if (((comp->flags & CameraFlags::EnableRender) != 0) && ((comp->layerMask & gridComp.layer) != 0)) {
          cameraComps.push_back(comp);
        }
      }
      else {
        auto * mcomp = context->cameraSystem->getMainCamera();
        if (((mcomp->flags & CameraFlags::EnableRender) != 0) && ((mcomp->layerMask & gridComp.layer) != 0)) {
          cameraComps.push_back(mcomp);
        }
      }
    }
 
    auto & gridData = getData(&gridComp);
    if (!gridData.initialized) {
      setupGridMesh(context, gridComp.tileResolutionLog2, gridData);
      setupProxyMesh(context, gridComp, gridData);
      setupDebugMesh(context, gridComp, gridData);
      gridData.tileEntitiesGroup = context->store->createChildEntity("Group", gridComp.getContainer(), "Tiles");
 
      auto renderComponent = gridComp.getContainer()->getComponent<RenderComponent>();
 
      if (renderComponent) {
        gridData.objectId = renderComponent->objectId;
      }
 
      gridData.initialized = true;
    }
 
    if (gridComp.hasChanged()) {
      setupGridMesh(context, gridComp.tileResolutionLog2, gridData);
    }
 
    auto trComp = gridComp.getComponent<TransformComponent>();
 
    context->transformSystem->updateTransformData(*trComp);
    const glm::mat4 worldToRoot = context->transformSystem->getLocalToWorld(trComp);
    const glm::mat4 rootToWorld = glm::inverse(worldToRoot);
    //const glm::mat4 worldToCam = cameraData.viewMatrix * worldToRoot;
    //const glm::mat4 camToWorld = rootToWorld * cameraData.inverseViewMatrix;
    //const glm::mat4 clipToWorld = rootToWorld * cameraData.inverseViewProjectionMatrix;
    //const glm::mat4 worldToClip = cameraData.viewProjection * worldToRoot;
 
    std::vector<CameraState> cameras;
    for (size_t i = 0; i < cameraComps.size(); i++) {
      auto * camComp = cameraComps[i];
      auto & camData = context->cameraSystem->getData(camComp);
      if (camData.viewportSize.x < 1.f || camData.viewportSize.y < 1.f) continue;
      if ((camComp->flags & CameraFlags::EnableRender) == 0) continue;
 
      // We run before the transform & camera systems (since we want to set position
      // on items that should be processed by these), and thus, the camera matrices
      // haven't been set up yet. We calc fresh view and use the projection from
      // last frame.
      auto * camTrComp = camComp->getComponent<TransformComponent>();
      context->transformSystem->updateTransformData(*camTrComp);
 
      const glm::mat4& inverseViewMatrix = context->transformSystem->getLocalToWorld(camTrComp);
      const glm::mat4 viewMatrix = glm::inverse(inverseViewMatrix);
 
 
      cameras.emplace_back(CameraState());
      auto & c = cameras.back();
 
 
      c.worldToView = viewMatrix * worldToRoot;
      c.worldToClip = camData.projectionMatrix * c.worldToView;
      c.viewToWorld = rootToWorld * inverseViewMatrix;
      c.worldUnitXY = normalize(vec3(c.viewToWorld[0])) + normalize(vec3(c.viewToWorld[1]));
      c.viewportSize = camData.viewportSize;
      c.nearDistance = camComp->nearDistance;
      c.farDistance = camComp->enableClippingPlaneAdjustment ? std::numeric_limits<float>::max() : camComp->farDistance;
    }
 
    if (gridComp.updateLoD) {
      if (!gridData.materialData.material) continue;
 
      if(!cameras.empty()) {
        auto lodComponent = gridComp.getComponent<LodComponent>();
 
        // Determine 2D orientation of camera z-axis, do determine patch draw order
        context->transformSystem->updateTransformData(*lodRefTrComp);
        const mat4 & viewToWorld = context->transformSystem->getLocalToWorld(lodRefTrComp);// cameraData.inverseViewMatrix;
        if (std::isfinite(viewToWorld[2].x) && std::isfinite(viewToWorld[2].y)) {
          const float theta = (float(gridData.tessellations) / (2.0f * glm::pi<float>())) * std::atan2(viewToWorld[2].y, viewToWorld[2].x);
          gridData.camYaw = max(0, min(gridData.tessellations - 1, int(round(theta + 0.5f * gridData.tessellations))));
        }
        else {
          gridData.camYaw = 0;
        }
 
        analyzeView(context, gridComp, gridData, cameras, false);
        refineTiles(gridComp, gridData, cameras, viewToWorld, lodComponent->geometricTolerance);
        updateTileTopology(gridData);
        updateTileEntities(context, gridComp, gridData, worldToRoot);
      }
 
      if (gridComp.debugGeometry) {
        // Clear any lingering debug geometry when LOD is enabled.
        auto meshComp = gridData.debugGraphics->getComponent<MeshComponent>();
        if (meshComp->meshHandle) {
          meshComp->meshHandle = MeshHandle::NoHandle;
        }
      }
    }
    else if (gridComp.debugGeometry) {
      // Create debug geometry, if requested, and if we currently don't update the
      // LoD, and it hasn't been done yet (meshHandle is null)
      auto meshComp = gridData.debugGraphics->getComponent<MeshComponent>();
 
      if (!meshComp->meshHandle) {
        meshComp->meshHandle = gridData.debugMesh;
 
        analyzeView(context, gridComp, gridData, cameras, true);
      }
    }
 
    updateTileMaterialInstances(gridComp, gridData);
 
    //const dvec2 extentMin = gridComp.extentMin - gridData.currentCoordinates + dvec2(gridComp.displaceMin);
    //const dvec2 extentMax = gridComp.extentMax - gridData.currentCoordinates + dvec2(gridComp.displaceMax);
 
    // Update bounding box. XY-bounds are defined via the transformation.
    // TODO: handle position.
    //gridData.proxyMesh->setBounds(Geometry::BoundingBox{ { extentMin.x, extentMin.y, gridComp.displaceMin.z }, { extentMax.x, extentMax.y, gridComp.displaceMax.z } });
  }
}