45 if (levelSet ==
nullptr) {
47 .addWarning(
"No level set was passed to CalculateNormalVectors.")
51 if (levelSet->getLevelSetWidth() < (maxValue * 4) + 1) {
53 .addWarning(
"CalculateNormalVectors: Level set width must be "
55 std::to_string((maxValue * 4) + 1) +
"!")
59 std::vector<std::vector<Vec3D<T>>> normalVectorsVector(
60 levelSet->getNumberOfSegments());
61 double pointsPerSegment =
62 double(2 * levelSet->getDomain().getNumberOfPoints()) /
63 double(levelSet->getLevelSetWidth());
65 auto grid = levelSet->getGrid();
68#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
72 p = omp_get_thread_num();
75 auto &normalVectors = normalVectorsVector[p];
76 normalVectors.reserve(pointsPerSegment);
78 viennahrle::Index<D>
const startVector =
79 (p == 0) ? grid.getMinGridPoint()
80 : levelSet->getDomain().getSegmentation()[p - 1];
82 viennahrle::Index<D>
const endVector =
83 (p !=
static_cast<int>(levelSet->getNumberOfSegments() - 1))
84 ? levelSet->getDomain().getSegmentation()[p]
85 : grid.incrementIndices(grid.getMaxGridPoint());
87 for (viennahrle::ConstSparseStarIterator<
89 neighborIt(levelSet->getDomain(), startVector);
90 neighborIt.getIndices() < endVector; neighborIt.next()) {
92 auto ¢er = neighborIt.getCenter();
93 if (!center.isDefined()) {
95 }
else if (std::abs(center.getValue()) > maxValue) {
98 normalVectors.push_back(tmp);
105 for (
int i = 0; i <
D; i++) {
106 T pos = neighborIt.getNeighbor(i).getValue() - center.getValue();
107 T neg = center.getValue() - neighborIt.getNeighbor(i +
D).getValue();
108 n[i] = (pos + neg) * 0.5;
109 denominator += n[i] * n[i];
112 denominator = std::sqrt(denominator);
113 if (std::abs(denominator) < 1e-12) {
114 std::ostringstream oss;
115 oss <<
"CalculateNormalVectors: Vector of length 0 at "
116 << neighborIt.getIndices();
117 Logger::getInstance().addWarning(oss.str()).print();
118 for (
unsigned i = 0; i <
D; ++i)
121 for (
unsigned i = 0; i <
D; ++i) {
126 normalVectors.push_back(n);
131 unsigned numberOfNormals = 0;
132 for (
unsigned i = 0; i < levelSet->getNumberOfSegments(); ++i) {
133 numberOfNormals += normalVectorsVector[i].size();
135 normalVectorsVector[0].reserve(numberOfNormals);
137 for (
unsigned i = 1; i < levelSet->getNumberOfSegments(); ++i) {
138 normalVectorsVector[0].insert(normalVectorsVector[0].end(),
139 normalVectorsVector[i].begin(),
140 normalVectorsVector[i].end());
144 auto &pointData = levelSet->getPointData();
147 if (vectorDataPointer ==
nullptr) {
148 pointData.insertNextVectorData(normalVectorsVector[0],
152 *vectorDataPointer = std::move(normalVectorsVector[0]);