190 lines
7.1 KiB
Python
190 lines
7.1 KiB
Python
|
||
import argparse
|
||
import json
|
||
from math import *
|
||
import os
|
||
import sys
|
||
|
||
import numpy as np
|
||
from plyfile import PlyData
|
||
from tqdm import tqdm, trange
|
||
|
||
|
||
EPS = sys.float_info.epsilon
|
||
|
||
|
||
class Triangle:
|
||
|
||
def __init__(self, v1, v2, v3):
|
||
# Vertices
|
||
self.v1 = v1
|
||
self.v2 = v2
|
||
self.v3 = v3
|
||
# Sort triangle vertices so that so that the normal points "outward"
|
||
if np.dot(self.v1, np.cross(self.v2 - self.v1, self.v3 - self.v1)) < 0:
|
||
self.v3, self.v2 = self.v2, self.v3
|
||
# Edges of the origin vertex and their L2 norms
|
||
self.e1 = self.v2 - self.v1
|
||
self.e2 = self.v3 - self.v1
|
||
self.e1n = np.linalg.norm(self.e1)
|
||
self.e2n = np.linalg.norm(self.e2)
|
||
# Normal vector of the triangle's plane
|
||
self.n = np.cross(self.e1 / self.e1n, self.e2 / self.e2n)
|
||
self.n = self.n / np.linalg.norm(self.n)
|
||
# The vertex coordinates in plane coordinates
|
||
self.v1p = self.plane_coords(self.v1)
|
||
self.v2p = self.plane_coords(self.v2)
|
||
self.v3p = self.plane_coords(self.v3)
|
||
|
||
def intersect(self, ray):
|
||
'''
|
||
:param ray: numpy.array(3) Ray cast from [0, 0, 0]
|
||
:return: A tuple of [bool, Optional[numpy.array(3)]]:
|
||
The first component is True if the ray intersects the triangle, False otherwise.
|
||
The second component marks the intersection of the triangle's plane (even if outside the triangle), or
|
||
None if the ray does not intersect the plane at all.
|
||
'''
|
||
# https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm
|
||
ray_cross_e2 = np.cross(ray, self.e2)
|
||
dot = np.dot(self.e1, ray_cross_e2)
|
||
if abs(dot) < EPS:
|
||
return False, None
|
||
s = -self.v1
|
||
u = np.dot(s, ray_cross_e2) / dot
|
||
s_cross_e1 = np.cross(s, self.e1)
|
||
v = np.dot(ray, s_cross_e1) / dot
|
||
t = np.dot(self.e2, s_cross_e1) / dot
|
||
if (u < 0 and abs(u) > EPS) or (u > 1 and abs(u - 1) > EPS):
|
||
return False, ray * t
|
||
if (v < 0 and abs(v) > EPS) or (u + v > 1 and abs(u + v - 1) > EPS):
|
||
return False, ray * t
|
||
if t <= EPS:
|
||
return False, None
|
||
return True, ray * t
|
||
|
||
def plane_coords(self, v):
|
||
u = np.dot(self.e1 / self.e1n, v - self.v1)
|
||
n = np.cross(self.n, self.e1 / self.e1n)
|
||
v = -np.dot(n, v - self.v1)
|
||
return np.array([u, v])
|
||
|
||
|
||
class Face:
|
||
|
||
def __init__(self, vs):
|
||
self.vertices = vs
|
||
self.tris = []
|
||
avg = np.average(self.vertices, axis=0)
|
||
self.center = convert_to_spherical(*avg)
|
||
# Sort vertices so that so that the normal points "outward"
|
||
if np.dot(self.vertices[0], np.cross(self.vertices[1] - self.vertices[0], self.vertices[2] - self.vertices[0])) < 0:
|
||
self.vertices = self.vertices[::-1]
|
||
for i in range(1, len(self.vertices)-1):
|
||
self.tris.append(Triangle(self.vertices[0], self.vertices[i], self.vertices[i+1]))
|
||
self.plane_coords = [self.tris[0].plane_coords(v) for v in self.vertices]
|
||
|
||
|
||
def convert_to_spherical(x, y, z):
|
||
xy = sqrt(x**2+y**2)
|
||
lon = atan2(z, xy)*180/pi
|
||
lat = atan2(y, x)*180/pi
|
||
mag = sqrt(x**2+y**2+z**2)
|
||
return np.array([lon, lat, mag])
|
||
|
||
|
||
def convert_to_cartesian(lon, lat):
|
||
x = cos(lon*pi/180)
|
||
y = sin(lon*pi/180)
|
||
z = tan(lat*pi/180)
|
||
mag = sqrt(x**2 + y**2 + z**2)
|
||
return np.array([x/mag, y/mag, z/mag])
|
||
|
||
|
||
def map_poly(tri, poly, ref=None):
|
||
if ref is None:
|
||
ref = tri
|
||
p0 = convert_to_cartesian(*poly[0])
|
||
p0b, p0i = tri.intersect(p0)
|
||
mapped = []
|
||
new = True
|
||
for i in range(1, len(poly)):
|
||
p1 = convert_to_cartesian(*poly[i])
|
||
p1b, p1i = tri.intersect(p1)
|
||
if p0b or p1b and p0i is not None and p1i is not None:
|
||
if new:
|
||
mapped.append([ref.plane_coords(p0i), ref.plane_coords(p1i)])
|
||
new = False
|
||
else:
|
||
mapped[-1].append(ref.plane_coords(p1i))
|
||
else:
|
||
new = True
|
||
p0 = p1
|
||
p0i = p1i
|
||
p0b = p1b
|
||
return mapped
|
||
|
||
|
||
def main(ns):
|
||
with open(ns.faces, 'rb') as f:
|
||
plydata = PlyData.read(f)
|
||
vs = plydata['vertex']
|
||
fs = plydata['face']['vertex_indices']
|
||
faces = [Face([np.array(vs[vi].tolist())*ns.scale for vi in fvs]) for fvs in fs]
|
||
with open(ns.geojson, 'r') as f:
|
||
geojson = json.load(f)
|
||
|
||
for i, face in tqdm(list(enumerate(faces)), desc='Faces '):
|
||
countries = {}
|
||
for t in face.tris:
|
||
for f in tqdm(geojson['features'], desc='Features', total=len(face.tris)*len(geojson['features'])):
|
||
cc = f['properties']['ADMIN']
|
||
if f['geometry']['type'] == 'MultiPolygon':
|
||
for poly in f['geometry']['coordinates']:
|
||
countries.setdefault(cc, []).extend(map_poly(t, poly[0], face.tris[0]))
|
||
elif f['geometry']['type'] == 'Polygon':
|
||
countries.setdefault(cc, []).extend(map_poly(t, f['geometry']['coordinates'][0], face.tris[0]))
|
||
minx = min(v[0] for v in face.plane_coords)
|
||
miny = min(v[1] for v in face.plane_coords)
|
||
maxx = max(v[0] for v in face.plane_coords)
|
||
maxy = max(v[1] for v in face.plane_coords)
|
||
w = maxx - minx
|
||
h = maxy - miny
|
||
formatters = dict(
|
||
faces=os.path.basename(ns.faces).rsplit('.', 1)[0],
|
||
geojson=ns.geojson,
|
||
face=i,
|
||
lon=face.center[0],
|
||
lat=face.center[1],
|
||
mag=face.center[2],
|
||
ilon=int(face.center[0]),
|
||
ilat=int(face.center[1]),
|
||
imag=int(face.center[2]),
|
||
)
|
||
outfile = ns.out.format(**formatters)
|
||
os.makedirs(os.path.dirname(outfile), exist_ok=True)
|
||
with open(outfile, 'w') as svg:
|
||
svg.write(f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="{minx} {miny} {w} {h}" width="{w}" height="{h}" id="face{i}">\n')
|
||
svg.write(f' <g id="cut">\n <path id="outline" fill="none" stroke="black" stroke-width="{0.01*ns.scale}" d="M {face.plane_coords[0][0]} {face.plane_coords[0][1]}')
|
||
for v in face.plane_coords[1:]:
|
||
svg.write(f' L {v[0]} {v[1]}')
|
||
svg.write(' Z" />\n </g>\n <g id="engrave">\n')
|
||
for c, l in countries.items():
|
||
if l:
|
||
svg.write(f' <g id="{c}">\n')
|
||
for ls in l:
|
||
svg.write(f' <path fill="none" stroke="red" stroke-width="{0.01*ns.scale}" d="M {ls[0][0]} {ls[0][1]}')
|
||
for x, y in ls[1:]:
|
||
svg.write(f' L {x} {y}')
|
||
svg.write('" />\n')
|
||
svg.write(' </g>\n')
|
||
svg.write(' </g>\n</svg>\n')
|
||
|
||
|
||
if __name__ == '__main__':
|
||
ap = argparse.ArgumentParser()
|
||
ap.add_argument('--geojson', '-g', default='countries.geojson')
|
||
ap.add_argument('--faces', '-f', default='shapes/cube.json')
|
||
ap.add_argument('--out', '-o', default='faces/{faces}/face{face}_{ilat}_{ilon}.svg')
|
||
ap.add_argument('--scale', '-s', type=float, default=1)
|
||
ns = ap.parse_args()
|
||
main(ns)
|