Analyze.cpp

 
#include "Services/Features.h"
#include "Platform/Instrumentation.h"
#include "Services/TaskManager.h"
#include "Services/Variables.h"
#include "Context.h"
 
#include "MarchingCubesTables.h"
#include "IsoSurfaces_internal.h"
 
#include "Foundation/BitTwiddling/PowerOfTwo.h"
#include "Foundation/Platform/Threads.h"
 
 
namespace {
  using namespace Cogs::Core;
 
  inline void buildOffsetTable(std::vector<int32_t>& table)
  {
    int32_t sum = 0;
    for (size_t i = 0; i < table.size(); i++) {
      int32_t v = table[i];
      table[i] = sum;
      sum += v;
    }
  }
 
 
  template<typename T>
  struct AnalyzeTileTask
  {
    Cogs::Core::IsoSurfaces::AnalyzeGlobalState * g = nullptr;
    glm::ivec3 id;
    void operator()();
  };
 
  template<> void AnalyzeTileTask<float>::operator()()
  {
#if !defined(EMSCRIPTEN) && !defined(__APPLE__)
    if (g->tryToUseAVX2 && g->context->features->supported(Cogs::Core::CPUFeature::AVX2)) {
      analyzeTile_f32_AVX2(g, id);
      return;
    }
    if (g->context->features->supported(Cogs::Core::CPUFeature::SSE2)) {
      analyzeTile_f32_SSE(g, id);
      return;
    }
#endif
    analyzeTile_f32(g, id);
  }
 
  template<> void AnalyzeTileTask<uint16_t>::operator()()
  {
#if !defined(EMSCRIPTEN) && !defined(__APPLE__)
    if (g->tryToUseAVX2 && g->context->features->supported(Cogs::Core::CPUFeature::AVX2)) {
      analyzeTile_u16_AVX2(g, id);
      return;
    }
    if (g->context->features->supported(Cogs::Core::CPUFeature::SSE2)) {
      analyzeTile_u16_SSE(g, id);
      return;
    }
#endif
    analyzeTile_u16(g, id);
  }
 
