180 if (!checkAndCalculateBounds()) {
182 sumSquaredDifferences = std::numeric_limits<T>::quiet_NaN();
183 sumDifferences = std::numeric_limits<T>::quiet_NaN();
188 const auto &gridTarget = levelSetTarget->getGrid();
189 const auto gridDelta = gridTarget.getGridDelta();
191 sumSquaredDifferences = 0.0;
192 sumDifferences = 0.0;
196 std::vector<Vec3D<T>> nodeCoordinates;
197 std::vector<std::array<unsigned, 1>> vertexIndices;
198 std::vector<T> differenceValues;
199 std::vector<T> squaredDifferenceValues;
202 const bool generateMesh = outputMesh !=
nullptr;
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();
213 nodeCoordinates.reserve(levelSetSample->getNumberOfPoints());
214 vertexIndices.reserve(levelSetSample->getNumberOfPoints());
215 differenceValues.reserve(levelSetSample->getNumberOfPoints());
216 squaredDifferenceValues.reserve(levelSetSample->getNumberOfPoints());
220 viennahrle::ConstSparseIterator<typename Domain<T, D>::DomainType> itSample(
221 levelSetSample->getDomain());
222 viennahrle::ConstSparseIterator<typename Domain<T, D>::DomainType> itTarget(
223 levelSetTarget->getDomain());
226 while (!itSample.isFinished()) {
227 if (!itSample.isDefined()) {
234 auto indices = itSample.getStartIndices();
237 T xCoord = indices[0] * gridDelta;
238 T yCoord = indices[1] * gridDelta;
242 if (useXRange && (xCoord < xRangeMin || xCoord > xRangeMax)) {
248 if (useYRange && (yCoord < yRangeMin || yCoord > yRangeMax)) {
254 T valueSample = itSample.getValue();
256 itTarget.goToIndicesSequential(indices);
257 T valueTarget = itTarget.getValue();
260 if (!itTarget.isDefined() || std::isinf(valueTarget) ||
261 std::isinf(valueSample)) {
267 T diff = std::abs(valueTarget - valueSample) * gridDelta;
268 T diffSquared = diff * diff;
269 sumDifferences += diff;
270 sumSquaredDifferences += diffSquared;
276 Vec3D<T> coords{xCoord, yCoord, zCoord};
279 nodeCoordinates.push_back(coords);
282 differenceValues.push_back(diff);
283 squaredDifferenceValues.push_back(diffSquared);
286 std::array<unsigned, 1> vertex = {
287 static_cast<unsigned>(nodeCoordinates.size() - 1)};
288 vertexIndices.push_back(vertex);
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]);
304 if (generateMesh && !nodeCoordinates.empty()) {
306 outputMesh->nodes = std::move(nodeCoordinates);
309 outputMesh->vertices = std::move(vertexIndices);
312 outputMesh->pointData.insertNextScalarData(std::move(differenceValues),
313 "Absolute differences");
314 outputMesh->pointData.insertNextScalarData(
315 std::move(squaredDifferenceValues),
"Squared differences");