001/** 002 * Copyright (c) 2011, The University of Southampton and the individual contributors. 003 * All rights reserved. 004 * 005 * Redistribution and use in source and binary forms, with or without modification, 006 * are permitted provided that the following conditions are met: 007 * 008 * * Redistributions of source code must retain the above copyright notice, 009 * this list of conditions and the following disclaimer. 010 * 011 * * Redistributions in binary form must reproduce the above copyright notice, 012 * this list of conditions and the following disclaimer in the documentation 013 * and/or other materials provided with the distribution. 014 * 015 * * Neither the name of the University of Southampton nor the names of its 016 * contributors may be used to endorse or promote products derived from this 017 * software without specific prior written permission. 018 * 019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 029 */ 030package org.openimaj.math.geometry.triangulation; 031 032import java.util.ArrayList; 033import java.util.Collections; 034import java.util.Comparator; 035import java.util.HashSet; 036import java.util.List; 037 038import org.openimaj.math.geometry.line.Line2d; 039import org.openimaj.math.geometry.point.Point2d; 040import org.openimaj.math.geometry.point.Point2dImpl; 041import org.openimaj.math.geometry.shape.Circle; 042import org.openimaj.math.geometry.shape.Triangle; 043 044/** 045 * The Delaunay Triangulation algorithm. Produces a triangulation of a set of 046 * points. 047 * <p> 048 * Originally ported from <a 049 * href="http://paulbourke.net/papers/triangulate/"Paul Bourke's 050 * triangulate.c</a>. 051 * <p> 052 * This OpenIMAJ version is based off <a 053 * href="http://www.florianjenett.de/">Florian Jenett's</a> Java port: fjenett, 054 * 20th february 2005, offenbach-germany. contact: http://www.florianjenett.de/ 055 */ 056public class DelaunayTriangulator { 057 /* 058 * From P Bourke's C prototype - 059 * 060 * qsort(p,nv,sizeof(XYZ),XYZCompare); 061 * 062 * int XYZCompare(void *v1,void *v2) { XYZ *p1,*p2; p1 = v1; p2 = v2; if 063 * (p1->x < p2->x) return(-1); else if (p1->x > p2->x) return(1); else 064 * return(0); } 065 */ 066 private static class XComparator implements Comparator<Point2d> { 067 068 @Override 069 public int compare(Point2d p1, Point2d p2) { 070 if (p1.getX() < p2.getX()) { 071 return -1; 072 } 073 else if (p1.getX() > p2.getX()) { 074 return 1; 075 } 076 else { 077 return 0; 078 } 079 } 080 } 081 082 static final float EPSILON = 0.00000001f; 083 084 /* 085 * Return TRUE if a point (xp,yp) is inside the circumcircle made up of the 086 * points (x1,y1), (x2,y2), (x3,y3) The circumcircle centre is returned in 087 * (xc,yc) and the radius r NOTE: A point on the edge is inside the 088 * circumcircle 089 */ 090 private static boolean circumCircle(Point2d p, Triangle t, Circle circle) { 091 float m1, m2, mx1, mx2, my1, my2; 092 float dx, dy, rsqr, drsqr; 093 094 /* Check for coincident points */ 095 if (Math.abs(t.firstVertex().getY() - t.secondVertex().getY()) < EPSILON 096 && Math.abs(t.secondVertex().getY() - t.thirdVertex().getY()) < EPSILON) 097 { 098 System.err.println("CircumCircle: Points are coincident."); 099 return false; 100 } 101 102 if (Math.abs(t.secondVertex().getY() - t.firstVertex().getY()) < EPSILON) { 103 m2 = -(t.thirdVertex().getX() - t.secondVertex().getX()) / (t.thirdVertex().getY() - t.secondVertex().getY()); 104 mx2 = (t.secondVertex().getX() + t.thirdVertex().getX()) / 2.0f; 105 my2 = (t.secondVertex().getY() + t.thirdVertex().getY()) / 2.0f; 106 circle.setX((t.secondVertex().getX() + t.firstVertex().getX()) / 2.0f); 107 circle.setY(m2 * (circle.getX() - mx2) + my2); 108 } 109 else if (Math.abs(t.thirdVertex().getY() - t.secondVertex().getY()) < EPSILON) { 110 m1 = -(t.secondVertex().getX() - t.firstVertex().getX()) / (t.secondVertex().getY() - t.firstVertex().getY()); 111 mx1 = (t.firstVertex().getX() + t.secondVertex().getX()) / 2.0f; 112 my1 = (t.firstVertex().getY() + t.secondVertex().getY()) / 2.0f; 113 circle.setX((t.thirdVertex().getX() + t.secondVertex().getX()) / 2.0f); 114 circle.setY(m1 * (circle.getX() - mx1) + my1); 115 } 116 else { 117 m1 = -(t.secondVertex().getX() - t.firstVertex().getX()) / (t.secondVertex().getY() - t.firstVertex().getY()); 118 m2 = -(t.thirdVertex().getX() - t.secondVertex().getX()) / (t.thirdVertex().getY() - t.secondVertex().getY()); 119 mx1 = (t.firstVertex().getX() + t.secondVertex().getX()) / 2.0f; 120 mx2 = (t.secondVertex().getX() + t.thirdVertex().getX()) / 2.0f; 121 my1 = (t.firstVertex().getY() + t.secondVertex().getY()) / 2.0f; 122 my2 = (t.secondVertex().getY() + t.thirdVertex().getY()) / 2.0f; 123 circle.setX((m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2)); 124 circle.setY(m1 * (circle.getX() - mx1) + my1); 125 } 126 127 dx = t.secondVertex().getX() - circle.getX(); 128 dy = t.secondVertex().getY() - circle.getY(); 129 rsqr = dx * dx + dy * dy; 130 circle.setRadius((float) Math.sqrt(rsqr)); 131 132 dx = p.getX() - circle.getX(); 133 dy = p.getY() - circle.getY(); 134 drsqr = dx * dx + dy * dy; 135 136 return drsqr <= rsqr; 137 } 138 139 /** 140 * Trianglate a set of vertices. 141 * 142 * @param pxyz 143 * vertices to triangulate. 144 * @return list of triangles arranged in a consistent clockwise order. 145 */ 146 public static List<Triangle> triangulate(List<? extends Point2d> pxyz) { 147 148 // sort vertex array in increasing x values 149 Collections.sort(pxyz, new XComparator()); 150 151 /* 152 * Find the maximum and minimum vertex bounds. This is to allow 153 * calculation of the bounding triangle 154 */ 155 float xmin = ((Point2d) pxyz.get(0)).getX(); 156 float ymin = ((Point2d) pxyz.get(0)).getY(); 157 float xmax = xmin; 158 float ymax = ymin; 159 160 for (final Point2d p : pxyz) { 161 if (p.getX() < xmin) 162 xmin = p.getX(); 163 if (p.getX() > xmax) 164 xmax = p.getX(); 165 if (p.getY() < ymin) 166 ymin = p.getY(); 167 if (p.getY() > ymax) 168 ymax = p.getY(); 169 } 170 171 final float dx = xmax - xmin; 172 final float dy = ymax - ymin; 173 final float dmax = (dx > dy) ? dx : dy; 174 final float xmid = (xmax + xmin) / 2.0f; 175 final float ymid = (ymax + ymin) / 2.0f; 176 177 final ArrayList<Triangle> triangles = new ArrayList<Triangle>(); // for 178 // the 179 // Triangles 180 final HashSet<Triangle> complete = new HashSet<Triangle>(); // for 181 // complete 182 // Triangles 183 184 /* 185 * Set up the supertriangle This is a triangle which encompasses all the 186 * sample points. The supertriangle coordinates are added to the end of 187 * the vertex list. The supertriangle is the first triangle in the 188 * triangle list. 189 */ 190 final Triangle superTriangle = new Triangle( 191 new Point2dImpl(xmid - 2.0f * dmax, ymid - dmax), 192 new Point2dImpl(xmid, ymid + 2.0f * dmax), 193 new Point2dImpl(xmid + 2.0f * dmax, ymid - dmax) 194 ); 195 triangles.add(superTriangle); 196 197 /* 198 * Include each point one at a time into the existing mesh 199 */ 200 final ArrayList<Line2d> edges = new ArrayList<Line2d>(); 201 for (final Point2d p : pxyz) { 202 edges.clear(); 203 204 /* 205 * Set up the edge buffer. If the point (xp,yp) lies inside the 206 * circumcircle then the three edges of that triangle are added to 207 * the edge buffer and that triangle is removed. 208 */ 209 final Circle circle = new Circle(0, 0, 0); 210 211 for (int j = triangles.size() - 1; j >= 0; j--) { 212 213 final Triangle t = triangles.get(j); 214 if (complete.contains(t)) { 215 continue; 216 } 217 218 final boolean inside = circumCircle(p, t, circle); 219 220 if (circle.getX() + circle.getRadius() < p.getX()) { 221 complete.add(t); 222 } 223 if (inside) { 224 edges.add(new Line2d(t.firstVertex(), t.secondVertex())); 225 edges.add(new Line2d(t.secondVertex(), t.thirdVertex())); 226 edges.add(new Line2d(t.thirdVertex(), t.firstVertex())); 227 triangles.remove(j); 228 } 229 230 } 231 232 /* 233 * Tag multiple edges Note: if all triangles are specified 234 * anticlockwise then all interior edges are opposite pointing in 235 * direction. 236 */ 237 for (int j = 0; j < edges.size() - 1; j++) { 238 final Line2d e1 = edges.get(j); 239 for (int k = j + 1; k < edges.size(); k++) { 240 final Line2d e2 = edges.get(k); 241 if (e1.begin == e2.end && e1.end == e2.begin) { 242 e1.begin = null; 243 e1.end = null; 244 e2.begin = null; 245 e2.end = null; 246 } 247 /* Shouldn't need the following, see note above */ 248 if (e1.begin == e2.begin && e1.end == e2.end) { 249 e1.begin = null; 250 e1.end = null; 251 e2.begin = null; 252 e2.end = null; 253 } 254 } 255 } 256 257 /* 258 * Form new triangles for the current point Skipping over any tagged 259 * edges. All edges are arranged in clockwise order. 260 */ 261 for (int j = 0; j < edges.size(); j++) { 262 final Line2d e = edges.get(j); 263 if (e.begin == null || e.end == null) { 264 continue; 265 } 266 triangles.add(new Triangle(e.begin, e.end, p)); 267 } 268 269 } 270 271 /* 272 * Remove triangles with supertriangle vertices 273 */ 274 for (int i = triangles.size() - 1; i >= 0; i--) { 275 final Triangle t = triangles.get(i); 276 if (t.sharesVertex(superTriangle)) { 277 triangles.remove(i); 278 } 279 } 280 281 return triangles; 282 } 283}