  template<typename Type>
  void analyzeGeneric(Context* context,
                      std::vector<int32_t>& vertexOffsets,  // [L]
                      std::vector<int32_t>& indexOffsets,   // [L]
                      std::vector<int32_t>& cellOffsets,    // [L]
                      Cogs::Memory::TypedBuffer<int32_t>& cellMap,
                      Cogs::Memory::TypedBuffer<uint8_t>& activeCellCases,
                      Cogs::Memory::TypedBuffer<int32_t>& activeCellVertexOffsets, // One element per active cell.
                      Cogs::Memory::TypedBuffer<int32_t>& activeCellIndexOffsets, // One element per active cell.
                      Cogs::Memory::TypedBuffer<int32_t>& activeCellIndices,            // One element per active cell.
                      Cogs::Memory::TypedBuffer<Cogs::Core::IsoSurfaces::Scratch::ijk_t>& activeCellIJK,
                      const std::vector<Type>& thresholds,
                      const bool exteriorIsLess,
                      const Type* field,
                      const glm::ivec3 fieldDim,
                      const glm::ivec3 gridA,
                      const glm::ivec3 gridB,
                      std::atomic<uint64_t>* elapsed_us)
  {
    CpuInstrumentationScope(SCOPE_ISOSURFACES, "IsoSurface::analyze");
 
    const auto * axesTable = MarchingCubes::axesTable().data();
    const auto * indexCountTable = MarchingCubes::indexCountTable().data();
    const size_t Nt = thresholds.size();
    const auto M = gridB - gridA;
 
    const size_t layerStride = (size_t)M.x * (size_t)M.y * (size_t)M.z;
 
    cellOffsets.clear();
    cellOffsets.resize((size_t)Nt + 1);
 
    Cogs::Memory::TypedBuffer<uint8_t> activeCellCasesTmp(Nt*layerStride);
    Cogs::Memory::TypedBuffer<int32_t> activeCellIndicesTmp(Nt*layerStride);
 
    cellMap.resize(Nt*layerStride);
    //auto * cellMapPtr = cellMap.data();
 
    auto * actCellCasesTmp = activeCellCasesTmp.data();
    auto * actCellIndicesTmp = activeCellIndicesTmp.data();
 
    auto Q = std::max(4u, static_cast<unsigned>(std::ceil(std::cbrt(8 * Cogs::Threads::hardwareConcurrency()))));
 
    glm::ivec3 taskSizeClassify;
    taskSizeClassify.x = Cogs::roundUpToPowerOfTwo((static_cast<unsigned>(M.x) + Q - 1) / Q);
    taskSizeClassify.y = Cogs::roundUpToPowerOfTwo((static_cast<unsigned>(M.y) + Q - 1) / Q);
    taskSizeClassify.z = Cogs::roundUpToPowerOfTwo((static_cast<unsigned>(M.z) + Q - 1) / Q);
 
    assert((0 < taskSizeClassify.x) && (0 < taskSizeClassify.y) && (0 < taskSizeClassify.z));
 
    const glm::ivec3 tiles((M.x + taskSizeClassify.x - 1) / taskSizeClassify.x,
                           (M.y + taskSizeClassify.y - 1) / taskSizeClassify.y,
                           (M.z + taskSizeClassify.z - 1) / taskSizeClassify.z);
 
    std::vector<std::atomic<int>> cellOffsetsAtomic(Nt);
    glm::ivec3 tile;
 
    Cogs::Core::IsoSurfaces::AnalyzeGlobalState aTileG;
    aTileG.context = context;
    aTileG.tileSize = taskSizeClassify;
    aTileG.tiles = tiles;
    aTileG.gridA = gridA;
    aTileG.fieldDim = fieldDim;
    aTileG.M = M;
    aTileG.field = field;
    aTileG.thresholds = thresholds.data();
    aTileG.Nt = (unsigned)Nt;
    aTileG.exteriorIsLess = exteriorIsLess;
    aTileG.cellOffsets = cellOffsetsAtomic.data();
    aTileG.cellMap = cellMap.data();
    aTileG.activeCellCases = activeCellCasesTmp.data();
    aTileG.activeCellIndices = activeCellIndicesTmp.data();
    aTileG.elapsed_us = elapsed_us;
    aTileG.tryToUseAVX2 = context->variables->get("Volumetric.IsoSurfaces.AVX2", false);
 
 
    auto analyzeGroup = context->taskManager->createGroup();
    for (tile.z = 0; tile.z < tiles.z; tile.z++) {
      for (tile.y = 0; tile.y < tiles.y; tile.y++) {
        for (tile.x = 0; tile.x < tiles.x; tile.x++) {
          AnalyzeTileTask<Type> task;
          task.g = &aTileG;
          task.id = tile;
          context->taskManager->enqueueChild(analyzeGroup, task);
        }
      }
    }
    context->taskManager->wait(analyzeGroup);
    context->taskManager->destroy(analyzeGroup);
    for (auto i : aTileG.scratchBuffers) {
      delete i;
    }
    aTileG.scratchBuffers.clear();
 
    for (size_t i = 0; i < Nt; i++) {
      cellOffsets[i] = cellOffsetsAtomic[i];
    }
 
    buildOffsetTable(cellOffsets);
 
    activeCellCases.resize(cellOffsets.back());
    activeCellVertexOffsets.resize((size_t)cellOffsets.back() + Nt);
    activeCellIndexOffsets.resize((size_t)cellOffsets.back() + Nt);
    activeCellIndices.resize(cellOffsets.back());
    activeCellIJK.resize(cellOffsets.back());
 
    auto eightConcurrency = std::max((size_t)1, (8 * static_cast<int>(Cogs::Threads::hardwareConcurrency()) + Nt - 1) / Nt);
    auto prefixSumGroup = context->taskManager->createGroup();
    std::vector<IsoSurfaces::AnalyzePopulateCounts::Global> analyzePopCntsG(Nt);
    for (size_t t = 0; t < Nt; t++) {
 
 
      const auto Nc = (size_t)(cellOffsets[t + 1] - cellOffsets[t]);
 
      analyzePopCntsG[t].context = context;
      analyzePopCntsG[t].axesTable = axesTable;
      analyzePopCntsG[t].indexCountTable = indexCountTable;
      analyzePopCntsG[t].actCellCasesIn = actCellCasesTmp + layerStride * t;
      analyzePopCntsG[t].actCellIndicesIn = actCellIndicesTmp + layerStride * t;
      analyzePopCntsG[t].M = M;
      analyzePopCntsG[t].gridA = gridA;
      analyzePopCntsG[t].gridB = gridB;
      analyzePopCntsG[t].actCellCasesOut = activeCellCases.data() + cellOffsets[t];
      analyzePopCntsG[t].actCellIndicesOut = activeCellIndices.data() + cellOffsets[t];
      analyzePopCntsG[t].actCellIJKOut = activeCellIJK.data() + cellOffsets[t];
      analyzePopCntsG[t].actCellVtxCntOut = activeCellVertexOffsets.data() + cellOffsets[t] + t;
      analyzePopCntsG[t].actCellIdxCntOut = activeCellIndexOffsets.data() + cellOffsets[t] + t;
      analyzePopCntsG[t].elapsed_us = elapsed_us;
 
      auto popCntGroup = context->taskManager->createGroup();
 
      auto taskSizePopCnt = std::max((size_t)1024, (Nc + eightConcurrency - 1) / eightConcurrency);
      for (size_t c = 0; c < Nc; c += taskSizePopCnt) {
        IsoSurfaces::AnalyzePopulateCounts popCntTask;
        popCntTask.g = &analyzePopCntsG[t];
        popCntTask.ca = c;
        popCntTask.cb = std::min(c + taskSizePopCnt, Nc);
        context->taskManager->enqueueChild(popCntGroup, popCntTask);
      }
 
      IsoSurfaces::AnalyzePrefixSum prefixSumTask;
      prefixSumTask.antecedent = popCntGroup;
      prefixSumTask.context = context;
      prefixSumTask.actCellVtxOff = activeCellVertexOffsets.data() + cellOffsets[t] + t;
      prefixSumTask.actCellIdxOff = activeCellIndexOffsets.data() + cellOffsets[t] + t;
      prefixSumTask.Nc = Nc;
      prefixSumTask.elapsed_us = elapsed_us;
      context->taskManager->enqueueChild(prefixSumGroup, prefixSumTask);
    }
    context->taskManager->wait(prefixSumGroup);
    context->taskManager->destroy(prefixSumGroup);
 
    vertexOffsets.resize(Nt + 1);
    indexOffsets.resize(Nt + 1);
    for (size_t t = 0; t < Nt; t++) {
      vertexOffsets[t] = activeCellVertexOffsets.data()[cellOffsets[t + 1] + t];
      indexOffsets[t] = activeCellIndexOffsets.data()[cellOffsets[t + 1] + t];
    }
    buildOffsetTable(vertexOffsets);
    buildOffsetTable(indexOffsets);
  }
}
 
