Counting.cpp

#if !defined(EMSCRIPTEN)
#include "IsoSurfaces_internal.h"
 
#include "Platform/Instrumentation.h"
#include "Services/Features.h"
#include "Services/TaskManager.h"
#include "Context.h"
 
#include "Foundation/Platform/Timer.h"
 
using namespace Cogs::Core;
 
void IsoSurfaces::AnalyzePrefixSum::operator()()
{
  context->taskManager->wait(antecedent);
  context->taskManager->destroy(antecedent);
  CpuInstrumentationScope(SCOPE_ISOSURFACES, "AnalyzePrefixSum");
  auto timer = Timer::startNew();
 
  actCellVtxOff[Nc] = 0;
  actCellIdxOff[Nc] = 0;
  int32_t vOff = 0;
  int32_t iOff = 0;
  for (int c = 0; c <= Nc; c++) {
    auto tv = actCellVtxOff[c];
    actCellVtxOff[c] = vOff;
    vOff += tv;
 
    auto ti = actCellIdxOff[c];
    actCellIdxOff[c] = iOff;
    iOff += ti;
  }
 
  if (elapsed_us != nullptr) {
    elapsed_us->fetch_add(timer.elapsedMicroseconds());
  }
}
 
void IsoSurfaces::AnalyzePopulateCounts::operator()()
{
  CpuInstrumentationScope(SCOPE_ISOSURFACES, "AnalyzePopCnt");
  auto timer = Timer::startNew();
 
  size_t c = ca;
  const auto Mx = static_cast<uint32_t>(g->M.x);
  const auto Mx_My = Mx*static_cast<uint32_t>(g->M.y);
  const auto Ux = g->gridB.x - g->gridA.x - 1;
  const auto Uy = g->gridB.y - g->gridA.y - 1;
  const auto Uz = g->gridB.z - g->gridA.z - 1;
 
#if !defined(__APPLE__)
  if (g->context->features->supported(CPUFeature::SSE42)) {
    __m128i ones4 = _mm_set1_epi32(1);
    for (; c + 4 < cb; c += 4) {
 
      __m128i r = _mm_loadu_si128(reinterpret_cast<const __m128i*>(g->actCellIndicesIn + c));
      _mm_storeu_si128(reinterpret_cast<__m128i*>(g->actCellIndicesOut + c), r);
 
      __m128i q = _mm_setzero_si128();
      q.m128i_u32[0] = *reinterpret_cast<const uint32_t*>(g->actCellCasesIn + c);
      *reinterpret_cast<uint32_t*>(g->actCellCasesOut + c) = q.m128i_u32[0];
 
      __m128i i_z;
      __m128i i_y;
      __m128i i_x;
      __m128i axes;
      __m128i indexCount;
      for (int i = 0; i < 4; i++) {
        
 
        auto rr = r.m128i_u32[i];
        i_z.m128i_u32[i] = rr / Mx_My;
        i_y.m128i_u32[i] = (rr % Mx_My) / Mx;
        i_x.m128i_u32[i] = ((rr % Mx_My) % Mx);
 
        Scratch::ijk_t tmp;
        tmp.i = (uint16_t)i_x.m128i_u32[i];
        tmp.j = (uint16_t)i_y.m128i_u32[i];
        tmp.k = (uint16_t)i_z.m128i_u32[i];
 
        //assert(tmp.i == i_x.m128i_u32[i]);
        //assert(tmp.j == i_y.m128i_u32[i]);
        //assert(tmp.k == i_z.m128i_u32[i]);
 
        tmp.code = q.m128i_u8[i];
        tmp.axes = g->axesTable[q.m128i_u8[i]];
        g->actCellIJKOut[c + i] = tmp;
 
        axes.m128i_u32[i] = g->axesTable[q.m128i_u8[i]];
        indexCount.m128i_u32[i] = g->indexCountTable[q.m128i_u8[i]];
      }
 
      __m128i axes_0 = _mm_and_si128(ones4, axes);
      __m128i axes_1 = _mm_and_si128(ones4, _mm_srli_epi32(axes, 1));
      __m128i axes_2 = _mm_and_si128(ones4, _mm_srli_epi32(axes, 2));
      __m128i Nv = _mm_add_epi32(axes_2, _mm_add_epi32(axes_1, axes_0));
      _mm_storeu_si128(reinterpret_cast<__m128i*>(g->actCellVtxCntOut + c), Nv);
 
      __m128i indexOut_x = _mm_cmplt_epi32(i_x, _mm_set1_epi32(Ux));
      __m128i indexOut_y = _mm_cmplt_epi32(i_y, _mm_set1_epi32(Uy));
      __m128i indexOut_z = _mm_cmplt_epi32(i_z, _mm_set1_epi32(Uz));
      __m128i indexOut = _mm_and_si128(_mm_and_si128(indexOut_x, indexOut_y), indexOut_z);
      __m128i Ni = _mm_and_si128(indexOut, indexCount);
      _mm_storeu_si128(reinterpret_cast<__m128i*>(g->actCellIdxCntOut + c), Ni);
    }
  }
#endif
 
  for (; c < cb; c++) {
    const auto code = g->actCellCasesIn[c];
    const auto uncompactedCellIndex = g->actCellIndicesIn[c];
 
    glm::ivec3 i;
    i.z = uncompactedCellIndex / Mx_My;
    i.y = (uncompactedCellIndex % Mx_My) / Mx;
    i.x = (uncompactedCellIndex % Mx_My) % Mx;
 
    Scratch::ijk_t tmp;
    tmp.i = (uint16_t)i.x;
    tmp.j = (uint16_t)i.y;
    tmp.k = (uint16_t)i.z;
 
    tmp.code = code;
    tmp.axes = g->axesTable[code];
    g->actCellIJKOut[c] = tmp;
 
    //const auto i = delinarize(uncompactedCellIndex, g->M, g->gridA);
 
    const auto axes = g->axesTable[code];
    const auto Nv = ((axes >> 2) & 1) + ((axes >> 1) & 1) + ((axes >> 0) & 1);
 
    const auto indexOut = (i.x < Ux) && (i.y < Uy) && (i.z < Uz);
    const auto Ni = indexOut ? g->indexCountTable[code] : 0;
 
    g->actCellCasesOut[c] = code;
    g->actCellIndicesOut[c] = uncompactedCellIndex;
    g->actCellVtxCntOut[c] = Nv;
    g->actCellIdxCntOut[c] = Ni;
  }
 
  if (g->elapsed_us != nullptr) {
    g->elapsed_us->fetch_add(timer.elapsedMicroseconds());
  }
}
 
#endif