UniformGridSystem_isosurf_ref.cpp

#include "Resources/Resources.h"
#include "Resources/MeshManager.h"
#include "Resources/VertexFormats.h"
 
#include "Context.h"
 
#include "../../../IsoSurfaces/MarchingCubesTables.h"
 
#include "Foundation/Logging/Logger.h"
#include "Foundation/Memory/MemoryBuffer.h"
#include "Foundation/Platform/Timer.h"
 
#include <glm/glm.hpp>
 
//#include "C:\utils\iaca-win64\iacaMarks.h"
namespace {
  using namespace Cogs::Core;
 
  Cogs::Logging::Log logger = Cogs::Logging::getLogger("UniformGridSystem");
 
  struct Vertex
  {
    glm::vec3 pos;
    glm::vec3 nrm;
    glm::vec2 tex;
  };
  typedef uint16_t Idx;
 
  void analyzeRef(uint8_t* tmp,
                  unsigned& vertexCount_,
                  unsigned& occupiedCells_,
                  unsigned& indexCount_,
                  const float* values,
                  const glm::uvec3 gridLayout,
                  const unsigned T_n,
                  const float* thresholds)
  {
    const auto cellLayout = gridLayout - glm::uvec3(1, 1, 1);
    const auto * vertexCountTable = MarchingCubes::vertexCountTable().data();
    const auto indexCountTable = MarchingCubes::indexCountTable();
 
    unsigned vertexCount = 0;
    unsigned occupiedCells = 0;
    unsigned indexCount = 0;
    for (unsigned l = 0; l < T_n; l++) {
      auto * offsets = (uint32_t*)(tmp + gridLayout.x*gridLayout.y*gridLayout.z * sizeof(uint32_t)*l);
      auto * codes = tmp + gridLayout.x*gridLayout.y*gridLayout.z*(sizeof(uint32_t)*T_n + l);
 
      const auto T = thresholds[l];
 
      auto * c = codes;
      auto * o = offsets;
      const auto * v = values;
      for (unsigned k = 0; k < gridLayout.z; k++) {
        for (unsigned j = 0; j < gridLayout.y; j++) {
          for (unsigned i = 0; i < gridLayout.x; i++) {
 
            auto ii = i < cellLayout.x ? 1 : 0;
            auto jj = j < cellLayout.y ? gridLayout.x : 0;
            auto kk = k < cellLayout.z ? gridLayout.x * gridLayout.y : 0;
 
            unsigned code =
              (v[0] < T ? 1 : 0) |
              (v[ii] < T ? 2 : 0) |
              (v[jj] < T ? 4 : 0) |
              (v[jj + ii] < T ? 8 : 0) |
              (v[kk + 0] < T ? 16 : 0) |
              (v[kk + ii] < T ? 32 : 0) |
              (v[kk + jj] < T ? 64 : 0) |
              (v[kk + jj + ii] < T ? 128 : 0);
            v++;
 
            if (code == 255) code = 0;
            *o++ = vertexCount;
            if (code) {
              vertexCount += vertexCountTable[code];
              if (i < cellLayout.x && j < cellLayout.y && k < cellLayout.z) {
                occupiedCells++;
                indexCount += indexCountTable[code];
              }
            }
            *c++ = uint8_t(code);
          }
        }
      }
    }
    vertexCount_ = vertexCount;
    occupiedCells_ = occupiedCells;
    indexCount_ = indexCount;
  }
 
