184 if (!checkAndCalculateBounds()) {
185 differentCellsCount =
186 std::numeric_limits<unsigned long int>::max();
187 customDifferentCellCount =
188 std::numeric_limits<unsigned long int>::max();
193 viennahrle::ConstDenseCellIterator<hrleDomainType> itTarget(
194 levelSetTarget->getDomain(), minIndex);
195 viennahrle::ConstDenseCellIterator<hrleDomainType> itSample(
196 levelSetSample->getDomain(), minIndex);
198 differentCellsCount = 0;
199 customDifferentCellCount = 0;
202 const bool generateMesh = outputMesh !=
nullptr;
206 for (
unsigned i = 0; i <
D; ++i) {
207 outputMesh->minimumExtent[i] = std::numeric_limits<T>::max();
208 outputMesh->maximumExtent[i] = std::numeric_limits<T>::lowest();
214 std::vector<T> cellDifference;
215 std::vector<T> incrementValues;
216 size_t currentPointId = 0;
217 std::unordered_map<viennahrle::Index<D>, size_t,
218 typename viennahrle::Index<D>::hash>
222 for (; itTarget.getIndices() < maxIndex; itTarget.next()) {
223 itSample.goToIndicesSequential(itTarget.getIndices());
226 T centerValueTarget = 0.;
227 T centerValueSample = 0.;
230 for (
int i = 0; i < (1 <<
D); ++i) {
231 centerValueTarget += itTarget.getCorner(i).getValue();
232 centerValueSample += itSample.getCorner(i).getValue();
234 centerValueTarget /= (1 <<
D);
235 centerValueSample /= (1 <<
D);
238 bool insideTarget = (centerValueTarget <= 0.);
239 bool insideSample = (centerValueSample <= 0.);
240 bool isDifferent = insideTarget != insideSample;
243 bool inXRange = useCustomXIncrement &&
244 (itTarget.getIndices()[0] * gridDelta >= xRangeMin &&
245 itTarget.getIndices()[0] * gridDelta <= xRangeMax);
246 bool inYRange = useCustomYIncrement &&
247 (itTarget.getIndices()[1] * gridDelta >= yRangeMin &&
248 itTarget.getIndices()[1] * gridDelta <= yRangeMax);
251 unsigned short int incrementToAdd = defaultIncrement;
252 if (inXRange && inYRange) {
254 incrementToAdd = customXIncrement + customYIncrement;
255 }
else if (inXRange) {
256 incrementToAdd = customXIncrement;
257 }
else if (inYRange) {
258 incrementToAdd = customYIncrement;
264 differentCellsCount += 1;
267 customDifferentCellCount += incrementToAdd;
272 if (generateMesh && (insideTarget || insideSample)) {
274 int difference = isDifferent ? 1 : 0;
276 std::array<unsigned, 1 << D> voxel;
277 bool addVoxel =
true;
281 for (
unsigned i = 0; i < (1 <<
D); ++i) {
282 viennahrle::Index<D> index;
283 for (
unsigned j = 0; j <
D; ++j) {
285 itTarget.getIndices(j) + itTarget.getCorner(i).getOffset()[j];
286 if (index[j] > maxIndex[j]) {
292 auto pointIdValue = std::make_pair(index, currentPointId);
293 auto pointIdPair = pointIdMapping.insert(pointIdValue);
294 voxel[i] = pointIdPair.first->second;
295 if (pointIdPair.second) {
304 if constexpr (
D == 3) {
305 std::array<unsigned, 8> hexa{voxel[0], voxel[1], voxel[3],
306 voxel[2], voxel[4], voxel[5],
308 outputMesh->hexas.push_back(hexa);
309 }
else if constexpr (
D == 2) {
310 std::array<unsigned, 4> quad{voxel[0], voxel[1], voxel[3],
312 outputMesh->tetras.push_back(quad);
315 cellDifference.push_back(difference);
318 incrementValues.push_back(isDifferent ?
static_cast<T>(incrementToAdd)
325 if (generateMesh && !pointIdMapping.empty()) {
327 outputMesh->nodes.resize(pointIdMapping.size());
328 for (
auto it = pointIdMapping.begin(); it != pointIdMapping.end(); ++it) {
330 for (
unsigned i = 0; i <
D; ++i) {
331 coords[i] = gridDelta * it->first[i];
333 if (coords[i] < outputMesh->minimumExtent[i]) {
334 outputMesh->minimumExtent[i] = coords[i];
335 }
else if (coords[i] > outputMesh->maximumExtent[i]) {
336 outputMesh->maximumExtent[i] = coords[i];
339 outputMesh->nodes[it->second] = coords;
343 assert(cellDifference.size() == incrementValues.size());
344 assert(incrementValues.size() ==
345 outputMesh->template getElements<1 <<
D>().size());
346 outputMesh->cellData.insertNextScalarData(cellDifference,
"Difference");
347 outputMesh->cellData.insertNextScalarData(incrementValues,