Class Triangle

1	/*
2	* Copyright (c) 2016 Vivid Solutions.
3	*
4	* All rights reserved. This program and the accompanying materials
5	* are made available under the terms of the Eclipse Public License 2.0
6	* and Eclipse Distribution License v. 1.0 which accompanies this distribution.
7	* The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
8	* and the Eclipse Distribution License is available at
9	*
10	* http://www.eclipse.org/org/documents/edl-v10.php.
11	*/
12	package org.locationtech.jts.geom;
13
14	import org.locationtech.jts.algorithm.Angle;
15	import org.locationtech.jts.algorithm.HCoordinate;
16	import org.locationtech.jts.algorithm.Orientation;
17	import org.locationtech.jts.math.DD;
18
19	/**
20	* Represents a planar triangle, and provides methods for calculating various
21	* properties of triangles.
22	*
23	* @version 1.7
24	*/
25	public class Triangle
26	{
27
28	/**
29	* Tests whether a triangle is acute. A triangle is acute iff all interior
30	* angles are acute. This is a strict test - right triangles will return
31	* <tt>false</tt> A triangle which is not acute is either right or obtuse.
32	* <p>
33	* Note: this implementation is not robust for angles very close to 90
34	* degrees.
35	*
36	* @param a
37	* a vertex of the triangle
38	* @param b
39	* a vertex of the triangle
40	* @param c
41	* a vertex of the triangle
42	* @return true if the triangle is acute
43	*/
44	public static boolean isAcute(Coordinate a, Coordinate b, Coordinate c)
45	{
46	if (!Angle.isAcute(a, b, c))
47	return false;
48	if (!Angle.isAcute(b, c, a))
49	return false;
50	if (!Angle.isAcute(c, a, b))
51	return false;
52	return true;
53	}
54
55	/**
56	* Computes the line which is the perpendicular bisector of the line segment
57	* a-b.
58	*
59	* @param a
60	* a point
61	* @param b
62	* another point
63	* @return the perpendicular bisector, as an HCoordinate
64	*/
65	public static HCoordinate perpendicularBisector(Coordinate a, Coordinate b)
66	{
67	// returns the perpendicular bisector of the line segment ab
68	double dx = b.x - a.x;
69	double dy = b.y - a.y;
70	HCoordinate l1 = new HCoordinate(a.x + dx / 2.0, a.y + dy / 2.0, 1.0);
71	HCoordinate l2 = new HCoordinate(a.x - dy + dx / 2.0, a.y + dx + dy / 2.0,
72	1.0);
73	return new HCoordinate(l1, l2);
74	}
75
76	/**
77	* Computes the circumcentre of a triangle. The circumcentre is the centre of
78	* the circumcircle, the smallest circle which encloses the triangle. It is
79	* also the common intersection point of the perpendicular bisectors of the
80	* sides of the triangle, and is the only point which has equal distance to
81	* all three vertices of the triangle.
82	*
83	* @param a
84	* a vertex of the triangle
85	* @param b
86	* a vertex of the triangle
87	* @param c
88	* a vertex of the triangle
89	* @return the circumcentre of the triangle
90	*/
91	/*
92	* // original non-robust algorithm public static Coordinate
93	* circumcentre(Coordinate a, Coordinate b, Coordinate c) { // compute the
94	* perpendicular bisector of chord ab HCoordinate cab =
95	* perpendicularBisector(a, b); // compute the perpendicular bisector of chord
96	* bc HCoordinate cbc = perpendicularBisector(b, c); // compute the
97	* intersection of the bisectors (circle radii) HCoordinate hcc = new
98	* HCoordinate(cab, cbc); Coordinate cc = null; try { cc = new
99	* Coordinate(hcc.getX(), hcc.getY()); } catch (NotRepresentableException ex)
100	* { // MD - not sure what we can do to prevent this (robustness problem) //
101	* Idea - can we condition which edges we choose? throw new
102	* IllegalStateException(ex.getMessage()); }
103	*
104	* //System.out.println("Acc = " + a.distance(cc) + ", Bcc = " +
105	* b.distance(cc) + ", Ccc = " + c.distance(cc) );
106	*
107	* return cc; }
108	*/
109
110	/**
111	* Computes the circumcentre of a triangle. The circumcentre is the centre of
112	* the circumcircle, the smallest circle which encloses the triangle. It is
113	* also the common intersection point of the perpendicular bisectors of the
114	* sides of the triangle, and is the only point which has equal distance to
115	* all three vertices of the triangle.
116	* <p>
117	* The circumcentre does not necessarily lie within the triangle. For example,
118	* the circumcentre of an obtuse isosceles triangle lies outside the triangle.
119	* <p>
120	* This method uses an algorithm due to J.R.Shewchuk which uses normalization
121	* to the origin to improve the accuracy of computation. (See <i>Lecture Notes
122	* on Geometric Robustness</i>, Jonathan Richard Shewchuk, 1999).
123	*
124	* @param a
125	* a vertex of the triangle
126	* @param b
127	* a vertex of the triangle
128	* @param c
129	* a vertex of the triangle
130	* @return the circumcentre of the triangle
131	*/
132	public static Coordinate circumcentre(Coordinate a, Coordinate b, Coordinate c)
133	{
134	double cx = c.x;
135	double cy = c.y;
136	double ax = a.x - cx;
137	double ay = a.y - cy;
138	double bx = b.x - cx;
139	double by = b.y - cy;
140
141	double denom = 2 * det(ax, ay, bx, by);
142	double numx = det(ay, ax * ax + ay * ay, by, bx * bx + by * by);
143	double numy = det(ax, ax * ax + ay * ay, bx, bx * bx + by * by);
144
145	double ccx = cx - numx / denom;
146	double ccy = cy + numy / denom;
147
148	return new Coordinate(ccx, ccy);
149	}
150
151	/**
152	* Computes the circumcentre of a triangle. The circumcentre is the centre of
153	* the circumcircle, the smallest circle which encloses the triangle. It is
154	* also the common intersection point of the perpendicular bisectors of the
155	* sides of the triangle, and is the only point which has equal distance to
156	* all three vertices of the triangle.
157	* <p>
158	* The circumcentre does not necessarily lie within the triangle. For example,
159	* the circumcentre of an obtuse isosceles triangle lies outside the triangle.
160	* <p>
161	* This method uses {@link DD} extended-precision arithmetic to
162	* provide more accurate results than {@link #circumcentre(Coordinate, Coordinate, Coordinate)}
163	*
164	* @param a
165	* a vertex of the triangle
166	* @param b
167	* a vertex of the triangle
168	* @param c
169	* a vertex of the triangle
170	* @return the circumcentre of the triangle
171	*/
172	public static Coordinate circumcentreDD(Coordinate a, Coordinate b, Coordinate c)
173	{
174	DD ax = DD.valueOf(a.x).subtract(c.x);
175	DD ay = DD.valueOf(a.y).subtract(c.y);
176	DD bx = DD.valueOf(b.x).subtract(c.x);
177	DD by = DD.valueOf(b.y).subtract(c.y);
178
179	DD denom = DD.determinant(ax, ay, bx, by).multiply(2);
180	DD asqr = ax.sqr().add( ay.sqr());
181	DD bsqr = bx.sqr().add( by.sqr());
182	DD numx = DD.determinant(ay, asqr, by, bsqr);
183	DD numy = DD.determinant(ax, asqr, bx, bsqr);
184
185	double ccx = DD.valueOf(c.x).subtract( numx.divide(denom) ).doubleValue();
186	double ccy = DD.valueOf(c.y).add( numy.divide(denom) ).doubleValue();
187
188	return new Coordinate(ccx, ccy);
189	}
190
191	/**
192	* Computes the determinant of a 2x2 matrix. Uses standard double-precision
193	* arithmetic, so is susceptible to round-off error.
194	*
195	* @param m00
196	* the [0,0] entry of the matrix
197	* @param m01
198	* the [0,1] entry of the matrix
199	* @param m10
200	* the [1,0] entry of the matrix
201	* @param m11
202	* the [1,1] entry of the matrix
203	* @return the determinant
204	*/
205	private static double det(double m00, double m01, double m10, double m11)
206	{
207	return m00 * m11 - m01 * m10;
208	}
209
210	/**
211	* Computes the incentre of a triangle. The <i>inCentre</i> of a triangle is
212	* the point which is equidistant from the sides of the triangle. It is also
213	* the point at which the bisectors of the triangle's angles meet. It is the
214	* centre of the triangle's <i>incircle</i>, which is the unique circle that
215	* is tangent to each of the triangle's three sides.
216	* <p>
217	* The incentre always lies within the triangle.
218	*
219	* @param a
220	* a vertex of the triangle
221	* @param b
222	* a vertex of the triangle
223	* @param c
224	* a vertex of the triangle
225	* @return the point which is the incentre of the triangle
226	*/
227	public static Coordinate inCentre(Coordinate a, Coordinate b, Coordinate c)
228	{
229	// the lengths of the sides, labelled by their opposite vertex
230	double len0 = b.distance(c);
231	double len1 = a.distance(c);
232	double len2 = a.distance(b);
233	double circum = len0 + len1 + len2;
234
235	double inCentreX = (len0 * a.x + len1 * b.x + len2 * c.x) / circum;
236	double inCentreY = (len0 * a.y + len1 * b.y + len2 * c.y) / circum;
237	return new Coordinate(inCentreX, inCentreY);
238	}
239
240	/**
241	* Computes the centroid (centre of mass) of a triangle. This is also the
242	* point at which the triangle's three medians intersect (a triangle median is
243	* the segment from a vertex of the triangle to the midpoint of the opposite
244	* side). The centroid divides each median in a ratio of 2:1.
245	* <p>
246	* The centroid always lies within the triangle.
247	*
248	*
249	* @param a
250	* a vertex of the triangle
251	* @param b
252	* a vertex of the triangle
253	* @param c
254	* a vertex of the triangle
255	* @return the centroid of the triangle
256	*/
257	public static Coordinate centroid(Coordinate a, Coordinate b, Coordinate c)
258	{
259	double x = (a.x + b.x + c.x) / 3;
260	double y = (a.y + b.y + c.y) / 3;
261	return new Coordinate(x, y);
262	}
263
264	/**
265	* Computes the length of the longest side of a triangle
266	*
267	* @param a
268	* a vertex of the triangle
269	* @param b
270	* a vertex of the triangle
271	* @param c
272	* a vertex of the triangle
273	* @return the length of the longest side of the triangle
274	*/
275	public static double longestSideLength(Coordinate a, Coordinate b,
276	Coordinate c)
277	{
278	double lenAB = a.distance(b);
279	double lenBC = b.distance(c);
280	double lenCA = c.distance(a);
281	double maxLen = lenAB;
282	if (lenBC > maxLen)
283	maxLen = lenBC;
284	if (lenCA > maxLen)
285	maxLen = lenCA;
286	return maxLen;
287	}
288
289	/**
290	* Computes the point at which the bisector of the angle ABC cuts the segment
291	* AC.
292	*
293	* @param a
294	* a vertex of the triangle
295	* @param b
296	* a vertex of the triangle
297	* @param c
298	* a vertex of the triangle
299	* @return the angle bisector cut point
300	*/
301	public static Coordinate angleBisector(Coordinate a, Coordinate b,
302	Coordinate c)
303	{
304	/**
305	* Uses the fact that the lengths of the parts of the split segment are
306	* proportional to the lengths of the adjacent triangle sides
307	*/
308	double len0 = b.distance(a);
309	double len2 = b.distance(c);
310	double frac = len0 / (len0 + len2);
311	double dx = c.x - a.x;
312	double dy = c.y - a.y;
313
314	Coordinate splitPt = new Coordinate(a.x + frac * dx, a.y + frac * dy);
315	return splitPt;
316	}
317
318	/**
319	* Computes the 2D area of a triangle. The area value is always non-negative.
320	*
321	* @param a
322	* a vertex of the triangle
323	* @param b
324	* a vertex of the triangle
325	* @param c
326	* a vertex of the triangle
327	* @return the area of the triangle
328	*
329	* @see #signedArea(Coordinate, Coordinate, Coordinate)
330	*/
331	public static double area(Coordinate a, Coordinate b, Coordinate c)
332	{
333	return Math
334	.abs(((c.x - a.x) * (b.y - a.y) - (b.x - a.x) * (c.y - a.y)) / 2);
335	}
336
337	/**
338	* Computes the signed 2D area of a triangle. The area value is positive if
339	* the triangle is oriented CW, and negative if it is oriented CCW.
340	* <p>
341	* The signed area value can be used to determine point orientation, but the
342	* implementation in this method is susceptible to round-off errors. Use
343	* {@link Orientation#index(Coordinate, Coordinate, Coordinate)}
344	* for robust orientation calculation.
345	*
346	* @param a
347	* a vertex of the triangle
348	* @param b
349	* a vertex of the triangle
350	* @param c
351	* a vertex of the triangle
352	* @return the signed 2D area of the triangle
353	*
354	* @see Orientation#index(Coordinate, Coordinate, Coordinate)
355	*/
356	public static double signedArea(Coordinate a, Coordinate b, Coordinate c)
357	{
358	/**
359	* Uses the formula 1/2 * \| u x v \| where u,v are the side vectors of the
360	* triangle x is the vector cross-product For 2D vectors, this formula
361	* simplifies to the expression below
362	*/
363	return ((c.x - a.x) * (b.y - a.y) - (b.x - a.x) * (c.y - a.y)) / 2;
364	}
365
366	/**
367	* Computes the 3D area of a triangle. The value computed is always
368	* non-negative.
369	*
370	* @param a
371	* a vertex of the triangle
372	* @param b
373	* a vertex of the triangle
374	* @param c
375	* a vertex of the triangle
376	* @return the 3D area of the triangle
377	*/
378	public static double area3D(Coordinate a, Coordinate b, Coordinate c)
379	{
380	/**
381	* Uses the formula 1/2 * \| u x v \| where u,v are the side vectors of the
382	* triangle x is the vector cross-product
383	*/
384	// side vectors u and v
385	double ux = b.x - a.x;
386	double uy = b.y - a.y;
387	double uz = b.getZ() - a.getZ();
388
389	double vx = c.x - a.x;
390	double vy = c.y - a.y;
391	double vz = c.getZ() - a.getZ();
392
393	// cross-product = u x v
394	double crossx = uy * vz - uz * vy;
395	double crossy = uz * vx - ux * vz;
396	double crossz = ux * vy - uy * vx;
397
398	// tri area = 1/2 * \| u x v \|
399	double absSq = crossx * crossx + crossy * crossy + crossz * crossz;
400	double area3D = Math.sqrt(absSq) / 2;
401
402	return area3D;
403	}
404
405	/**
406	* Computes the Z-value (elevation) of an XY point on a three-dimensional
407	* plane defined by a triangle whose vertices have Z-values. The defining
408	* triangle must not be degenerate (in other words, the triangle must enclose
409	* a non-zero area), and must not be parallel to the Z-axis.
410	* <p>
411	* This method can be used to interpolate the Z-value of a point inside a
412	* triangle (for example, of a TIN facet with elevations on the vertices).
413	*
414	* @param p
415	* the point to compute the Z-value of
416	* @param v0
417	* a vertex of a triangle, with a Z ordinate
418	* @param v1
419	* a vertex of a triangle, with a Z ordinate
420	* @param v2
421	* a vertex of a triangle, with a Z ordinate
422	* @return the computed Z-value (elevation) of the point
423	*/
424	public static double interpolateZ(Coordinate p, Coordinate v0, Coordinate v1,
425	Coordinate v2)
426	{
427	double x0 = v0.x;
428	double y0 = v0.y;
429	double a = v1.x - x0;
430	double b = v2.x - x0;
431	double c = v1.y - y0;
432	double d = v2.y - y0;
433	double det = a * d - b * c;
434	double dx = p.x - x0;
435	double dy = p.y - y0;
436	double t = (d * dx - b * dy) / det;
437	double u = (-c * dx + a * dy) / det;
438	double z = v0.getZ() + t * (v1.getZ() - v0.getZ()) + u * (v2.getZ() - v0.getZ());
439	return z;
440	}
441
442	/**
443	* The coordinates of the vertices of the triangle
444	*/
445	public Coordinate p0, p1, p2;
446
447	/**
448	* Creates a new triangle with the given vertices.
449	*
450	* @param p0
451	* a vertex
452	* @param p1
453	* a vertex
454	* @param p2
455	* a vertex
456	*/
457	public Triangle(Coordinate p0, Coordinate p1, Coordinate p2)
458	{
459	this.p0 = p0;
460	this.p1 = p1;
461	this.p2 = p2;
462	}
463
464	/**
465	* Computes the incentre of this triangle. The <i>incentre</i> of a triangle
466	* is the point which is equidistant from the sides of the triangle. It is
467	* also the point at which the bisectors of the triangle's angles meet. It is
468	* the centre of the triangle's <i>incircle</i>, which is the unique circle
469	* that is tangent to each of the triangle's three sides.
470	*
471	* @return the point which is the inCentre of this triangle
472	*/
473	public Coordinate inCentre()
474	{
475	return inCentre(p0, p1, p2);
476	}
477
478	/**
479	* Tests whether this triangle is acute. A triangle is acute iff all interior
480	* angles are acute. This is a strict test - right triangles will return
481	* <tt>false</tt> A triangle which is not acute is either right or obtuse.
482	* <p>
483	* Note: this implementation is not robust for angles very close to 90
484	* degrees.
485	*
486	* @return true if this triangle is acute
487	*/
488	public boolean isAcute()
489	{
490	return isAcute(this.p0, this.p1, this.p2);
491	}
492
493	/**
494	* Computes the circumcentre of this triangle. The circumcentre is the centre
495	* of the circumcircle, the smallest circle which encloses the triangle. It is
496	* also the common intersection point of the perpendicular bisectors of the
497	* sides of the triangle, and is the only point which has equal distance to
498	* all three vertices of the triangle.
499	* <p>
500	* The circumcentre does not necessarily lie within the triangle.
501	* <p>
502	* This method uses an algorithm due to J.R.Shewchuk which uses normalization
503	* to the origin to improve the accuracy of computation. (See <i>Lecture Notes
504	* on Geometric Robustness</i>, Jonathan Richard Shewchuk, 1999).
505	*
506	* @return the circumcentre of this triangle
507	*/
508	public Coordinate circumcentre()
509	{
510	return circumcentre(this.p0, this.p1, this.p2);
511	}
512
513	/**
514	* Computes the centroid (centre of mass) of this triangle. This is also the
515	* point at which the triangle's three medians intersect (a triangle median is
516	* the segment from a vertex of the triangle to the midpoint of the opposite
517	* side). The centroid divides each median in a ratio of 2:1.
518	* <p>
519	* The centroid always lies within the triangle.
520	*
521	* @return the centroid of this triangle
522	*/
523	public Coordinate centroid()
524	{
525	return centroid(this.p0, this.p1, this.p2);
526	}
527
528	/**
529	* Computes the length of the longest side of this triangle
530	*
531	* @return the length of the longest side of this triangle
532	*/
533	public double longestSideLength()
534	{
535	return longestSideLength(this.p0, this.p1, this.p2);
536	}
537
538	/**
539	* Computes the 2D area of this triangle. The area value is always
540	* non-negative.
541	*
542	* @return the area of this triangle
543	*
544	* @see #signedArea()
545	*/
546	public double area()
547	{
548	return area(this.p0, this.p1, this.p2);
549	}
550
551	/**
552	* Computes the signed 2D area of this triangle. The area value is positive if
553	* the triangle is oriented CW, and negative if it is oriented CCW.
554	* <p>
555	* The signed area value can be used to determine point orientation, but the
556	* implementation in this method is susceptible to round-off errors. Use
557	* {@link Orientation#index(Coordinate, Coordinate, Coordinate)}
558	* for robust orientation calculation.
559	*
560	* @return the signed 2D area of this triangle
561	*
562	* @see Orientation#index(Coordinate, Coordinate, Coordinate)
563	*/
564	public double signedArea()
565	{
566	return signedArea(this.p0, this.p1, this.p2);
567	}
568
569	/**
570	* Computes the 3D area of this triangle. The value computed is always
571	* non-negative.
572	*
573	* @return the 3D area of this triangle
574	*/
575	public double area3D()
576	{
577	return area3D(this.p0, this.p1, this.p2);
578	}
579
580	/**
581	* Computes the Z-value (elevation) of an XY point on a three-dimensional
582	* plane defined by this triangle (whose vertices must have Z-values). This
583	* triangle must not be degenerate (in other words, the triangle must enclose
584	* a non-zero area), and must not be parallel to the Z-axis.
585	* <p>
586	* This method can be used to interpolate the Z-value of a point inside this
587	* triangle (for example, of a TIN facet with elevations on the vertices).
588	*
589	* @param p
590	* the point to compute the Z-value of
591	* @return the computed Z-value (elevation) of the point
592	*/
593	public double interpolateZ(Coordinate p)
594	{
595	if (p == null)
596	throw new IllegalArgumentException("Supplied point is null.");
597	return interpolateZ(p, this.p0, this.p1, this.p2);
598	}
599
600	}
601