Fichier:Bowyer-Watson_0.png
De Wikipedia, l'encyclopédie encyclopedia
Fichier d’origine (1 200 × 900 pixels, taille du fichier : 62 kio, type MIME : image/png)
Ce fichier et sa description proviennent de Wikimedia Commons.
Description
DescriptionBowyer-Watson 0.png |
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 :
|
Date | |
Source | Travail personnel (see source code below) |
Auteur | Nojhan |
Conditions d’utilisation
- Vous êtes libre :
- de partager – de copier, distribuer et transmettre cette œuvre
- d’adapter – de modifier cette œuvre
- Sous les conditions suivantes :
- paternité – Vous devez donner les informations appropriées concernant l'auteur, fournir un lien vers la licence et indiquer si des modifications ont été faites. Vous pouvez faire cela par tout moyen raisonnable, mais en aucune façon suggérant que l’auteur vous soutient ou approuve l’utilisation que vous en faites.
- partage à l’identique – Si vous modifiez, transformez, ou vous basez sur cette œuvre, vous devez distribuer votre contribution sous la même licence ou une licence compatible avec celle de l’original.
Vous avez la permission de copier, distribuer et modifier ce document selon les termes de la GNU Free Documentation License version 1.2 ou toute version ultérieure publiée par la Free Software Foundation, sans sections inaltérables, sans texte de première page de couverture et sans texte de dernière page de couverture. Un exemplaire de la licence est inclus dans la section intitulée GNU Free Documentation License.http://www.gnu.org/copyleft/fdl.htmlGFDLGNU Free Documentation Licensetruetrue |
Cette œuvre est un logiciel libre ; elle peut être redistribuée ou modifiée selon les termes de la licence CeCILL.
Le texte de cette licence est consultable sur www.cecill.info.CEA CNRS INRIA Logiciel LibreCeCILLhttp://www.cecill.info/licences/Licence_CeCILL_V2-en.htmltrue |
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()
Éléments décrits dans ce fichier
dépeint
Valeur sans élément de Wikidata
26 mars 2014
image/png
7d7e09b23f0c7ec7c67fc2b2f0ecc33215f95f1b
63 926 octet
900 pixel
1 200 pixel
Historique du fichier
Cliquer sur une date et heure pour voir le fichier tel qu'il était à ce moment-là.
Date et heure | Vignette | Dimensions | Utilisateur | Commentaire | |
---|---|---|---|---|---|
actuel | 26 mars 2014 à 16:37 | 1 200 × 900 (62 kio) | Nojhan | User created page with UploadWizard |
Utilisation du fichier
La page suivante utilise ce fichier :
Usage global du fichier
Les autres wikis suivants utilisent ce fichier :
- Utilisation sur en.wikipedia.org
- Utilisation sur es.wikipedia.org
- Utilisation sur ru.wikipedia.org
Métadonnées
Ce fichier contient des informations supplémentaires, probablement ajoutées par l'appareil photo numérique ou le numériseur utilisé pour le créer.
Si le fichier a été modifié depuis son état original, certains détails peuvent ne pas refléter entièrement l'image modifiée.
Résolution horizontale | 59,05 pt/cm |
---|---|
Résolution verticale | 59,05 pt/cm |