  void analyzeSomewhatOptimized(uint8_t* tmp,
                                unsigned& vertexCount_,
                                unsigned& occupiedCells_,
                                unsigned& indexCount_,
                                const float* values,
                                const glm::uvec3 gridLayout,
                                const unsigned T_n,
                                const float* thresholds)
  {
    assert((gridLayout.x & 7) == 0 && (gridLayout.y & 7) == 0 && (gridLayout.z & 7) == 0);
    assert((gridLayout.y*gridLayout.x & 31) == 0);
 
    const auto cellLayout = gridLayout - glm::uvec3(1, 1, 1);
    const auto * vertexCountTable = MarchingCubes::vertexCountTable().data();
    const auto indexCountTable = MarchingCubes::indexCountTable();
 
    unsigned vertexCount = 0;
    unsigned occupiedCells = 0;
    unsigned indexCount = 0;
    for (unsigned l = 0; l < T_n; l++) {
      auto * offsets = (uint32_t*)(tmp + gridLayout.x*gridLayout.y*gridLayout.z * sizeof(uint32_t)*l);
      auto * codes = (tmp + gridLayout.x*gridLayout.y*gridLayout.z*(sizeof(uint32_t)*T_n + l));
      const auto T = thresholds[l];
      for (unsigned k = 0; k < gridLayout.z*gridLayout.y*gridLayout.x; k++) {
        codes[k] = values[k] < T ? 1 : 0;
      }
 
      auto * c = codes;
      const auto m = gridLayout.y*gridLayout.x;
      for (unsigned k = 0; k < cellLayout.z; k++) {
        for (unsigned j = 0; j < gridLayout.y*gridLayout.x; j++) {
          unsigned code = (c[0] << 0) | (c[m] << 4);
          *c++ = uint8_t(code);
        }
      }
      for (unsigned j = 0; j < gridLayout.y* gridLayout.x; j++) {
        unsigned code = (c[0] << 0) | (c[0] << 4);
        *c++ = uint8_t(code);
      }
 
      c = codes;
      for (unsigned k = 0; k < gridLayout.z; k++) {
        for (unsigned i = 0; i < cellLayout.y*gridLayout.x; i++) {
          c[0] = c[0] | (c[gridLayout.x] << 2);
          c++;
        }
        for (unsigned i = 0; i < gridLayout.x; i++) {
          c[0] = c[0] | (c[0] << 2);
          c++;
        }
      }
 
      c = codes;
      auto * o = offsets;
      for (unsigned k = 0; k < gridLayout.z; k++) {
        for (unsigned j = 0; j < gridLayout.y; j++) {
          for (unsigned i = 0; i < cellLayout.x; i++) {
            unsigned code = (c[0] << 0) | (c[1] << 1);
            if (code == 255) code = 0;
            *o++ = vertexCount;
            *c++ = uint8_t(code);
            if (code == 0) continue;
 
            vertexCount += vertexCountTable[code];
            if (i < cellLayout.x && j < cellLayout.y && k < cellLayout.z) {
              occupiedCells++;
              indexCount += indexCountTable[code];
            }
          }
          unsigned code = (c[0] << 0) | (c[0] << 1);
          if (code == 255) code = 0;
          *o++ = vertexCount;
          if (code) {
            vertexCount += vertexCountTable[code];
          }
          *c++ = uint8_t(code);
        }
      }
    }
 
    vertexCount_ = vertexCount;
    occupiedCells_ = occupiedCells;
    indexCount_ = indexCount;
  }
 
