ViennaLS
Loading...
Searching...
No Matches
lsCompareSparseField.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <hrleSparseIterator.hpp>
4#include <lsDomain.hpp>
5#include <lsExpand.hpp>
6#include <lsMesh.hpp>
8#include <lsReduce.hpp>
9
10#include <cmath>
11#include <limits>
12
13namespace viennals {
14
15using namespace viennacore;
16
23template <class T, int D = 2> class CompareSparseField {
24 using hrleIndexType = viennahrle::IndexType;
25
26 SmartPointer<Domain<T, D>> levelSetTarget = nullptr;
27 SmartPointer<Domain<T, D>> levelSetSample = nullptr;
28
29 // Variables for x and y range restrictions
30 T xRangeMin = std::numeric_limits<T>::lowest();
31 T xRangeMax = std::numeric_limits<T>::max();
32 T yRangeMin = std::numeric_limits<T>::lowest();
33 T yRangeMax = std::numeric_limits<T>::max();
34 bool useXRange = false;
35 bool useYRange = false;
36
37 // Fields to store the calculation results
38 T sumSquaredDifferences = 0.0;
39 T sumDifferences = 0.0;
40 unsigned numPoints = 0;
41
42 // Add mesh output capability
43 SmartPointer<Mesh<T>> outputMesh = nullptr;
44
45 bool checkAndCalculateBounds() {
46 if (levelSetTarget == nullptr || levelSetSample == nullptr) {
47 Logger::getInstance()
48 .addWarning("Missing level set in CompareSparseField.")
49 .print();
50 return false;
51 }
52
53 // Check if the grids are compatible
54 const auto &gridTarget = levelSetTarget->getGrid();
55 const auto &gridSample = levelSetSample->getGrid();
56
57 if (gridTarget.getGridDelta() != gridSample.getGridDelta()) {
58 Logger::getInstance()
59 .addWarning("Grid delta mismatch in CompareSparseField. The grid "
60 "deltas of the two level sets must be equal.")
61 .print();
62 return false;
63 }
64
65 // // Check if the x extents of both level sets are equal
66 // const auto &domainTarget = levelSetTarget->getDomain();
67 // const auto &domainSample = levelSetSample->getDomain();
68
69 // hrleIndexType targetMinX = gridTarget.isNegBoundaryInfinite(0)
70 // ? domainTarget.getMinRunBreak(0)
71 // : gridTarget.getMinIndex(0);
72 // hrleIndexType targetMaxX = gridTarget.isPosBoundaryInfinite(0)
73 // ? domainTarget.getMaxRunBreak(0)
74 // : gridTarget.getMaxIndex(0);
75 // hrleIndexType sampleMinX = gridSample.isNegBoundaryInfinite(0)
76 // ? domainSample.getMinRunBreak(0)
77 // : gridSample.getMinIndex(0);
78 // hrleIndexType sampleMaxX = gridSample.isPosBoundaryInfinite(0)
79 // ? domainSample.getMaxRunBreak(0)
80 // : gridSample.getMaxIndex(0);
81
82 // if (targetMinX != sampleMinX || targetMaxX != sampleMaxX) {
83 // Logger::getInstance()
84 // .addWarning("X extent mismatch in CompareSparseField. The x extents
85 // "
86 // "of both level sets must be equal.")
87 // .print();
88 // return false;
89 // }
90
91 // Check if target level set width is sufficient
92 if (levelSetTarget->getLevelSetWidth() < 50) {
93 Logger::getInstance()
94 .addWarning(
95 "Target level set width is insufficient. It must have a width of "
96 "at least 50. \n"
97 " CORRECTION: The expansion was performed. \n"
98 "ALTERNATIVE: Alternatively, please expand the target yourself "
99 "using lsExpand before passing it to this function. \n")
100 .print();
101 Expand<T, D>(levelSetTarget, 50).apply();
102 }
103
104 // Reduce the sample level set to a sparse field if necessary
105 if (levelSetSample->getLevelSetWidth() > 1) {
106 Logger::getInstance()
107 .addWarning(
108 "Sample level set width is too large. It must be reduced to a "
109 "sparse field. \n"
110 " CORRECTION: The reduction was performed. \n"
111 "ALTERNATIVE: Alternatively, please reduce the sample yourself "
112 "using lsReduce before passing it to this function. \n")
113 .print();
114 Reduce<T, D>(levelSetSample, 1).apply();
115 }
116
117 return true;
118 }
119
120public:
122 static_assert(
123 D == 2 &&
124 "CompareSparseField is currently only implemented for 2D level sets.");
125 }
126
127 CompareSparseField(SmartPointer<Domain<T, D>> passedLevelSetTarget,
128 SmartPointer<Domain<T, D>> passedLevelSetSample)
129 : levelSetTarget(passedLevelSetTarget),
130 levelSetSample(passedLevelSetSample) {
131 static_assert(
132 D == 2 &&
133 "CompareSparseField is currently only implemented for 2D level sets.");
134 }
135
136 void setLevelSetTarget(SmartPointer<Domain<T, D>> passedLevelSet) {
137 levelSetTarget = passedLevelSet;
138 }
139
140 void setLevelSetSample(SmartPointer<Domain<T, D>> passedLevelSet) {
141 levelSetSample = passedLevelSet;
142 }
143
145 void setXRange(T minXRange, T maxXRange) {
146 xRangeMin = minXRange;
147 xRangeMax = maxXRange;
148 useXRange = true;
149 }
150
152 void setYRange(T minYRange, T maxYRange) {
153 yRangeMin = minYRange;
154 yRangeMax = maxYRange;
155 useYRange = true;
156 }
157
159 void clearXRange() {
160 useXRange = false;
161 xRangeMin = std::numeric_limits<T>::lowest();
162 xRangeMax = std::numeric_limits<T>::max();
163 }
164
166 void clearYRange() {
167 useYRange = false;
168 yRangeMin = std::numeric_limits<T>::lowest();
169 yRangeMax = std::numeric_limits<T>::max();
170 }
171
173 void setOutputMesh(SmartPointer<Mesh<T>> passedMesh) {
174 outputMesh = passedMesh;
175 }
176
178 void apply() {
179 // Perform compatibility checks
180 if (!checkAndCalculateBounds()) {
181 // If checks fail, return NaN
182 sumSquaredDifferences = std::numeric_limits<T>::quiet_NaN();
183 sumDifferences = std::numeric_limits<T>::quiet_NaN();
184 numPoints = 0;
185 return;
186 }
187
188 const auto &gridTarget = levelSetTarget->getGrid();
189 const auto gridDelta = gridTarget.getGridDelta();
190
191 sumSquaredDifferences = 0.0;
192 sumDifferences = 0.0;
193 numPoints = 0;
194
195 // Direct storage for points and differences
196 std::vector<Vec3D<T>> nodeCoordinates; // 3D necessary for lsMesh
197 std::vector<std::array<unsigned, 1>> vertexIndices;
198 std::vector<T> differenceValues;
199 std::vector<T> squaredDifferenceValues;
200
201 // Prepare mesh output if needed
202 const bool generateMesh = outputMesh != nullptr;
203 if (generateMesh) {
204 outputMesh->clear();
205
206 // Initialize mesh extent
207 for (unsigned i = 0; i < D; ++i) {
208 outputMesh->minimumExtent[i] = std::numeric_limits<T>::max();
209 outputMesh->maximumExtent[i] = std::numeric_limits<T>::lowest();
210 }
211
212 // Reserve space for mesh data
213 nodeCoordinates.reserve(levelSetSample->getNumberOfPoints());
214 vertexIndices.reserve(levelSetSample->getNumberOfPoints());
215 differenceValues.reserve(levelSetSample->getNumberOfPoints());
216 squaredDifferenceValues.reserve(levelSetSample->getNumberOfPoints());
217 }
218
219 // Create sparse iterators for the level sets
220 viennahrle::ConstSparseIterator<typename Domain<T, D>::DomainType> itSample(
221 levelSetSample->getDomain());
222 viennahrle::ConstSparseIterator<typename Domain<T, D>::DomainType> itTarget(
223 levelSetTarget->getDomain());
224
225 // Iterate over all defined points in the sample level set
226 while (!itSample.isFinished()) {
227 if (!itSample.isDefined()) {
228 // this block is necessary to skip undefined points, I tested it:
229 // std::cout << "Skipping undefined point" << std::endl;
230 itSample.next();
231 continue;
232 }
233
234 auto indices = itSample.getStartIndices();
235
236 // Calculate coordinates
237 T xCoord = indices[0] * gridDelta;
238 T yCoord = indices[1] * gridDelta;
239 T zCoord = 0.0; // Always use 0 for z-coordinate in 2D
240
241 // Skip if outside the specified x-range
242 if (useXRange && (xCoord < xRangeMin || xCoord > xRangeMax)) {
243 itSample.next();
244 continue;
245 }
246
247 // Skip if outside the specified y-range
248 if (useYRange && (yCoord < yRangeMin || yCoord > yRangeMax)) {
249 itSample.next();
250 continue;
251 }
252
253 // Get sample value
254 T valueSample = itSample.getValue();
255
256 itTarget.goToIndicesSequential(indices);
257 T valueTarget = itTarget.getValue();
258
259 // Check for infinite or extreme values that might cause numerical issues
260 if (!itTarget.isDefined() || std::isinf(valueTarget) ||
261 std::isinf(valueSample)) {
262 itSample.next();
263 continue;
264 }
265
266 // Calculate difference and add to sum
267 T diff = std::abs(valueTarget - valueSample) * gridDelta;
268 T diffSquared = diff * diff;
269 sumDifferences += diff;
270 sumSquaredDifferences += diffSquared;
271 numPoints++;
272
273 // Store difference in mesh if required
274 if (generateMesh) {
275 // Create a new point with the coordinates of the sample level set point
276 Vec3D<T> coords{xCoord, yCoord, zCoord}; // lsMesh needs 3D
277
278 // Store the coordinates
279 nodeCoordinates.push_back(coords);
280
281 // Store the difference value (squared and absolute)
282 differenceValues.push_back(diff);
283 squaredDifferenceValues.push_back(diffSquared);
284
285 // Create a vertex for this point
286 std::array<unsigned, 1> vertex = {
287 static_cast<unsigned>(nodeCoordinates.size() - 1)};
288 vertexIndices.push_back(vertex);
289
290 // Update the mesh extent
291 for (unsigned i = 0; i < D; ++i) {
292 outputMesh->minimumExtent[i] =
293 std::min(outputMesh->minimumExtent[i], coords[i]);
294 outputMesh->maximumExtent[i] =
295 std::max(outputMesh->maximumExtent[i], coords[i]);
296 }
297 }
298
299 // Move to next point
300 itSample.next();
301 }
302
303 // Finalize mesh output
304 if (generateMesh && !nodeCoordinates.empty()) {
305 // Store the node coordinates directly
306 outputMesh->nodes = std::move(nodeCoordinates);
307
308 // Store the vertices for point visualization
309 outputMesh->vertices = std::move(vertexIndices);
310
311 // Store the difference values as point data
312 outputMesh->pointData.insertNextScalarData(std::move(differenceValues),
313 "Absolute differences");
314 outputMesh->pointData.insertNextScalarData(
315 std::move(squaredDifferenceValues), "Squared differences");
316 }
317 }
318
320 T getSumSquaredDifferences() const { return sumSquaredDifferences; }
321
323 T getSumDifferences() const { return sumDifferences; }
324
326 unsigned getNumPoints() const { return numPoints; }
327
329 T getRMSE() const {
330 return (numPoints > 0) ? std::sqrt(sumSquaredDifferences / numPoints)
331 : std::numeric_limits<T>::infinity();
332 }
333};
334
335// Add template specializations for this class
336PRECOMPILE_PRECISION_DIMENSION(CompareSparseField)
337
338} // namespace viennals
unsigned getNumPoints() const
Return the number of points used in the comparison.
Definition lsCompareSparseField.hpp:326
T getRMSE() const
Calculate the root mean square error from previously computed values.
Definition lsCompareSparseField.hpp:329
CompareSparseField()
Definition lsCompareSparseField.hpp:121
void setLevelSetSample(SmartPointer< Domain< T, D > > passedLevelSet)
Definition lsCompareSparseField.hpp:140
void setXRange(T minXRange, T maxXRange)
Set the x-coordinate range to restrict the comparison area.
Definition lsCompareSparseField.hpp:145
void setOutputMesh(SmartPointer< Mesh< T > > passedMesh)
Set the output mesh where difference values will be stored.
Definition lsCompareSparseField.hpp:173
void setYRange(T minYRange, T maxYRange)
Set the y-coordinate range to restrict the comparison area.
Definition lsCompareSparseField.hpp:152
CompareSparseField(SmartPointer< Domain< T, D > > passedLevelSetTarget, SmartPointer< Domain< T, D > > passedLevelSetSample)
Definition lsCompareSparseField.hpp:127
void setLevelSetTarget(SmartPointer< Domain< T, D > > passedLevelSet)
Definition lsCompareSparseField.hpp:136
T getSumSquaredDifferences() const
Return the sum of squared differences calculated by apply().
Definition lsCompareSparseField.hpp:320
T getSumDifferences() const
Return the sum of differences calculated by apply().
Definition lsCompareSparseField.hpp:323
void clearYRange()
Clear the y-range restriction.
Definition lsCompareSparseField.hpp:166
void apply()
Apply the comparison and calculate the sum of squared differences.
Definition lsCompareSparseField.hpp:178
void clearXRange()
Clear the x-range restriction.
Definition lsCompareSparseField.hpp:159
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
Reduce the level set size to the specified width. This means all level set points with value <= 0....
Definition lsReduce.hpp:14
void apply()
Reduces the leveleSet to the specified number of layers. The largest value in the levelset is thus wi...
Definition lsReduce.hpp:55
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
Definition lsAdvect.hpp:36
constexpr int D
Definition pyWrap.cpp:71
double T
Definition pyWrap.cpp:69