lsBooleanOperation.hpp

Go to the documentation of this file.

#pragma once
 
#include <lsPreCompileMacros.hpp>
 
#include <hrleSparseStarIterator.hpp>
 
#include <lsDomain.hpp>
#include <lsPrune.hpp>
 
#include <vcLogger.hpp>
#include <vcSmartPointer.hpp>
#include <vcVectorType.hpp>
 
#include <functional>
 
namespace viennals {
 
using namespace viennacore;

enum struct BooleanOperationEnum : unsigned {
  INTERSECT = 0,
  UNION = 1,
  RELATIVE_COMPLEMENT = 2,
  INVERT = 3,
  CUSTOM = 4
};

template <class T, int D> class BooleanOperation {
public:
  using ComparatorType =
      std::function<std::pair<T, bool>(const T &, const T &)>;
 
private:
  typedef typename Domain<T, D>::DomainType hrleDomainType;
  SmartPointer<Domain<T, D>> levelSetA = nullptr;
  SmartPointer<Domain<T, D>> levelSetB = nullptr;
  BooleanOperationEnum operation = BooleanOperationEnum::INTERSECT;
  ComparatorType operationComp = nullptr;
  bool updatePointData = true;
  bool pruneResult = true;
 
  void booleanOpInternal(ComparatorType comp) {
    auto &grid = levelSetA->getGrid();
    auto newlsDomain = SmartPointer<Domain<T, D>>::New(grid);
    typename Domain<T, D>::DomainType &newDomain = newlsDomain->getDomain();
    typename Domain<T, D>::DomainType &domain = levelSetA->getDomain();
 
    newDomain.initialize(domain.getNewSegmentation(), domain.getAllocation());
 
    const bool updateData = updatePointData;
    // save how data should be transferred to new level set
    // list of indices into the old pointData vector
    std::vector<std::vector<unsigned>> newDataSourceIds;
    std::vector<std::vector<bool>> newDataLS;
    if (updateData) {
      newDataSourceIds.resize(newDomain.getNumberOfSegments());
      newDataLS.resize(newDataSourceIds.size());
    }
 
#pragma omp parallel for
    for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
 
      auto &domainSegment = newDomain.getDomainSegment(p);
 
      viennahrle::Index<D> currentVector =
          (p == 0) ? grid.getMinGridPoint()
                   : newDomain.getSegmentation()[p - 1];
 
      viennahrle::Index<D> const endVector =
          (p != static_cast<int>(newDomain.getNumberOfSegments() - 1))
              ? newDomain.getSegmentation()[p]
              : grid.incrementIndices(grid.getMaxGridPoint());
 
      viennahrle::ConstSparseIterator<hrleDomainType> itA(
          levelSetA->getDomain(), currentVector);
      viennahrle::ConstSparseIterator<hrleDomainType> itB(
          levelSetB->getDomain(), currentVector);
 
      while (currentVector < endVector) {
        const auto &comparison = comp(itA.getValue(), itB.getValue());
        const auto &currentValue = comparison.first;
 
        if (currentValue != Domain<T, D>::NEG_VALUE &&
            currentValue != Domain<T, D>::POS_VALUE) {
          domainSegment.insertNextDefinedPoint(currentVector, currentValue);
          if (updateData) {
            // if taken from A, set to true
            const bool originLS = comparison.second;
            newDataLS[p].push_back(originLS);
            const auto originPointId =
                (originLS) ? itA.getPointId() : itB.getPointId();
            newDataSourceIds[p].push_back(originPointId);
          }
        } else {
          domainSegment.insertNextUndefinedPoint(
              currentVector, (currentValue < 0) ? Domain<T, D>::NEG_VALUE
                                                : Domain<T, D>::POS_VALUE);
        }
 
        switch (Compare(itA.getEndIndices(), itB.getEndIndices())) {
        case -1:
          itA.next();
          break;
        case 0:
          itA.next();
        default:
          itB.next();
        }
        currentVector =
            Max(itA.getStartIndices().get(), itB.getStartIndices().get());
      }
    }
 
    // merge data
    for (unsigned i = 1; i < newDataLS.size(); ++i) {
      newDataLS[0].insert(newDataLS[0].end(), newDataLS[i].begin(),
                          newDataLS[i].end());
 
      newDataSourceIds[0].insert(newDataSourceIds[0].end(),
                                 newDataSourceIds[i].begin(),
                                 newDataSourceIds[i].end());
    }
 
    // transfer data from the old LSs to new LS
    // Only do so if the same data exists in both LSs
    // If this is not the case, the data is invalid
    // and therefore not needed anyway.
    if (updateData) {
      const auto &AData = levelSetA->getPointData();
      const auto &BData = levelSetB->getPointData();
      auto &newData = newlsDomain->getPointData();
 
      // scalars — preserve A's fields even when B doesn't have them.
      // Previously, fields present only on A (e.g. OxVelocity on the oxide
      // in a RELATIVE_COMPLEMENT with the mask) were silently dropped because
      // the mask has no matching field.  For A-sourced result points we always
      // use A's value; for B-sourced points without a matching B field we use
      // 0.
      for (unsigned i = 0; i < AData.getScalarDataSize(); ++i) {
        auto scalarDataLabel = AData.getScalarDataLabel(i);
        auto APointer = AData.getScalarData(i);
        auto BPointer =
            BData.getScalarData(scalarDataLabel, true); // may be nullptr
        typename Domain<T, D>::PointDataType::ScalarDataType scalars;
        scalars.resize(newlsDomain->getNumberOfPoints(), T(0));
        for (unsigned j = 0; j < newlsDomain->getNumberOfPoints(); ++j) {
          if (newDataLS[0][j])
            scalars[j] = APointer->at(newDataSourceIds[0][j]);
          else if (BPointer != nullptr)
            scalars[j] = BPointer->at(newDataSourceIds[0][j]);
          // else: B-sourced point with no matching B field → stays 0
        }
        newData.insertNextScalarData(scalars, scalarDataLabel);
      }
 
      // vectors — same asymmetric treatment as scalars above.
      for (unsigned i = 0; i < AData.getVectorDataSize(); ++i) {
        auto vectorDataLabel = AData.getVectorDataLabel(i);
        auto APointer = AData.getVectorData(i);
        auto BPointer =
            BData.getVectorData(vectorDataLabel, true); // may be nullptr
        using Vec =
            typename Domain<T, D>::PointDataType::VectorDataType::value_type;
        typename Domain<T, D>::PointDataType::VectorDataType vectors;
        vectors.resize(newlsDomain->getNumberOfPoints(), Vec{});
        for (unsigned j = 0; j < newlsDomain->getNumberOfPoints(); ++j) {
          if (newDataLS[0][j])
            vectors[j] = APointer->at(newDataSourceIds[0][j]);
          else if (BPointer != nullptr)
            vectors[j] = BPointer->at(newDataSourceIds[0][j]);
          // else: B-sourced point with no matching B field → stays
          // zero-initialised
        }
        newData.insertNextVectorData(vectors, vectorDataLabel);
      }
    }
 
    newDomain.finalize();
    newDomain.segment();
    newlsDomain->setLevelSetWidth(levelSetA->getLevelSetWidth());
 
    if (pruneResult) {
      auto pruner = Prune<T, D>(newlsDomain);
      pruner.setRemoveStrayZeros(true);
      pruner.apply();
 
      // now we need to prune, to remove stray defined points
      Prune<T, D>(newlsDomain).apply();
    }
 
    levelSetA->deepCopy(newlsDomain);
  }
 
