Bugs in first version of bpypolyskel

[polarkernel]

I just have commited a first version of the bpypolyskel modules. The demo-module requires the mathutils library to be installed in your interpreter using 'pip install mathutils'. A demo running in Blender will follow later.

Please report and discuss eventual bugs or otpimizations here.

[vvoovv]

A quite simple polygon. What's wrong with it?

bpypolyskel.py.txt

[polarkernel]

A quite simple polygon. What's wrong with it?

Oh no !! This issue is produced by me. I implemented polygonize() expecting a clockwise order of vertices and not counter-clockwise, as defined. I will fix that as soon as possible, sorry.

[vvoovv]

Yes, the order of vertices is counter-clockwise for the outer ring
The order of vertices is clockwise for the holes.

[polarkernel]

I will fix that as soon as possible,

I am sorry, but this will take a while. The only proper way is to change the vertex order in Botffy's code. In the original paper, the vertex order is anti-clockwise in a coordinate system as we like to use it. But I assume that he used the coordinate system of the images he introduced for debugging. In many image libraries, the y-axis goes from top to bottom.

Reversing vertices lists and indices in polygonize() leads immediately to an unreadable code and I like to avoid this. I think it is better if I take the time to change skeletonize(). It will even help me to understand the algorithm for the straight skeleton a bit better, I assume.

[vvoovv]

Ok. No problem.

[polarkernel]

Yes, the order of vertices is counter-clockwise for the outer ring
The order of vertices is clockwise for the holes.

Fixed and commited. This was really an intensive exercise in computational geometry, but I hope it works now.

[vvoovv]

I was away for a few days.
I'll try it today or tomorrow.

[vvoovv]

Got a weird result.

(1) bpypolyskel_test1.py.txt

---

The same footprint, but imported relative to another reference point.
(2)
bpypolyskel_test2.py.txt

[polarkernel]

Got a weird result.

The same footprint, but imported relative to another reference point.

Fixed. Was damned diffcult to find. Finally, it required only a '>' instead of a '<'.

[vvoovv]

Great! Thank you so much!
I've fully integrated bpypolyskel with the addon. Hipped roofs with textures are now generated for polygons without a hole.

[polarkernel]

Fantastic! The buildings look really very natural. I hope, the package will work also for buildings with holes. I tested only one example of them.

[vvoovv]

I notices the faces are returned in the arbitrary order rather than follow the order of vertices in the verts input parameter,

Would it be possible to return the faces in the order of the verts without sorting? If not possible, no problem at all. I can cope with it.

[polarkernel]

notices the faces are returned in the arbitrary order rather than follow the order of vertices in the verts input parameter,

I can't find this issue. The vertices in faces should always start with an edge of the original polygon (or hole). Thus, to follow the order of vertices in the verts, their number should always increase, except for the edge from the last vertice to the first. If I take your example bpypolyskel_test1.py above and print faces, I get the following values:

[[10, 11, 21, 20], [8, 9, 17, 16], [15, 8, 16], [13, 14, 18, 19], [14, 15, 16, 17, 18], 
 [9, 10, 20, 19, 18, 17], [11, 12, 21], [12, 13, 19, 20, 21]]

As you can see, the numbers are always increasing, except for the last (15) to the first (8) vertice. In my example with the Blender addon in the repo, the face normals are all looking outside, as it should be:

Do you have an example that shows this issue?

[vvoovv]

So the vertices of the polygon have the indices [8, 9, 10, 11, 12, 13, 14, 15]

Under the same order I meant the following result:

[
[8, 9, 17, 16],
[9, 10, 20, 19, 18, 17],
[10, 11, 21, 20],
[11, 12, 21],
[12, 13, 19, 20, 21],
[13, 14, 18, 19],
[14, 15, 16, 17, 18],
[15, 8, 16]
]

Note the first vertex in each of the resulting face. It is the vertex of the original polygon for the same index.

If sorting is required to get that kind of result, then keep the resulting faces as they are now.

