Jump to content

File:Bowyer-Watson 0.png

Page contents not supported in other languages.
This is a file from the Wikimedia Commons
From Wikipedia, the free encyclopedia

Original file (1,200 × 900 pixels, file size: 62 KB, MIME type: image/png)

Summary

Description
Français : 1e étape intermédiaire de l'algorithme de Bowyer-Watson, construisant une triangulation de Delaunay sur 5 points dans le plan euclidien.

Légende :

  • points rouge : nœud ajouté à cette étape,
  • points noirs : nœuds déjà présents,
  • points gris : sommets du triangle englobant,
  • traits bleus : triangulation existante,
  • traits rouges : enveloppe englobante du nœud inséré dans la triangulation existante,
  • traits bleus pointillés : arrêtes supprimées à cette itération,
  • traits verts : arrêtes crées à cette itération,
  • cercles gris : cercles circonscrits ne contenant pas le nœud inséré,
  • cercles jaunes : cercles circonscrits contenant le nœud inséré,
  • rectangle magenta : boite englobante des nœuds à considérer.
Date
Source Own work (see source code below)
Author Nojhan

Licensing

w:en:Creative Commons
attribution share alike
This file is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license.
Attribution: © Johann Dréo — CC-BY-SA
You are free:
  • to share – to copy, distribute and transmit the work
  • to remix – to adapt the work
Under the following conditions:
  • attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
  • share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
GNU head Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled GNU Free Documentation License.
CeCILL This work is free software; you can redistribute it or modify it under the terms of the CeCILL. The terms of the CeCILL license are available at www.cecill.info.

Source code

Python script to produce the graphics. Copy-paste in a file named "bw.py" and use it as "python bw.py" to output triangulate the same points, or "python bw.py 11" to triangulate 11 random points.

This code is released under the GPL license (version 3 or later).

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

import sys
import math
from itertools import ifilterfalse as filter_if_not

import matplotlib.pyplot as plot
from matplotlib.path import Path
import matplotlib.patches as patches



def LOG( *args ):
    """Print something on stderr and flush"""
    for msg in args:
        sys.stderr.write( str(msg) )
        sys.stderr.write(" ")
    sys.stderr.flush()

def LOGN( *args ):
    """Print something on stderr, with a trailing new line, and flush"""
    LOG( *args )
    LOG("\n")

def tour(lst):
    # consecutive pairs in lst  + last-to-first element
    for a,b in zip(lst, lst[1:] + [lst[0]]):
        yield (a,b)

def plot_segments( ax, segments, **kwargs ):
    for start,end in segments:
        verts = [start,end,(0,0)]
        codes = [Path.MOVETO,Path.LINETO,Path.STOP]
        path = Path(verts, codes)
        patch = patches.PathPatch(path, facecolor='none', **kwargs )
        ax.add_patch(patch)


def scatter_segments( ax, segments, **kwargs  ):
    xy = [ ((i[0],j[0]),(i[1],j[1])) for (i,j) in segments]
    x = [i[0] for i in xy]
    y = [i[1] for i in xy]
    ax.scatter( x,y, s=20, marker='o', **kwargs)



# Based on http://paulbourke.net/papers/triangulate/
# Efficient Triangulation Algorithm Suitable for Terrain Modelling
#     An Algorithm for Interpolating Irregularly-Spaced Data
#     with Applications in Terrain Modelling
# Written by Paul Bourke
# Presented at Pan Pacific Computer Conference, Beijing, China.
# January 1989

def x( point ):
    return point[0]

def y( point ):
    return point[1]

def mid( xy, pa, pb ):
    return ( xy(pa) + xy(pb) ) / 2.0

def middle( pa, pb ):
    return mid(x,pa,pb),mid(y,pa,pb)

def mtan( pa, pb ):
    return -1 * ( x(pa) - x(pb) ) / ( y(pa) - y(pb) )


class CoincidentPointsError(Exception):
    """Coincident points"""
    pass

