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]

If I somehow can help in localizing the bug, please let me know.

[polarkernel]

If I somehow can help in localizing the bug, please let me know.

It's weird. I constructed a reduced footprint of your example, which contains only the right part of the building, where the issue occurs in the original (see bldg_nikolskaya_street_reduced.py.txt). Although the border which closes the polygon now seems not to contribute to the issue, skeletonize() work correctly for this example. Now, I try to find out, where the difference occurs and why. Don't know if you can contribute anything.

[polarkernel]

It's weird.

I am sorry, but I have to take a break. This is the most boring task I ever did! Two days now and I still haven't the slightest idea on where this issue occurs. Everything depends on what has been computed before. Just as an example, after the first time, the _nextevents of all Vertex objects are computed, the result looks like below for your example (blue: EdgeEvent, red: SplitEvent, number: distance): The same for the reduced footprint looks like this: The first produces the error, but the latter not. Mainly the split events are different. But where does the issue start? I have no idea.

[vvoovv]

I am sorry, but I have to take a break. This is the most boring task I ever did! Two days now and I still haven't the slightest idea on where this issue occurs.

I understand it.

I noticed the distances for the lower edge events are different: 2.83 vs 2.84. But I don't think it's related to the problem with that polygon.

[polarkernel]

It's weird.

After some more hours, I have two new findings:

1. The difference of events between the reduced polygon and the original depends mainly on one reflex vertice (red circle). With my new example, I get (almost) the same events and positions as the original, and then, the crossing of skeleton edges reappears. Somehow, the issue must occur later, during the processing of the event queue
2. There is a difference between the original paper and the actual implementation. Felkel and Obdržálek write there (page 7, point 1c): Store the nearest intersection point I of these three ones into the priority queue, while in the implementation all valid intersection points of reflex vertices get enqueued as SplitEvent. I tried the paper's proposition once, but there was no difference. It seems that the priority queue solved this by presenting the nearest events first.

[vvoovv]

Another weirdness. It looks like the same issue as the latest one.

hole.py.txt

Selected is the single weird face generated by bpypolyskel:

[polarkernel]

Another weirdness. It looks like the same issue as the latest one.

This one helped me to locate a (strange) issue. About three weeks ago , I changed the '>' in line 131 to '<', which was also logical to me (the description in my comment at the time was the wrong way round). The smaller the absolute value of the dot-product is, the less parallel are the vectors, in my understanding. At this time, this solved the issues that we had. Now I changed it back to the original, which does the contrary of its intention, but the results for the skeleton are better.

Its weird, this algorithm is still irritating to me.

[polarkernel]

Its weird, this algorithm is still irritating to me.

I updated all examples to the non-contiguous contours version. Doing this I found the counter-example for the last fix. The symmetrical case in _holesymmetrical.py, where, by this usage of the dot-product, the selected egdes are exactly parallel. And this leads to a wrong result, as expected.

[vvoovv]

Is there a way to differentiate between those 2 cases: > and < ?

[polarkernel]

Is there a way to differentiate between those 2 cases: > and < ?

The simplest way I have found until now is to use something like this

if prevdot < EPSILON:
    selfedge = self.edge_prev
elif nextdot < EPSILON:
    selfedge = self.edge_next
else:
    selfedge = self.edge_prev if prevdot > nextdot else self.edge_next 

This decision works for all examples we have at time. However, I do not know if it can be generalized. Actually I try to find out, why '<' does not work.

[vvoovv]

Some good knews for now. I've integrated the hipped roofs for polygons with holes into the addon:

[polarkernel]

Some good knews for now.

Yeah, that looks really great! Congratulations! Again another great and encouraging step.

[vvoovv]

Another test case: c_shape.py.txt

If there is > in the line 131 of bpypolyskep.py:

Since the edges of the middle part are parallel, the related edge of the straight skeleton must be also parallel to them. However it is not the case. The straight skeleton doesn't look to be the correct one.

---

If there is < in the line 131 of bpypolyskep.py. The result looks ok.

[vvoovv]

c_shape_2.py.txt

If there is > in the line 131 of bpypolyskep.py:

---

If there is < in the line 131 of bpypolyskep.py: Endless or at least very long lasting cycle in the line

faces3D = graph.faces(embedding, firstSkelIndex)

