156 std::vector<std::vector<std::vector<IndexType>>> componentList(
157 domain->getNumberOfSegments());
158 for (
unsigned segmentId = 0; segmentId < domain->getNumberOfSegments();
160 componentList[segmentId].resize(
D + 1);
161 for (
int dim = -1; dim <
D; ++dim) {
162 componentList[segmentId][dim + 1].resize(
163 domain->getDomain().getNumberOfRuns(segmentId, dim), -1);
167 std::size_t numberOfComponents = 0;
168 std::vector<IndexType> pointsPerComponent;
173 neighborIt(domain->getDomain());
174 !neighborIt.isFinished(); neighborIt.next()) {
175 auto ¢er = neighborIt.getCenter();
178 auto ¤tComponentId =
179 componentList[center.getSegmentId()][center.getLevel()]
180 [center.getRunTypePosition()];
183 if (currentComponentId == -1) {
184 for (
int k = 0; k < 2 *
D; ++k) {
185 auto &neighbor = neighborIt.getNeighbor(k);
186 const auto neighborComponentId =
187 componentList[neighbor.getSegmentId()][neighbor.getLevel()]
188 [neighbor.getRunTypePosition()];
192 if (neighborComponentId != -1) {
193 if (areConnected(center.getValue(), neighbor.getValue())) {
194 currentComponentId = neighborComponentId;
195 if (center.getValue() >= 0.)
196 ++pointsPerComponent[currentComponentId];
204 if (currentComponentId == -1) {
205 currentComponentId = numberOfComponents;
206 pointsPerComponent.push_back((center.getValue() > 0.) ? 1 : 0);
208 ++numberOfComponents;
212 for (
int k = 0; k < 2 *
D; ++k) {
213 auto &neighbor = neighborIt.getNeighbor(k);
214 auto &neighborComponentId =
215 componentList[neighbor.getSegmentId()][neighbor.getLevel()]
216 [neighbor.getRunTypePosition()];
217 if (areConnected(center.getValue(), neighbor.getValue())) {
219 if (neighborComponentId != -1) {
221 if (currentComponentId != neighborComponentId) {
225 neighborComponentId = currentComponentId;
226 if (neighbor.getValue() >= 0.)
227 ++pointsPerComponent[neighborComponentId];
236 (reverseVoidDetection) ? components[0] : components.back();
240 auto pointsPerConnected =
241 mergeComponentCounts(components, pointsPerComponent);
244 if (detectLargestSurface) {
245 topComponent = calculateTopID(components, pointsPerConnected);
247 while (pointsPerConnected[topComponent] == 0) {
248 if (reverseVoidDetection)
256 std::vector<T> voidPointMarkers;
257 voidPointMarkers.resize(domain->getNumberOfPoints());
259 std::vector<T> componentMarkers;
261 componentMarkers.resize(domain->getNumberOfPoints());
266 neighborIt(domain->getDomain());
267 !neighborIt.isFinished(); neighborIt.next()) {
268 auto center = neighborIt.getCenter();
270 if (!center.isDefined())
274 if (center.getValue() >= 0) {
275 const int &oldComponentId = componentList[center.getSegmentId()][0]
276 [center.getRunTypePosition()];
277 voidPointMarkers[center.getPointId()] =
278 (components[oldComponentId] != topComponent);
284 for (k = 0; k < 2 *
D; ++k) {
285 auto &neighbor = neighborIt.getNeighbor(k);
286 if (std::signbit(neighbor.getValue()) ==
287 std::signbit(center.getValue()))
289 const int &oldneighborComponentId =
290 componentList[neighbor.getSegmentId()][neighbor.getLevel()]
291 [neighbor.getRunTypePosition()];
292 if (components[oldneighborComponentId] == topComponent) {
296 voidPointMarkers[center.getPointId()] = (k == 2 *
D);
299 if (saveComponents) {
300 const int &oldComponentId = componentList[center.getSegmentId()][0]
301 [center.getRunTypePosition()];
302 componentMarkers[center.getPointId()] = components[oldComponentId];
306 auto &pointData = domain->getPointData();
307 auto voidMarkersPointer = pointData.getScalarData(
voidPointLabel,
true);
309 if (voidMarkersPointer ==
nullptr) {
312 *voidMarkersPointer = std::move(voidPointMarkers);
315 if (saveComponents) {
316 auto componentMarkersPointer =
317 pointData.getScalarData(
"ConnectedComponentId",
true);
319 if (componentMarkersPointer ==
nullptr) {
320 pointData.insertNextScalarData(componentMarkers,
321 "ConnectedComponentId");
323 *componentMarkersPointer = std::move(componentMarkers);