2010-02-27 52 views
6

我使用Fortune算法來查找一組點的Voronoi圖。 我回來的是一個線段列表,但我需要知道哪些線段構成了封閉多邊形,並將它們放在一個由它們圍繞的原點散列的對象中。獲得由Voronoi線段形成的凸多邊形集合的最快方法

什麼可能是最快的方式來找到這些? 我應該從算法中保存一些關鍵信息嗎?如果是這樣?

這是我實現在Java中財富的算法從C++實現here

class Voronoi { 

// The set of points that control the centers of the cells 
private LinkedList<Point> pts; 
// A list of line segments that defines where the cells are divided 
private LinkedList<Edge> output; 
// The sites that have not yet been processed, in acending order of X coordinate 
private PriorityQueue sites; 
// Possible upcoming cirlce events in acending order of X coordinate 
private PriorityQueue events; 
// The root of the binary search tree of the parabolic wave front 
private Arc root; 

void runFortune(LinkedList pts) { 

    sites.clear(); 
    events.clear(); 
    output.clear(); 
    root = null; 

    Point p; 
    ListIterator i = pts.listIterator(0); 
    while (i.hasNext()) { 
     sites.offer(i.next()); 
    } 

    // Process the queues; select the top element with smaller x coordinate. 
    while (sites.size() > 0) { 
     if ((events.size() > 0) && ((((CircleEvent) events.peek()).xpos) <= (((Point) sites.peek()).x))) { 
      processCircleEvent((CircleEvent) events.poll()); 
     } else { 
      //process a site event by adding a curve to the parabolic front 
      frontInsert((Point) sites.poll()); 
     } 
    } 

    // After all points are processed, do the remaining circle events. 
    while (events.size() > 0) { 
     processCircleEvent((CircleEvent) events.poll()); 
    } 

    // Clean up dangling edges. 
    finishEdges(); 

} 

private void processCircleEvent(CircleEvent event) { 
    if (event.valid) { 
     //start a new edge 
     Edge edgy = new Edge(event.p); 

     // Remove the associated arc from the front. 
     Arc parc = event.a; 
     if (parc.prev != null) { 
      parc.prev.next = parc.next; 
      parc.prev.edge1 = edgy; 
     } 
     if (parc.next != null) { 
      parc.next.prev = parc.prev; 
      parc.next.edge0 = edgy; 
     } 

     // Finish the edges before and after this arc. 
     if (parc.edge0 != null) { 
      parc.edge0.finish(event.p); 
     } 
     if (parc.edge1 != null) { 
      parc.edge1.finish(event.p); 
     } 

     // Recheck circle events on either side of p: 
     if (parc.prev != null) { 
      checkCircleEvent(parc.prev, event.xpos); 
     } 
     if (parc.next != null) { 
      checkCircleEvent(parc.next, event.xpos); 
     } 

    } 
} 

void frontInsert(Point focus) { 
    if (root == null) { 
     root = new Arc(focus); 
     return; 
    } 

    Arc parc = root; 
    while (parc != null) { 
     CircleResultPack rez = intersect(focus, parc); 
     if (rez.valid) { 
      // New parabola intersects parc. If necessary, duplicate parc. 

      if (parc.next != null) { 
       CircleResultPack rezz = intersect(focus, parc.next); 
       if (!rezz.valid){ 
        Arc bla = new Arc(parc.focus); 
        bla.prev = parc; 
        bla.next = parc.next; 
        parc.next.prev = bla; 
        parc.next = bla; 
       } 
      } else { 
       parc.next = new Arc(parc.focus); 
       parc.next.prev = parc; 
      } 
      parc.next.edge1 = parc.edge1; 

      // Add new arc between parc and parc.next. 
      Arc bla = new Arc(focus); 
      bla.prev = parc; 
      bla.next = parc.next; 
      parc.next.prev = bla; 
      parc.next = bla; 

      parc = parc.next; // Now parc points to the new arc. 

      // Add new half-edges connected to parc's endpoints. 
      parc.edge0 = new Edge(rez.center); 
      parc.prev.edge1 = parc.edge0; 
      parc.edge1 = new Edge(rez.center); 
      parc.next.edge0 = parc.edge1; 

      // Check for new circle events around the new arc: 
      checkCircleEvent(parc, focus.x); 
      checkCircleEvent(parc.prev, focus.x); 
      checkCircleEvent(parc.next, focus.x); 

      return; 
     } 

     //proceed to next arc 
     parc = parc.next; 
    } 

    // Special case: If p never intersects an arc, append it to the list. 
    parc = root; 
    while (parc.next != null) { 
     parc = parc.next; // Find the last node. 
    } 
    parc.next = new Arc(focus); 
    parc.next.prev = parc; 
    Point start = new Point(0, (parc.next.focus.y + parc.focus.y)/2); 
    parc.next.edge0 = new Edge(start); 
    parc.edge1 = parc.next.edge0; 

} 

void checkCircleEvent(Arc parc, double xpos) { 
    // Invalidate any old event. 
    if ((parc.event != null) && (parc.event.xpos != xpos)) { 
     parc.event.valid = false; 
    } 
    parc.event = null; 

    if ((parc.prev == null) || (parc.next == null)) { 
     return; 
    } 

    CircleResultPack result = circle(parc.prev.focus, parc.focus, parc.next.focus); 
    if (result.valid && result.rightmostX > xpos) { 
     // Create new event. 
     parc.event = new CircleEvent(result.rightmostX, result.center, parc); 
     events.offer(parc.event); 
    } 

} 

// Find the rightmost point on the circle through a,b,c. 
CircleResultPack circle(Point a, Point b, Point c) { 
    CircleResultPack result = new CircleResultPack(); 

    // Check that bc is a "right turn" from ab. 
    if ((b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y) > 0) { 
     result.valid = false; 
     return result; 
    } 

    // Algorithm from O'Rourke 2ed p. 189. 
    double A = b.x - a.x; 
    double B = b.y - a.y; 
    double C = c.x - a.x; 
    double D = c.y - a.y; 
    double E = A * (a.x + b.x) + B * (a.y + b.y); 
    double F = C * (a.x + c.x) + D * (a.y + c.y); 
    double G = 2 * (A * (c.y - b.y) - B * (c.x - b.x)); 

    if (G == 0) { // Points are co-linear. 
     result.valid = false; 
     return result; 
    } 

    // centerpoint of the circle. 
    Point o = new Point((D * E - B * F)/G, (A * F - C * E)/G); 
    result.center = o; 

    // o.x plus radius equals max x coordinate. 
    result.rightmostX = o.x + Math.sqrt(Math.pow(a.x - o.x, 2.0) + Math.pow(a.y - o.y, 2.0)); 

    result.valid = true; 
    return result; 
} 

// Will a new parabola at point p intersect with arc i? 
CircleResultPack intersect(Point p, Arc i) { 
    CircleResultPack res = new CircleResultPack(); 
    res.valid = false; 
    if (i.focus.x == p.x) { 
     return res; 
    } 

    double a = 0.0; 
    double b = 0.0; 
    if (i.prev != null) // Get the intersection of i->prev, i. 
    { 
     a = intersection(i.prev.focus, i.focus, p.x).y; 
    } 
    if (i.next != null) // Get the intersection of i->next, i. 
    { 
     b = intersection(i.focus, i.next.focus, p.x).y; 
    } 

    if ((i.prev == null || a <= p.y) && (i.next == null || p.y <= b)) { 
     res.center = new Point(0, p.y); 

     // Plug it back into the parabola equation to get the x coordinate 
     res.center.x = (i.focus.x * i.focus.x + (i.focus.y - res.center.y) * (i.focus.y - res.center.y) - p.x * p.x)/(2 * i.focus.x - 2 * p.x); 

     res.valid = true; 
     return res; 
    } 
    return res; 
} 

// Where do two parabolas intersect? 
Point intersection(Point p0, Point p1, double l) { 
    Point res = new Point(0, 0); 
    Point p = p0; 

    if (p0.x == p1.x) { 
     res.y = (p0.y + p1.y)/2; 
    } else if (p1.x == l) { 
     res.y = p1.y; 
    } else if (p0.x == l) { 
     res.y = p0.y; 
     p = p1; 
    } else { 
     // Use the quadratic formula. 
     double z0 = 2 * (p0.x - l); 
     double z1 = 2 * (p1.x - l); 

     double a = 1/z0 - 1/z1; 
     double b = -2 * (p0.y/z0 - p1.y/z1); 
     double c = (p0.y * p0.y + p0.x * p0.x - l * l)/z0 - (p1.y * p1.y + p1.x * p1.x - l * l)/z1; 

     res.y = (-b - Math.sqrt((b * b - 4 * a * c)))/(2 * a); 
    } 
    // Plug back into one of the parabola equations. 
    res.x = (p.x * p.x + (p.y - res.y) * (p.y - res.y) - l * l)/(2 * p.x - 2 * l); 
    return res; 
} 

void finishEdges() { 
    // Advance the sweep line so no parabolas can cross the bounding box. 
    double l = gfx.width * 2 + gfx.height; 

    // Extend each remaining segment to the new parabola intersections. 
    Arc i = root; 
    while (i != null) { 
     if (i.edge1 != null) { 
      i.edge1.finish(intersection(i.focus, i.next.focus, l * 2)); 
     } 
     i = i.next; 
    } 
} 

class Point implements Comparable<Point> { 