[vvoovv]

c_shape_2_duplication.py.txt

If there is < in the line 131 of bpypolyskep.py: I am getting vertex index duplications

[polarkernel]

Actually I try to find out, why '<' does not work.

I think I am one step further. Currently, I try to verify step by step whether the algorithm delivers correct results, without only looking at the final result. Felkel and Obdržálek write in their paper about this step:

We then compute the coordinates of the candidate point Bi as the intersection between the bisector at V and the axis of the angle between one of the edges starting at V and the tested line segment ei (see fig. 4). Simple check should be performed to exclude the case when one of the line segments starting at V is parallel to ei.

This means that either selfedge should deliver the same result, as long as it is not parallel to the line segment. For the actual implementation, this is not the case. Here are the results for '<' (on the left) and '>' (on the right) for your example hole.py.txt:

Let me first discuss the case '<' on the left, which results in the less parallel edge. The red edge is selfedge with a black elongation to the intersection (red dot) with the line segment e (green line). The dotted ray to the right is the bisector between these two lines currently used by the code. The dotted ray to the left is the bisector of the Vertex under consideration. As can be easily seen, these rays do not intersect!

On the right is the result for the case '>', the more parallel selfedge. Here, the intersection between selfedge and e occurs very far away from the polygon. But in this case, the bisectors intersect in the polygon (almost no more visible in this scale). However, with these long lines, this intersection cannot be very precise. So, either selfedge do not deliver the same result, something mus be wrong!

If we go back to the image at the left, we see that the bisector at the intersection is on the wrong side. It should be on the left of the intersection at the same side as the bisector of the Vertex. This condition can be forced by replacing the condition in line 138 of bypolyskel.py by

             if self.bisector.v.cross(linvec) < 0:

Using this condition and both cases '<' (on the left) and '>' (on the right) , the results for both selfedges get:

Now we get an (almost) equal intersection in both cases, which fulfills the description of Felkel and Obdržálek. The computed split-events (red dots) look reasonable to me: Note the blue double-next-events, but this is to solve later, if required. The final skeleton then becomes:

I was very happy when I looked at this result, but ... ... for your example c_shape.py.txt I got also reasonable positions for the split-events (red) ...

... but the skeleton edges again have intersections:

