For posterity and for those interested, I finally put all the pieces together on the thing I was working on. This script does the following:
1. For a given polygon shape, create a convex hull.
2. From the convex hull, determine the minimum-area oriented bounding box (obb)
3. Set the chunk region to be this oriented bounding box (conform_region_to_polygon).
This was all needed because we want to use known GEOJSON shapes of property lines to shrink the processing region automatically.
NOTE: the script assumes that:
1. You've already run "Align Photos..." on a set of georeferenced photos
2. The chunk coordinate system is WGS 84.
3. The chunk has a single polygon shape with at least 3 vertices.
__author__ = 'Colin Basnett (https://github.com/cmbasnett)'
import math
from functools import cmp_to_key
from PhotoScan import *
# Three points are a counter-clockwise turn if ccw > 0, clockwise if
# ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
# gives twice the signed area of the triangle formed by p1, p2 and p3.
def ccw(p1, p2, p3):
return (p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x)
def angle(p1, p2):
x = (p2.x - p1.x)
y = (p2.y - p1.y)
a = math.atan2(y, x)
return a
def points_cmp(a, b):
if a.y == b.y:
return a.x - b.x
else:
return a.y - b.y
# Removes duplicates from a sorted list (duplicate items are guaranteed to be adjacent)
def remove_duplicates(points):
old = points[-1]
for i in range(len(points) - 2, 0, -1):
if points[i] == old:
points.pop(i)
else:
old = points[i]
# Graham Scan to get a convex hull from a set of points
def get_convex_hull(points):
# Sort by y, then x
points = sorted(points, key=cmp_to_key(points_cmp))
# Remove duplicates from sorted list
remove_duplicates(points)
# Now sort by polar angle from the first point
points[1:] = sorted(points[1:], key=lambda p: angle(points[0], p))
# Remove adjacent collinear line segments.
for i in range(len(points) - 2, 0, -1):
a = angle(points[i - 1], points[i])
b = angle(points[i], points[i + 1])
if a == b:
points.pop(i)
stack = points[0:3]
for i in range(3, len(points)):
while ccw(stack[-2], stack[-1], points[i]) <= 0:
stack.pop()
stack.append(points[i])
return stack
def distance(p1, p2):
return math.hypot(p1.x - p2.x, p1.y - p2.y)
# Numerically stable version of Heron's formula
def area_of_triangle(points):
assert len(points) == 3
a, b, c = points
a, b, c = sorted([distance(a, b), distance(b, c), distance(c, a)], reverse=True)
return 0.25 * math.sqrt((a + (b + c)) * (c - (a - b)) * (c + (a -b)) * (a + (b - c)))
# Area of polygon as a sum of triangle areas
def area_of_convex_hull(convex_hull):
poly_area = 0
for i in range(1, len(convex_hull) - 1, 1):
points = [convex_hull[0], convex_hull[i], convex_hull[i + 1]]
poly_area += area_of_triangle(points)
return poly_area
def rotate_points(points, origin, angle):
new_points = list(map(lambda x: x - origin, points))
for p in new_points:
s = math.sin(angle)
c = math.cos(angle)
px = (p.x * c) - (p.y * s)
py = (p.x * s) + (p.y * c)
p.x = px
p.y = py
return list(map(lambda x: x + origin, new_points))
def fit_obb(points):
min_area = None
obb = None
for i in range(len(points) - 1):
theta = -angle(points[i], points[i + 1])
origin = points[i]
rotated_polygon = rotate_points(points, origin, theta)
aabb = aabb_from_points(rotated_polygon)
area = aabb.area()
if min_area is None or area < min_area:
min_area = area
size = aabb.size()
rot = -theta
obb_vertices = rotate_points(aabb.vertices(), origin, rot)
center = (obb_vertices[0] + obb_vertices[2]) * 0.5
obb = OBB(center, rot, size)
return obb
class OBB(object):
def __init__(self, center, rot, size):
self.center = center
self.rot = rot
self.size = size
def vertices(self):
points = [
Vector([-self.size.x / 2, -self.size.y / 2, 0]),
Vector([ self.size.x / 2, -self.size.y / 2, 0]),
Vector([ self.size.x / 2, self.size.y / 2, 0]),
Vector([-self.size.x / 2, self.size.y / 2, 0]),
]
points = rotate_points(points, Vector([0, 0, 0]), self.rot)
return list(map(lambda x: x + self.center, points))
def aabb_from_points(points):
min_x = min(points, key=lambda p: p.x).x
min_y = min(points, key=lambda p: p.y).y
max_x = max(points, key=lambda p: p.x).x
max_y = max(points, key=lambda p: p.y).y
return AABB(Vector([min_x, min_y, 0]), Vector([max_x, max_y, 0]))
class AABB(object):
def __init__(self, min_, max_):
self.min = min_
self.max = max_
def center(self):
return (self.max + self.min) * 0.5
def area(self):
a = self.max - self.min
return a.x * a.y
def vertices(self):
return [
self.min,
Vector([self.max.x, self.min.y, 0]),
self.max,
Vector([self.min.x, self.max.y, 0])
]
def size(self):
return self.max - self.min
def set_rot_of_matrix(transform, rot):
transform[0, 0] = rot[0, 0]
transform[0, 1] = rot[0, 1]
transform[0, 2] = rot[0, 2]
transform[1, 0] = rot[1, 0]
transform[1, 1] = rot[1, 1]
transform[1, 2] = rot[1, 2]
transform[2, 0] = rot[2, 0]
transform[2, 1] = rot[2, 1]
transform[2, 2] = rot[2, 2]
def region_transform(chunk):
transform = Matrix.Translation(chunk.region.center)
set_rot_of_matrix(transform, chunk.region.rot)
return transform
def local_to_region(chunk, points):
transform = region_transform(chunk)
points = [transform.inv().mulp(x) for x in points]
return points
def convert_to_region(chunk, points):
T = chunk.transform.matrix
crs = chunk.shapes.crs
points = [T.inv().mulp(crs.unproject(x)) for x in points]
return local_to_region(chunk, points)
# Shrinks the region to a minimum bounding box around a polygon.
def conform_region_to_polygon(chunk, shape):
assert chunk is not None
assert shape is not None
assert shape.type is Shape.Polygon
assert len(shape.vertices) >= 3
chunk.resetRegion()
# Save the old region rotation as we will use it in the final region translation.
old_region_rot = chunk.region.rot
# Get the WSG vertices in region space
vertices = convert_to_region(chunk, chunk.shapes[0].vertices)
# Generate a convex hull around the shape vertices
convex_hull = get_convex_hull(vertices)
# Fit a minimal bounding box around the convex hull
obb = fit_obb(convex_hull)
# Maintain the chunk region height
obb.size.z = chunk.region.size.z
# Resize the chunk to the size of the OBB
chunk.region.size = obb.size
# Rotate the region by the rotation of the OBB
chunk.region.rot = chunk.region.rot * Utils.ypr2mat(Vector([math.degrees(obb.rot), 0, 0])).inv()
# Translate the region center to the OBB's center
transform = Matrix.Diag([1, 1, 1, 1])
set_rot_of_matrix(transform, old_region_rot)
chunk.region.center += transform.mulp(obb.center)
# Do the deed!
chunk = PhotoScan.app.document.chunk
conform_region_to_polygon(chunk, chunk.shapes[0])
Topic: Translating WGS84 coordinates to local coordinates (Read 3274 times)