Cogs::Memory::MemoryBuffer* IsoSurfaces::AnalyzeGlobalState::scratchAcquire(uint32_t byteCount)
{
  Cogs::Memory::MemoryBuffer* scratch_ = nullptr;
  {
    std::lock_guard<std::mutex> guard(scratchLock);
    if (!scratchBuffers.empty()) {
      scratch_ = scratchBuffers.back();
      scratchBuffers.pop_back();
    }
  }
  if (!scratch_) {
    const glm::ivec3 scratchSize = tileSize + glm::ivec3(1);
    scratch_ = new Cogs::Memory::MemoryBuffer(byteCount);
  }
  else {
    assert(scratch_->size() == byteCount);
  }
  return scratch_;
}
 
void IsoSurfaces::AnalyzeGlobalState::scratchRelease(Cogs::Memory::MemoryBuffer* scratch_)
{
  std::lock_guard<std::mutex> guard(scratchLock);
  scratchBuffers.push_back(scratch_);
}
 
 
 
void Cogs::Core::IsoSurfaces::analyze(Context* context,
                                      std::vector<int32_t>& vertexOffsets,  // [L]
                                      std::vector<int32_t>& indexOffsets,   // [L]
                                      std::vector<int32_t>& cellOffsets,    // [L]
                                      Cogs::Memory::TypedBuffer<int32_t>& cellMap,
                                      Cogs::Memory::TypedBuffer<uint8_t>& activeCellCases,
                                      Cogs::Memory::TypedBuffer<int32_t>& activeCellVertexOffsets, // One element per active cell.
                                      Cogs::Memory::TypedBuffer<int32_t>& activeCellIndexOffsets, // One element per active cell.
                                      Cogs::Memory::TypedBuffer<int32_t>& activeCellIndices,            // One element per active cell.
                                      Cogs::Memory::TypedBuffer<Scratch::ijk_t>& activeCellIJK,
                                      const std::vector<float>& thresholds,
                                      const bool exteriorIsLess,
                                      const float* field,
                                      const glm::ivec3 fieldDim,
                                      const glm::ivec3 gridA,
                                      const glm::ivec3 gridB,
                                      std::atomic<uint64_t>* elapsed_us)
{
  analyzeGeneric(context,
                 vertexOffsets, indexOffsets, cellOffsets, cellMap,
                 activeCellCases, activeCellVertexOffsets, activeCellIndexOffsets, activeCellIndices, activeCellIJK,
                 thresholds, exteriorIsLess,
                 field, fieldDim, gridA, gridB,
                 elapsed_us);
}