ViennaLS
Loading...
Searching...
No Matches
lsCompareNarrowBand.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <hrleDenseCellIterator.hpp>
4#include <lsDomain.hpp>
5#include <lsExpand.hpp>
6#include <lsMesh.hpp>
8
9#include <unordered_map>
10
11namespace viennals {
12
13using namespace viennacore;
14
19template <class T, int D = 2> class CompareNarrowBand {
20 using hrleIndexType = viennahrle::IndexType;
21 SmartPointer<Domain<T, D>> levelSetTarget = nullptr;
22 SmartPointer<Domain<T, D>> levelSetSample = nullptr;
23 viennahrle::Index<D> minIndex, maxIndex;
24
25 // Variables for x and y range restrictions
26 T xRangeMin = std::numeric_limits<T>::lowest();
27 T xRangeMax = std::numeric_limits<T>::max();
28 T yRangeMin = std::numeric_limits<T>::lowest();
29 T yRangeMax = std::numeric_limits<T>::max();
30 bool useXRange = false;
31 bool useYRange = false;
32
33 // Fields to store the calculation results
34 T sumSquaredDifferences = 0.0;
35 T sumDifferences = 0.0;
36 unsigned numPoints = 0;
37
38 // Add mesh output capability
39 SmartPointer<Mesh<T>> outputMesh = nullptr;
40 bool outputMeshSquaredDifferences = true;
41
42 bool checkAndCalculateBounds() {
43 if (levelSetTarget == nullptr || levelSetSample == nullptr) {
44 Logger::getInstance()
45 .addWarning("Missing level set in CompareNarrowBand.")
46 .print();
47 return false;
48 }
49
50 // Check if the grids are compatible
51 const auto &gridTarget = levelSetTarget->getGrid();
52 const auto &gridSample = levelSetSample->getGrid();
53
54 if (gridTarget.getGridDelta() != gridSample.getGridDelta()) {
55 Logger::getInstance()
56 .addWarning("Grid delta mismatch in CompareNarrowBand. The grid "
57 "deltas of the two level sets must be equal.")
58 .print();
59 return false;
60 }
61
62 // Check if the x extents of both level sets are equal
63 const auto &domainTarget = levelSetTarget->getDomain();
64 const auto &domainSample = levelSetSample->getDomain();
65
66 hrleIndexType targetMinX = gridTarget.isNegBoundaryInfinite(0)
67 ? domainTarget.getMinRunBreak(0)
68 : gridTarget.getMinIndex(0);
69 hrleIndexType targetMaxX = gridTarget.isPosBoundaryInfinite(0)
70 ? domainTarget.getMaxRunBreak(0)
71 : gridTarget.getMaxIndex(0);
72 hrleIndexType sampleMinX = gridSample.isNegBoundaryInfinite(0)
73 ? domainSample.getMinRunBreak(0)
74 : gridSample.getMinIndex(0);
75 hrleIndexType sampleMaxX = gridSample.isPosBoundaryInfinite(0)
76 ? domainSample.getMaxRunBreak(0)
77 : gridSample.getMaxIndex(0);
78
79 if (targetMinX != sampleMinX || targetMaxX != sampleMaxX) {
80 Logger::getInstance()
81 .addWarning("X extent mismatch in CompareNarrowBand. The x extents "
82 "of both level sets must be equal.")
83 .print();
84 return false;
85 }
86
87 // Expand the sample level set using lsExpand to a default width of 5
88 if (levelSetSample->getLevelSetWidth() < 5) {
89 Logger::getInstance()
90 .addWarning("Sample level set width is insufficient. Expanding it to "
91 "a width of 5.")
92 .print();
93 Expand<T, D>(levelSetSample, 5).apply();
94 }
95
96 // Check if target level set width is sufficient
97 if (levelSetTarget->getLevelSetWidth() <
98 levelSetSample->getLevelSetWidth() + 50) {
99 Logger::getInstance()
100 .addWarning(
101 "Target level set width is insufficient. It must exceed sample "
102 "width by least 50. \n"
103 " CORRECTION: The expansion was performed. \n"
104 "ALTERNATIVE: Alternatively, please expand the target yourself "
105 "using lsExpand before passing it to this function. \n")
106 .print();
107 Expand<T, D>(levelSetTarget, levelSetSample->getLevelSetWidth() + 50)
108 .apply();
109 }
110
111 // Initialize min and max indices
112 for (unsigned i = 0; i < D; ++i) {
113 minIndex[i] = std::numeric_limits<hrleIndexType>::max();
114 maxIndex[i] = std::numeric_limits<hrleIndexType>::lowest();
115 }
116
117 // Calculate actual bounds
118 for (unsigned i = 0; i < D; ++i) {
119 minIndex[i] = std::min({minIndex[i],
120 (gridTarget.isNegBoundaryInfinite(i))
121 ? domainTarget.getMinRunBreak(i)
122 : gridTarget.getMinIndex(i),
123 (gridSample.isNegBoundaryInfinite(i))
124 ? domainSample.getMinRunBreak(i)
125 : gridSample.getMinIndex(i)});
126
127 maxIndex[i] = std::max({maxIndex[i],
128 (gridTarget.isPosBoundaryInfinite(i))
129 ? domainTarget.getMaxRunBreak(i)
130 : gridTarget.getMaxIndex(i),
131 (gridSample.isPosBoundaryInfinite(i))
132 ? domainSample.getMaxRunBreak(i)
133 : gridSample.getMaxIndex(i)});
134 }
135
136 return true;
137 }
138
139public:
141 static_assert(
142 D == 2 &&
143 "CompareNarrowBand is currently only implemented for 2D level sets.");
144 }
145
146 CompareNarrowBand(SmartPointer<Domain<T, D>> passedLevelSetTarget,
147 SmartPointer<Domain<T, D>> passedlevelSetSample)
148 : levelSetTarget(passedLevelSetTarget),
149 levelSetSample(passedlevelSetSample) {
150 static_assert(
151 D == 2 &&
152 "CompareNarrowBand is currently only implemented for 2D level sets.");
153 }
154
155 void setLevelSetTarget(SmartPointer<Domain<T, D>> passedLevelSet) {
156 levelSetTarget = passedLevelSet;
157 }
158
159 void setLevelSetSample(SmartPointer<Domain<T, D>> passedLevelSet) {
160 levelSetSample = passedLevelSet;
161 }
162
164 void setXRange(T minXRange, T maxXRange) {
165 xRangeMin = minXRange;
166 xRangeMax = maxXRange;
167 useXRange = true;
168 }
169
171 void setYRange(T minYRange, T maxYRange) {
172 yRangeMin = minYRange;
173 yRangeMax = maxYRange;
174 useYRange = true;
175 }
176
178 void clearXRange() {
179 useXRange = false;
180 xRangeMin = std::numeric_limits<T>::lowest();
181 xRangeMax = std::numeric_limits<T>::max();
182 }
183
185 void clearYRange() {
186 useYRange = false;
187 yRangeMin = std::numeric_limits<T>::lowest();
188 yRangeMax = std::numeric_limits<T>::max();
189 }
190
192 void setOutputMesh(SmartPointer<Mesh<T>> passedMesh,
193 bool outputMeshSquaredDiffs = true) {
194 outputMesh = passedMesh;
195 outputMeshSquaredDifferences = outputMeshSquaredDiffs;
196 }
197
201 outputMeshSquaredDifferences = value;
202 }
203
205 void apply() {
206 // Perform compatibility checks and calculate bounds
207 if (!checkAndCalculateBounds()) {
208 // If checks fail, return NaN
209 sumSquaredDifferences = std::numeric_limits<T>::quiet_NaN();
210 numPoints = 0;
211 return;
212 }
213
214 const auto &gridTarget = levelSetTarget->getGrid();
215 double gridDelta = gridTarget.getGridDelta();
216
217 // Set up iterators for both level sets
218 viennahrle::ConstDenseCellIterator<typename Domain<T, D>::DomainType>
219 itSample(levelSetSample->getDomain(), minIndex);
220 viennahrle::ConstDenseCellIterator<typename Domain<T, D>::DomainType>
221 itTarget(levelSetTarget->getDomain(), minIndex);
222
223 sumSquaredDifferences = 0.0;
224 numPoints = 0;
225
226 // Prepare mesh output if needed
227 std::unordered_map<viennahrle::Index<D>, size_t,
228 typename viennahrle::Index<D>::hash>
229 pointIdMapping;
230 std::vector<T> differenceValues;
231 size_t currentPointId = 0;
232
233 const bool generateMesh = outputMesh != nullptr;
234 if (generateMesh) {
235 outputMesh->clear();
236
237 // Initialize mesh extent
238 for (unsigned i = 0; i < D; ++i) {
239 outputMesh->minimumExtent[i] = std::numeric_limits<T>::max();
240 outputMesh->maximumExtent[i] = std::numeric_limits<T>::lowest();
241 }
242 }
243
244 // Iterate through the domain defined by the bounding box
245 for (; itSample.getIndices() < maxIndex; itSample.next()) {
246 // Check if current point is within specified x and y ranges
247 T xCoord = itSample.getIndices()[0] * gridDelta;
248 T yCoord = (D > 1) ? itSample.getIndices()[1] * gridDelta : 0;
249
250 // Skip if outside the specified x-range
251 if (useXRange && (xCoord < xRangeMin || xCoord > xRangeMax)) {
252 continue;
253 }
254
255 // Skip if outside the specified y-range (only check in 2D and 3D)
256 if (D > 1 && useYRange && (yCoord < yRangeMin || yCoord > yRangeMax)) {
257 continue;
258 }
259
260 // Move the second iterator to the same position
261 itTarget.goToIndicesSequential(itSample.getIndices());
262
263 // Get values at current position
264 T valueTarget = 0.0;
265 T valueSample = 0.0;
266
267 // Calculate average value at cell center
268 for (int i = 0; i < (1 << D); ++i) {
269 valueSample += itSample.getCorner(i).getValue();
270 valueTarget += itTarget.getCorner(i).getValue();
271 }
272 valueTarget /= (1 << D);
273 valueSample /= (1 << D);
274
275 // Check for infinite or extreme values that might cause numerical issues
276 if (std::isinf(valueTarget) || std::isinf(valueSample) ||
277 std::abs(valueTarget) > 1000 || std::abs(valueSample) > 1000) {
278 continue;
279 }
280
281 // Calculate difference and add to sum
282 T diff = std::abs(valueTarget - valueSample) * gridDelta;
283 T diffSquared = diff * diff;
284 sumSquaredDifferences += diffSquared;
285 sumDifferences += diff;
286 numPoints++;
287
288 // Store difference in mesh if required
289 if (generateMesh) {
290 std::array<unsigned, 1 << D> voxel;
291 bool addVoxel = true;
292 // TODO: possibly remove this addVoxel check
293 // Insert all points of voxel into pointList
294 for (unsigned i = 0; i < (1 << D); ++i) {
295 viennahrle::Index<D> index;
296 for (unsigned j = 0; j < D; ++j) {
297 index[j] =
298 itSample.getIndices(j) + itSample.getCorner(i).getOffset()[j];
299 if (index[j] > maxIndex[j]) {
300 addVoxel = false;
301 break;
302 }
303 }
304 if (addVoxel) {
305 auto pointIdValue = std::make_pair(index, currentPointId);
306 auto pointIdPair = pointIdMapping.insert(pointIdValue);
307 voxel[i] = pointIdPair.first->second;
308 if (pointIdPair.second) {
309 ++currentPointId;
310 }
311 } else {
312 break;
313 }
314 }
315
316 if (addVoxel) {
317 if constexpr (D == 3) {
318 std::array<unsigned, 8> hexa{voxel[0], voxel[1], voxel[3],
319 voxel[2], voxel[4], voxel[5],
320 voxel[7], voxel[6]};
321 outputMesh->hexas.push_back(hexa);
322 } else if constexpr (D == 2) {
323 std::array<unsigned, 4> quad{voxel[0], voxel[1], voxel[3],
324 voxel[2]};
325 outputMesh->tetras.push_back(quad);
326 }
327
328 // Add difference value to cell data depending on whether squared
329 // differences are requested
330 differenceValues.push_back(outputMeshSquaredDifferences ? diffSquared
331 : diff);
332 }
333 }
334 }
335
336 // Finalize mesh output
337 if (generateMesh && !pointIdMapping.empty()) {
338 outputMesh->nodes.resize(pointIdMapping.size());
339 for (auto it = pointIdMapping.begin(); it != pointIdMapping.end(); ++it) {
340 std::array<T, 3> coords{};
341 for (unsigned i = 0; i < D; ++i) {
342 coords[i] = gridDelta * it->first[i];
343
344 if (coords[i] < outputMesh->minimumExtent[i]) {
345 outputMesh->minimumExtent[i] = coords[i];
346 } else if (coords[i] > outputMesh->maximumExtent[i]) {
347 outputMesh->maximumExtent[i] = coords[i];
348 }
349 }
350 outputMesh->nodes[it->second] = coords;
351 }
352
353 assert(differenceValues.size() ==
354 outputMesh->template getElements<1 << D>().size());
355 outputMesh->cellData.insertNextScalarData(std::move(differenceValues),
356 outputMeshSquaredDifferences
357 ? "Squared differences"
358 : "Absolute differences");
359 }
360 }
361
363 T getSumSquaredDifferences() const { return sumSquaredDifferences; }
364
365 // Return the sum of differences calculated by apply().
366 T getSumDifferences() const { return sumDifferences; }
367
369 unsigned getNumPoints() const { return numPoints; }
370
372 T getRMSE() const {
373 return (numPoints > 0) ? std::sqrt(sumSquaredDifferences / numPoints)
374 : std::numeric_limits<T>::infinity();
375 }
376};
377
378// Add template specializations for this class
379PRECOMPILE_PRECISION_DIMENSION(CompareNarrowBand)
380
381} // namespace viennals
void setYRange(T minYRange, T maxYRange)
Set the y-coordinate range to restrict the comparison area.
Definition lsCompareNarrowBand.hpp:171
void setOutputMeshSquaredDifferences(bool value)
Set whether to output squared differences (true) or absolute differences (false)
Definition lsCompareNarrowBand.hpp:200
void setXRange(T minXRange, T maxXRange)
Set the x-coordinate range to restrict the comparison area.
Definition lsCompareNarrowBand.hpp:164
void setOutputMesh(SmartPointer< Mesh< T > > passedMesh, bool outputMeshSquaredDiffs=true)
Set the output mesh where difference values will be stored.
Definition lsCompareNarrowBand.hpp:192
void setLevelSetTarget(SmartPointer< Domain< T, D > > passedLevelSet)
Definition lsCompareNarrowBand.hpp:155
T getRMSE() const
Calculate the root mean square error from previously computed values.
Definition lsCompareNarrowBand.hpp:372
void setLevelSetSample(SmartPointer< Domain< T, D > > passedLevelSet)
Definition lsCompareNarrowBand.hpp:159
unsigned getNumPoints() const
Return the number of points used in the comparison.
Definition lsCompareNarrowBand.hpp:369
void clearYRange()
Clear the y-range restriction.
Definition lsCompareNarrowBand.hpp:185
T getSumSquaredDifferences() const
Return the sum of squared differences calculated by apply().
Definition lsCompareNarrowBand.hpp:363
CompareNarrowBand(SmartPointer< Domain< T, D > > passedLevelSetTarget, SmartPointer< Domain< T, D > > passedlevelSetSample)
Definition lsCompareNarrowBand.hpp:146
void apply()
Apply the comparison and calculate the sum of squared differences.
Definition lsCompareNarrowBand.hpp:205
void clearXRange()
Clear the x-range restriction.
Definition lsCompareNarrowBand.hpp:178
T getSumDifferences() const
Definition lsCompareNarrowBand.hpp:366
CompareNarrowBand()
Definition lsCompareNarrowBand.hpp:140
Class containing all information about the level set, including the dimensions of the domain,...
Definition lsDomain.hpp:27
Expands the levelSet to the specified number of layers. The largest value in the levelset is thus wid...
Definition lsExpand.hpp:17
void apply()
Apply the expansion to the specified width.
Definition lsExpand.hpp:44
This class holds an explicit mesh, which is always given in 3 dimensions. If it describes a 2D mesh,...
Definition lsMesh.hpp:21
#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