def circumcircle( triangle, epsilon = sys.float_info.epsilon ):
    """Compute the circumscribed circle of a triangle and 
    Return a 2-tuple: ( (center_x, center_y), radius )"""

    assert( len(triangle) == 3 )
    p0,p1,p2 = triangle
    assert( len(p0) == 2 )
    assert( len(p1) == 2 )
    assert( len(p2) == 2 )

    dy01 = abs( y(p0) - y(p1) )
    dy12 = abs( y(p1) - y(p2) )

    if dy01 < epsilon and dy12 < epsilon:
        # coincident points
        raise CoincidentPointsError

    elif dy01 < epsilon:
        m12 = mtan( p2,p1 )
        mx12,my12 = middle( p1, p2 )
        cx = mid( x, p1, p0 )
        cy = m12 * (cx - mx12) + my12

    elif dy12 < epsilon:
        m01 = mtan( p1, p0 )
        mx01,my01 = middle( p0, p1 )
        cx = mid( x, p2, p1 )
        cy = m01 * ( cx - mx01 ) + my01

    else:
        m01 =  mtan( p1, p0 )
        m12 =  mtan( p2, p1 )
        mx01,my01 = middle( p0, p1 )
        mx12,my12 = middle( p1, p2 )
        cx = ( m01 * mx01 - m12 * mx12 + my12 - my01 ) / ( m01 - m12 )
        if dy01 > dy12:
            cy = m01 * ( cx - mx01 ) + my01
        else:
            cy = m12 * ( cx - mx12 ) + my12

    dx1 = x(p1) - cx
    dy1 = y(p1) - cy
    r = math.sqrt(dx1**2 + dy1**2)

    return (cx,cy),r


def in_circle( p, center, radius, epsilon  = sys.float_info.epsilon ):
    """Return True if the given point p is in the given circle"""

    assert( len(p) == 2 )
    cx,cy = center

    dxp = x(p) - cx
    dyp = y(p) - cy
    dr = math.sqrt(dxp**2 + dyp**2)

    if (dr - radius) <= epsilon:
        return True
    else:
        return False


def in_circumcircle( p, triangle, epsilon  = sys.float_info.epsilon ):
    """Return True if the given point p is in the circumscribe circle of the given triangle"""

    assert( len(p) == 2 )
    (cx,cy),r = circumcircle( triangle, epsilon )
    return in_circle( p, (cx,cy), r, epsilon )


def bounds( vertices ):
    """Return the iso-axis rectangle enclosing the given points"""
    # find vertices set bounds
    xmin = x(vertices[0])
    ymin = y(vertices[0])
    xmax = xmin
    ymax = ymin

    # we do not use min(vertices,key=x) because it would iterate 4 times over the list, instead of just one
    for v in vertices:
        xmin = min(x(v),xmin)
        xmax = max(x(v),xmax)
        ymin = min(y(v),ymin)
        ymax = max(y(v),ymax)
    return (xmin,ymin),(xmax,ymax)


def edges_of( triangulation ):
    """Return a list containing the edges of the given list of 3-tuples of points"""
    edges = []
    for t in triangulation:
        for e in utils.tour(list(t)):
            edges.append( e )
    return edges


def supertriangle( vertices, delta = 0.1 ):
    """Return a super-triangle that encloses all given points.
    The super-triangle has its base at the bottom and encloses the bounding box at a distance given by:
        delta*max(width,height)
    """

    # Iso-rectangle bounding box.
    (xmin,ymin),(xmax,ymax) = bounds( vertices )

    dx = xmax - xmin
    dy = ymax - ymin
    dmax = max( dx, dy )
    xmid = (xmax + xmin) / 2.0

    supertri = ( ( xmin-dy-dmax*delta, ymin-dmax*delta ),
                 ( xmax+dy+dmax*delta, ymin-dmax*delta ),
                 ( xmid              , ymax+(xmax-xmid)+dmax*delta ) )

    return supertri


