Python: найти область многоугольника из координат XYZ
Я пытаюсь использовать shapely.geometry.Polygon модуль, чтобы найти площадь полигонов, но он выполняет все вычисления на xy самолет. Это хорошо для некоторых моих полигонов, но у других есть z измерение тоже, поэтому он не совсем делает то, что я хотел бы.
Есть ли пакет, который либо даст мне площадь плоского многоугольника из xyz координаты или, альтернативно, пакет или алгоритм для поворота многоугольника в xy самолет, чтобы я мог использовать shapely.geometry.Polygon().area?
Полигоны представлены в виде списка кортежей в виде [(x1,y1,z1),(x2,y2,z3),...(xn,yn,zn)],
8 ответов
Вот вывод формулы для расчета площади трехмерного плоского многоугольника
Вот код Python, который его реализует:
#determinant of matrix a
def det(a):
return a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + a[0][2]*a[1][0]*a[2][1] - a[0][2]*a[1][1]*a[2][0] - a[0][1]*a[1][0]*a[2][2] - a[0][0]*a[1][2]*a[2][1]
#unit normal vector of plane defined by points a, b, and c
def unit_normal(a, b, c):
x = det([[1,a[1],a[2]],
[1,b[1],b[2]],
[1,c[1],c[2]]])
y = det([[a[0],1,a[2]],
[b[0],1,b[2]],
[c[0],1,c[2]]])
z = det([[a[0],a[1],1],
[b[0],b[1],1],
[c[0],c[1],1]])
magnitude = (x**2 + y**2 + z**2)**.5
return (x/magnitude, y/magnitude, z/magnitude)
#dot product of vectors a and b
def dot(a, b):
return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
#cross product of vectors a and b
def cross(a, b):
x = a[1] * b[2] - a[2] * b[1]
y = a[2] * b[0] - a[0] * b[2]
z = a[0] * b[1] - a[1] * b[0]
return (x, y, z)
#area of polygon poly
def area(poly):
if len(poly) < 3: # not a plane - no area
return 0
total = [0, 0, 0]
for i in range(len(poly)):
vi1 = poly[i]
if i is len(poly)-1:
vi2 = poly[0]
else:
vi2 = poly[i+1]
prod = cross(vi1, vi2)
total[0] += prod[0]
total[1] += prod[1]
total[2] += prod[2]
result = dot(total, unit_normal(poly[0], poly[1], poly[2]))
return abs(result/2)
И чтобы проверить это, вот квадрат 10х5, который наклоняется:
>>> poly = [[0, 0, 0], [10, 0, 0], [10, 3, 4], [0, 3, 4]]
>>> poly_translated = [[0+5, 0+5, 0+5], [10+5, 0+5, 0+5], [10+5, 3+5, 4+5], [0+5, 3+5, 4+5]]
>>> area(poly)
50.0
>>> area(poly_translated)
50.0
>>> area([[0,0,0],[1,1,1]])
0
Изначально проблема заключалась в том, что я упростил. Необходимо рассчитать единичный вектор, нормальный к плоскости. Площадь составляет половину точечного произведения этого и суммы всех перекрестных произведений, а не половину суммы всех величин перекрестных произведений.
Это можно немного исправить (классы матриц и векторов сделали бы это лучше, если они у вас есть, или стандартные реализации детерминанта / кросс-продукта / точечного продукта), но это должно быть концептуально обоснованно.
Это последний код, который я использовал. Он не использует фигурную форму, но реализует теорему Стокса для непосредственного расчета площади. Он основан на ответе @Tom Smilack, который показывает, как это сделать без всяких ошибок.
import numpy as np
#unit normal vector of plane defined by points a, b, and c
def unit_normal(a, b, c):
x = np.linalg.det([[1,a[1],a[2]],
[1,b[1],b[2]],
[1,c[1],c[2]]])
y = np.linalg.det([[a[0],1,a[2]],
[b[0],1,b[2]],
[c[0],1,c[2]]])
z = np.linalg.det([[a[0],a[1],1],
[b[0],b[1],1],
[c[0],c[1],1]])
magnitude = (x**2 + y**2 + z**2)**.5
return (x/magnitude, y/magnitude, z/magnitude)
#area of polygon poly
def poly_area(poly):
if len(poly) < 3: # not a plane - no area
return 0
total = [0, 0, 0]
N = len(poly)
for i in range(N):
vi1 = poly[i]
vi2 = poly[(i+1) % N]
prod = np.cross(vi1, vi2)
total[0] += prod[0]
total[1] += prod[1]
total[2] += prod[2]
result = np.dot(total, unit_normal(poly[0], poly[1], poly[2]))
return abs(result/2)
Кстати, вот такой же алгоритм в Mathematica, с детским юнит-тестом
ClearAll[vertexPairs, testPoly, area3D, planeUnitNormal, pairwise];
pairwise[list_, fn_] := MapThread[fn, {Drop[list, -1], Drop[list, 1]}];
vertexPairs[Polygon[{points___}]] := Append[{points}, First[{points}]];
testPoly = Polygon[{{20, -30, 0}, {40, -30, 0}, {40, -30, 20}, {20, -30, 20}}];
planeUnitNormal[Polygon[{points___}]] :=
With[{ps = Take[{points}, 3]},
With[{p0 = First[ps]},
With[{qs = (# - p0) & /@ Rest[ps]},
Normalize[Cross @@ qs]]]];
area3D[p : Polygon[{polys___}]] :=
With[{n = planeUnitNormal[p], vs = vertexPairs[p]},
With[{areas = (Dot[n, #]) & /@ pairwise[vs, Cross]},
Plus @@ areas/2]];
area3D[testPoly]
#pythonn код для полигональной области в 3D (оптимизированная версия)
def polygon_area(poly):
#shape (N, 3)
if isinstance(poly, list):
poly = np.array(poly)
#all edges
edges = poly[1:] - poly[0:1]
# row wise cross product
cross_product = np.cross(edges[:-1],edges[1:], axis=1)
#area of all triangles
area = np.linalg.norm(cross_product, axis=1)/2
return sum(area)
if __name__ == "__main__":
poly = [[0+5, 0+5, 0+5], [10+5, 0+5, 0+5], [10+5, 3+5, 4+5], [0+5, 3+5, 4+5]]
print(polygon_area(poly))
Спасибо за подробные ответы, но я немного удивлен, что нет простого ответа, чтобы получить площадь.
Итак, я просто публикую упрощенный подход для расчета площади с использованием трехмерных координат многоугольника или поверхности с использованием pyny3d.
#Install pyny3d as:
pip install pyny3d
#Calculate area
import numpy as np
import pyny3d.geoms as pyny
coords_3d = np.array([[0, 0, 0],
[7, 0, 0],
[7, 10, 2],
[0, 10, 2]])
polygon = pyny.Polygon(coords_3d)
print(f'Area is : {polygon.get_area()}')
Площадь 2D многоугольника может быть рассчитана с использованием Numpy в качестве однострочного...
poly_Area(vertices) = np.sum( [0.5, -0.5] * vertices * np.roll( np.roll(vertices, 1, axis=0), 1, axis=1) )
То же, что и ответ @Tom Smilack, но в javascript
//determinant of matrix a
function det(a) {
return a[0][0] * a[1][1] * a[2][2] + a[0][1] * a[1][2] * a[2][0] + a[0][2] * a[1][0] * a[2][1] - a[0][2] * a[1][1] * a[2][0] - a[0][1] * a[1][0] * a[2][2] - a[0][0] * a[1][2] * a[2][1];
}
//unit normal vector of plane defined by points a, b, and c
function unit_normal(a, b, c) {
let x = math.det([
[1, a[1], a[2]],
[1, b[1], b[2]],
[1, c[1], c[2]]
]);
let y = math.det([
[a[0], 1, a[2]],
[b[0], 1, b[2]],
[c[0], 1, c[2]]
]);
let z = math.det([
[a[0], a[1], 1],
[b[0], b[1], 1],
[c[0], c[1], 1]
]);
let magnitude = Math.pow(Math.pow(x, 2) + Math.pow(y, 2) + Math.pow(z, 2), 0.5);
return [x / magnitude, y / magnitude, z / magnitude];
}
// dot product of vectors a and b
function dot(a, b) {
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
// cross product of vectors a and b
function cross(a, b) {
let x = (a[1] * b[2]) - (a[2] * b[1]);
let y = (a[2] * b[0]) - (a[0] * b[2]);
let z = (a[0] * b[1]) - (a[1] * b[0]);
return [x, y, z];
}
// area of polygon poly
function area(poly) {
if (poly.length < 3) {
console.log("not a plane - no area");
return 0;
} else {
let total = [0, 0, 0]
for (let i = 0; i < poly.length; i++) {
var vi1 = poly[i];
if (i === poly.length - 1) {
var vi2 = poly[0];
} else {
var vi2 = poly[i + 1];
}
let prod = cross(vi1, vi2);
total[0] = total[0] + prod[0];
total[1] = total[1] + prod[1];
total[2] = total[2] + prod[2];
}
let result = dot(total, unit_normal(poly[0], poly[1], poly[2]));
return Math.abs(result/2);
}
}Этот код совпадает с ответом @Michael на JavaScript, но на C# (требуется System.Numerics) - проверен с помощью решения @Tom Smilack и получил те же результаты.
using System.Numerics;
public static float Area( List<Vector3> poly) {
if (poly.Count < 3)
{
throw new InvalidOperationException("Planar area require at least 3 vertices");
}
Vector3 total = Vector3.Zero;
for (var i = 0; i < poly.Count; i++) {
Vector3 vi1 = poly[i];
Vector3 vi2;
if (i == poly.Count - 1) {
vi2 = poly[0];
} else {
vi2 = poly[i + 1];
}
var prod = Vector3.Cross( vi1, vi2 );
total.X += prod.X;
total.Y += prod.Y;
total.Z += prod.Z;
}
// Unit normal vector can be calculated using this answer as well
// https://stackoverflow.com/a/1966605/15408268
// but since Plane class is available in System.Numerics - we can use it
var plane = Plane.CreateFromVertices( poly[0], poly[1], poly[2] );
var result = Plane.DotNormal( plane, total );
return Math.Abs(result/2f);
}
var l = new List<Vector3>();
l.Add( new Vector3( 5,5,5) );
l.Add( new Vector3( 15,5,5) );
l.Add( new Vector3( 15,8,9) );
l.Add( new Vector3( 5,8,9) );
Area(l); // <-- 50
l = new List<Vector3>();
l.Add( new Vector3( 0,0,0) );
l.Add( new Vector3( 10,0,0) );
l.Add( new Vector3( 10,3,4) );
l.Add( new Vector3( 0,3,4) );
Area(l); // <-- 50
Вам нужно отсортировать вершины по часовой стрелке, чтобы это решение работало правильно. См. /questions/7906130/sortirovka-tochek-po-chasovoj-strelke/7906141#7906141 для идей.
Или вы можете просто использовать этот код: (Эта функция использует функцию Tuple из C # 7.0 https://docs.microsoft.com/en-us/dotnet/csharp/language-reference/builtin-types/value-tuples )
/// <summary>
/// Sort the list of points clockwise
/// </summary>
/// <param name="toSort"></param>
/// <returns></returns>
public IEnumerable<T> SortPointClockwise<T>( IEnumerable<T> pointList, Func<T, (double x, double y)> pointSelector)
{
(double x, double y) center = (0, 0);
var input = pointList.Select( e => ( pointSelector(e), e ) ).ToList();
// Center is the average of all points in each direction
center.x = input.Average(pt => pt.Item1.x);
center.y = input.Average(pt => pt.Item1.y);
//https://stackoverflow.com/a/46635372
// a comes before b in counterclockwise
Func<(double x, double y), (double x, double y), bool> abCounterClockwise = (a, b) =>
{
// Computes the quadrant for a and b (0-3):
// ^
// 1 | 0
// ---+-->
// 2 | 3
var dax = ((a.x - center.x) > 0) ? 1 : 0;
var day = ((a.y - center.y) > 0) ? 1 : 0;
var qa = (1 - dax) + (1 - day) + ((dax & (1 - day)) << 1);
/* The previous computes the following:
const int qa =
( (a.x() > center.x())
? ((a.y() > center.y())
? 0 : 3)
: ((a.y() > center.y())
? 1 : 2)); */
var dbx = ((b.x - center.x) > 0) ? 1 : 0;
var dby = ((b.y - center.y) > 0) ? 1 : 0;
var qb = (1 - dbx) + (1 - dby) + ((dbx & (1 - dby)) << 1);
if (qa == qb)
{
return (b.x - center.x) * (a.y - center.y) < (b.y - center.y) * (a.x - center.x);
}
else
{
return qa < qb;
}
};
// sort clockwise (for counter-clockwise switch 1 and -1)
input.Sort( (a, b) => abCounterClockwise(a.Item1, b.Item1) ? 1 : -1 );
return input.Select( sorted => sorted.e );
}
// these points need to be sorted
var l = new List<Vector3>();
l.Add( new Vector3( 5,8,9) );
l.Add( new Vector3( 5,5,5) );
l.Add( new Vector3( 15,8,9) );
l.Add( new Vector3( 15,5,5) );
var sorted = SortPointClockwise( l, v3 => (v3.X, v3.Y) ).ToList();
Area(sorted); // <-- 50
Спасибо всем в StackOverflow за эти решения.
Полный код: https://dotnetfiddle.net/7jG4Wy