There must be more bugs in the processing after the initialization. :-(
I will investigate them further, but this will take some time, at the moment I'm very busy elsewhere. Maybe you find the time to investigate which examples work now correctly with these changes and which not. I did not yet commit this solution, but I thik it's easy to introduce it manually. Sorry for the very long comment.

[vvoovv]

Maybe you find the time to investigate which examples work now correctly with these changes and which not. I did not yet commit this solution, but I thik it's easy to introduce it manually. Sorry for the very long comment.

Thank you for the detailed description. I'll check it out.

[vvoovv]

             if self.bisector.v.cross(linvec) < 0:

I commited that line.

The problems are exhibited in:

• c_shape.py.txt
• c_shape_2.py.txt

[polarkernel]

I commited that line.

Thanks.

The problems are exhibited in:

Thanks for searching. Meanwhile I assume there is again a floating point precision issue, but I was not yet able to catch it.

[polarkernel]

I commited a new update. This time, the fixes are more based on intuition than on real knowledge, but it seems to help. For the c_shape examples, I observed that after initialization, the two top events have the same distance. But their position in the priority queue is random. If I processed one of them first, the result was garbage, while the other got discarded during the following handling process. The same occured, if I started the process with the other one.

As a first idea, I tried to extract and process both of them before new events get extracted from the queue. I had to use a robust floating point comparison of their distances to make it work correctly. The result was ways better, but not for all examples. After some more experiments, I found that this idea is not correct, if the top events share not only the same distance but also the same position. Then, the second event has to be pushed back in the queue.

I can't explain why it works and it is also not part of the original paper by Felkel and Obdržálek. But it seems to work for all examples. Now let me know if you get happy with it.

[vvoovv]

Many discoveries are made through the intuition.

Yes, it works for all our cases so far! Including bldg_nikolskaya_street.py.txt.

[vvoovv]

Now let me know if you get happy with it.

Happiness probably doesn't exist.

Encountered a new weirdness.

berlin_karlshorst_zwieseler_strasse_40.py.txt

That one returns the index -1 for one of the resulting faces. And it happens twice for that face:

[107, 108, 151, 158, -1, 157, 164, 163, 166, 165, 65, 66, 130, 132, 134, 157, -1, 158, 159, 160, 174, 172]

[polarkernel]

Happiness probably doesn't exist.

Seems to be an endless story. This bug got fixed. At the end of skeletonize(), singular nodes (those without sinks) of the skeleton got removed, but not eventual sinks form other nodes, that pointed to these singular nodes. Like this, vertices remained that haven't been contour vertices nor skeleton nodes. So maybe: happiness += 1

[vvoovv]

Seems to be an endless story.

Exactly. Duplicated vertices in a resulting face: berlin_karlshorst_roemerweg_117.py.txt

[14, 15, 19, 15, 8, 20, 18, 22, 23, 22, 21] The indices 15 and 22 are repeated twice.

[polarkernel]

Duplicated vertices in a resulting face:

After my examinations I have to state that we are moving slowly towards a delicate terrain. It seems that the straight skeleton algorithm produces a skeleton that can't be cleaned by the current measures. Floating point precision begins to play a dominant role. The example _berlin_karlshorst_roemerweg117.py uses quite high coordinates in the range of 800 and 300. The raw skeleton algorithm delivers the following skeleton before cleaning:

After moving it near to the origin of the coordinate system, this skeleton looks like this:

The result depends on the size of the coordinates, which indicates that there is a floating point precision problem. It looks very similar to the many footprints of this type we already have in this bug list, where the ghost vertice indicated by the red arrow could successfully get eliminated. However, this time this is not the case. The following plot shows a magnified view of this position:

There are two vertices, with a distance of 0.0006729 between them. If the coordinates are in meters, the distance is less than one millimeter, which cannot make any architectural sense. When these vertices get merged to one, the result looks fine:

I have now commited an experimental version, which includes an internal shift of the footprint's center of gravity to the center of the coordinate system. Additionally, I replaced the original function _mergesources(), which merged only exact duplicates, by a new function mergeNearbyNodes() which merges skeleton nodes, as the one described above, when they are closer to each other than 0.001. This distance is currently a (default) parameter of this function, but it could also be introduced to skeletonize() or even to polygonize().

Maybe you could test it now, as you have many more footprint examples available than I do. The main questions remaining are:

• What is a reasonable parameter for mergeNearbyNodes()?
• One of the nearby nodes is merged into the other. Does it play a role which node is the target?

[vvoovv]

.. an internal shift of the footprint's center of gravity to the center of the coordinate system.

That's a good idea!

combs = combinations(arc.sinks,2)

There are N*(N-1) combinations where N is the number of sinks. It would be good to optimize it.

d = (arc1.source-arc2.source).magnitude

Here _.lengthsquared can be used for optimization purposes instead of .magnitude,

• One of the nearby nodes is merged into the other. Does it play a role which node is the target?

Do you mean: which of the nodes should be preserved?

[vvoovv]

Vertex duplication again: berlin_karlshorst_robert_siewert_strasse_52.py.txt

[40, 41, 63, 69, 68, 65, 64, 36, 37, 61, 67, 68, 69, 71, 66, 62, 33, 34, 65, 68, 67] The node 69 is encountered twice. Thenode 68 in encountered thrice.

[polarkernel]

There are N*(N-1) combinations where N is the number of sinks. It would be good to optimize it.

I agree. However, I'll wait until we have confidence that merging nearby nodes is really useful correction and doesn't produce unwanted side effects.

Here .length_squared can be used for optimization purposes instead of .magnitude,

I agree. Maybe even using the Manhattan distance could work.

Do you mean: which of the nodes should be preserved?

Yes.

[polarkernel]

Vertex duplication again: berlin_karlshorst_robert_siewert_strasse_52.py.txt

This example leads us to a new principle problem. The straight skeleton algorithm produces the following skeleton:

If we take a closer look at the center of the apse (red arrow) we see, that a multitude of nodes has been constructed around this center:

For debugging purposes I use another view, where I add an increasing height to all nodes and their arcs. In this view we see that a total of 15 nodes has been constructed near this position, where 4 of them even dont have sinks:

Obviosly, the unavoidable inaccuracies of the vertices of the polygon, produced by the user that entered this building to OSM, get amplified when computing the intersections of the bisectors. I don't see a way to properly clean up them using any postprocessing.

I did not yet a deep search, but maybe it would be more helpful to "rectify" somehow the footprint by a preprocessing, before computing the skeleton. Do you perhaps have more knowledge from your work with architectural algorithms?

[vvoovv]

• What is a reasonable parameter for mergeNearbyNodes()?

I suggest to add an optional parameter for the function polygonize(..).

• One of the nearby nodes is merged into the other. Does it play a role which node is the target?

At the moment I can only suggest to pickup the first node out the merged nodes.

[vvoovv]

If we take a closer look at the center of the apse (red arrow) we see, that a multitude of nodes has been constructed around this center:

How about merging those nodes? I tried to increase the merge distance in mergeNearbyNodes(..), but got errors.

An alternative way of merging would to group the nodes with nearly the save time of the edge event and then perform another grouping of those nodes by neighboring polygon edges.

[polarkernel]

How about merging those nodes?

I commited a new version using a function called mergeNodeClusters() to merge nearby nodes, much more sophisticated than the previous one. It starts with a sort, running in O(NlogN), and continues, replacing node clusters by merged nodes, all this running in O(N). First it finds pairs of nearby nodes where their source vertices lie in a square of size mergeRange. This range value has also been introduced as additional parameter at the end of polygonize(), using a default value of 0.15, which was the minimum working value for the example _berlin_karlshorst_robert_siewert_strasse52.py.

These pairs are then collected to clusters (there may be more than one in a skeleton) within an increased range of currently 5*mergeRange and then merged into one node per cluster. Maybe this range has to be adapted, depending on the experience on other examples. As center of this node, the center of gravity of the merged nodes is used and its height becomes the average of the heights of the merged nodes. All the sinks of the clustered nodes, that leave the cluster, become sinks of the cluster node and all sinks pointing to the clustered nodes get redirected to the cluster node. Then all clustered nodes get removed.

This seems to work for the example _berlin_karlshorst_robert_siewert_strasse52.py and all the other examples I have:

I hope it works also for future examples.

[vvoovv]

The roof now looks mostly ok:

But there are a couple of issues: (1) How are these edges formed?

Below is the magnified version of the above image:

Is it the result of the closed nodes removal?

---

(2) There is a gap near the vertex where the selected edges from the previous case meet:

Is it also the result of the closed nodes removal?

[vvoovv]

Another weirdness from Karlshorst district in Berlin: berlin_karlshorst_rheinpfalzallee_82.py.txt

-1 (twice) and duplicated index 114 in the face [58, 59, 112, -1, 114, 102, 66, 67, 115, 77, 78, 103, 114, -1, 111, 92, 54, 55, 94]

-1 and duplicated index 97 in the face [86, 87, 97, 113, 97, 44, 45, 91, 48, 49, 89, 52, 53, 92, 111, -1, 112, 95]

In general it looks quite weird:

[polarkernel]

But there are a couple of issues: (1) How are these edges formed?

The issue can be found at the two nodes marked by the rectangle:

From this coarse view, I would expect that the short edge between them would be considered as ghost edge and be removed, as we encountered in earlier examples. Let us look at the details: There are more than two nodes there. The two nodes in the upper right have a distance of 0.012, are therefore considered as cluster and get merged: But the edges to the left from the merged node are no more considered as parallel, as it should be the case for ghost edges. Therefore, this spike remains in the skeleton.

(2) There is a gap near the vertex where the selected edges from the previous case meet:

I couldn't find this gap within the faces. But I'm afraid that this issue is produced by height errors. Both merged nodes, the one at top of the aspe and this one, get the average height of the merged nodes as new height. For faces with more than 3 nodes, as there are some around the merged nodes, this could move the node out of the plane of the face. I don't know how Blender treats faces, when their vertices are not in one plane.

My conclusion: I'm afraid that we are getting more and more lost in a maze of special cases. I don't know how to solve that, but the main issue still remains in this imprecise skeleton algorithm.

[vvoovv]

A narrow gap would be ok.

However the spike doesn't look nice.

What is the difference between the examples _berlin_karlshorst_roemerweg117.py and _berlin_karlshorst_robert_siewert_strasse52.py?

The spike gets eliminated in the former case and it's preserved in the latter one?

[polarkernel]

What is the difference between the examples berlin_karlshorst_roemerweg_117.py and berlin_karlshorst_robert_siewert_strasse_52.py?

The latter is less parallel than the former. The limit for this correction is set in line 446 of bpypolyskel.py and is currently EPSILON. You may change it to a maximum angle alpha, that you like to accept for such spikes, by setting it to abs(1-cos(alpha)). For _berlin_karlshorst_robert_siewert_strasse52.py, a value of 1.0e-4 removes the spike, but you may encounter other examples. Please set it as you like and commit the change, once it is defined.

I am still working on the issue in example _berlin_karlshorst_rheinpfalzallee82.py but it is difficult to debug. I hope I find a solution today.

[vvoovv]

Please set it as you like and commit the change, once it is defined.

Ok, I'll do it.

Here is a slightly modified berlin_karlshorst_zwieseler_strasse_40.py.txt

I set the roof height to 3. But the resulting roof height after the call of polygonize(..) is more than 100 meters.

[vvoovv]

I set the roof height to 3. But the resulting roof height after the call of polygonize(..) is more than 100 meters.

The same problem is in the other examples. Somehow the supplied roof height isn't respected.

[polarkernel]

I set the roof height to 3. But the resulting roof height after the call of polygonize(..) is more than 100 meters.

Fixed. In the old version, the maximum height was found as height of the last skeleton node. As the nodes get sorted in mergeNodeClusters() now, this was no more the case. No more unwanted towers in Berlin!

[polarkernel]

Another weirdness from Karlshorst district in Berlin: berlin_karlshorst_rheinpfalzallee_82.py.txt

Fixed. There were several bugs. I had to set the 'spike parameter' to 1.e-2 to make this example work. One bug was irritating me for quite a long time and I still do not under stand, why this happened. The follwing code does change skeleton:

        for arc in skeleton:
            if arc.source not in mergedSources:
                for i,sink in enumerate(arc.sinks):
                    if sink in mergedSources:
                        arc.sinks[i] = new_source

while this solution does not:

        for arc in skeleton:
            if arc.source not in mergedSources:
                for sink in arc.sinks:
                    if sink in mergedSources:
                        sink = new_source

I can't understand, why not. But never mind, I am waiting for the next special case.

[polarkernel]

Fixed another bug in _cleanskeleton().

[vvoovv]

Example _berlin_karlshorst_rheinpfalzallee82.py.txt How are the selected nodes formed? They aren't connected to the vertices of the original polygon:

[vvoovv]

But never mind, I am waiting for the next special case.

They will never end!

berlin_karlshorst_roemerweg_62.py.txt

Got an error.

[polarkernel]

They will never end!

You know that my lifetime is limited! However, by end of next year, all cleaning and corrections will create a skeleton without skeletonize(). :-)

