207 if (!checkAndCalculateBounds()) {
209 sumSquaredDifferences = std::numeric_limits<T>::quiet_NaN();
214 const auto &gridTarget = levelSetTarget->getGrid();
215 double gridDelta = gridTarget.getGridDelta();
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);
223 sumSquaredDifferences = 0.0;
227 std::unordered_map<viennahrle::Index<D>, size_t,
228 typename viennahrle::Index<D>::hash>
230 std::vector<T> differenceValues;
231 size_t currentPointId = 0;
233 const bool generateMesh = outputMesh !=
nullptr;
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();
245 for (; itSample.getIndices() < maxIndex; itSample.next()) {
247 T xCoord = itSample.getIndices()[0] * gridDelta;
248 T yCoord = (
D > 1) ? itSample.getIndices()[1] * gridDelta : 0;
251 if (useXRange && (xCoord < xRangeMin || xCoord > xRangeMax)) {
256 if (
D > 1 && useYRange && (yCoord < yRangeMin || yCoord > yRangeMax)) {
261 itTarget.goToIndicesSequential(itSample.getIndices());
268 for (
int i = 0; i < (1 <<
D); ++i) {
269 valueSample += itSample.getCorner(i).getValue();
270 valueTarget += itTarget.getCorner(i).getValue();
272 valueTarget /= (1 <<
D);
273 valueSample /= (1 <<
D);
276 if (std::isinf(valueTarget) || std::isinf(valueSample) ||
277 std::abs(valueTarget) > 1000 || std::abs(valueSample) > 1000) {
282 T diff = std::abs(valueTarget - valueSample) * gridDelta;
283 T diffSquared = diff * diff;
284 sumSquaredDifferences += diffSquared;
285 sumDifferences += diff;
290 std::array<unsigned, 1 << D> voxel;
291 bool addVoxel =
true;
294 for (
unsigned i = 0; i < (1 <<
D); ++i) {
295 viennahrle::Index<D> index;
296 for (
unsigned j = 0; j <
D; ++j) {
298 itSample.getIndices(j) + itSample.getCorner(i).getOffset()[j];
299 if (index[j] > maxIndex[j]) {
305 auto pointIdValue = std::make_pair(index, currentPointId);
306 auto pointIdPair = pointIdMapping.insert(pointIdValue);
307 voxel[i] = pointIdPair.first->second;
308 if (pointIdPair.second) {
317 if constexpr (
D == 3) {
318 std::array<unsigned, 8> hexa{voxel[0], voxel[1], voxel[3],
319 voxel[2], voxel[4], voxel[5],
321 outputMesh->hexas.push_back(hexa);
322 }
else if constexpr (
D == 2) {
323 std::array<unsigned, 4> quad{voxel[0], voxel[1], voxel[3],
325 outputMesh->tetras.push_back(quad);
330 differenceValues.push_back(outputMeshSquaredDifferences ? diffSquared
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];
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];
350 outputMesh->nodes[it->second] = coords;
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");