ViennaLS
Loading...
Searching...
No Matches
lsExpand.hpp
Go to the documentation of this file.
1#pragma once
2
4
5#include <hrleSparseStarIterator.hpp>
6#include <lsDomain.hpp>
7
8#include <vcVectorType.hpp>
9
10namespace viennals {
11
12using namespace viennacore;
13
17template <class T, int D> class Expand {
18 SmartPointer<Domain<T, D>> levelSet = nullptr;
19 int width = 0;
20 bool updatePointData = true;
21
22public:
23 Expand() = default;
24
25 Expand(SmartPointer<Domain<T, D>> passedlsDomain)
26 : levelSet(passedlsDomain) {}
27
28 Expand(SmartPointer<Domain<T, D>> passedlsDomain, int passedWidth)
29 : levelSet(passedlsDomain), width(passedWidth) {}
30
31 void setLevelSet(SmartPointer<Domain<T, D>> passedlsDomain) {
32 levelSet = passedlsDomain;
33 }
34
37 void setWidth(int passedWidth) { width = passedWidth; }
38
41 void setUpdatePointData(bool update) { updatePointData = update; }
42
44 void apply() {
45 if (levelSet == nullptr) {
46 Logger::getInstance()
47 .addWarning("No level set passed to Expand. Not expanding.")
48 .print();
49 return;
50 }
51
52 if (width <= levelSet->getLevelSetWidth())
53 return;
54
55 if (levelSet->getNumberOfPoints() == 0)
56 return;
57
58 const T totalLimit = width * 0.5;
59 const int startWidth = levelSet->getLevelSetWidth();
60 const int numberOfRequiredCycles = width - startWidth;
61
62 for (int currentCycle = 0; currentCycle < numberOfRequiredCycles;
63 ++currentCycle) {
64
65 const int allocationFactor =
66 1 + 1.0 / static_cast<double>(startWidth + currentCycle);
67 const T limit = (startWidth + currentCycle + 1) * T(0.5);
68
69 auto &grid = levelSet->getGrid();
70 auto newlsDomain = SmartPointer<Domain<T, D>>::New(grid);
71 auto &newDomain = newlsDomain->getDomain();
72 auto &domain = levelSet->getDomain();
73
74 newDomain.initialize(domain.getNewSegmentation(),
75 domain.getAllocation() * allocationFactor);
76
77 const bool updateData = updatePointData;
78 // save how data should be transferred to new level set
79 // list of indices into the old pointData vector
80 std::vector<std::vector<unsigned>> newDataSourceIds;
81 if (updateData)
82 newDataSourceIds.resize(newDomain.getNumberOfSegments());
83
84#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
85 {
86 int p = 0;
87#ifdef _OPENMP
88 p = omp_get_thread_num();
89#endif
90
91 auto &domainSegment = newDomain.getDomainSegment(p);
92
93 viennahrle::Index<D> const startVector =
94 (p == 0) ? grid.getMinGridPoint()
95 : newDomain.getSegmentation()[p - 1];
96
97 viennahrle::Index<D> const endVector =
98 (p != static_cast<int>(newDomain.getNumberOfSegments() - 1))
99 ? newDomain.getSegmentation()[p]
100 : grid.incrementIndices(grid.getMaxGridPoint());
101
102 for (viennahrle::SparseStarIterator<typename Domain<T, D>::DomainType,
103 1>
104 neighborIt(domain, startVector);
105 neighborIt.getIndices() < endVector; neighborIt.next()) {
106
107 auto &centerIt = neighborIt.getCenter();
108 if (std::abs(centerIt.getValue()) <= totalLimit) {
109 domainSegment.insertNextDefinedPoint(neighborIt.getIndices(),
110 centerIt.getValue());
111 if (updateData)
112 newDataSourceIds[p].push_back(centerIt.getPointId());
113 } else {
114 if (centerIt.getValue() > -std::numeric_limits<T>::epsilon()) {
115 T distance = Domain<T, D>::POS_VALUE;
116 int neighbor = -1;
117 for (int i = 0; i < 2 * D; i++) {
118 T newValue = neighborIt.getNeighbor(i).getValue() + T(1);
119 if (distance > newValue) {
120 distance = newValue;
121 neighbor = i;
122 }
123 }
124 if (distance <= limit) {
125 domainSegment.insertNextDefinedPoint(neighborIt.getIndices(),
126 distance);
127 if (updateData)
128 newDataSourceIds[p].push_back(
129 neighborIt.getNeighbor(neighbor).getPointId());
130 } else {
131 // TODO: use insertNextUndefinedRunType
132 domainSegment.insertNextUndefinedPoint(neighborIt.getIndices(),
134 }
135 } else {
136 T distance = Domain<T, D>::NEG_VALUE;
137 int neighbor = -1;
138 for (int i = 0; i < 2 * D; i++) {
139 T newValue = neighborIt.getNeighbor(i).getValue() - T(1);
140 if (distance < newValue) {
141 distance = newValue;
142 neighbor = i;
143 }
144 }
145 if (distance >= -limit) {
146 domainSegment.insertNextDefinedPoint(neighborIt.getIndices(),
147 distance);
148 if (updateData)
149 newDataSourceIds[p].push_back(
150 neighborIt.getNeighbor(neighbor).getPointId());
151 } else {
152 // TODO: use insertNextUndefinedRunType
153 domainSegment.insertNextUndefinedPoint(neighborIt.getIndices(),
155 }
156 }
157 }
158 }
159 }
160
161 // now copy old data into new level set
162 if (updateData) {
163 newlsDomain->getPointData().translateFromMultiData(
164 levelSet->getPointData(), newDataSourceIds);
165 }
166
167 newDomain.finalize();
168 levelSet->deepCopy(newlsDomain);
169 }
170 levelSet->getDomain().segment();
171 levelSet->finalize(width);
172 }
173};
174
175// add all template specialisations for this class
177
178} // namespace viennals
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
static constexpr T NEG_VALUE
Definition lsDomain.hpp:52
static constexpr T POS_VALUE
Definition lsDomain.hpp:51
Expand(SmartPointer< Domain< T, D > > passedlsDomain, int passedWidth)
Definition lsExpand.hpp:28
void setUpdatePointData(bool update)
Set whether to update the point data stored in the LS during this algorithm. Defaults to true.
Definition lsExpand.hpp:41
void apply()
Apply the expansion to the specified width.
Definition lsExpand.hpp:44
Expand()=default
Expand(SmartPointer< Domain< T, D > > passedlsDomain)
Definition lsExpand.hpp:25
void setWidth(int passedWidth)
Set how far the level set should be extended. Points with value width*0.5 will be added by this algor...
Definition lsExpand.hpp:37
void setLevelSet(SmartPointer< Domain< T, D > > passedlsDomain)
Definition lsExpand.hpp:31
#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