208 if (!checkAndCalculateBounds()) {
210 sumSquaredDifferences = std::numeric_limits<T>::quiet_NaN();
215 const auto &gridTarget = levelSetTarget->getGrid();
216 double gridDelta = gridTarget.getGridDelta();
219 viennahrle::ConstDenseCellIterator<typename Domain<T, D>::DomainType>
220 itSample(levelSetSample->getDomain(), minIndex);
221 viennahrle::ConstDenseCellIterator<typename Domain<T, D>::DomainType>
222 itTarget(levelSetTarget->getDomain(), minIndex);
224 sumSquaredDifferences = 0.0;
228 std::unordered_map<viennahrle::Index<D>, size_t,
229 typename viennahrle::Index<D>::hash>
231 std::vector<T> differenceValues;
232 size_t currentPointId = 0;
234 const bool generateMesh = outputMesh !=
nullptr;
239 for (
unsigned i = 0; i <
D; ++i) {
240 outputMesh->minimumExtent[i] = std::numeric_limits<T>::max();
241 outputMesh->maximumExtent[i] = std::numeric_limits<T>::lowest();
246 for (; itSample.getIndices() < maxIndex; itSample.next()) {
248 T xCoord = itSample.getIndices()[0] * gridDelta;
249 T yCoord = (
D > 1) ? itSample.getIndices()[1] * gridDelta : 0;
252 if (useXRange && (xCoord < xRangeMin || xCoord > xRangeMax)) {
257 if (
D > 1 && useYRange && (yCoord < yRangeMin || yCoord > yRangeMax)) {
262 itTarget.goToIndicesSequential(itSample.getIndices());
269 for (
int i = 0; i < (1 <<
D); ++i) {
270 valueSample += itSample.getCorner(i).getValue();
271 valueTarget += itTarget.getCorner(i).getValue();
273 valueTarget /= (1 <<
D);
274 valueSample /= (1 <<
D);
277 if (std::isinf(valueTarget) || std::isinf(valueSample) ||
278 std::abs(valueTarget) > 1000 || std::abs(valueSample) > 1000) {
283 T diff = std::abs(valueTarget - valueSample) * gridDelta;
284 T diffSquared = diff * diff;
285 sumSquaredDifferences += diffSquared;
286 sumDifferences += diff;
291 std::array<unsigned, 1 << D> voxel;
292 bool addVoxel =
true;
295 for (
unsigned i = 0; i < (1 <<
D); ++i) {
296 viennahrle::Index<D> index;
297 for (
unsigned j = 0; j <
D; ++j) {
299 itSample.getIndices(j) + itSample.getCorner(i).getOffset()[j];
300 if (index[j] > maxIndex[j]) {
306 auto pointIdValue = std::make_pair(index, currentPointId);
307 auto pointIdPair = pointIdMapping.insert(pointIdValue);
308 voxel[i] = pointIdPair.first->second;
309 if (pointIdPair.second) {
318 if constexpr (
D == 3) {
319 std::array<unsigned, 8> hexa{voxel[0], voxel[1], voxel[3],
320 voxel[2], voxel[4], voxel[5],
322 outputMesh->hexas.push_back(hexa);
323 }
else if constexpr (
D == 2) {
324 std::array<unsigned, 4> quad{voxel[0], voxel[1], voxel[3],
326 outputMesh->tetras.push_back(quad);
331 differenceValues.push_back(outputMeshSquaredDifferences ? diffSquared
338 if (generateMesh && !pointIdMapping.empty()) {
339 outputMesh->nodes.resize(pointIdMapping.size());
340 for (
auto it = pointIdMapping.begin(); it != pointIdMapping.end(); ++it) {
341 std::array<T, 3> coords{};
342 for (
unsigned i = 0; i <
D; ++i) {
343 coords[i] = gridDelta * it->first[i];
345 if (coords[i] < outputMesh->minimumExtent[i]) {
346 outputMesh->minimumExtent[i] = coords[i];
347 }
else if (coords[i] > outputMesh->maximumExtent[i]) {
348 outputMesh->maximumExtent[i] = coords[i];
351 outputMesh->nodes[it->second] = coords;
354 assert(differenceValues.size() ==
355 outputMesh->template getElements<1 <<
D>().size());
356 outputMesh->cellData.insertNextScalarData(std::move(differenceValues),
357 outputMeshSquaredDifferences
358 ?
"Squared differences"
359 :
"Absolute differences");