[polarkernel]

Example berlin_karlshorst_rheinpfalzallee_82.py.txt How are the selected nodes formed? They aren't connected to the vertices of the original polygon:

We have a similar case of parallelism for adjacent parallel edges that become connected in this case. Due to merging of nodes, they haven't been exactly parallel any more. I assume that both cases may share the same limit and I have introduced a parameter PARALLEL on top of bpypolyskel.py that is used for both. Let's see if my assumption persists.

[polarkernel]

They will never end! berlin_karlshorst_roemerweg_62.py.txt Got an error.

The value for PARALLEL of 1.0e-2 was too agressive for this example. If set to 1.0e-3, it works fine. I have made some minor changes to _cleanskeleton() in order to avoid a program crash in case of destroyed skeletons at this place. But this would not help.

[vvoovv]

You know that my lifetime is limited! However, by end of next year, all cleaning and corrections will create a skeleton without skeletonize(). :-)

That's an opportunity to publish a paper about the new straight skeleton algorithm!

[vvoovv]

I've fixed the indentation problem in the line 460 of bpypolyskel.py.

[vvoovv]

I've made some changes to have a more concise code and removed an excessive Python list inside the function max(..).

It is more concise and slightly more efficient to have

if candidates:
    ...

instead of

if len(candidates) > 0:
    ...

and

if not candidates:
    ...

instead of

if len(candidates) == 0:
    ...