From 4f0d3ce83d6c46d2867bc052d5979c76303db176 Mon Sep 17 00:00:00 2001 From: s3lph Date: Tue, 17 Dec 2024 03:22:16 +0100 Subject: [PATCH] feat: tessellation of multiple multi-triangle faces --- map.py | 114 ++++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 76 insertions(+), 38 deletions(-) diff --git a/map.py b/map.py index 4adee05..fad2927 100644 --- a/map.py +++ b/map.py @@ -4,6 +4,8 @@ import numpy as np import sys import json +from tqdm import tqdm, trange + EPS = sys.float_info.epsilon @@ -32,23 +34,30 @@ class Triangle: 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 None + return False, None s = -self.v1 u = np.dot(s, ray_cross_e2) / dot - if (u < 0 and abs(u) > EPS) or (u > 1 and abs(u - 1) > EPS): - return None s_cross_e1 = np.cross(s, self.e1) v = np.dot(ray, s_cross_e1) / dot - if (v < 0 and abs(v) > EPS) or (u + v > 1 and abs(u + v - 1) > EPS): - return None 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 None - return ray * t + return False, None + return True, ray * t def plane_coords(self, v): u = np.dot(self.e1 / self.e1n, v - self.v1) @@ -57,6 +66,16 @@ class Triangle: return np.array([u, v]) +class Face: + + def __init__(self, vs): + self.vertices = vs + self.tris = [] + for i in range(1, len(vs)-1): + self.tris.append(Triangle(vs[0], vs[i], vs[i+1])) + self.plane_coords = [self.tris[0].plane_coords(v) for v in self.vertices] + + def convert_to_cartesian(lon, lat): x = cos(lon*pi/180) y = sin(lon*pi/180) @@ -65,52 +84,71 @@ def convert_to_cartesian(lon, lat): return np.array([x/mag, y/mag, z/mag]) -def map_poly(tri, poly): +def map_poly(tri, poly, ref=None): + if ref is None: + ref = tri p0 = convert_to_cartesian(*poly[0]) - p0i = tri.intersect(p0) + p0b, p0i = tri.intersect(p0) mapped = [] new = True for i in range(1, len(poly)): p1 = convert_to_cartesian(*poly[i]) - p1i = tri.intersect(p1) - if p0i is not None and p1i is not None: + p1b, p1i = tri.intersect(p1) + if p0b or p1b and p0i is not None and p1i is not None: if new: - mapped.append([tri.plane_coords(p0i), tri.plane_coords(p1i)]) + mapped.append([ref.plane_coords(p0i), ref.plane_coords(p1i)]) new = False else: - mapped[-1].append(tri.plane_coords(p1i)) + mapped[-1].append(ref.plane_coords(p1i)) else: new = True p0 = p1 p0i = p1i + p0b = p1b return mapped if __name__ == '__main__': - t = Triangle(np.array([5,0,0]), np.array([0,5,0]), np.array([0,0,5])) - countries = {} + fs = [Face([a, b, c, d]) for a, b, c, d in [ + (np.array([1,1,1]), np.array([-1,1,1]), np.array([-1,-1,1]), np.array([1,-1,1])), # top + (np.array([1,1,-1]), np.array([-1,1,-1]), np.array([-1,-1,-1]), np.array([1,-1,-1])), # bottom + (np.array([-1,-1,-1]), np.array([-1,1,-1]), np.array([-1,1,1]), np.array([-1,-1,1])), # left + (np.array([1,-1,-1]), np.array([1,1,-1]), np.array([1,1,1]), np.array([1,-1,1])), # right + (np.array([-1,-1,-1]), np.array([1,-1,-1]), np.array([1,-1,1]), np.array([-1,-1,1])), # front + (np.array([-1,1,-1]), np.array([1,1,-1]), np.array([1,1,1]), np.array([-1,1,1])), # back + ]] with open('countries.geojson', 'r') as f: j = json.load(f) - for f in j['features']: - cc = f['properties']['ADMIN'] - if f['geometry']['type'] == 'MultiPolygon': - countries[cc] = [] - for poly in f['geometry']['coordinates']: - countries[cc].extend(map_poly(t, poly[0])) - elif f['geometry']['type'] == 'Polygon': - countries[cc] = map_poly(tri, f['geometry']['coordinates'][0]) - minx = min(t.v1p[0], t.v2p[0], t.v3p[0]) - miny = min(t.v1p[1], t.v2p[1], t.v3p[1]) - maxx = max(t.v1p[0], t.v2p[0], t.v3p[0]) - maxy = max(t.v1p[1], t.v2p[1], t.v3p[1]) - w = maxx - minx - h = maxy - miny - print(f'') - print(f'') - for c, l in countries.items(): - for ls in l: - print(f'') - print('') + + for i, face in tqdm(list(enumerate(fs)), desc='Faces '): + countries = {} + for t in face.tris: + for f in tqdm(j['features'], desc='Features', total=len(face.tris)*len(j['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 + with open(f'face{i}.svg', 'w') as svg: + svg.write(f'\n') + svg.write(f' \n \n \n \n') + for c, l in countries.items(): + if l: + svg.write(f' \n') + for ls in l: + svg.write(f' \n') + svg.write(' \n') + svg.write(' \n\n')