BasicAPI.cpp

#include <algorithm>
#include "../Extensions/GeometryProcessing/GeometryProcessing.h"
#include "Resources/MeshManager.h"
#include "Context.h"
#include "IsoSurfaces.h"
 
using namespace Cogs::Core;
 
namespace {
 
  struct VertexEmitStore
  {
    glm::vec3* P = nullptr;
    glm::vec3* N = nullptr;
    glm::vec2* T = nullptr;
    const float* F = nullptr;
    const float* DTA = nullptr;
    glm::ivec3 maxFieldIndex;
    glm::uvec3 fieldSize;
    glm::vec3 scale;
    glm::vec3 offset;
    float threshold = 0.f;
    float normalizedLayer = 0.f;
    unsigned vertexOffset = 0;
  };
 
  struct VertexEmitNoData {
    VertexEmitStore s;
    VertexEmitNoData(VertexEmitStore& s) : s(s) {}
    VertexEmitNoData(const VertexEmitNoData& o) : s(o.s) {}
 
    void operator()(int32_t ixIndex, glm::ivec4* samples, const uint32_t ixN)
    {
      auto Fx = s.fieldSize.x;
      auto Fy = s.fieldSize.y;
 
      auto I0 = glm::clamp(glm::ivec3(samples[0][0],
                                      samples[1][0],
                                      samples[2][0]), glm::ivec3(0), s.maxFieldIndex);
      auto value0 = s.F[((I0.z*Fy) + I0.y)*Fx + I0.x] - s.threshold;
      auto pos0 = s.scale*glm::vec3(I0) + s.offset;
 
      for (uint32_t k = 0; k < ixN; k++) {
        auto I1 = glm::clamp(glm::ivec3(samples[0][1 + k],
                                        samples[1][1 + k],
                                        samples[2][1 + k]), glm::ivec3(0), s.maxFieldIndex);
        auto value1 = s.F[((I1.z*Fy) + I1.y)*Fx + I1.x] - s.threshold;
        auto pos1 = s.scale*glm::vec3(I1) + s.offset;
 
        float den = (value1 - value0);
        float t = abs(den) < std::numeric_limits<float>::min() ? 0.5f : static_cast<float>(-value0 / den);
 
        int vo = s.vertexOffset + ixIndex + k;
        s.P[vo] = (1.f - t)*pos0 + t*pos1;
        s.T[vo] = glm::vec2(s.normalizedLayer, 0.5f);
      }
    }
  };
 
  struct VertexEmitWithData {
    VertexEmitStore s;
    VertexEmitWithData(VertexEmitStore& s) : s(s) {}
    VertexEmitWithData(const VertexEmitWithData& o) : s(o.s) {}
 
    void operator()(int32_t ixIndex, glm::ivec4* samples, const uint32_t ixN)
    {
      auto Fx = s.fieldSize.x;
      auto Fy = s.fieldSize.y;
 
      auto I0 = glm::clamp(glm::ivec3(samples[0][0],
                                      samples[1][0],
                                      samples[2][0]), glm::ivec3(0), s.maxFieldIndex);
      auto value0_ = s.F[((I0.z*Fy) + I0.y)*Fx + I0.x];
      auto value0 = value0_ - s.threshold;
      auto pos0 = s.scale*glm::vec3(I0) + s.offset;
      auto data0 = s.DTA[((I0.z*Fy) + I0.y)*Fx + I0.x];
 
      for (uint32_t k = 0; k < ixN; k++) {
        auto I1 = glm::clamp(glm::ivec3(samples[0][1 + k],
                                        samples[1][1 + k],
                                        samples[2][1 + k]), glm::ivec3(0), s.maxFieldIndex);
        auto value1_ = s.F[((I1.z*Fy) + I1.y)*Fx + I1.x];
        auto value1 = value1_ - s.threshold;
        auto pos1 = s.scale*glm::vec3(I1) + s.offset;
        auto data1 = s.DTA[((I1.z*Fy) + I1.y)*Fx + I1.x];
 
        float den = value1 - value0;
        float t = abs(den) < std::numeric_limits<float>::min() ? 0.5f : static_cast<float>(-value0 / den);
 
        float dta = ((1.f - t)*data0 + t*data1) / s.threshold;
 
        int vo = s.vertexOffset + ixIndex + k;
        s.P[vo] = (1.f - t)*pos0 + t*pos1;
        s.T[vo] = glm::vec2(s.normalizedLayer, dta);
      }
    }
  };
 
