ExtractIndices.cpp

#include "MarchingCubesTables.h"
#include "IsoSurfaces_internal.h"
 
#include "Context.h"
#include "Platform/Instrumentation.h"
#include "Services/Features.h"
 
#include "Foundation/Platform/Timer.h"
#include "Foundation/Platform/Threads.h"
 
namespace {
 
  struct ExtIdxTask
  {
    Cogs::Core::Context* context;
    const uint8_t* axesTable;
    const uint8_t* indexCountTable;
    const uint8_t* indexTable;
    const int32_t* cellMap;
    const uint8_t* activeCellCases;
    const int32_t* activeCellVertexOffsets;
    const int32_t* activeCellIndexOffsets;
    const int32_t* activeCellIndices;
    const Cogs::Core::IsoSurfaces::Scratch::ijk_t* activeCellIJK;
    glm::ivec3 M;
    glm::ivec3 gridA;
    int ca, cb;
    std::atomic<uint64_t>* elapsed_us;
    Cogs::Core::IsoSurfaces::IndexEmitFunc emit;
 
    void operator()()
    {
      CpuInstrumentationScope(SCOPE_ISOSURFACES, "ExtIdx");
      auto timer = Cogs::Timer::startNew();
 
      uint32_t Mx = M.x;
      uint32_t MxMy = M.x*M.y;
 
#if !defined(EMSCRIPTEN) && !defined(__APPLE__)
      if (context->features->supported(Cogs::Core::CPUFeature::SSE42)) {
        for (int c = ca; c < cb; c++) {
          const auto oc = activeCellIndexOffsets[c];
          const auto on = activeCellIndexOffsets[c + 1];
          if (oc == on) continue;
 
          uint32_t index = oc;
 
          auto ijk = activeCellIJK[c];
 
          int cellIndex = activeCellIndices[c];
 
          const auto code = activeCellCases[c];
          const auto n = static_cast<uint32_t>(indexCountTable[code]);
 
          uint32_t indices[16];
          for (uint32_t k = 0; k < n; k += 3) {
            const auto axisShift0 = indexTable[16 * code + k + 0];
            const auto axisShift1 = indexTable[16 * code + k + 1];
            const auto axisShift2 = indexTable[16 * code + k + 2];
 
            __m128i axisShift = _mm_set_epi32(0, axisShift2, axisShift1, axisShift0);
 
            __m128i t0 = _mm_andnot_si128(_mm_cmpeq_epi32(_mm_and_si128(_mm_set1_epi32(8), axisShift), _mm_setzero_si128()), _mm_set1_epi32(1));
            __m128i t1 = _mm_andnot_si128(_mm_cmpeq_epi32(_mm_and_si128(_mm_set1_epi32(16), axisShift), _mm_setzero_si128()), _mm_set1_epi32(Mx));
            __m128i t2 = _mm_andnot_si128(_mm_cmpeq_epi32(_mm_and_si128(_mm_set1_epi32(32), axisShift), _mm_setzero_si128()), _mm_set1_epi32(MxMy));
 
            __m128i cellOfVertex = _mm_add_epi32(_mm_add_epi32(_mm_set1_epi32(cellIndex), t0),
                                                 _mm_add_epi32(t1, t2));
 
            uint32_t cellOfVertex0 = cellMap[cellOfVertex.m128i_i32[0]];
            uint32_t cellOfVertex1 = cellMap[cellOfVertex.m128i_i32[1]];
            uint32_t cellOfVertex2 = cellMap[cellOfVertex.m128i_i32[2]];
 
            auto codeShifted0 = activeCellCases[cellOfVertex0];
            auto codeShifted1 = activeCellCases[cellOfVertex1];
            auto codeShifted2 = activeCellCases[cellOfVertex2];
            __m128i axes = _mm_and_si128(_mm_set_epi32(0,
                                                       axesTable[codeShifted2],
                                                       axesTable[codeShifted1],
                                                       axesTable[codeShifted0]),
                                         axisShift);
            __m128i t3 = _mm_and_si128(_mm_srli_epi32(axes, 1), _mm_set1_epi32(1));
            __m128i t4 = _mm_and_si128(_mm_srli_epi32(axes, 2), _mm_set1_epi32(1));
 
            auto offsetShifted0 = activeCellVertexOffsets[cellOfVertex0];
            auto offsetShifted1 = activeCellVertexOffsets[cellOfVertex1];
            auto offsetShifted2 = activeCellVertexOffsets[cellOfVertex2];
            __m128i offsetShifted = _mm_add_epi32(_mm_set_epi32(0, offsetShifted2, offsetShifted1, offsetShifted0),
                                                  _mm_add_epi32(t3, t4));
            _mm_storeu_si128(reinterpret_cast<__m128i*>(indices + k), offsetShifted);
          }
          emit(index, glm::ivec3(ijk.i, ijk.j, ijk.k) + gridA, indices, n);
        }
      }
      else {
#else
      {
#endif
        for (int c = ca; c < cb; c++) {
          const auto oc = activeCellIndexOffsets[c];
          const auto on = activeCellIndexOffsets[c + 1];
          if (oc == on) continue;
 
          uint32_t index = oc;
 
          auto ijk = activeCellIJK[c];
 
          int cellIndex = activeCellIndices[c];
 
          //auto i = delinarize(cellIndex, M, gridA);
 
          const auto code = activeCellCases[c];
          const auto n = static_cast<uint32_t>(indexCountTable[code]);
 
          uint32_t indices[16];
          for (uint32_t k = 0; k < n; k++) {
            const auto axisShift = indexTable[16 * code + k];
            assert(axisShift != 255);
 
            uint32_t cellOfVertex = cellMap[cellIndex +
              ((axisShift & 8) ? 1 : 0) +
              ((axisShift & 16) ? Mx : 0) +
              ((axisShift & 32) ? MxMy : 0)];
            assert(0 <= cellOfVertex);
 
            auto codeShifted = activeCellCases[cellOfVertex];
            auto offsetShifted = activeCellVertexOffsets[cellOfVertex];
            auto axes = axesTable[codeShifted] & axisShift;
 
            offsetShifted +=
              ((axes & 4) ? 1 : 0) +
              ((axes & 2) ? 1 : 0);
 
            indices[k] = offsetShifted;
          }
          emit(index, glm::ivec3(ijk.i, ijk.j, ijk.k) + gridA, indices, n);
        }
      }
      if (elapsed_us != nullptr) {
        elapsed_us->fetch_add(timer.elapsedMicroseconds());
      }
    }
 
  };
 
}
 
 
void Cogs::Core::IsoSurfaces::extractIndices(Context* context,
                                             TaskId group,
                                             IndexEmitFunc&   emitFunc,
                                             const int32_t* cellMap,
                                             const uint8_t* activeCellCases,
                                             const int32_t* activeCellVertexOffsets,
                                             const int32_t* activeCellIndexOffsets,
                                             const int32_t* activeCellIndices,
                                             const Scratch::ijk_t* activeCellIJK,
                                             const int32_t    activeCellCount,
                                             const glm::ivec3 gridA,
                                             const glm::ivec3 gridB,
                                             std::atomic<uint64_t>* elapsed_us)
{
  const auto & indexCountTable = MarchingCubes::indexCountTable();
  const auto & indexTable = MarchingCubes::indexTable();
  const auto & axesTable = MarchingCubes::axesTable();
  const glm::ivec3 M = gridB - gridA;
 
  const auto Q = static_cast<int>(Threads::hardwareConcurrency());
  const auto taskSize = std::max(1024, (activeCellCount + Q - 1) / Q);
  for (int c = 0; c < activeCellCount; c += taskSize) {
 
    ExtIdxTask task;
    task.context = context;
    task.axesTable = axesTable.data();
    task.indexCountTable = indexCountTable.data();
    task.indexTable = indexTable.data();
    task.cellMap = cellMap;
    task.activeCellCases = activeCellCases;
    task.activeCellVertexOffsets = activeCellVertexOffsets;
    task.activeCellIndexOffsets = activeCellIndexOffsets;
    task.activeCellIndices = activeCellIndices;
    task.activeCellIJK = activeCellIJK;
    task.M = M;
    task.gridA = gridA;
    task.ca = c;
    task.cb = std::min(c + taskSize, activeCellCount);
    task.elapsed_us = elapsed_us;
    task.emit = emitFunc;
 
    context->taskManager->enqueueChild(group, task);
 
  }
}