GSoC2010 Sweepline Algorithm
Hi, Boost Community! I am interested in participation in GSoC 2010. The project I liked the most is the generic implementation of Sweepline Algorithm. Basically fast and robust Sweepline Algorithm requires implementation of next data structure: 1) Event Queue (priority queue). 2) Sweep Line (self-balancing binary search tree). 3) Output datastructures (voronoi vertices, triangulation edges, etc.). Based on that I got next questions: - Which kind of self-balancing binary trees is better for that kind of approach? I am thinking about Red-Black Trees in prefer to AVL trees, as sweepline algorithm requires more adding and deletion operations. And AVL trees are better for look up operations. - STL library includes two generic implementations of self-balancing binary trees (red-black): map, set. I am quite new with Boost library, but probably there are some more specific self-balancing binary tree implementations there? - For Voronoi diagram and delaunay triangulation input data is set of points. What about input of the Medial Axis problem? -- Best, Andrii Sydorchuk
Andriy Sydorchuk wrote:
Basically fast and robust Sweepline Algorithm requires implementation of next data structure: 1) Event Queue (priority queue). 2) Sweep Line (self-balancing binary search tree). 3) Output datastructures (voronoi vertices, triangulation edges, etc.).
These are intersting technical questions. I guess the most important thing to do is "writing the proposal for GSoC". But you are certainly right that a good understanding of the technical details helps writing a good proposal (and helps deciding whether you want to work on this problem at all).
Based on that I got next questions:
- Which kind of self-balancing binary trees is better for that kind of approach? I am thinking about Red-Black Trees in prefer to AVL trees, as sweepline algorithm requires more adding and deletion operations. And AVL trees are better for look up operations.
Most implementations I have seen so far simply tried to keep the implementation effort manageable, and made use of existing stl containers, using suitable custom comparison functors.
- STL library includes two generic implementations of self-balancing binary trees (red-black): map, set.
Indeed.
I am quite new with Boost library, but probably there are some more specific self-balancing binary tree implementations there?
Boost.Intrusive could also be intersting to look at.
- For Voronoi diagram and delaunay triangulation input data is set of points. What about input of the Medial Axis problem?
Strictly speaking, a single polygon (with potential holes). However, the actual input is often a polygonal domain (= set of non-overlapping polygons), and the output is normally equivalent to computing the medial axis for each polygon of the polygonal domain separately. Regards, Thomas
Andriy Sydorchuk wrote:
Hi, Boost Community!
I am interested in participation in GSoC 2010. The project I liked the most is the generic implementation of Sweepline Algorithm.
Basically fast and robust Sweepline Algorithm requires implementation of next data structure: 1) Event Queue (priority queue). 2) Sweep Line (self-balancing binary search tree). 3) Output datastructures (voronoi vertices, triangulation edges, etc.).
Based on that I got next questions:
- Which kind of self-balancing binary trees is better for that kind of approach? I am thinking about Red-Black Trees in prefer to AVL trees, as sweepline algorithm requires more adding and deletion operations. And AVL trees are better for look up operations. - STL library includes two generic implementations of self-balancing binary trees (red-black): map, set. I am quite new with Boost library, but probably there are some more specific self-balancing binary tree implementations there? - For Voronoi diagram and delaunay triangulation input data is set of points. What about input of the Medial Axis problem?
You should absolutely use stl map or set. See: Ruud, B. Building a Mutable Set, 2003. Retrieved March 3, 2009, from Dr. Dobb's Portal: http://www.ddj.com/cpp/184401664 For an explanation of how to use map or set for a sweepline. You don't need to use the optimization of casting away const of the key. Just remove and reinsert with a new key. That's what I do. I have already implemented something like 12 sweepline algorithms in Boost.Polygon using the STL data structures. You can use a map or set for the event queue, otherwise use a stl vector with the stl heap operations if you are performance conscious. This is an optimization that can be performed after you get the algorithm working. It is better to take all of your original input events and sort them in a vector then merge them with events generated by the sweepline itself on the fly. Since a sorted vector satisfies the properties of a heap you could use it with the stl heap operations, however, it requires that your sweep algorithm know the data structure its input is coming from, and I have always abstracted that away with iterator interfaces. It is reasonable to assume that the input is a vector, however. If you implement Medial Axis with sweepline that works on non-convex polygons with holes you have implmented an algorithm that can also handle multiple non-overlapping polygons. Regards, Luke
For an explanation of how to use map or set for a sweepline. You don't need to use the optimization of casting away const of the key. Just remove and reinsert with a new key. That's what I do.
That was the first idea I was thinking: creating generic wrapper around set datastructure (Adapter pattern) with definition of node comparer. So it will look smth like this: template <class T> struct NodeComparer { public: } template <class T> class BeachLine { public: ... private: ... typedef std::set<T, NodeComparer<T> > bl_set; ... bl_set beach_line_; } You can use a map or set for the event queue, otherwise use a stl vector
with the stl heap operations if you are performance conscious. This is an optimization that can be performed after you get the algorithm working. It is better to take all of your original input events and sort them in a vector then merge them with events generated by the sweepline itself on the fly. Since a sorted vector satisfies the properties of a heap you could use it with the stl heap operations, however, it requires that your sweep algorithm know the data structure its input is coming from, and I have always abstracted that away with iterator interfaces. It is reasonable to assume that the input is a vector, however.
Yes, I've read about this kind of approach at this article .
For an explanation of how to use map or set for a sweepline. You don't need to use the optimization of casting away const of the key. Just remove and reinsert with a new key. That's what I do.
That was the first idea I was thinking: creating generic wrapper around set datastructure (Adapter pattern) with definition of node comparer. So it will look smth like this: template <class T> struct NodeComparer { public: bool operator() (const T& node1, const T& node2) const { ... } } template <class T, class Comparer = NodeComparer<T> > class BeachLine { public: ... private: ... typedef std::set<T, Comparer > bl_set; ... bl_set beach_line_; } You can use a map or set for the event queue, otherwise use a stl vector
with the stl heap operations if you are performance conscious. This is an optimization that can be performed after you get the algorithm working. It is better to take all of your original input events and sort them in a vector then merge them with events generated by the sweepline itself on the fly. Since a sorted vector satisfies the properties of a heap you could use it with the stl heap operations, however, it requires that your sweep algorithm know the data structure its input is coming from, and I have always abstracted that away with iterator interfaces. It is reasonable to assume that the input is a vector, however.
Yes, I've read about this kind of approach at this article: "An Efficient Implementation of Fortune’s Plane-Sweep Algorithm for Voronoi Diagrams" KennyWong, Hausi A. Muller. There are also few more nice ideas there. -- Best regards, Andrii Sydorchuk
Andriy Sydorchuk wrote:
For an explanation of how to use map or set for a sweepline. You don't need to use the optimization of casting away const of the key. Just remove and reinsert with a new key. That's what I do.
That was the first idea I was thinking: creating generic wrapper around set datastructure (Adapter pattern) with definition of node comparer. So it will look smth like this: template <class T> struct NodeComparer { public:
} template <class T> class BeachLine { public: ... private: ... typedef std::set<T, NodeComparer<T> > bl_set; ... bl_set beach_line_; }
Yep, from my code for the arbitrary angle polygon formation sweepline algorithm: template <typename Unit> class polygon_arbitrary_formation : public scanline_base<Unit> { ... protected: typedef std::map<vertex_half_edge, active_tail_arbitrary*, less_vertex_half_edge> scanline_data; typedef typename scanline_data::iterator iterator; typedef typename scanline_data::const_iterator const_iterator; //data scanline_data scanData_; Unit x_; //current x location bool justBefore_; //how to break ties using slope comparison ... }; In my code I call it scanline instead of sweepline. You want to use a map instead of a set because each element in the sweepline will need to have a data structure associated with it storing whatever data the problem you are solving requires. My sweeplines are often parameterized on the element type for the map and I use the same sweepline template for polygon clipping and connectivity extraction, for example, which require different data be stored on the sweepline associated with the line segments crossing the sweepline. Another tip that will help you with the sweepline algorithm is to use my approach of comparing keys in the sweepline by breaking ties based on slope comparion. Typically sweepline considers keys to be equal if they evaluate to the same y value at the current x. However, this makes it hard to manage the data structure if you are using a std::map without epsilon hackery. Instead you can remove all elements from the sweepline data structure that are going to cross at the current x location and reinsert them in the reverse order by having a pointer to justBefore_ in the comparison functor that breaks ties in y value by slope in reverse order when removing elements and forward order when reinserting elements. That way the data structure maintains its full proper ordering (and all important invariant) at all points during the execution of the algorithm.
Yes, I've read about this kind of approach at this article: "An Efficient Implementation of Fortune's Plane-Sweep Algorithm for Voronoi Diagrams" KennyWong, Hausi A. Muller. There are also few more nice ideas there.
I'll have to have a look at this paper. Thanks for posting the reference. Regards, Luke
You want to use a map instead of a set because each element in the sweepline will need to have a data structure associated with it storing whatever data the problem you are solving requires.
Well, we still can use set for this kind of a problem and store all the information in a key. But I will definitely agree that it is better to split key and value data using map. My sweeplines are often parameterized on the element type for the map and I
use the same sweepline template for polygon clipping and connectivity extraction, for example, which require different data be stored on the sweepline associated with the line segments crossing the sweepline.
I will have a look at your implementations. Thanks for the reference. Instead you can remove all elements from the sweepline data structure that
are going to cross at the current x location and reinsert them in the reverse order by having a pointer to justBefore_ in the comparison functor that breaks ties in y value by slope in reverse order when removing elements and forward order when reinserting elements.
I am a little bit confused. Lets say we have poitns A (0, 1), B (0, -1), C (1, 0). Our sweep line is moving along OX axis. Lets assume that points A and B are already processed by algorithm. So the next event is site C. Our map in this case will contain just one node with references to A and B points. So when we are adding point C we compare it's y-coord with y-coord of intersection of arcs made by A and B points. As you said in this case keys will be considered equal, as they have the same y-coordinate. Did I understood everything correctly by this moment? After that you proposed to remove current node (A, B) and reinsert new nodes (A, C), (C, B), am I right? Could you explain also "a pointer to justBefore_" meaning? Thank you forward. -- Best, Andrii Sydorchuk
Andriy Sydorchuk wrote:
Instead you can remove all elements from the sweepline data structure that
are going to cross at the current x location and reinsert them in the reverse order by having a pointer to justBefore_ in the comparison functor that breaks ties in y value by slope in reverse order when removing elements and forward order when reinserting elements.
I am a little bit confused. Lets say we have poitns A (0, 1), B (0, -1), C (1, 0). Our sweep line is moving along OX axis. Lets assume that points A and B are already processed by algorithm. So the next event is site C. Our map in this case will contain just one node with references to A and B points. So when we are adding point C we compare it's y-coord with y-coord of intersection of arcs made by A and B points. As you said in this case keys will be considered equal, as they have the same y-coordinate. Did I understood everything correctly by this moment? After that you proposed to remove current node (A, B) and reinsert new nodes (A, C), (C, B), am I right? Could you explain also "a pointer to justBefore_" meaning? Thank you forward.
What I mean by a pointer to justBefore_ has to do with what you have to do in the eventuality that you need to remove and insert many elements in the sweepline data structure at more than one colinear point. Imagine you have a 10x10 matrix of sites in the set of points that are arranged in a rectangular grid and a sweeline that is vertical and moving along the X axis. The sweepline will reach 10 points simultaneously at the same x stopping point. Each will require that elements need to be removed from the sweepline and new elements inserted. If you try to remove elements and then insert elements for each of the ten points one at a time how do you manage the sweepline data structure such that the invariant that all its elements are in strict sorting order at all stages of execution? The way I do it is by breaking ties in y intercept on the sweepline with the slope of the intercept on the sweepline. I remove all elements first using a reverse sorting order on slope and then change the sorting order on slope to forward and then insert all new elements. The justBefore_ value tells the comparison functor to sort slopes in reverse order if it is true, otherwise forward order. An alternative approach would be to make the comparison depend on the current Y. You scan event points from low to high on the sweepline. Your comparison functor can compare slopes in reverse order if their y intercept is larger than the current y and forward order if it is less than or equal to the current y. template <element_type> struct less_functor_by_x_and_y { coordinate *x_; coordinate *y_; less_functor_by_x_and_y (coordinate x, coordinate y) : x_(x), y_(y) {} bool operator()(const element_type& e1, const element_type& e2) { coordinate y1 = eval_at_x(e1, x_); coordinate y2 = eval_at_x(e2, x_); if(y1 < y2) return true; if(y2 < y1) return false; slope s1 = eval_slope_at_x(e1, x_); slope s2 = eval_slope_at_x(e2, x_); if(y1 < y_) return s2 < s1; //compare slopes in reverse order return s1 < s1; //y_ >= y1 and y2, compare slopes in forward order } }; The order of elements that are going to cross at a certain x and y are in the correct sorting order when the sweepline reaches that point if their slopes are compared in reverse order. To remove and reinsert those elements such that they are effectively swapped you remove them then update y_ to be equal to y and reinsert them. If you ever try to use the tree when its invariant that all elements are in strict sorting order is violated you will get undefined behavior and either crash or infinitely loop. It is managing the sweepline data structure and getting the comparison functor right that is the tricky part of implementing sweepline. If your eval_at_x and eval_slope_at_x are not numerically robust then you can't ensure a strict ordering of elements in the sweepline and the execution of the code will fail catastrophically. Regards, Luke
Maybe it is better to write proposal at first and after that discuss technical questions. But still I didn't get everything.
Each will require that elements need to be removed from the sweepline and new elements inserted.
What is your definition of element? Does leaf and node elements contain the same data? I am thinking about element like (leaves and nodes are the same): struct SweepLineElement { private: const Point2d* point_left; const Point2d* point_right; double getComparerValue(double line_x); }; So each node actually contains information about intersection point of two neighbour arcs from the sweep line, created by point_left and point_right at current X. Probably it's smth similar to the second approach. If you try to remove elements and then insert elements for each of the ten
points one at a time how do you manage the sweepline data structure such that the invariant that all its elements are in strict sorting order at all stages of execution?
This is the main idea of sweep line data structure, isn't it? We insert new items based on comparison functor. Or do you mean precision problems as mentioned further ("If you ever try to use the tree when its invariant that all elements are in strict sorting order is violated you will get undefined behavior and either crash or infinitely loop")? slope s1 = eval_slope_at_x(e1, x_); What does "slope" exactly mean? An alternative approach would be to make the comparison depend on the
current Y.
It should depend on the current Y in anyway, or am I wrong? Site events with different y-coordinates change sweep line data in different ways. Thanks, Andrii Sydorchuk
Andriy Sydorchuk wrote:
What is your definition of element? Does leaf and node elements contain the same data? I am thinking about element like (leaves and nodes are the same): struct SweepLineElement { private: const Point2d* point_left; const Point2d* point_right;
double getComparerValue(double line_x); }; So each node actually contains information about intersection point of two neighbour arcs from the sweep line, created by point_left and point_right at current X. Probably it's smth similar to the second approach.
Using double for the comparer value won't be robust. If the coordinate type is a built-in it is better to store points by value in you data structures because you save one level of pointer indirection and get better cache utilization. What happens when a new input event comes in with a point between point_left and point_right?
This is the main idea of sweep line data structure, isn't it? We insert new items based on comparison functor. Or do you mean precision problems as mentioned further ("If you ever try to use the tree when its invariant that all elements are in strict sorting order is violated you will get undefined behavior and either crash or infinitely loop")?
slope s1 = eval_slope_at_x(e1, x_);
What does "slope" exactly mean?
When working with line segments as I do I use comparison of the "rise over run" slope of the line segments to break ties when two line segments touch the sweepline at the same y value. When working with arcs you could use the tangent slope of arcs. Regards, Luke
Using double for the comparer value won't be robust. If the coordinate type is a built-in it is better to store points by value in you data structures because you save one level of pointer indirection and get better cache utilization. What happens when a new input event comes in with a point between point_left and point_right?
Well it depends on the type of input event: 1) Site event: a) Find, position where new node should be inserted, using lower_bound function with temporary node in the sweep line. For using lower_bound function we implement our NodeComparer: template <class T> struct NodeComparer { const double& line_y_; NodeComparer(double line_y) : line_y_(line_y) {} // Compares two nodes based on the parabolas interection points. bool operator() (const T& t1, const T& t2) const { return t1.getComparerValue(line_y) < t2.getComparerValue(line_y_); } }; getComparerValue may be implemented in such a way (but as you said there probably might be problems with precisions): // Get intersection point(x-coord) of parabolic arcs made by points from current node and sweep line at line_y. double SweepLineElement::getComparerValue(double line_y) const { // If right and left points have the same y-coordinate. if (DOUBLE_EQ(point_left->y, point_right->y)) return (point_left->x + point_right->x) / 2.0; // If left point and sweep line have the same y-coord return x-coord. if (DOUBLE_EQ(point_left->y, line_y)) return point_left->x; // If right point and sweep line have the same y-coord return x-coord. if (DOUBLE_EQ(point_right->y, line_y)) return point_right->x; // Parabola representation: y = a*X*X + b*X + c. // We get our parabola parameters based on: // (x - point->x)^2 + (y - point->y)^2 = (y - line_y)^2. // Left parabola parameters. double a1 = 0.5 / (point_left->y - line_y); double b1 = -a1 * 2.0 * point_left->x; double c1 = a1 * SQR(point_left->x) + point_left->y + line_y; // Right parabola parameters. double a2 = 0.5 / (point_right->y - line_y); double b2 = -a2 * 2.0 * point_right->x; double c2 = a2 * SQR(point_right->x) + point_right->y + line_y; // As left and right points don't have the same y-coords, there // will be two intersection points. We can find them from next condition: // a1*X^2 + b1*X + c1 = a2*X^2 + b2*X + c2 or // (a1 - a2)*X^2 + (b1 - b2)*X + c1 - c2 = 0. double a = a1 - a2; double b = b1 - b2; double c = c1 - c2; double sqrt_D = sqrt(SQR(b) - 4.0 * a * c); double x1 = (-b - sqrt_D) * 0.5 / a; double x2 = (-b + sqrt_D) * 0.5 / a; // We are interested only in intersection point x-coord that is created // by left parabola on the left and right parabola on the right. double deriv1 = 2.0 * a1 * x1 + b1; double deriv2 = 2.0 * a2 * x1 + b2; if(deriv1 > deriv2) return x1; return x2; } b) Insert two new nodes in our sweep line structure at appropriate postions, update circle events queue, update output data structures. 2) Circle event: a) erase nodes that made current circle event; b) insert new node, update circle events queue, update output data structures. What happens when a new input event comes in with a point between
point_left and point_right?
Basically we find a node that corresponds to (point_left, point_right). Find to which point intersection with new arc corresponds and insert to new arcs, like: (point_left, new_point), (new_point, point_left). When working with line segments as I do I use comparison of the "rise over
run" slope of the line segments to break ties when two line segments touch the sweepline at the same y value. When working with arcs you could use the tangent slope of arcs.
Probably the reason I misunderstood you at first, is that voronoi sweep line functionality differs a little bit from finding intersections of set of segment. -- Best regards, Andrii Sydorchuk
Andriy Sydorchuk wrote:
Using double for the comparer value won't be robust. If the coordinate type is a built-in it is better to store points by value in you data structures because you save one level of pointer indirection and get better cache utilization. What happens when a new input event comes in with a point between point_left and point_right?
...
Basically we find a node that corresponds to (point_left, point_right). Find to which point intersection with new arc corresponds and insert to new arcs, like: (point_left, new_point), (new_point, point_left).
That's what I thought. Why not use a point plus implicit parabola above and below as your sweepline element instead of point pair and implicit parabolas between? That way instead of two points your element contains only one and inserting a new point doesn't require the removal of any elements. Regards, Luke
That's what I thought. Why not use a point plus implicit parabola above and below as your sweepline element instead of point pair and implicit parabolas between? That way instead of two points your element contains only one and inserting a new point doesn't require the removal of any elements.
The reason is I read such kind of approach in this book: "Computational Geometry: Algorithms and Applicatoins" by* *Mark de Berg. Quota (page 156): "The beach line is represented by a balanced binary search tree T; it is the status structure. Its leaves correspond to the arcs of the beach line - which is x-monotone - in an ordered manner; the leftmost leaf represents the leftmost arc, the next leaf represents the second leftmost arc, and so on. Each leaf L stores the site that defines the arc it represents. The internal nodes of T represent the breakpoint on the beach line. A breakpoint is stored at an internal node by an ordered tuple of sites (pi, pj), where pi defines the parabola left of the breakpoint and pj defines the parabola to the right". I simplified this approach to be used without leaves, in this case all elements of the map will be the same structures. Why not use a point... Do you mean intersection point of above and below parabolas, or just site point? If the answer is intersection point, you may skip next questions. In case our element contains only one point, how would you implement insertion operation? For example we insert new element, how do we compare new element with elements from sweep line (we only know site point corresponding to the element, it seems to me that this information is not enough)? -- Best regards, Andrii Sydorchuk
Andriy Sydorchuk wrote:
That's what I thought. Why not use a point plus implicit parabola above and below as your sweepline element instead of point pair and implicit parabolas between? That way instead of two points your element contains only one and inserting a new point doesn't require the removal of any elements.
The reason is I read such kind of approach in this book: "Computational Geometry: Algorithms and Applicatoins" by* *Mark de Berg.
Quota (page 156): "The beach line is represented by a balanced binary search tree T; it is the status structure. Its leaves correspond to the arcs of the beach line - which is x-monotone - in an ordered manner; the leftmost leaf represents the leftmost arc, the next leaf represents the second leftmost arc, and so on. Each leaf L stores the site that defines the arc it represents. The internal nodes of T represent the breakpoint on the beach line. A breakpoint is stored at an internal node by an ordered tuple of sites (pi, pj), where pi defines the parabola left of the breakpoint and pj defines the parabola to the right". I simplified this approach to be used without leaves, in this case all elements of the map will be the same structures.
Why not use a point...
Do you mean intersection point of above and below parabolas, or just site point? If the answer is intersection point, you may skip next questions. In case our element contains only one point, how would you implement insertion operation? For example we insert new element, how do we compare new element with elements from sweep line (we only know site point corresponding to the element, it seems to me that this information is not enough)?
I see now that your element is actually a bisector of two sites, the moving point at which two parabolas intersect that traces line segments of the veronoi diagram as the sweepline advances. Regards, Luke
I see now that your element is actually a bisector of two sites, the moving point at which two parabolas intersect that traces line segments of the veronoi diagram as the sweepline advances.
Yes, you are right. Thanks for all the tips. I've already written my proposal. If you implement Medial Axis with sweepline that works on non-convex
polygons with holes you have implmented an algorithm that can also handle multiple non-overlapping polygons.
Ok, I see. At the moment I am looking for docs about Medial Axis sweepline algo implementation. Maybe you have some references? Best, Andrii Sydorchuk
Andriy Sydorchuk wrote:
I see now that your element is actually a bisector of two sites, the moving point at which two parabolas intersect that traces line segments of the veronoi diagram as the sweepline advances.
Yes, you are right. Thanks for all the tips. I've already written my proposal.
If you implement Medial Axis with sweepline that works on non-convex
polygons with holes you have implmented an algorithm that can also handle multiple non-overlapping polygons.
Ok, I see. At the moment I am looking for docs about Medial Axis sweepline algo implementation. Maybe you have some references?
I'm surprised that you have written your proposal without detailed discussion of medial axis. Here is a reference I was able to quickly find: http://www.springerlink.com/content/796w76637460t384/ Medial Axis, also known as straight skeleton, is similar to the veronoi diagram which considers line segments as sites. It is solved by a sweepline over bisectors and the algorithm is similar to Fortune's algorithm. Regards, Luke
Simonson, Lucanus J wrote:
Andriy Sydorchuk wrote:
Ok, I see. At the moment I am looking for docs about Medial Axis sweepline algo implementation. Maybe you have some references?
I'm surprised that you have written your proposal without detailed discussion of medial axis.
Here is a reference I was able to quickly find: http://www.springerlink.com/content/796w76637460t384/
Some more references: [CGAL] http://www.cgal.org/Manual/last/doc_html/cgal_manual/Straight_skeleton_2/Cha... [AA95] O. Aichholzer and F. Aurenhammer. Straight skeletons for general polygonal figures. Technical Report 432, Inst. for Theor. Comput. Sci., Graz Univ. of Technology, Graz, Austria, 1995. [AAAG95] O. Aichholzer, D. Alberts, F. Aurenhammer, and B. Gärtner. A novel type of skeleton for polygons. J. Universal Comput. Sci., 1(12):752-761, 1995. [LD03] R. G. Laycock and A. M. Day. Automatically generating roof models from building footprints. In The 11-th International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision'2003. Journal of WSCG - FULL Papers, volume 11, 2003. The given link is actually the reference [AA95].
Medial Axis, also known as straight skeleton, is similar to the veronoi diagram which considers line segments as sites. It is solved by a sweepline over bisectors and the algorithm is similar to Fortune's algorithm.
There are actually subtle differences between "Medial Axis", "straight skeleton" and the "voronoi diagram" of a polygon. A more detailed explanation of the differences can be found in the reference [CGAL], for example. However, from the comment "Medial Axis, also known as straight skeleton", I conclude that the project Luke has in mind is to compute the "straight skeleton", which he prefers to call "Medial Axis". This mixing of names is not uncommon, so you better get used to it. Regards, Thomas
I'm surprised that you have written your proposal without detailed discussion of medial axis.
It is my fault, but I've just returned from my internship in another country. That's why I had not enough time to discuss that and student application deadline is 9th of April (actually today). Probably we still can discuss it before end of the Interim Period (April 21). Some more references: [CGAL]
http://www.cgal.org/Manual/last/doc_html/cgal_manual/Straight_skeleton_2/Cha...
... The given link is actually the reference [AA95]. Thanks for the links. I actually started to look at http://www.cgal.org/Manual/last/doc_html/cgal_manual/Segment_Delaunay_graph_... .
Medial Axis, also known as straight skeleton, is similar to the veronoi
diagram which considers line segments as sites. It is solved by a
sweepline over bisectors and the algorithm is similar to Fortune's
algorithm.
There are actually subtle differences between "Medial Axis", "straight skeleton" and the "voronoi diagram" of a polygon. A more detailed explanation of the differences can be found in the reference [CGAL], for example. However, from the comment "Medial Axis, also known as straight skeleton", I conclude that the project Luke has in mind is to compute the "straight skeleton", which he prefers to call "Medial Axis". This mixing of names is not uncommon, so you better get used to it.
So do we need to compute straight skeleton (it consists only of straight line segments) or medial axis (may also include parabolic arcs)? Medial axis seems to be more similar to Voronoi diagrams (instead of sites input it uses edges as Luke said). So I am quite not sure if we should compute straight skeleton there.
Andriy Sydorchuk wrote:
I'm surprised that you have written your proposal without detailed discussion of medial axis.
It is my fault, but I've just returned from my internship in another country. That's why I had not enough time to discuss that and student application deadline is 9th of April (actually today). Probably we still can discuss it before end of the Interim Period (April 21).
Some more references:
[CGAL]
http://www.cgal.org/Manual/last/doc_html/cgal_manual/Straight_skeleton_2/Cha...
...
The given link is actually the reference [AA95].
Thanks for the links. I actually started to look at http://www.cgal.org/Manual/last/doc_html/cgal_manual/Segment_Delaunay_graph_... .
Medial Axis, also known as straight skeleton, is similar to the veronoi
diagram which considers line segments as sites. It is solved by a
sweepline over bisectors and the algorithm is similar to Fortune's
algorithm.
There are actually subtle differences between "Medial Axis", "straight skeleton" and the "voronoi diagram" of a polygon. A more detailed explanation of the differences can be found in the reference [CGAL], for example. However, from the comment "Medial Axis, also known as straight skeleton", I conclude that the project Luke has in mind is to compute the "straight skeleton", which he prefers to call "Medial Axis". This mixing of names is not uncommon, so you better get used to it.
So do we need to compute straight skeleton (it consists only of straight line segments) or medial axis (may also include parabolic arcs)? Medial axis seems to be more similar to Voronoi diagrams (instead of sites input it uses edges as Luke said). So I am quite not sure if we should compute straight skeleton there.
So I am learning along with the students. I'd like to thank Thomas for pointing out that my terminology has been imprecise. In answer to Andriy's question, it is actually straight skeleton that I intended to suggest for the GSOC project. I think it may be easier to focus on straight skeleton because if the output does not contain arcs then we don't need to provide interfaces and data structures for non-piecewise-linear geometry. That is a bigger issue than it may seem at first because of numerical issues inherent in the way arcs are chosen to be represented. On the other hand, if the student has a preference for implementing medial axis I would support that. Sorry for the confusion, Luke
So I am learning along with the students. I'd like to thank Thomas for pointing out that my terminology has been imprecise. In answer to Andriy's question, it is actually straight skeleton that I intended to suggest for the GSOC project. I think it may be easier to focus on straight skeleton because if the output does not contain arcs then we don't need to provide interfaces and data structures for non-piecewise-linear geometry. That is a bigger issue than it may seem at first because of numerical issues inherent in the way arcs are chosen to be represented. On the other hand, if the student has a preference for implementing medial axis I would support that
In any case during last few days I was trying to find how to solve both of this problems (medial axis and straight skeleton) with sweepline algorithm and haven't found anything (theory, articles). Are you sure that any of them can be solved with sweepline? Best, Andrii
Andriy Sydorchuk wrote:
In any case during last few days I was trying to find how to solve both of this problems (medial axis and straight skeleton) with sweepline algorithm and haven't found anything (theory, articles). Are you sure that any of them can be solved with sweepline?
Of course. It's the same argument as for the voronoi diagram. You sweep the line over the polygon, and nothing right of the sweepline enters into the "current virtual construction". So for a specific sweepline position, the virtual polygon that you have is the one that is the original polygon clipped by the sweepline, and your "current virtual construction" is probably the straight line skeletton (or the medial axis) of this virtual polygon. (When "virtually" moving the sweepline now, you know how the "current virtual construction" will change, if all the clipped sides just continue straight. But this is the easy (and fun) part of the problem, and the sweepline paradigm is probably clear anyway. It get's difficult when you start fighting with floating point arithmetic or integer coordinates rounding logic. Regards, Thomas
So for a specific sweepline position, the virtual polygon that you have is the one that is the original polygon clipped by the sweepline, and your "current virtual construction" is probably the straight line skeletton (or the medial axis) of this virtual polygon
Well this is general approach used to solve problem with sweepline algorithm. What I am asking about are algorithm details. For voronoi diagram there is exact prototype of algorithm (as example "An Efficient Implementation of Fortune’s Plane-Sweep Algorithm for Voronoi Diagrams" KennyWong, Hausi A. Muller, page 18). This prototype is definitely not the same for medial axis and straight skeleton problems. If it is common approach to solve this problems with sweepline, so there should be some articles and theory about that, but I haven't found them. Maybe I am expected to develop this prototype on my own, am I? (When "virtually" moving the sweepline now, you know how the "current
virtual construction" will change, if all the clipped sides just continue straight. But this is the easy (and fun) part of the problem, and the sweepline paradigm is probably clear anyway.
Well lets consider that I know "current virtual construction", but it is not so clear how to update it with moving further. For example how we should deal with edges that are not visible yet but have influence on the "current virtual construction". In the sweepline implementation for voronoi diagram our "current virtual construction" lies above the sweepline datastructure (union of arcs) and it is mathematically proven that it won't change with further algorithm execution. It get's difficult when you start fighting with floating point arithmetic
or integer coordinates rounding logic.
Yes, I understand that. But at first I'd like to be clear with all the algorithm details. Best regards, Andrii
Andriy Sydorchuk wrote:
Well lets consider that I know "current virtual construction", but it is not so clear how to update it with moving further. For example how we should deal with edges that are not visible yet but have influence on the "current virtual construction". In the sweepline implementation for voronoi diagram our "current virtual construction" lies above the sweepline datastructure (union of arcs) and it is mathematically proven that it won't change with further algorithm execution.
I didn't mean to deliver a complete algorithm description, I just gave a detailed answer to the question "Are you sure that any of them can be solved with sweepline?". The short answer would have just been "Yes", because I can be sure that they can be solved with sweepline, even if this would not be true.
For example how we should deal with edges that are not visible yet but have influence on the "current virtual construction".
I know that you wrote "for example", so you are probably not especially interested in an answer to this specific question, because the question only serves to illustrate your point that a detailed algorithm description would be handy, but I will answer anyway. The sweepline algorithm for the Voronoi diagram is called Fortune's algorithm. In Fortune's algorithm, there is both a "beach line" and a "sweep line". (The "beach line" is the set of points that have the same distance to the "sweep line" as the minimum distance to the points left of the "sweep line".) 1) The "current virtual construction" left of the "beach line" is already finished, so no points or edges left of the "sweep line" have any influence on it, even when the "sweep line" will move further to the right. 2) The "current virtual construction" right of the "beach line" and left of the "sweep line" is still empty, but this property isn't too important. But I guess it will also be satisfied by the sweepline algorithm for "straight skeleton" and "Medial Axis". 3) No points right of the "sweep line" have any influence on the "current virtual construction". This is the important property of any sweepline algorithm. So, the answer to the question is that edges that are not visible yet means edges that are right of the current sweepline position. By These edges don't have any influence on the "current virtual construction", by definition of the "current virtual construction". The analogon to the "beach line" for the "straight skeleton" or the "Medial Axis" should be fairly obvious. Other details of the algorithm might be less obvious, but none of these details will pose any fundamental difficulty for the sweepline paradigm. Regards, Thomas
For example how we should deal with edges that are not visible yet but have influence on the "current virtual construction".
I agree, this statement was wrong by the definition of "beach line" and "sweep line". Basically I was trying to find analog of "beach line" in case of "straight skeleton" and "Medial Axis". In case of "Medial Axis" it will be something similar to Fortune's "beach line". However it is not so clear what it is in case of "straight skeleton" problem. I'll do further research in the given area. Thanks, Andrii Sydorchuk
Andriy Sydorchuk wrote:
For example how we should deal with edges that are not visible yet but have influence on the "current virtual construction".
I agree, this statement was wrong by the definition of "beach line" and "sweep line". Basically I was trying to find analog of "beach line" in case of "straight skeleton" and "Medial Axis". In case of "Medial Axis" it will be something similar to Fortune's "beach line". However it is not so clear what it is in case of "straight skeleton" problem. I'll do further research in the given area.
You are probably right that there doesn't exist a "beach line" in case of the "straight skeleton" problem for an arbitrary non-convex polygon. For a reflective vertex with angle alpha, the actual distance of the spliting point can actually be up to "1 / sin (alpha/2)" further away from the "sweep line" than the direct distance to the "sweep line". Since the angle alpha can be arbitrary small, the influence range of these reflective vertices cannot be bounded, hence there doesn't exist a real "beach line" as required for the straightforward application of the sweepline paradigm. So whatever the references Luke found described, it probably wasn't the straightforward application of the sweepline paradigm to the "straight skeleton" problem, as this would not be directly applicable to polygons without a lower bound for the minimum reflective vertex angle. However, convex polygons, manhattan polygons and 45 degree angle polygons all have such a lower bound, so a "straightforward" sweepline algorithm for the "straight skeleton" problem might still make sense to have for Boost.Polygon. Let's hope I finally managed to create some confusion. Regards, Thomas
Thomas Klimpel wrote:
Andriy Sydorchuk wrote:
For example how we should deal with edges that are not visible yet but have influence on the "current virtual construction".
I agree, this statement was wrong by the definition of "beach line" and "sweep line". Basically I was trying to find analog of "beach line" in case of "straight skeleton" and "Medial Axis". In case of "Medial Axis" it will be something similar to Fortune's "beach line". However it is not so clear what it is in case of "straight skeleton" problem. I'll do further research in the given area.
You are probably right that there doesn't exist a "beach line" in case of the "straight skeleton" problem for an arbitrary non-convex polygon. For a reflective vertex with angle alpha, the actual distance of the spliting point can actually be up to "1 / sin (alpha/2)" further away from the "sweep line" than the direct distance to the "sweep line". Since the angle alpha can be arbitrary small, the influence range of these reflective vertices cannot be bounded, hence there doesn't exist a real "beach line" as required for the straightforward application of the sweepline paradigm.
So whatever the references Luke found described, it probably wasn't the straightforward application of the sweepline paradigm to the "straight skeleton" problem, as this would not be directly applicable to polygons without a lower bound for the minimum reflective vertex angle. However, convex polygons, manhattan polygons and 45 degree angle polygons all have such a lower bound, so a "straightforward" sweepline algorithm for the "straight skeleton" problem might still make sense to have for Boost.Polygon.
Let's hope I finally managed to create some confusion.
There has been some confusion from the beginning. What you have finally managed to create is doubt, which has lead me to recheck my assumptions. This is what I get for not doing the research *before* suggesting the project idea. In any case the Felkel algorithm for straight skeleton ( http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.40.5505&rep=rep1&ty... ) appears to be the one to implement according to Fernando, who ought to know. Felkel provides a step by step breakdown of the algorithm, which is more of a spiraling inward on closed polygons and shows how it can be applied to non-convex polygons with holes. This isn't a sweepline, it is a wavefront algorithm. If you imagine the roofline analogy the algorithm is rising water. I for some reason had in my mind a veronoi diagram -> medial axis -> straight skeleton progression of algorithm. It turns out that straight skeleton, despite being superficially similar to medial axis, is not solved the same way. So while I would like to see a straight skeleton implementation because it leads to a number of very nice polygon operations, straight skeleton may have to wait until next year or I may end up implementing it myself. Thanks Thomas, Luke
In any case the Felkel algorithm for straight skeleton ( http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.40.5505&rep=rep1&ty... ) appears to be the one to implement according to Fernando, who ought to know. Felkel provides a step by step breakdown of the algorithm, which is more of a spiraling inward on closed polygons and shows how it can be applied to non-convex polygons with holes. This isn't a sweepline, it is a wavefront algorithm. If you imagine the roofline analogy the algorithm is rising water.
This is exactly what I was looking for. Thanks a lot. I've already explored part connected with convex polygons, everything seems to be clear. I for some reason had in my mind a veronoi diagram -> medial axis ->
straight skeleton progression of algorithm. It turns out that straight skeleton, despite being superficially similar to medial axis, is not solved the same way. So while I would like to see a straight skeleton implementation because it leads to a number of very nice polygon operations, straight skeleton may have to wait until next year or I may end up implementing it myself
At the moment my proposal only includes Voronoi and Delaunay triangulation problems, as I removed everything about medial axis or straight skeleton before things would become clear. As I am not able to edit my proposal after 9th of April I'll add changes there: During GSoC I am ready to implement everything that I mentioned in my proposal ( http://socghop.appspot.com/gsoc/student_proposal/private/google/gsoc2010/asy...) + straight skeleton generic implementation (with visualization tool). After GSoC I am ready to continue our collaboration and implement medial axis algorithm also. Proposed Milestones and Schedule will look like this: Proposed Milestones and Schedule Present - end of May: design basic architecture, write test cases, read articles connected with the topic, learn Boost; start of June - end of June: finish Voronoi and Delaunay triangulation implementation; start of July - end of July: finish straight skeleton generic implementation; start of August - 20th of August: improve logic, make refactoring, work on performance, testing, finishing documentation, benchmark tests. If you expect something else, please let me know. Thanks Luke, Thomas! Andrii
Andriy Sydorchuk wrote:
During GSoC I am ready to implement everything that I mentioned in my proposal ( http://socghop.appspot.com/gsoc/student_proposal/private/google/gsoc201 0/asydorchuk/t127068181731)
I get a "You do not have the required role." Can you make your proposal public, so that also non-mentors can view it? Well, on the other hand, I shouldn't spend too much time into this, so it's probably good that I can't see it. Regards, Thomas
Hi Fernando, In any case the Felkel algorithm for straight skeleton (
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.40.5505&rep=rep1&ty... ) appears to be the one to implement according to Fernando, who ought to know.
Thanks for the reference to the Felkel's algorithm. Best, Andrii
Hi Thomas,
I get a "You do not have the required role." Can you make your proposal public, so that also non-mentors can view it? Well, on the other hand, I shouldn't spend too much time into this, so it's probably good that I can't see it.
It seems to me that I am not able to make it public, that's why I'll add it to this message. Personal Details - Name: Andrii Sydorchuk - College/University: Kyiv Taras Schevchenko National Universtiy - Course/Major: Cybernetics/Applied Mathematics - Degree Program (BS, MS, PhD, etc.): Master - Email: sydorchuk.andriy@gmail.com - Homepage: - - Availability: I plan to spend 6-8 hours per day (more if required). I've already started (reading articles), 20th of August. I don't think that any of this factors with affect my availability. Background Information Educational background: I've earned Bachelors Degree in Applied Mathematics (GPA 4.95). At the moment I am on the first course of my Master's program in Applied Mathematics. Courses taken: Algebra and Geometry, Mathematical Analysis, Discrete Mathematics, Programming, Differential Equations, Numeric Methods, Mathematical Logic and Algorithms Theory, Mathematica Modeling, etc. Programming background: 1) January 2010 - April 2010 - Google Company, Software Engineer Intern (Zurich, Switzerland). Responsible for managing workflow pipelines. 2) April 2009 - December 2009 - Scopicsoftware Company, C++ Algorithm Software Developer. Responsible for 2D/3D applications development. Programming interests: Algorithms, Data Structures, Computational Geometry; Software competitions: 2010 - Internal Google Code Jam 2010 - 7th place. 2009 - All Ukrainian Collegiate Programming Contest 2009 - 3d diploma. 2008 - Advanced to Google Code Jam 2008 semifinal (London). Languages, technologies and tools: C++, C#, Python, STL, MFC, Qt, OpenGL; 3D rendering; Interest in contribution to the Boost C++ libraries: become more familiar with Boost libraries; make contribution to the library used by lots of people; Interest in the project: I am interested in algorithms, data structures, analytical geometry. And all of them are included as part of Sweepline project. But the main idea is to become more familiar with Boost library (not only part I'll work on). Previous work in this area (Scopicsoftware): 1) Making improvement to the algorithm of finding center lines of standard fonts, however this project used wave algorithm approach. 2) Developing 3D CAD system for creating office cubicles. Implementation of Collada[1] 3D models importer. Adding support of YafaRay raytracing engine[2]. Development and improvement of algorithms for cubicles clustering and objects arrangement. Plans beyond Summer of Code: 1) continue implementation of computational geometry algorithms in Boost. 2) become more familiar with other parts of Boost library and work on them also. C++ (4.2) C++ Standard Library (4.3) Boost C++ Library (0.0, starting to learn) Subversion (2.0, only SVN) Software development environments: Visual Studio, emacs. Documentation tool: I haven't used any of them (will study one of them during April - May). Project Proposal Implement generic and robust sweepline algorithm for solving 2D Voronoi diagram and Delaunay triangulation problems. I will also implement Sweep Line algorithm visualization tool using OpenGL. This will slightly help during testing. Basically it will visualize all the execution process: site events, circle events, adding new sites to the beach line, removing sites and so on. During last month I'll implement benchmark tests. It is also a good idea to start with writing tests and add new cases during development. Basically implementation of sweepline requires next data structures: 1) Event queue (priority queue). We can split site events and circle events into two data structures[3]: a) site events may be stored in stl vector. Which will be sorted after reading input data. Complexity of sorting n*log(n), where n - is number of input sites. b) circle events may be stored in stl priority_queue. Adding new element, lookup, deletion will require log(n) time, where n is number of elements in the priority queue. 2) Sweepline (self-balanced binary tree). For sweepline it is good choice to use wrapper class (adapter pattern) of stl map or set. It is better to use map, because we can split keys and associated data (storing whatever our problem requires). For elements we can choose bisectores of two neihbour sites[4]. We will also need to define robust key comparer[5], as there might be precision problems. 3) Output data structures (half edges, double linked lists, etc.). I'll do further research in this area. Next degeneracies will require careful handling: 1) 2 or more events on a common vertical line; 2) event points coincide, e.g. 4 or more co-circular points; 3) a site coincides with event; Proposed Milestones and Schedule Present - end of May: design basic architecture, write test cases, read articles connected with the topic, learn Boost; start of June - mid of July: finish Voronoi and Delaunay triangulation implementation; mid of July - 20th of August: improve logic, make refactoring, work on performance, testing, finishing documentation, benchmark tests. *References* [1] - http://code.google.com/p/opencollada/ [2] - http://www.yafaray.org/ [3] - Kenny Wong, Hausi A. Muller. An Efficient Implementation of Fortune's Plane-Sweep Algorithm for Voronoi Diagrams. [4] - Mark de Berg. Computational Geometry: Algorithms and Applicatoins. Page.156. [5] - Ruud, B. Building a Mutable Set, 2003
Hi Luke,
I for some reason had in my mind a veronoi diagram -> medial axis -> straight skeleton progression of algorithm.
In fact, Aichholzer's original paper specially argues how the straight skeleton cosntruction cannot be derived as a voronoi-diagram-like algorithm because of the non-locality of the interactions caused by reflex vertices. OTOH...
It turns out that straight skeleton, despite being superficially similar to medial axis, is not solved the same way.
For a convex polygon, where there are no reflex vertices, the procedure is the same, and in fact both structures are exactly equivalent. Best -- Fernando Cacciola SciSoft Consulting, Founder http://www.scisoft-consulting.com
Fernando Cacciola wrote:
Hi Luke,
I for some reason had in my mind a veronoi diagram -> medial axis -> straight skeleton progression of algorithm.
In fact, Aichholzer's original paper specially argues how the straight skeleton cosntruction cannot be derived as a voronoi-diagram-like algorithm because of the non-locality of the interactions caused by reflex vertices.
I'm not sure that it is impossible to use a modification of Fortune's algorithm for the straight skeleton construction, all I admitted was that a straightforward application of the sweepline paradigm would not work (and that there would be no "beach line" that clearly marks the parts of the construction that are already finished and won't change anymore when the "sweep line" advances). Especially when using Fortune's algorithm, the reflex vertices could be brought in with a larger angle (such that one of the sides would be parallel to the sweepline), and then shrunk down to the real angle. Of course, this process would make the "beach line" go backwards, and discart some part of the already existing construction. This is not nice, and probably also destroys the prof O(n log(n)) complexity of the algorithm. However, it would still use the sweepline paradigm, and would be identical to Fortune's algorithm in case the polygon is convex. Regards, Thomas
Andriy Sydorchuk wrote:
Well this is general approach used to solve problem with sweepline algorithm. What I am asking about are algorithm details. For voronoi diagram there is exact prototype of algorithm (as example "An Efficient Implementation of Fortune's Plane-Sweep Algorithm for Voronoi Diagrams" KennyWong, Hausi A. Muller, page 18). This prototype is definitely not the same for medial axis and straight skeleton problems. If it is common approach to solve this problems with sweepline, so there should be some articles and theory about that, but I haven't found them. Maybe I am expected to develop this prototype on my own, am I?
There are papers that you can read. At a certain point you may need take a leap of intuition from the description of the algorithm to an understanding of the algorithm, and a well written paper should at least get you to the point where you are capable of making that leap. I like to use google scholar for literature searches: http://scholar.google.com/scholar?hl=en&q=straight+skeleton+sweep+line&btnG=... The first two hits appear spot on to me after skimming them: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.40.5505&rep=rep1&ty... http://www.springerlink.com/index/796W76637460T384.pdf Look at the references of these papers and papers citing them (using features of google scholar or equivalent literature search engines) to find similar papers. Often it is a paper that is written later that contributes very little itself that will explain things clearest or help you understand something that was confusing in the original publication of an algorithm. Early papers tend to focus on proofs while later papers focus on implementation details and practical considerations. Regards, Luke
Hi Andriy,
So do we need to compute straight skeleton (it consists only of straight line segments) or medial axis (may also include parabolic arcs)? Medial axis seems to be more similar to Voronoi diagrams (instead of sites input it uses edges as Luke said). So I am quite not sure if we should compute straight skeleton there.
There is one major thing that is not evident in most drawn samples of a straight skeleton: it might very well be anything but "medial" if the polygon is non-convex. That is, in a non-convex polygon, the inner bisectors might be significantly pushed toward one side of the polygon. Thus, a medial axis is a subset of a straight skeleton ONLY if the polygon is convex. If not, the medial axis is a subset of the voronoi diagram instead. Best -- Fernando Cacciola SciSoft Consulting, Founder http://www.scisoft-consulting.com
participants (4)
-
Andriy Sydorchuk -
Fernando Cacciola -
Simonson, Lucanus J -
Thomas Klimpel