ViennaLS
Loading...
Searching...
No Matches
lsCalculateNormalVectors.hpp
Go to the documentation of this file.
1#pragma once
2
4
5#include <algorithm>
6
7#include <hrleSparseStarIterator.hpp>
8
9#include <lsDomain.hpp>
10
11#include <vcLogger.hpp>
12#include <vcSmartPointer.hpp>
13#include <vcVectorType.hpp>
14
15namespace viennals {
16
17using namespace viennacore;
18
25template <class T, int D> class CalculateNormalVectors {
26 SmartPointer<Domain<T, D>> levelSet = nullptr;
27 T maxValue = 0.5;
28
29public:
30 static constexpr char normalVectorsLabel[] = "Normals";
31
33
34 CalculateNormalVectors(SmartPointer<Domain<T, D>> passedLevelSet,
35 T passedMaxValue = 0.5)
36 : levelSet(passedLevelSet), maxValue(passedMaxValue) {}
37
38 void setLevelSet(SmartPointer<Domain<T, D>> passedLevelSet) {
39 levelSet = passedLevelSet;
40 }
41
42 void setMaxValue(const T passedMaxValue) { maxValue = passedMaxValue; }
43
44 void apply() {
45 if (levelSet == nullptr) {
46 Logger::getInstance()
47 .addWarning("No level set was passed to CalculateNormalVectors.")
48 .print();
49 }
50
51 if (levelSet->getLevelSetWidth() < (maxValue * 4) + 1) {
52 Logger::getInstance()
53 .addWarning("CalculateNormalVectors: Level set width must be "
54 "greater than " +
55 std::to_string((maxValue * 4) + 1) + "!")
56 .print();
57 }
58
59 std::vector<std::vector<Vec3D<T>>> normalVectorsVector(
60 levelSet->getNumberOfSegments());
61 double pointsPerSegment =
62 double(2 * levelSet->getDomain().getNumberOfPoints()) /
63 double(levelSet->getLevelSetWidth());
64
65 auto grid = levelSet->getGrid();
66
68#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
69 {
70 int p = 0;
71#ifdef _OPENMP
72 p = omp_get_thread_num();
73#endif
74
75 auto &normalVectors = normalVectorsVector[p];
76 normalVectors.reserve(pointsPerSegment);
77
78 viennahrle::Index<D> const startVector =
79 (p == 0) ? grid.getMinGridPoint()
80 : levelSet->getDomain().getSegmentation()[p - 1];
81
82 viennahrle::Index<D> const endVector =
83 (p != static_cast<int>(levelSet->getNumberOfSegments() - 1))
84 ? levelSet->getDomain().getSegmentation()[p]
85 : grid.incrementIndices(grid.getMaxGridPoint());
86
87 for (viennahrle::ConstSparseStarIterator<
88 typename Domain<T, D>::DomainType, 1>
89 neighborIt(levelSet->getDomain(), startVector);
90 neighborIt.getIndices() < endVector; neighborIt.next()) {
91
92 auto &center = neighborIt.getCenter();
93 if (!center.isDefined()) {
94 continue;
95 } else if (std::abs(center.getValue()) > maxValue) {
96 // push an empty vector to keep ordering correct
97 Vec3D<T> tmp = {};
98 normalVectors.push_back(tmp);
99 continue;
100 }
101
102 Vec3D<T> n;
103
104 T denominator = 0;
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];
110 }
111
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)
119 n[i] = 0.;
120 } else {
121 for (unsigned i = 0; i < D; ++i) {
122 n[i] /= denominator;
123 }
124 }
125
126 normalVectors.push_back(n);
127 }
128 }
129
130 // copy all normals
131 unsigned numberOfNormals = 0;
132 for (unsigned i = 0; i < levelSet->getNumberOfSegments(); ++i) {
133 numberOfNormals += normalVectorsVector[i].size();
134 }
135 normalVectorsVector[0].reserve(numberOfNormals);
136
137 for (unsigned i = 1; i < levelSet->getNumberOfSegments(); ++i) {
138 normalVectorsVector[0].insert(normalVectorsVector[0].end(),
139 normalVectorsVector[i].begin(),
140 normalVectorsVector[i].end());
141 }
142
143 // insert into pointData of levelSet
144 auto &pointData = levelSet->getPointData();
145 auto vectorDataPointer = pointData.getVectorData(normalVectorsLabel, true);
146 // if it does not exist, insert new normals vector
147 if (vectorDataPointer == nullptr) {
148 pointData.insertNextVectorData(normalVectorsVector[0],
150 } else {
151 // if it does exist, just swap the old with the new values
152 *vectorDataPointer = std::move(normalVectorsVector[0]);
153 }
154 }
155};
156
157// add all template specialisations for this class
158PRECOMPILE_PRECISION_DIMENSION(CalculateNormalVectors)
159
160} // namespace viennals
CalculateNormalVectors(SmartPointer< Domain< T, D > > passedLevelSet, T passedMaxValue=0.5)
Definition lsCalculateNormalVectors.hpp:34
static constexpr char normalVectorsLabel[]
Definition lsCalculateNormalVectors.hpp:30
void setLevelSet(SmartPointer< Domain< T, D > > passedLevelSet)
Definition lsCalculateNormalVectors.hpp:38
void apply()
Definition lsCalculateNormalVectors.hpp:44
void setMaxValue(const T passedMaxValue)
Definition lsCalculateNormalVectors.hpp:42
Class containing all information about the level set, including the dimensions of the domain,...
Definition lsDomain.hpp:27
viennahrle::Domain< T, D > DomainType
Definition lsDomain.hpp:32
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
Definition lsAdvect.hpp:36
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68