  unsigned vertexExtractRef(Cogs::Core::MappedStream<Vertex>& vertices,
                            uint8_t* tmp,
                            const float* values,
                            const glm::uvec3 gridLayout,
                            const unsigned T_n,
                            const float* thresholds)
  {
    unsigned vix = 0;
    const auto cellLayout = gridLayout - glm::uvec3(1, 1, 1);
    const auto * axesTable = MarchingCubes::axesTable().data();
    for (unsigned l = 0; l < T_n; l++) {
      //auto * offsets = (uint32_t*)(tmp + gridLayout.x*gridLayout.y*gridLayout.z * sizeof(uint32_t)*l);
      auto * codes = tmp + gridLayout.x*gridLayout.y*gridLayout.z*(sizeof(uint32_t)*T_n + l);
      auto * c = codes;
      const auto tp = (l + 0.5f) * (1.0f / 4.f);
      const auto T = thresholds[l];
      for (unsigned k = 0; k < gridLayout.z; k++) {
        for (unsigned j = 0; j < gridLayout.y; j++) {
          for (unsigned i = 0; i < gridLayout.x; i++) {
            auto code = *c++;
            if (code == 0) continue;
 
            auto axes = axesTable[code];
            unsigned ll = 0;
            if (axes) {
              auto ip = std::min(cellLayout.x, i + 1);
              auto jp = std::min(cellLayout.y, j + 1);
              auto kp = std::min(cellLayout.z, k + 1);
 
              float a = values[(k*gridLayout.y + j)*gridLayout.x + i];
 
              float ax = values[(k*gridLayout.y + j)*gridLayout.x + ip] - a;  // gradient via forward differences
              float ay = values[(k*gridLayout.y + jp)*gridLayout.x + i] - a;
              float az = values[(kp*gridLayout.y + j)*gridLayout.x + i] - a;
 
              if (axes & 4) {
                ll++;
                auto kk = k + 1;
                auto kkp = std::min(cellLayout.z, kk + 1);
                //assert(kk < gridLayout.z);
 
                float b = values[(kk*gridLayout.y + j)*gridLayout.x + i];
                float bx = values[(kk*gridLayout.y + j)*gridLayout.x + ip] - b;
                float by = values[(kk*gridLayout.y + jp)*gridLayout.x + i] - b;
                float bz = values[(kkp*gridLayout.y + j)*gridLayout.x + i] - b;
 
                float t = (T - a) / (b - a);  // solve linear eq
                //assert(0 <= t && t <= 1.f);
                float s = 1.f - t;
                vertices[vix++] = Vertex{ glm::vec3(i, j, k + t),
                                          -glm::vec3(s*ax + t * bx, s*ay + t * by, s*az + t * bz),
                                          glm::vec2(tp, 0.5) };
              }
 
              if (axes & 2) {
                ll++;
                auto jj = j + 1;
                auto jjp = std::min(cellLayout.y, jj + 1);
                //assert(jj < gridLayout.y);
 
                float b = values[(k*gridLayout.y + jj)*gridLayout.x + i];
                float bx = values[(k*gridLayout.y + jj)*gridLayout.x + ip] - b;
                float by = values[(k*gridLayout.y + jjp)*gridLayout.x + i] - b;
                float bz = values[(kp*gridLayout.y + jj)*gridLayout.x + i] - b;
 
                float t = (T - a) / (b - a);
                //assert(0 <= t && t <= 1.f);
                float s = 1.f - t;
                vertices[vix++] = Vertex{ glm::vec3(i, j + t, k),
                                          -glm::vec3(s*ax + t * bx, s*ay + t * by, s*az + t * bz),
                                          glm::vec2(tp, 0.5) };
              }
              if (axes & 1) {
                ll++;
                auto ii = i + 1;
                auto iip = std::min(cellLayout.y, ii + 1);
                //assert(ii < gridLayout.x);
 
                float b = values[(k*gridLayout.y + j)*gridLayout.x + ii];
                float bx = values[(k*gridLayout.y + j)*gridLayout.x + iip] - b;
                float by = values[(k*gridLayout.y + jp)*gridLayout.x + ii] - b;
                float bz = values[(kp*gridLayout.y + j)*gridLayout.x + ii] - b;
 
                float t = (T - a) / (b - a);
                //assert(0 <= t && t <= 1.f);
                float s = 1.f - t;
                vertices[vix++] = Vertex{ glm::vec3(i, j + t, k),
                                          -glm::vec3(s*ax + t * bx, s*ay + t * by, s*az + t * bz),
                                          glm::vec2(tp, 0.5) };
              }
              //auto q = vertexCountTable[code];
              //assert(l == q);
            }
 
          }
        }
      }
    }
 
    return vix;
  }
 
}
 
