lsRemoveStrayPoints.hpp

Go to the documentation of this file.

#pragma once
 
#include <lsDomain.hpp>
#include <lsMarkVoidPoints.hpp>
 
namespace viennals {
 
using namespace viennacore;
 
template <class T, int D> class RemoveStrayPoints {
  SmartPointer<Domain<T, D>> levelSet = nullptr;
  VoidTopSurfaceEnum voidTopSurface = VoidTopSurfaceEnum::LARGEST;
 
public:
  RemoveStrayPoints() {}
 
  RemoveStrayPoints(SmartPointer<Domain<T, D>> passedLevelSet)
      : levelSet(passedLevelSet) {}
 
  void setLevelSet(SmartPointer<Domain<T, D>> passedLevelSet) {
    levelSet = passedLevelSet;
  }
 
  void setVoidTopSurface(VoidTopSurfaceEnum topSurface) {
    voidTopSurface = topSurface;
  }
 
  void apply() {
    if (levelSet == nullptr) {
      Logger::getInstance()
          .addWarning("No level set was passed to lsPrune.")
          .print();
      return;
    }
    if (levelSet->getNumberOfPoints() == 0) {
      return;
    }
 
    // Mark which points are voids
    {
      MarkVoidPoints<T, D> marker;
      marker.setLevelSet(levelSet);
      marker.setVoidTopSurface(voidTopSurface);
      marker.apply();
    }
 
    auto voidMarkers =
        levelSet->getPointData().getScalarData("VoidPointMarkers", true);
    if (voidMarkers == nullptr) {
      Logger::getInstance()
          .addWarning("RemoveStrayPoints: No scalar data for void point "
                      "markers found. Cannot remove stray points.")
          .print();
    }
 
    // now iterate through the domain and remove points which are void points
    auto &grid = levelSet->getGrid();
    auto newlsDomain = SmartPointer<Domain<T, D>>::New(grid);
    typename Domain<T, D>::DomainType &newDomain = newlsDomain->getDomain();
    typename Domain<T, D>::DomainType &domain = levelSet->getDomain();
 
    newDomain.initialize(domain.getNewSegmentation(), domain.getAllocation());
 
    std::vector<typename Domain<T, D>::PointValueVectorType> newPoints;
    newPoints.resize(newDomain.getNumberOfSegments());
 
#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
    {
      int p = 0;
#ifdef _OPENMP
      p = omp_get_thread_num();
#endif
 
      auto &domainSegment = newDomain.getDomainSegment(p);
 
      hrleVectorType<hrleIndexType, D> startVector =
          (p == 0) ? grid.getMinGridPoint()
                   : newDomain.getSegmentation()[p - 1];
 
      hrleVectorType<hrleIndexType, D> endVector =
          (p != static_cast<int>(newDomain.getNumberOfSegments() - 1))
              ? newDomain.getSegmentation()[p]
              : grid.incrementIndices(grid.getMaxGridPoint());
 
      for (hrleConstSparseIterator<typename Domain<T, D>::DomainType> it(
               domain, startVector);
           it.getStartIndices() < endVector; it.next()) {
        if (it.isDefined() && !voidMarkers->at(it.getPointId())) {
          newPoints[p].push_back(
              std::make_pair(it.getStartIndices(), it.getValue()));
        }
      }
    }
 
    // merge points
    for (unsigned i = 1; i < newDomain.getNumberOfSegments(); ++i) {
      newPoints[0].insert(newPoints[0].end(), newPoints[i].begin(),
                          newPoints[i].end());
    }
 
    newlsDomain->insertPoints(newPoints[0]);
 
    // distribute evenly across segments and copy
    newDomain.finalize();
    newDomain.segment();
    levelSet->deepCopy(newlsDomain);
    levelSet->finalize(2);
  }
};
 
} // namespace viennals