  void invert() {
    auto &hrleDomain = levelSetA->getDomain();
#pragma omp parallel for
    for (unsigned p = 0; p < hrleDomain.getNumberOfSegments(); ++p) {
      auto &domainSegment = hrleDomain.getDomainSegment(p);
 
      // change all defined values
      for (unsigned i = 0; i < domainSegment.definedValues.size(); ++i) {
        domainSegment.definedValues[i] = -domainSegment.definedValues[i];
      }
 
      // add undefined values if missing
      if (domainSegment.undefinedValues.size() < 1) {
        domainSegment.undefinedValues.push_back(T(Domain<T, D>::NEG_VALUE));
      }
      if (domainSegment.undefinedValues.size() < 2) {
        if (domainSegment.undefinedValues[0] == Domain<T, D>::NEG_VALUE) {
          domainSegment.undefinedValues.push_back(T(Domain<T, D>::POS_VALUE));
        } else {
          domainSegment.undefinedValues.push_back(T(Domain<T, D>::NEG_VALUE));
        }
      }
 
      // change all runTypes
      // there are only two undefined runs: negative undefined (UNDEF_PT)
      // and positive undefined (UNDEF_PT+1)
      for (unsigned dim = 0; dim < D; ++dim) {
        for (unsigned c = 0; c < domainSegment.runTypes[dim].size(); ++c) {
          if (domainSegment.runTypes[dim][c] ==
              viennahrle::RunTypeValues::UNDEF_PT) {
            domainSegment.runTypes[dim][c] =
                viennahrle::RunTypeValues::UNDEF_PT + 1;
          } else if (domainSegment.runTypes[dim][c] ==
                     viennahrle::RunTypeValues::UNDEF_PT + 1) {
            domainSegment.runTypes[dim][c] =
                viennahrle::RunTypeValues::UNDEF_PT;
          }
        }
      }
    }
    levelSetA->finalize();
  }
 