    public double x, y; 
    //public Point goal; 

    public Point(double X, double Y) { 
     x = X; 
     y = Y; 
    } 

    public int compareTo(Point foo) { 
     return ((Double) this.x).compareTo((Double) foo.x); 
    } 
} 

class CircleEvent implements Comparable<CircleEvent> { 

    public double xpos; 
    public Point p; 
    public Arc a; 
    public boolean valid; 

    public CircleEvent(double X, Point P, Arc A) { 
     xpos = X; 
     a = A; 
     p = P; 
     valid = true; 
    } 

    public int compareTo(CircleEvent foo) { 
     return ((Double) this.xpos).compareTo((Double) foo.xpos); 
    } 
} 

class Edge { 

    public Point start, end; 
    public boolean done; 

    public Edge(Point p) { 
     start = p; 
     end = new Point(0, 0); 
     done = false; 
     output.add(this); 
    } 

    public void finish(Point p) { 
     if (done) { 
      return; 
     } 
     end = p; 
     done = true; 
    } 
} 

class Arc { 
    //parabolic arc is the set of points eqadistant from a focus point and the beach line 

    public Point focus; 
    //these object exsit in a linked list 
    public Arc next, prev; 
    // 
    public CircleEvent event; 
    // 
    public Edge edge0, edge1; 