[polarkernel]

Oh my god. what a stupid misunderstanding. I wondered because the order of the vertices is given as counterclockwise, but I didn't get the self-evident meaning, sorry.

No, it is not possible to create this order without sorting. It gets porduced by the construction of a so called counterclockwise embedding in the graph class. It is based on a set of edges, every edge with both directions. A set in Python is by defintion unsorted, may be due to hash values used to identify the items in it. But I don't know the underlying algorithm. However, it makes finding faces in counterclockwise order very fast. It could be changed, but I think this would result in even more computation time than sorting.

[vvoovv]

Ok, thank you for the explanation. Let's preserve the current realisation.

[vvoovv]

Unit vectors as the input parameter don't work for me:
bpypolyskel_test.py.txt

[polarkernel]

Unit vectors as the input parameter don't work for me:

Fixed. Accidentally, a unit vector for all vertices in verts was expected , which was clearly wrong.

[vvoovv]

Fixed. Accidentally, a unit vector for all vertices in verts was expected , which was clearly wrong.

Thank you! It works now.

[vvoovv]

Is it possible to solve the problem of the extra edge created for a t-shaped polygon?
t_shape_extra_edge.py.txt

[polarkernel]

Is it possible to solve the problem of the extra edge created for a t-shaped polygon?

Difficult. I wouldn't touch skeltonize(), because its task is to deliver one straight skeleton and not two partial skeletons. Maybe some post-processing could solve this issue: Find adjacent coplanar faces and merge them. Could maybe Blender solve this? As far as I know, there exists a built-in operator

bpy.ops.mesh.select_similar(type='COPLANAR', threshold=???)

I don't know if it could be applied already to the mesh of polygon faces. Maybe Blender could then merge (dissolve egde(s)?) these. However, my skills in Blender are very limited, if this idea is not feasible, I could try to code some simple (=fast) post-processing, but in Python this could be slower than in Blender.

[polarkernel]

Is it possible to solve the problem of the extra edge created for a t-shaped polygon?

Wait, before you waist time on the Blender solution proposal. I got an idea how this could quite simply be solved in polygonize(), but I have to elaborate it. This will take some time.

[vvoovv]

I think that extra edge is erroneously delivered by skeletonize(..).

[polarkernel]

I know and try to debug. In the skeleton there is more than one erroneous edge. Even one as loop from a node to itself.

[polarkernel]

Is it possible to solve the problem of the extra edge created for a t-shaped polygon?

I give up!! I can't find a real bug in skeletonize(). It seems that the algorithm in the corresponding paper is realized without bugs. Maybe the original algorithm even didn't solve this issue. Whatever, I couldn't solve it. :-(

As an intermediate solution, I have introduced a function clean_skeleton(), which solves issues with extra edges at the end of skeltonize() by cleaning up the skeleton. Let me explain its idea using your example t_shape_extra_edge.py.txt:

One of the skeleton nodes in the result of skeletonize() is node (60,50). One of its sinks points to itself and forms a graph loop. Two other sinks, one is (40,50) and the other is (60,70), form two parallel skeleton edges, the first to another skeleton node and the second to a polygon vertex. Cleaning means now removing the loop, moving the sink (60,70) of the polygon vertex to the sinks list of the node (40,50), the first edge points to, and remove it from the sinks list of (60,50).

Like this, there are no more extra edges in the skeleton. Its not for free, but I think it works quite fast.

[vvoovv]

I tried t_shape_extra_edge.py.txt again. The extra edge was still there.

[polarkernel]

I tried t_shape_extra_edge.py.txt again. The extra edge was still there.

I removed two extra edges and you are still not happy?

Today, I studied some literature on straight skeletons and have found the expression ghost edge for this type of extra edge. Unfortunately, I couldn't find any solution to remove them properly, except by a completely new algorithm (see weighted straight skeleton).

As a workaround, I have removed this edge in my solution. However, I don't know, if this can be generalized for all types of ghost edges. Let us see if you encounter any new example, where this type of issue (or a side effect of my solution) appears.