  static std::pair<T, bool> minComp(const T &a, const T &b) {
    bool AIsSmaller = a < b;
    if (AIsSmaller)
      return std::make_pair(a, true);
    else
      return std::make_pair(b, false);
  }
 
  static std::pair<T, bool> maxComp(const T &a, const T &b) {
    bool AIsLarger = a > b;
    if (AIsLarger)
      return std::make_pair(a, true);
    else
      return std::make_pair(b, false);
  }
 
  static std::pair<T, bool> relativeComplementComp(const T &a, const T &b) {
    return maxComp(a, -b);
  }
 
public:
  BooleanOperation() = default;
 
  BooleanOperation(
      SmartPointer<Domain<T, D>> passedlsDomain,
      BooleanOperationEnum passedOperation = BooleanOperationEnum::INVERT)
      : levelSetA(passedlsDomain), operation(passedOperation) {}
 
  BooleanOperation(
      SmartPointer<Domain<T, D>> passedlsDomainA,
      SmartPointer<Domain<T, D>> passedlsDomainB,
      BooleanOperationEnum passedOperation = BooleanOperationEnum::INTERSECT)
      : levelSetA(passedlsDomainA), levelSetB(passedlsDomainB),
        operation(passedOperation){};

  void setLevelSet(SmartPointer<Domain<T, D>> passedlsDomain) {
    levelSetA = passedlsDomain;
  }

  void setSecondLevelSet(SmartPointer<Domain<T, D>> passedlsDomain) {
    levelSetB = passedlsDomain;
  }

  void setBooleanOperation(BooleanOperationEnum passedOperation) {
    operation = passedOperation;
  }

  void setBooleanOperationComparator(ComparatorType passedOperationComp) {
    operationComp = passedOperationComp;
  }

  void setUpdatePointData(bool update) { updatePointData = update; }

  void setPruneResult(bool pR) { pruneResult = pR; }

  void apply() {
    if (levelSetA == nullptr) {
      VIENNACORE_LOG_ERROR("No level set was passed to BooleanOperation.");
      return;
    }
 
    if (static_cast<unsigned>(operation) < 3) {
      if (levelSetB == nullptr) {
        VIENNACORE_LOG_ERROR(
            "Only one level set was passed to BooleanOperation, "
            "although two were required.");
        return;
      }
    }
 
    switch (operation) {
    case BooleanOperationEnum::INTERSECT:
      booleanOpInternal(&BooleanOperation::maxComp);
      break;
    case BooleanOperationEnum::UNION:
      booleanOpInternal(&BooleanOperation::minComp);
      break;
    case BooleanOperationEnum::RELATIVE_COMPLEMENT:
      booleanOpInternal(&BooleanOperation::relativeComplementComp);
      break;
    case BooleanOperationEnum::INVERT:
      invert();
      break;
    case BooleanOperationEnum::CUSTOM:
      if (operationComp == nullptr) {
        VIENNACORE_LOG_ERROR(
            "No comparator supplied to custom BooleanOperation.");
        return;
      }
      booleanOpInternal(operationComp);
    }
  }
};
 
// add all template specialisations for this class
PRECOMPILE_PRECISION_DIMENSION(BooleanOperation)
 
} // namespace viennals