Botffy / polyskel

Straight skeleton implementation in Python
GNU Lesser General Public License v3.0
79 stars 21 forks source link

Question about handling a split event #9

Closed warrenbocphet closed 4 years ago

warrenbocphet commented 4 years ago

Hi, I was reading through your code to understand the paper better. There is one point that I'm not quite sure about.

In step 2e (for non convex polygon), dot point 4 (I'm referring back to the paper), the author define point X and Y as end point and start point of the opposite line segment. The way I understand this is literally the 2 end points of one line segment.

But after reading your code (and to be fair I didn't spend enough time to process them all), I think the way you implemented it is the 2 vertexes that is currently pointing to that line segment. Which is not similar to my understanding, since the vertexes pointing to one edge can change during the computation...

Another thing I wasn't quite understand (from the paper) is that on step 2c (for non convex polygon), the authors mentioned outputing the "local peak" of the roof. I can see you did the same for edge event, but not with split event. Why is that?

Thanks!

Botffy commented 4 years ago

Hi!

2e/4: X and Y are the endpoints of the line segment opposite to the vertex causing the split event, yes. Now, while conceptually the vertices do change as the polygon shrinks, the algorithm never actually recalculates the positions of the vertices. The algorithm works with the angles between the edges, and those never change during the shrinking. Vertex positions are calculated only if a new vertex appears, and then we calculate its bisector angles using the original line segments.

2c: split events may split a LAV in two so that one of the resulting LAVs contains only two vertices -- that is, the polygon collapses into a line. In this case, the two endpoints of the line are peaks. This case is handled on line 298 (in the code peaks are gathered in the sinks list). However, off the top of my head I don't know why I'm outputting only one of the points. It MAY have something to do with the fact the l.head there refers to the v1 vertex created at the intersection point, at the same position as the v2 that goes into the other LAV, so the v2 still takes part in the calculation. But I'm not sure, actually. I may look into it in the near future.

warrenbocphet commented 4 years ago

Wow, that clears things up a lot! Another question I have is, after the initialisation, you will have edge events and split events. So is it possible that another new split event to occur during the execution of the second step? In the paper, the author said it is, but I have seen some papers claim it is not...

Botffy commented 4 years ago

Yes, split events may happen later during the computation. A split event means a concave vertex hits an edge, thus splitting the polygon in two. But the edge it hits is not necessarily an edge of the original polygon, it might have been created only after a series of other events.

Also, consider that the two polygons born after a split events are not necessarily convex either, so another split event may very well occur, splitting them further. But split events get rarer and rarer, as the polygons get "more and more convex" during the shrinking. An interesting point here is that concave vertices may be only created when two split events meet at the same intersection point. These are what Huber calls vertex events (p5-6). The algorithm of Felkel and Obdrzálek doesn't deal with those.

Also, I thought a bit more about the problem of the split event that creates a polygon with only two vertices. Well, of course we record only one of them as a peak, as the two vertices weren't created at the same time, so one of them is "higher" up than the other. I also understand why I chose the already existing vertex, but now I'm unsure if that case is even possible (we should do some geometry to see that).

warrenbocphet commented 4 years ago

Thanks a lot. I have been working towards an implementation of this library in Matlab (since most of my work is tied to Matlab). I just realised there is a case where the intersection link to two parallel edges (such as in here). So how would you calculate the bisector in this case? Also, would two same intersection be acceptable? Like in here, the letter V is symmetric, so it is possible to have 2 intersections with the same coordinates, but pointing to different set of vertexes.

Botffy commented 4 years ago

The bisectors of A and B collide in C. We can calculate the bisector of C by taking the bisector angle of the left edge of A and the right edge of B, which we can get by simply adding the vectors defining those edges (since we aren't interesting in the length, we don't even need to normalize them or anything). It doesn't matter that they're parallel, we'll just get the same angle. No problem there.

But the symmetry can be problematic if there's a reflex vertex involved, yes. It's possible that two split events will happen at the same place, and this implementation (IIRC) will cause the second split event to fail (produce incorrect results). Luckily, this is fairly rare in most use-cases.

(If you decide to make your Matlab implementation public, I'll happily add a link it in the Readme)

warrenbocphet commented 4 years ago

Definitely. Thanks for the help so far! Quick question, I was trying out Yongha's fork and I realize it seems like polyskel only accept vertexes with positive coordinate? Is that from your branch or is it Yongha's modification? If it is yours, is there any reason why?

I also tried a couple of different shape, other than in the demo and the result wasn't quite as expected. For example: polygon = np.array([[200,600], [200,200], [500,200], [500,600], [400,600], [400,300], [300,300], [300,600]]) which is a letter U will give a skeleton that is outside the polygon.

polygon = np.array([[160,100], [160,60], [180,40], [240,40], [260,60], [260,100], [240,120], [180,120]]) will give no skeleton at all. I'm pretty sure I got the correct orientation (CCW) and for both case, the hole is [].


import polyskel
from polyskel import skeletonize

import matplotlib.pyplot as plt
import numpy as np

# Rectangular
polygon = np.array([[7,3], [7,7], [13,7], [13,3]])

# # Letter U
polygon = np.array([[200,600], [200,200], [500,200], [500,600], [400,600], [400,300], [300,300], [300,600]])

# Modified letter U
polygon = np.array([[2,10], [2,4], [4,2], [8,2], [10,4], [10,10], [8,10], [8,6], [6,4], [4,6], [4,10]])

# Diamond
polygon = np.array([[160,100], [160,60], [180,40], [240,40], [260,60], [260,100], [240,120], [180,120]])

skeleton = skeletonize(polygon, [])

toPlot = np.append(polygon, [polygon[0]], axis=0)
print(skeleton)

plt.gca().set_aspect('equal', adjustable='box')

plt.plot(toPlot[:,0], toPlot[:,1])

for arc in skeleton:
    for sink in arc.sinks:
        plt.plot([arc.source.x, sink.x], [arc.source.y, sink.y], c="black")

plt.show()      
Botffy commented 4 years ago

Wow, that's pretty horrible! (the letter U)

I don't know. I may look into it in the weekend. I do know symmetry, especially with reflex vertices, may very well cause bad behaviour. But this thing with the U looks extreme, are you sure it's not the orientation?

Botffy commented 4 years ago

Hi, "in the weekend", heh :(

So, I looked into the letter U, and it's about the order as you've suspected. A-and it's not ccw for polyskel because the coordinate system I use has the y axis pointing down. Because I don't know, because that seemed natural to me because computer graphics and everything.

But this is a bug indeed, because that's not obvious or intuitive at all. At the very least I should document it.