[vvoovv]

Still not happy. Sorry.
There is an extra vertex. Would it possible to get rid of it as well?

[vvoovv]

Another interesting case:
t_shape_extra_edge_2.py.txt
Everything is ok here, except of the fact that I'd prefer to have a quadrangle instead of the pentagon for the polygon with the selected vertices.

[vvoovv]

A brute force way to solve the problem of the extra vertices would be to iterate through the resulting faces that are polygons with at least 5 vertices and check if there are vertices with collinear adjacent edges.

Since those candidates for extra vertices are known in advance, would it be possible to optimize the task?

[vvoovv]

In order to be flexible for future I would extend it to two global variables.

I like your idea. debugRequests should also be a Python dictionary, since it's a slightly faster to check if a key is in the dictionary than to check if an element in the small list.

I'll regenerate the automated test to check for edge intersections.

[polarkernel]

I like your idea. debugRequests should also be a Python dictionary, since it's a slightly faster to check if a key is in the dictionary than to check if an element in the small list.

Implemented for 'skeleton'. I just changed the name debugOutput to debugOutputs, because there my be several.

[polarkernel]

Apse: I removed this feature for the moment and will try to find a better solution.

I like to demonstrate another solution that fails to create proper apses, using the file test_231994916_batticaloa_pillaiyar_kovil_street.py. Although this roof should not be tagged as 'hipped' in my opinion, it shows well the issues of a straight skeleton for apses, because it is extreme.

Using the current code, the skeleton becomes (click on the images to enlarge them):

I thought that the short edges produce imprecise bisectors, so that the edges do not end in one point. Therefore I detected the apse using the method already described and adjusted bisectors and edges to a circle. The red dots in the following image are the detected vertices belonging to the circle. The dotted lines are the bisectors after correction, now pointing all to the center.

But surprisingly, the skeleton does not accept this center, as the vertices not belonging to the circle produce events at lower heights that propagate to all bisectors (detail of the center in right image):

This is also the case for all apses that have structures near to their ends. For example see the result of the file test_95211505_berlin_oehlertring_55.py:

It is questionable, whether a straight skeleton is able to build the roof of an apse.

[vvoovv]

I regenerated all script to include testing of edge crossing. Edge crossings are encountered for all cases with vertex duplications and faces with less than 3 vertices.

I ended up with a single Python dictionary debugOutputs that has also the functions of debugRequests

[vvoovv]

I moved the script 102735465_gdansk_aleja_grunwaldzka_186.py to debug/edge_crossing

[polarkernel]

Again a step forward! In the remaining test files with issues, I found that dormers produce often weird skeleton edges, that lead to crossed edges some processing steps later. Therefore, I introduced a dormer-event, that is able to clean up these problems. But let me first explain the reason for these issues.

The following graph shows on the left the shape of what I named an dormer. Four vertices with a sequence of right - left - left - right turns (RLLR) in counterclockwise order by angles of almost 90°. In the middle, the progress of the moving vertices along the bisectors is depicted. At the end of the dormer-step, the green lines collapse to one line, which forms one of the edges of the skeleton. The red lines get unified to one line, which moves upwards from then on. This is what Tom Kelly called a parallel consecutive edge (PCE) degeneracy. The resulting skeleton part is what we like to have for roofs.

However, this is not what the algorithm by Felkel and Obdržálek delivers. Ghost edges leave the dormer and interact in an unpredictable way with other edges and events, as for example shown in the following images:

My solution starts by a pattern matching of the outer polygon contour, similar to what I tried for apses, but this time with success. In a first step, patterns of edge rotations with directions RLLR get detected. The shape of the pattern has then to meet some additional restrictions, as for instance:

• The length of the edges into and out of the dormer must have about the same length.
• The edges on both sides of the dormer (those that get unifyed later) must have a minimal length, depending on the directions of the turns at their far ends. This is required to avoid conflicts with other skeleton edges.
• The base width of the dormer is currently limited to 10 meters.
• Two dormers must not share vertices as for example here:

Once such a pattern is found, a DormerEvent is created at the position of the intersection of the outer bisectors, while the position of the intersection of the inner bisectors is stored within the event. The event distance is set as half the base width of the dormer.

A new event handler finally processes this event, once the event distance gets reached by the event loop. It removes the four vertices from the LAV, so that the red edges get unifyed to one edge that crosses the peak of the dormer. The restrictions for the dormer-event prevent any crossing edges. Two subtrees get created, one with a node at the position of the inner bisector intersection with sinks to the inner vertices and another at the position of the outer bisector intersection with sinks to the outer vertices and finally a third sink to the first node. And then, and this is important, no new event gets created, that could disturb other edges of the skeleton.

The result is remarkable. For instance the test file 4804904_london_temperate_house (you remember that I once refused to debug such a large footprint) passes all the tests without problems, while it produced crossing edges, faces with less than 3 vertices and duplicate vertices before (left before and right after introduction of dormer-events):

7 more test files could be moved from the debug folder to the tests folder, while 5 files still remain in debug. After about 240000 tests you made, this corresponds to an error rate of about 2.1 ppm. All of the remaining files in debug start with edge crossings, therefore I moved them to edge_crossing. However. these edge crossing also produce duplicated indices in polygonize(), so that they get detected without the edge-crossing test.

What a good state for the end of this year!

[vvoovv]

OSM footprints are never perfect. So dormers are unlikely to be perfectly rectangular and to have the sides of the equal length.

[vvoovv]

What a good state for the end of this year!

Sorry for being an evil hero who spoiled the good state for the end of the year.

I've just finished executing bpypolyskel against OSM footprints with the hipped roof. Nearly 320 thousand hipped roofs were tested totally.

The last 40 thousand footprints contain more problems than usual:

• 1 "float division by zero"
• 1 "list index out of range"
• 13 footprints with duplicated indices

The increase of the relative number of the problems can be partially explained by the polygons with holes. All of them are contained at the very end of the dataset.

I'll add the related scripts in the coming days.

[polarkernel]

OSM footprints are never perfect.

I know. However, I was able to allow quite some tolerances and it still worked. The remaining files in debug have other reasons for crossing edges.

Sorry for being an evil hero who spoiled the good state for the end of the year.

Never mind. Maybe a good reason to perfect the library, there were few holes until now. Did you already use the newest version of the code?

[vvoovv]

Did you already use the newest version of the code?

Yes, the newest one.

• 1 "float division by zero"

That one was likely caused by a degenerate hole that touched itself. I'll fix it right in the OSM database.

[vvoovv]

• 1 "float division by zero"

Fixed the multipolygon. Below is the result. Looks nice.

[vvoovv]

Added a few new scripts to debug/edge_crossing and one script to debug/repeated_indices (it passed the test edge_crossing).

[vvoovv]

Added the remaining tests.

Bpypolyskel fails for 30 footprints from the whole dataset (99.99% are successful):

• 2 footprints: only the test for duplications fails
• 28 footprints: the test for duplications and the one for edge crossing fail

[polarkernel]

Added the remaining tests.

Today a small progress. The center of apses becomes a cloud of skeleton nodes, even when the generationg edges are almost perfect. Computing the intersections of the bisectors from a sequence of short edges amplifies floating point errors. Parts of this cloud may then be subject of merging by the clustering algorithm. Often, this produces edge crossings, as for example here:

To avoid this effect, I revived the apse detector. But it is used only to prevent cluster merging in a region around the detected center of the apse arc. This allowed to move the following files from edge_crossing to tests:

test_530431938_denver_basilica
test_336079890_kropivnitsky_yaroshenka_1

The reason of issues for the remaining files stays strange.

[vvoovv]

The apse of the Basilica in Denver:

[vvoovv]

test_336079890_kropivnitsky_yaroshenka_1: similar apses