def delaunay_bowyer_watson( points, supertri = None, superdelta = 0.1, epsilon = sys.float_info.epsilon,
        do_plot = None, plot_filename = "Bowyer-Watson_%i.png" ):
    """Return the Delaunay triangulation of the given points

    epsilon: used for floating point comparisons, two points are considered equals if their distance is < epsilon.
    do_plot: if not None, plot intermediate steps on this matplotlib object and save them as images named: plot_filename % i
    """

    if do_plot and len(points) > 10:
        print "WARNING it is a bad idea to plot each steps of a triangulation of many points"

    # Sort points first on the x-axis, then on the y-axis.
    vertices = sorted( points )

    # LOGN( "super-triangle",supertri )
    if not supertri:
        supertri = supertriangle( vertices, superdelta )

    # It is the first triangle of the list.
    triangles = [ supertri ]

    completed = { supertri: False }

    # The predicate returns true if at least one of the vertices
    # is also found in the supertriangle.
    def match_supertriangle( tri ):
        if tri[0] in supertri or \
           tri[1] in supertri or \
           tri[2] in supertri:
            return True

    # Returns the base of each plots, with points, current triangulation, super-triangle and bounding box.
    def plot_base(ax,vi = len(vertices), vertex = None):
        ax.set_aspect('equal')
        # regular points
        scatter_x = [ p[0] for p in vertices[:vi]]
        scatter_y = [ p[1] for p in vertices[:vi]]
        ax.scatter( scatter_x,scatter_y, s=30, marker='o', facecolor="black")
        # super-triangle vertices
        scatter_x = [ p[0] for p in list(supertri)]
        scatter_y = [ p[1] for p in list(supertri)]
        ax.scatter( scatter_x,scatter_y, s=30, marker='o', facecolor="lightgrey", edgecolor="lightgrey")
        # current vertex
        if vertex:
            ax.scatter( vertex[0],vertex[1], s=30, marker='o', facecolor="red", edgecolor="red")
        # current triangulation
        plot_segments( ax, edges_of(triangles), edgecolor = "blue", alpha=0.5, linestyle='solid' )
        # bounding box
        (xmin,ymin),(xmax,ymax) = bounds(vertices)
        plot_segments( ax, tour([(xmin,ymin),(xmin,ymax),(xmax,ymax),(xmax,ymin)]), edgecolor = "magenta", alpha=0.2, linestyle='dotted' )


    # Insert vertices one by one.
    LOG("Insert vertices: ")
    if do_plot:
        it=0
    for vi,vertex in enumerate(vertices):
        # LOGN( "\tvertex",vertex )
        assert( len(vertex) == 2 )

        if do_plot:
            ax = do_plot.add_subplot(111)
            plot_base(ax,vi,vertex)

        # All the triangles whose circumcircle encloses the point to be added are identified,
        # the outside edges of those triangles form an enclosing polygon.

        # Forget previous candidate polygon's edges.
        enclosing = []

        removed = []
        for triangle in triangles:
            # LOGN( "\t\ttriangle",triangle )
            assert( len(triangle) == 3 )

            # Do not consider triangles already tested.
            # If completed has a key, test it, else return False.
            if completed.get( triangle, False ):
                # LOGN( "\t\t\tAlready completed" )
                # if do_plot:
                    # plot_segments( ax, tour(list(triangle)), edgecolor = "magenta", alpha=1, lw=1, linestyle='dotted' )
                continue

            # LOGN( "\t\t\tCircumcircle" ) 
            assert( triangle[0] != triangle[1] and triangle[1] != triangle [2] and triangle[2] != triangle[0] )
            center,radius = circumcircle( triangle, epsilon )

            # If it match Delaunay's conditions.
            if x(center) < x(vertex) and math.sqrt((x(vertex)-x(center))**2) > radius:
                # LOGN( "\t\t\tMatch Delaunay, mark as completed" ) 
                completed[triangle] = True

            # If the current vertex is inside the circumscribe circle of the current triangle,
            # add the current triangle's edges to the candidate polygon.
            if in_circle( vertex, center, radius, epsilon ):
                # LOGN( "\t\t\tIn circumcircle, add to enclosing polygon",triangle )
                if do_plot:
                    circ = plot.Circle(center, radius, facecolor='yellow', edgecolor="orange", alpha=0.2, clip_on=False)
                    ax.add_patch(circ)

                for p0,p1 in tour(list(triangle)):
                    # Then add this edge to the polygon enclosing the vertex,
                    enclosing.append( (p0,p1) )
                # and remove the corresponding triangle from the current triangulation.
                removed.append( triangle )
                completed.pop(triangle,None)

            elif do_plot:
                circ = plot.Circle(center, radius, facecolor='lightgrey', edgecolor="grey", alpha=0.2, clip_on=False)
                ax.add_patch(circ)

        # end for triangle in triangles

        # The triangles in the enclosing polygon are deleted and
        # new triangles are formed between the point to be added and
        # each outside edge of the enclosing polygon. 

        # Actually remove triangles.
        for triangle in removed:
            triangles.remove(triangle)


        # Remove duplicated edges.
        # This leaves the edges of the enclosing polygon only,
        # because enclosing edges are only in a single triangle,
        # but edges inside the polygon are at least in two triangles.
        hull = []
        for i,(p0,p1) in enumerate(enclosing):
            # Clockwise edges can only be in the remaining part of the list.
            # Search for counter-clockwise edges as well.
            if (p0,p1) not in enclosing[i+1:] and (p1,p0) not in enclosing:
                hull.append((p0,p1))
            elif do_plot:
                plot_segments( ax, [(p0,p1)], edgecolor = "white", alpha=1, lw=1, linestyle='dotted' )



        if do_plot:
            plot_segments( ax, hull, edgecolor = "red", alpha=1, lw=1, linestyle='solid' )


        # Create new triangles using the current vertex and the enclosing hull.
        # LOGN( "\t\tCreate new triangles" )
        for p0,p1 in hull:
            assert( p0 != p1 )
            triangle = tuple([p0,p1,vertex])
            # LOGN("\t\t\tNew triangle",triangle)
            triangles.append( triangle )
            completed[triangle] = False

            if do_plot:
                plot_segments( ax, [(p0,vertex),(p1,vertex)], edgecolor = "green", alpha=1, linestyle='solid' )

        if do_plot:
            plot.savefig( plot_filename % it, dpi=150)
            plot.clf()

            it+=1
        LOG(".")

    # end for vertex in vertices
    LOGN(" done")


    # Remove triangles that have at least one of the supertriangle vertices.
    # LOGN( "\tRemove super-triangles" ) 

    # Filter out elements for which the predicate is False,
    # here: *keep* elements that *do not* have a common vertex.
    # The filter is a generator, so we must make a list with it to actually get the data.
    triangulation = list(filter_if_not( match_supertriangle, triangles ))

    if do_plot:
            ax = do_plot.add_subplot(111)
            plot_base(ax)
            plot_segments( ax, edges_of(triangles), edgecolor = "red", alpha=0.5, linestyle='solid' )
            plot_segments( ax, edges_of(triangulation), edgecolor = "blue", alpha=1, linestyle='solid' )
            plot.savefig( plot_filename % it, dpi=150)
            plot.clf()

    return triangulation


if __name__ == "__main__":
    import random
    import utils
    import matplotlib.pyplot as plot

    if len(sys.argv) > 1:
        scale = 100
        nb = int(sys.argv[1])
        points = [ (scale*random.random(),scale*random.random()) for i in range(nb)]
    else:
        points = [
                (0,40),
                (100,60),
                (40,0),
                (50,100),
                (90,10),
                # (50,50),
                ]

    fig = plot.figure()

    triangles = delaunay_bowyer_watson( points, do_plot = fig )

    edges = edges_of( triangles )

    ax = fig.add_subplot(111)
    ax.set_aspect('equal')
    scatter_segments( ax, edges, facecolor = "red" )
    plot_segments( ax, edges, edgecolor = "blue" )
    plot.show()

Captions

Add a one-line explanation of what this file represents

Items portrayed in this file

depicts

26 March 2014

image/png

7d7e09b23f0c7ec7c67fc2b2f0ecc33215f95f1b

63,926 byte

900 pixel

1,200 pixel

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current14:37, 26 March 2014Thumbnail for version as of 14:37, 26 March 20141,200 × 900 (62 KB)NojhanUser created page with UploadWizard

The following page uses this file:

Global file usage

The following other wikis use this file:

Metadata