namespace Cogs::Core::EchoSounder {
 
 
//#pragma optimize( "", off )
  MeshHandle createIsoSurfacesReference(Context* context,
                                        uint64_t& analyze,
                                        uint64_t& vtx,
                                        uint64_t& idx,
                                        Cogs::Memory::MemoryBuffer& scratch,
                                        const float* values,
                                        const glm::uvec3 gridLayout,
                                        const glm::vec3 /*minCorner*/,
                                        const glm::vec3 /*maxCorner*/,
                                        const float *thresholds, size_t count)
  {
    assert(gridLayout.x != 0 && gridLayout.y != 0 && gridLayout.z != 0);
 
    const uint32_t T_n = (uint32_t)count;
 
    // grid layout is the layout of sample positions
    // a cell must be surrounded by samples, so cell layout is one less along each dimension
    //const auto indexCountTable = MarchingCubes::indexCountTable();
 
    //const auto * vertexCountTable = MarchingCubes::vertexCountTable().data();
    const auto * axesTable = MarchingCubes::axesTable().data();
    const auto * indexTable = MarchingCubes::indexTable().data();
 
    const auto cellLayout = gridLayout - glm::uvec3(1, 1, 1);
 
    scratch.resize(T_n*(sizeof(uint32_t) + sizeof(uint8_t) + 3*sizeof(uint32_t))*gridLayout.x*gridLayout.y*gridLayout.z, false);
    auto * tmp = (uint8_t*)scratch.data();
 
 
    unsigned occupiedCells = 0;
    unsigned indexCount = 0;
    unsigned vertexCount = 0;
 
    auto timer = Timer::startNew();
    //analyzeRef(tmp, vertexCount, occupiedCells, indexCount,  values, gridLayout, T_n, thresholds);
    analyzeSomewhatOptimized(tmp, vertexCount, occupiedCells, indexCount, values, gridLayout, 
                             static_cast<const unsigned>(T_n), thresholds);
    analyze = timer.elapsedMicroseconds();
 
    if (occupiedCells == 0) {
      vtx = 0;
      idx = 0;
      return MeshHandle::NoHandle;
    }
 
    auto mesh = context->meshManager->createLocked();
    //mesh->setBounds(Cogs::Geometry::BoundingBox{ glm::vec3(0,0,0), glm::vec3(gridLayout.x, gridLayout.y, gridLayout.z) });
    mesh->setMeshFlag(MeshFlags::ClockwiseWinding);
    mesh->primitiveType = PrimitiveType::TriangleList;
 
    auto vertices = mesh->map<Vertex>(VertexDataType::Interleaved0, VertexFormats::Pos3fNorm3fTex2f, vertexCount);
 
    // Extract vertices
    timer = Timer::startNew();
    auto vix = vertexExtractRef(vertices, tmp, values, gridLayout,
                                static_cast<const unsigned>(T_n), thresholds);
    vtx = timer.elapsedMicroseconds();
 
    // extract indices
    std::vector<uint32_t> sub_mesh(T_n+1, 0);
    mesh->clearIndexes();
    Idx *indices = (Idx*)mesh->mapStream(VertexDataType::Indexes, 0, indexCount, sizeof(Idx), true);
 
    unsigned iix = 0;
    timer = Timer::startNew();
    for (unsigned l = 0; l < T_n; l++) {
      auto * offsets = (uint32_t*)(tmp + gridLayout.x*gridLayout.y*gridLayout.z * sizeof(uint32_t)*l);
      auto * codes = tmp + gridLayout.x*gridLayout.y*gridLayout.z*(sizeof(uint32_t)*T_n + l);
      auto * c = codes;
      for (unsigned k = 0; k < cellLayout.z; k++) {
        for (unsigned j = 0; j < cellLayout.y; j++) {
          for (unsigned i = 0; i < cellLayout.x; i++) {
            auto code = *c++;
            if (code == 0) continue;
 
            auto * axisShifts = indexTable + 16 * code;
            while(true)
            {
              auto axisShift = *axisShifts++;
              if (axisShift == 255) break;
 
              auto ii = i + (axisShift & 8 ? 1 : 0);
              auto jj = j + (axisShift & 16 ? 1 : 0);
              auto kk = k + (axisShift & 32 ? 1 : 0);
 
              auto shiftedCell = (kk*gridLayout.y + jj)*gridLayout.x + ii;
 
              const auto shiftedCode = codes[shiftedCell];
              const auto axes = axesTable[shiftedCode] & axisShift;
 
              const auto ix = offsets[shiftedCell] + ((axes & 4) ? 1 : 0) + ((axes & 2) ? 1 : 0);
              assert(ix < vertexCount);
 
              assert(sizeof(Idx)!=2 || ix < 0xffff);
              indices[iix++] = (Idx)ix;
            }
          }
 
          c++;
        }
        c += gridLayout.x;
      }
      sub_mesh[l+1] = iix;
    }
    assert(vix == vertexCount);
    assert(iix == indexCount);
 
    mesh->setMeshFlag(MeshFlags::Indexed);
    mesh->setMeshFlag(MeshFlags::IndexesChanged);
 
    auto subMeshes = mesh->mapSubMeshes(T_n);
    for (unsigned l = 0; l < T_n; l++) {
      uint32_t start = sub_mesh[l];
      uint32_t size = sub_mesh[l+1]-sub_mesh[l];
      subMeshes[l] = {start, size, PrimitiveType::TriangleList};
    }
    mesh->setCount(indexCount);
 
    auto box = Geometry::makeEmptyBoundingBox<Geometry::BoundingBox>();
    for(unsigned i=0; i<vertexCount; i++){
      glm::vec3 pos = vertices[i].pos;
      box.min = glm::min(box.min, pos);
      box.max = glm::max(box.max, pos);
    }
    mesh->setBounds(box);
 
    idx = timer.elapsedMicroseconds();
    return mesh.getHandle();
  }
//#pragma optimize( "", on )
 
}