[vvoovv]

I added the code to the Blender-OSM addon to generate flat roofs as a fallback solution if there is a vertex duplication.

Do you think it's possible to improve anything? Or we could leave it as is. I think it is really an outstanding result in its current state.

[polarkernel]

The apse of the Basilica in Denver:

I know, but I can't solve it, see my comment here (I didn't commit this solution then). The current solution only prevents edge crossings due to merging.

Do you think it's possible to improve anything?

No, I don't think so. Since several days, all ideas that improve files in debug lead to about the same number of errors in tests. I think at some point is should be good enough. Let us mature the library at the users.

[vvoovv]

Just an idea.

Apses never have a sharp end in the real world.

Suppose an apse was detected. Would it help somehow if a small polygon was generated for the end of the apse instead of the sharp end?

[polarkernel]

Just an idea.

I will think about. Maybe this should be done at the end of polygonize(), for instance by cutting the peak around the center of the apse. It can't be done in skeletonize(), as the algorithm is 2D and it would be very difficult to add height to the small polygon.

[vvoovv]

Do you think it would be possible to include the feature #9 to bpypolyskel?

[polarkernel]

Just an idea.
Apses never have a sharp end in the real world.
Suppose an apse was detected. Would it help somehow if a small polygon was generated for the end of the apse instead of the sharp end?

No chance! I did some more experiments towards this goal. In principle, one could try all or some of the following steps to get a better peak of an apse:

1. Detect vertices that belong to an apse.
2. Fit a circle to these vertices.
3. Adjust these vertices and their edges and edge normals to be on the circle.
4. Compute the vertexes and the LAV that correspond to this new contour.
5. Construct the edge events positioned at the center of the circle for these vertexes.
6. Initialize all other events as usual, but not those for the vertexes that belong to the circle.

I did that once for the file test_336079890_kropivnitsky_yaroshenka_1 you mentioned above. On the left you see the all events for this example and on the right an enlarged detail of the skeleton, constructed by the algorithm, at the center of the right tower:

There are edge crossings produced during the run by events outside of the apse. In anext experiment, I changed step 5, so that the events got created on a small circle with a radius of 5% of the original radius and an adjusted height to still fit to the contour. On the left a detail of the events and on the right a magnified detail of the skeleton around the center of the circle:

The constructed events get continued to the center of the circle. But at a closer look at this new center we see that the skeleton got detached, producing impossible faces:

If one omits the steps 5 and 6 and just corrects the contour, the events get distributed due to floating point errors. Again, the result is a weird skeleton:

Conclusion:

This algorithm HATES symmetries, rectangles and precise edges! It bases only on local information and propagates any change, even if it's only a millimeter, uncompromisingly to the bitter end, often onky visible in the last processing step.

My idea of cutting the peak around the center of the apse will also not work. In contrast to ropivnitsky_yaroshenka, the peaks have often the same height as the ridge of the roof (see left for the Basilica in Denver) or even have no peak at all (right):

[vvoovv]

Ok. Thank you for the explanation!

[vvoovv]

I noticed there are imports

import numpy as np
from numpy.linalg import inv
import cmath

in bpypolyskel/bpyeuclid.py that are never used.

Can they be removed?

[polarkernel]

Can they be removed?

Sorry, these are relicts from some trials.

import numpy as np
from numpy.linalg import inv

can be removed.

import cmath

is required by fitCircle3Points(). Complex numbers provide a very fast method to fit a circle through three points.

[vvoovv]

import cmath

is required by fitCircle3Points(). Complex numbers provide a very fast method to fit a circle through three points.

No function from cmath is used in the code. The Python object complex to create a complex number is available without any import. I'd suggest to remove that import as well.

[polarkernel]

No function from cmath is used in the code

Didn't know. Then it's clear, this import may be removed.

[vvoovv]

I removed those imports and tagged the latest version as v1.1.1

[polarkernel]

removed those imports and tagged the latest version as v1.1.1

Thanks!!