  struct IndexEmitStore
  {
    uint32_t* indices;
    uint32_t vertexOffset;
  };
 
  struct IndexEmit
  {
    IndexEmitStore s;
    IndexEmit(IndexEmitStore& s) : s(s) {}
    IndexEmit(const IndexEmit& o) : s(o.s) {}
 
    void operator()(uint32_t ixIndex, glm::ivec3 /*cell*/, const uint32_t* cellIndices, const uint32_t ixN)
    {
      for (uint32_t k = 0; k < ixN; k += 3) {
        s.indices[ixIndex + k + 0] = cellIndices[k + 0] + s.vertexOffset;
        s.indices[ixIndex + k + 1] = cellIndices[k + 1] + s.vertexOffset;
        s.indices[ixIndex + k + 2] = cellIndices[k + 2] + s.vertexOffset;
      }
    }
  };
 
  struct IndexEmitFlip
  {
    IndexEmitStore s;
    IndexEmitFlip(IndexEmitStore& s) : s(s) {}
    IndexEmitFlip(const IndexEmitFlip& o) : s(o.s) {}
 
    void operator()(uint32_t ixIndex, glm::ivec3 /*cell*/, const uint32_t* cellIndices, const uint32_t ixN)
    {
      for (uint32_t k = 0; k < ixN; k += 3) {
        s.indices[ixIndex + k + 0] = cellIndices[k + 0] + s.vertexOffset;
        s.indices[ixIndex + k + 1] = cellIndices[k + 2] + s.vertexOffset;
        s.indices[ixIndex + k + 2] = cellIndices[k + 1] + s.vertexOffset;
      }
    }
  };
 
}
 
 
MeshHandle IsoSurfaces::createIsoSurfacesMesh(Context* context,
                                              IsoSurfaces::Scratch& scratch,
                                              const std::vector<float>& thresholds,
                                              const float* densityField,
                                              const float* dataField,
                                              const glm::uvec3 fieldSize,
                                              const glm::vec3 minCorner,
                                              const glm::vec3 maxCorner,
                                              const bool /*exteriorIsInside*/,
                                              const bool flipOrientation,
                                              const bool closeZNeg,
                                              const bool closeZPos,
                                              const bool closeYNeg,
                                              const bool closeYPos,
                                              const bool closeXNeg,
                                              const bool closeXPos,
                                              MeshHandle mesh)
{
  const glm::ivec3 gridA(closeXNeg ? -1 : 0,
                         closeYNeg ? -1 : 0,
                         closeZNeg ? -1 : 0);
  const glm::ivec3 gridB(fieldSize.x + (closeXPos ? 1 : 0),
                         fieldSize.y + (closeYPos ? 1 : 0),
                         fieldSize.z + (closeZPos ? 1 : 0));
  const auto M = gridB - gridA;
 
  const uint32_t layers = static_cast<uint32_t>(std::min(size_t(16),thresholds.size()));
 
  std::vector<int32_t> vertexOffsets;
  std::vector<int32_t> indexOffsets;
  std::vector<int32_t> cellOffsets;
 
  analyze(context,
          vertexOffsets, 
          indexOffsets,
          cellOffsets,
          scratch.cellMap,
          scratch.activeCellCases,
          scratch.activeCellVertexOffsets,
          scratch.activeCellIndexOffsets,
          scratch.activeCellIndices,
          scratch.activeCellIJK,
          thresholds, true, densityField, glm::ivec3(fieldSize), gridA, gridB, nullptr);
 
  const auto vertexCount = vertexOffsets.back();
  const auto indexCount = indexOffsets.back();
 
  if (0 < vertexCount) {
    if (!HandleIsValid(mesh)) {
      mesh = context->meshManager->create();
    }
 
    auto proxy = context->meshManager->lock(mesh);
 
    auto P = proxy->mapPositions(0, vertexCount);
    auto N = proxy->map<glm::vec3>(VertexDataType::Normals, VertexFormats::Norm3f, 0, vertexCount);
    auto T = proxy->map<glm::vec2>(VertexDataType::TexCoords0, VertexFormats::Tex2f, 0, vertexCount);
    scratch.indices.resize(indexCount);
 
    auto extractGroup = context->taskManager->createGroup();
    for (size_t i = 0; i < layers; i++) {
 
      VertexEmitFunc vertexEmit;
 
      VertexEmitStore vtxStore;
      vtxStore.P = P.data;
      vtxStore.N = N.data;
      vtxStore.T = T.data;
      vtxStore.F = densityField;
      vtxStore.DTA = dataField;
      vtxStore.maxFieldIndex = glm::max(glm::uvec3(1), fieldSize) - glm::uvec3(1);
      vtxStore.fieldSize = fieldSize;
      vtxStore.scale = (maxCorner - minCorner) / glm::vec3(fieldSize - glm::uvec3(1));
      vtxStore.offset = minCorner;
      vtxStore.threshold = thresholds[i];
      vtxStore.normalizedLayer = (i + 0.5f) / layers;
      vtxStore.vertexOffset = vertexOffsets[i];
      if (dataField == nullptr) {
        vertexEmit = VertexEmitNoData(vtxStore);
      }
      else {
        vertexEmit = VertexEmitWithData(vtxStore);
      }
 
      extractVertices(context,
                      extractGroup,
                      vertexEmit,
                      scratch.activeCellCases.data() + cellOffsets[i],
                      scratch.activeCellVertexOffsets.data() + cellOffsets[i] + i,
                      scratch.activeCellIndices.data() + cellOffsets[i],
                      scratch.activeCellIJK.data() + cellOffsets[i],
                      cellOffsets[i + 1] - cellOffsets[i],
                      gridA, gridB);
 
      IndexEmitFunc indexEmit;
 
      IndexEmitStore idxStore;
      idxStore.indices = scratch.indices.data() + indexOffsets[i];
      idxStore.vertexOffset = vertexOffsets[i];
      if (flipOrientation == false) {
        indexEmit = IndexEmit(idxStore);
      }
      else {
        indexEmit = IndexEmitFlip(idxStore);
      }
 
      extractIndices(context,
                     extractGroup,
                     indexEmit,
                     scratch.cellMap.data() + i * M.x*M.y*M.z,
                     scratch.activeCellCases.data() + cellOffsets[i],
                     scratch.activeCellVertexOffsets.data() + cellOffsets[i] + i,
                     scratch.activeCellIndexOffsets.data() + cellOffsets[i] + i,
                     scratch.activeCellIndices.data() + cellOffsets[i],
                     scratch.activeCellIJK.data() + cellOffsets[i],
                     cellOffsets[i + 1] - cellOffsets[i],
                     gridA, gridB);
 
    }
    context->taskManager->wait(extractGroup);
    context->taskManager->destroy(extractGroup);
 
    TaskId normalVecGroup = context->taskManager->createGroup();
    context->taskManager->enqueueChild(normalVecGroup, [&]
    {
      auto * indicesPtr = scratch.indices.data();
          GeometryProcessing::normalsFromIndexedTriangles(context,
                                                      reinterpret_cast<float*>(N.data), 3,
                                                      reinterpret_cast<const float*>(P.data), 3,
                                                      vertexOffsets.back(),
                                                      indicesPtr,
                                                      indexOffsets.back(),
                                                      10000);
        });
    context->taskManager->wait(normalVecGroup);
    context->taskManager->destroy(normalVecGroup);
 
    auto bbox = Cogs::Geometry::makeEmptyBoundingBox<Cogs::Geometry::BoundingBox>();
    bbox.min = minCorner;
    bbox.max = maxCorner;
    proxy->setBounds(bbox);
    proxy->setIndexData(scratch.indices.data(), indexOffsets.back());
    proxy->mapSubMeshes(layers);
    for (size_t l = 0; l < layers; l++) {
      auto & subMesh = proxy->getSubMeshes()[l];
      subMesh.startIndex = indexOffsets[l];
      subMesh.numIndexes = indexOffsets[l + 1] - indexOffsets[l];
      subMesh.primitiveType = Cogs::PrimitiveType::TriangleList;
    }
    proxy->setMeshFlag(MeshFlags::ClockwiseWinding);
  }
 
  return mesh;
}