Export mesh instead of polygon soup

[raphaelsulzer]

Hi and thanks for making this nice work available!

I would like to use your library to export the cell complex and final surface. If I am not mistaken both can only be exported as polygon soups, ie with each face being a single component and many duplicate vertices.

Do you know if there is an easy way to access common vertices between adjacent cells and even adjacent facets through sage?

Kind regards

[chenzhaiyu]

Hi @raphaelsulzer, you're right that currently polygon soups are exported.

I didn't manage to locate common vertices across cells in the first implementation. That was difficult with native Sage functions (e.g., with Polyhedron.intersect(Polyhedron) we lose track of the vertex indices as in the original cells). Common vertices between adjacent facets shouldn't be a concern as the Polyhedron class already handles it well - now they're more of a by-product of those between cells. I think it would nevertheless be possible to record the vertices globally and sort out the common vertices only when exporting, which would involve some out-of-the-box implementations that are much slower. This is also a reminder for me to take a fresh look at what Sage offers at this point. In the meantime, if you find anything useful, please share it here or shoot a PR :)

[chenzhaiyu]

Also as a quick note: I'm not sure to what extent these functions in Mapple could help you as a post-processing, in case an immediate fix is needed. CC @LiangliangNan.

[raphaelsulzer]

Thanks for the reply. I really have no experience with sage. I guess finding adjacencies has to be done via the graph, because the cell complex is only a list of independent cells.

I implemented some quick fix for getting a surface, by simply collecting all vertices, finding unique once, and linking them to the interface facets. Unfortunately, I cannot manage to orient the facets by forcing the orientation to be towards the reachable cells. Not sure why the orientation is still not correct, but the rest works.

You can find the function here:
https://github.com/raphaelsulzer/abspy/blob/6ca3a5ad66ba85e9fa32865ffcb425214f1c9df2/abspy/graph.py#L433

[chenzhaiyu]

Thanks for sharing your solution here! For re-orientation, I have written some lines before (with rough edges) which may be integrated as a fix: https://gist.github.com/chenzhaiyu/179fde5aedd1f909bb6ab8145f4366d0

As mentioned before it would be optimal to include the info on common vertices already during partitioning. I will look into this when having time, otherwise if impossible there, your solution will be the workaround that we need to take :)

[raphaelsulzer]

In fact, there seems to be a bigger problem with the cell complex.
When running the solution I present above I cannot manage to merge all the duplicate vertices. Even with rational number types from SAGE some vertices that should probably be the same are not the same in practice. I attached an example output surface from abspy that has two disconnected components. Maybe the problem lies already somewhere in the cell complex construct method?

anchor.zip

[chenzhaiyu]

Sage itself should handle well its rational number types in the box - I have never bumped into any arithmetically degenerate cases using pure Sage. However, for operations outside Sage, whether the exactness is retained really depends on implementations. I guess most operators and predicates in Numpy, for example the ones in your function (https://github.com/raphaelsulzer/abspy/blob/6ca3a5ad66ba85e9fa32865ffcb425214f1c9df2/abspy/graph.py#L485 and https://github.com/raphaelsulzer/abspy/blob/6ca3a5ad66ba85e9fa32865ffcb425214f1c9df2/abspy/graph.py#L497), might not support rational numbers. I have no idea how Numpy internally coerces Sage's rational numbers, but I think there should be a breaking point. There are other libs (including the native fractions) that you could use for exact operations out of the box.

[raphaelsulzer]

For the example I attached above I kept the exact type from SAGE by using dtype=object in the numpy arrays.
The function for producing that output is here: https://github.com/raphaelsulzer/abspy/blob/64d41e3c2d541d0c648da54afccf235831ee4947/abspy/graph.py#L222

[chenzhaiyu]

https://github.com/raphaelsulzer/abspy/blob/64d41e3c2d541d0c648da54afccf235831ee4947/abspy/graph.py#L257
It was difficult to read this line undocumented but are you sure this is the expected behavior with np.isin()?

>>> pset = np.array([(1, 2, 3), (4, 5, 6)])
>>> p = np.array([1, 3, 2])
>>> np.isin(pset, p)
array([[ True,  True,  True],
       [False, False, False]])

In addition, as I already mentioned, please make sure inside all the Numpy functions you used there are no type coercing happening, though you specified dtype=object at some point.

[raphaelsulzer]

You are right that it is probably better to write pset == p to respect the order of dimensions (or even np.equal(pset, p, dtype=object) for that matter), but it doesn't solve the problem.

I still believe the issue lies somewhere else. With exhaustive partioning, I can export a watertight surface. I compare below adaptive and exhaustive cell complex at one of the vertices in question. Should't the adaptive partition be included in the exhaustive partition, ie the two lines close to the vertex should not exist?

Non-watertight surface. Component 1 grey, component 2 red.

Adaptive cell complex.

Exhaustive cell complex.

[chenzhaiyu]

Now clear and interesting finding! By exhaustive do you mean Sage's native hyperplane arrangement or simply switching on the following exhaustive option?

abspy/abspy/complex.py

Line 419 in 307d857

def construct(self, exhaustive=False, num_workers=0):

[raphaelsulzer]

switching on your exhaustive option.

[chenzhaiyu]

Sorry for the very late reply. I couldn't produce your result but I checked your results anchored here:

In fact, there seems to be a bigger problem with the cell complex... Maybe the problem lies already somewhere in the cell complex construct method?
anchor.zip

I can't find the discrepancies you mentioned:

Should't the adaptive partition be included in the exhaustive partition, ie the two lines close to the vertex should not exist?
Non-watertight surface. Component 1 grey, component 2 red.

Adaptive cell complex.

Exhaustive cell complex.

From my end, the two lines do not exist and the adaptive partition is indeed included in the exhaustive partition:

Adaptive

Exhaustive

[raphaelsulzer]

You have to zoom in more on the lower vertex. You can also verify it in Meshlab by right clicking on the mesh and then "Split in Connected Components".

[chenzhaiyu]

I have currently no better way to dig deeper on this but one quick idea: the "components" might also be the result of non-manifoldness. For example, if two polytopes touch each other at one vertex they are not regarded as adjacent in the graph; the duplicate vertices in this case are not properly removed within your code. You can easily change this behavior here:

abspy/abspy/complex.py

Line 395 in 6eb26ac

def _intersect_neighbour(self, kwargs):