    public Arc(Point p) { 
     focus = p; 
     next = null; 
     prev = null; 
     event = null; 
     edge0 = null; 
     edge1 = null; 
    } 
} 

class CircleResultPack { 
    // stupid Java doesnt let me return multiple variables without doing this 

    public boolean valid; 
    public Point center; 
    public double rightmostX; 
} 
} 

(我知道這不會編譯,數據結構需要被初始化,它的失蹤進口)

移植我想要的是:

LinkedList<Poly> polys; 
//contains all polygons created by Voronoi edges 

class Poly { 
    //defines a single polygon 
    public Point locus; 
    public LinkedList<Points> verts; 
} 

我能想到的最直接的蠻力方法是創建一個無向圖中點的圖形(邊的端點),每個點有單個條目,點之間每個邊的單個連接(不重複),然後去查找此圖中的所有循環,然後爲每個共享3個或更多點的一組循環,扔掉一切,但最短的循環。然而,這會太慢。

回答

4

Voronoi圖的對偶是Delaunay三角剖分。這意味着Voroni圖上的每個頂點都連接到三個邊 - 意味着每個頂點都屬於三個區域。因爲只有3每個頂點段

for each vertex in Voronoi Diagram 
    for each segment next to this point 
     "walk around the perimeter" (just keep going counter-clockwise) 
        until you get back to the starting vertex 

這應該是O(N)

我的算法來使用,這將是。你還必須做一些記賬工作,以確保你不會做同樣的地區兩次(一個簡單的方法是保持每一個傳出邊緣的bool,當你走路,標記它),並記住點在無窮遠處,但這個想法應該足夠了。

+0

這看起來像我需要做的找到多邊形, 但是,我怎麼把多邊形與它們圍繞的點(原始輸入集合中的那些)相關聯 – Nathan

+0

沒有在你的代碼中看得太深,你應該保持在拋物線中追蹤它們,當它們創建硬邊時,將邊與點關聯起來。 – Larry