684 lines
27 KiB
Python
684 lines
27 KiB
Python
|
#!/usr/bin/env python3
|
|||
|
|
# -*- coding: utf-8 -*-
|
|||
|
|
"""
|
|||
|
|
三维模型体积计算器
|
|||
|
|
基于指定地理范围内的三维模型数据，使用三角构网方法计算体积
|
|||
|
|
"""
|
|||
|
|
|
|||
|
|
import json
|
|||
|
|
import struct
|
|||
|
|
import os
|
|||
|
|
import sys
|
|||
|
|
from pathlib import Path
|
|||
|
|
import numpy as np
|
|||
|
|
from scipy.spatial import Delaunay
|
|||
|
|
import math
|
|||
|
|
from tqdm import tqdm
|
|||
|
|
|
|||
|
|
try:
|
|||
|
|
import DracoPy
|
|||
|
|
except ImportError:
|
|||
|
|
print("警告: DracoPy库未安装，无法处理Draco压缩的数据")
|
|||
|
|
print("请运行: pip install DracoPy")
|
|||
|
|
DracoPy = None
|
|||
|
|
|
|||
|
|
class VolumeCalculator:
|
|||
|
|
def __init__(self, location_file):
|
|||
|
|
self.location_bounds = self.load_location_bounds(location_file)
|
|||
|
|
self.all_vertices = []
|
|||
|
|
self.filtered_vertices = []
|
|||
|
|
|
|||
|
|
def load_location_bounds(self, location_file):
|
|||
|
|
"""加载地理范围边界"""
|
|||
|
|
try:
|
|||
|
|
with open(location_file, 'r', encoding='utf-8') as f:
|
|||
|
|
coords = json.load(f)
|
|||
|
|
|
|||
|
|
# 提取经纬度范围
|
|||
|
|
lons = [coord[0] for coord in coords]
|
|||
|
|
lats = [coord[1] for coord in coords]
|
|||
|
|
elevs = [coord[2] for coord in coords]
|
|||
|
|
|
|||
|
|
bounds = {
|
|||
|
|
'min_lon': min(lons),
|
|||
|
|
'max_lon': max(lons),
|
|||
|
|
'min_lat': min(lats),
|
|||
|
|
'max_lat': max(lats),
|
|||
|
|
'min_elev': min(elevs),
|
|||
|
|
'max_elev': max(elevs),
|
|||
|
|
'coords': coords
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
print(f"地理范围边界:")
|
|||
|
|
print(f" 经度: {bounds['min_lon']:.8f} ~ {bounds['max_lon']:.8f}")
|
|||
|
|
print(f" 纬度: {bounds['min_lat']:.8f} ~ {bounds['max_lat']:.8f}")
|
|||
|
|
print(f" 高程: {bounds['min_elev']:.2f} ~ {bounds['max_elev']:.2f}")
|
|||
|
|
|
|||
|
|
return bounds
|
|||
|
|
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"加载地理范围文件失败: {e}")
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
def wgs84_to_cartesian(self, lon, lat, elev):
|
|||
|
|
"""WGS84坐标转换为笛卡尔坐标（高精度算法）"""
|
|||
|
|
# WGS84椭球参数（EPSG:4326）
|
|||
|
|
a = 6378137.0 # 长半轴 (米)
|
|||
|
|
f = 1/298.257223563 # 扁率
|
|||
|
|
e2 = 2*f - f*f # 第一偏心率平方
|
|||
|
|
|
|||
|
|
# 角度转弧度
|
|||
|
|
lon_rad = math.radians(lon)
|
|||
|
|
lat_rad = math.radians(lat)
|
|||
|
|
|
|||
|
|
# 计算卯酉圈曲率半径
|
|||
|
|
sin_lat = math.sin(lat_rad)
|
|||
|
|
cos_lat = math.cos(lat_rad)
|
|||
|
|
sin_lon = math.sin(lon_rad)
|
|||
|
|
cos_lon = math.cos(lon_rad)
|
|||
|
|
|
|||
|
|
N = a / math.sqrt(1 - e2 * sin_lat * sin_lat)
|
|||
|
|
|
|||
|
|
# 高精度笛卡尔坐标计算
|
|||
|
|
x = (N + elev) * cos_lat * cos_lon
|
|||
|
|
y = (N + elev) * cos_lat * sin_lon
|
|||
|
|
z = (N * (1 - e2) + elev) * sin_lat
|
|||
|
|
|
|||
|
|
return [x, y, z]
|
|||
|
|
|
|||
|
|
def cartesian_to_wgs84(self, x, y, z):
|
|||
|
|
"""笛卡尔坐标转换为WGS84坐标（高精度迭代法）"""
|
|||
|
|
# WGS84椭球参数
|
|||
|
|
a = 6378137.0 # 长半轴
|
|||
|
|
f = 1/298.257223563 # 扁率
|
|||
|
|
e2 = 2*f - f*f # 第一偏心率平方
|
|||
|
|
ep2 = e2 / (1 - e2) # 第二偏心率平方
|
|||
|
|
|
|||
|
|
# 计算经度（精确值）
|
|||
|
|
lon = math.atan2(y, x)
|
|||
|
|
|
|||
|
|
# 计算纬度和高程（使用Bowring迭代法）
|
|||
|
|
p = math.sqrt(x*x + y*y)
|
|||
|
|
|
|||
|
|
if p == 0:
|
|||
|
|
# 极点情况
|
|||
|
|
lat = math.pi/2 if z > 0 else -math.pi/2
|
|||
|
|
elev = abs(z) - a * math.sqrt(1 - e2)
|
|||
|
|
return [math.degrees(lon), math.degrees(lat), elev]
|
|||
|
|
|
|||
|
|
# 初始估计
|
|||
|
|
theta = math.atan2(z, p * (1 - f))
|
|||
|
|
lat_prev = math.atan2(z + ep2 * a * (1 - f) * math.sin(theta)**3,
|
|||
|
|
p - e2 * a * math.cos(theta)**3)
|
|||
|
|
|
|||
|
|
# 迭代求解纬度
|
|||
|
|
max_iterations = 10
|
|||
|
|
tolerance = 1e-12
|
|||
|
|
|
|||
|
|
for i in range(max_iterations):
|
|||
|
|
N = a / math.sqrt(1 - e2 * math.sin(lat_prev)**2)
|
|||
|
|
elev = p / math.cos(lat_prev) - N
|
|||
|
|
|
|||
|
|
# 更新纬度估计
|
|||
|
|
lat_new = math.atan2(z + e2 * N * math.sin(lat_prev), p)
|
|||
|
|
|
|||
|
|
# 检查收敛性
|
|||
|
|
if abs(lat_new - lat_prev) < tolerance:
|
|||
|
|
break
|
|||
|
|
|
|||
|
|
lat_prev = lat_new
|
|||
|
|
|
|||
|
|
# 最终计算高程
|
|||
|
|
N = a / math.sqrt(1 - e2 * math.sin(lat_prev)**2)
|
|||
|
|
elev = p / math.cos(lat_prev) - N
|
|||
|
|
|
|||
|
|
return [math.degrees(lon), math.degrees(lat_prev), elev]
|
|||
|
|
|
|||
|
|
def is_point_in_bounds(self, vertex):
|
|||
|
|
"""检查点是否在指定的地理范围内"""
|
|||
|
|
if not self.location_bounds:
|
|||
|
|
return True
|
|||
|
|
|
|||
|
|
# 将笛卡尔坐标转换为WGS84
|
|||
|
|
try:
|
|||
|
|
lon, lat, elev = self.cartesian_to_wgs84(vertex[0], vertex[1], vertex[2])
|
|||
|
|
|
|||
|
|
# 检查是否在边界范围内
|
|||
|
|
return (self.location_bounds['min_lon'] <= lon <= self.location_bounds['max_lon'] and
|
|||
|
|
self.location_bounds['min_lat'] <= lat <= self.location_bounds['max_lat'])
|
|||
|
|
except:
|
|||
|
|
return False
|
|||
|
|
|
|||
|
|
def multiply_matrix_vector(self, matrix, vector):
|
|||
|
|
"""4x4矩阵与4D向量相乘"""
|
|||
|
|
m = [
|
|||
|
|
[matrix[0], matrix[4], matrix[8], matrix[12]],
|
|||
|
|
[matrix[1], matrix[5], matrix[9], matrix[13]],
|
|||
|
|
[matrix[2], matrix[6], matrix[10], matrix[14]],
|
|||
|
|
[matrix[3], matrix[7], matrix[11], matrix[15]]
|
|||
|
|
]
|
|||
|
|
|
|||
|
|
v = [vector[0], vector[1], vector[2], 1.0]
|
|||
|
|
|
|||
|
|
result = [
|
|||
|
|
m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2] + m[0][3]*v[3],
|
|||
|
|
m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2] + m[1][3]*v[3],
|
|||
|
|
m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2] + m[2][3]*v[3]
|
|||
|
|
]
|
|||
|
|
|
|||
|
|
return result
|
|||
|
|
|
|||
|
|
def multiply_matrices(self, m1, m2):
|
|||
|
|
"""两个4x4矩阵相乘"""
|
|||
|
|
def list_to_matrix(lst):
|
|||
|
|
return [
|
|||
|
|
[lst[0], lst[4], lst[8], lst[12]],
|
|||
|
|
[lst[1], lst[5], lst[9], lst[13]],
|
|||
|
|
[lst[2], lst[6], lst[10], lst[14]],
|
|||
|
|
[lst[3], lst[7], lst[11], lst[15]]
|
|||
|
|
]
|
|||
|
|
|
|||
|
|
def matrix_to_list(mat):
|
|||
|
|
return [
|
|||
|
|
mat[0][0], mat[1][0], mat[2][0], mat[3][0],
|
|||
|
|
mat[0][1], mat[1][1], mat[2][1], mat[3][1],
|
|||
|
|
mat[0][2], mat[1][2], mat[2][2], mat[3][2],
|
|||
|
|
mat[0][3], mat[1][3], mat[2][3], mat[3][3]
|
|||
|
|
]
|
|||
|
|
|
|||
|
|
a = list_to_matrix(m1)
|
|||
|
|
b = list_to_matrix(m2)
|
|||
|
|
|
|||
|
|
result = [[0 for _ in range(4)] for _ in range(4)]
|
|||
|
|
for i in range(4):
|
|||
|
|
for j in range(4):
|
|||
|
|
for k in range(4):
|
|||
|
|
result[i][j] += a[i][k] * b[k][j]
|
|||
|
|
|
|||
|
|
return matrix_to_list(result)
|
|||
|
|
|
|||
|
|
def apply_transform_to_vertices(self, vertices, transform_matrix):
|
|||
|
|
"""对顶点应用变换矩阵"""
|
|||
|
|
if not transform_matrix:
|
|||
|
|
return vertices
|
|||
|
|
|
|||
|
|
transformed_vertices = []
|
|||
|
|
for vertex in vertices:
|
|||
|
|
transformed = self.multiply_matrix_vector(transform_matrix, vertex)
|
|||
|
|
transformed_vertices.append(transformed)
|
|||
|
|
|
|||
|
|
return transformed_vertices
|
|||
|
|
|
|||
|
|
def parse_tileset_json(self, tileset_path, parent_transform=None):
|
|||
|
|
"""解析tileset.json文件，收集B3DM文件和变换矩阵"""
|
|||
|
|
try:
|
|||
|
|
with open(tileset_path, 'r', encoding='utf-8') as f:
|
|||
|
|
tileset_data = json.load(f)
|
|||
|
|
|
|||
|
|
b3dm_files = []
|
|||
|
|
|
|||
|
|
def process_node(node, base_path, accumulated_transform):
|
|||
|
|
current_transform = node.get('transform')
|
|||
|
|
|
|||
|
|
if current_transform and accumulated_transform:
|
|||
|
|
final_transform = self.multiply_matrices(accumulated_transform, current_transform)
|
|||
|
|
elif current_transform:
|
|||
|
|
final_transform = current_transform
|
|||
|
|
else:
|
|||
|
|
final_transform = accumulated_transform
|
|||
|
|
|
|||
|
|
if 'content' in node and 'uri' in node['content']:
|
|||
|
|
uri = node['content']['uri']
|
|||
|
|
if uri.endswith('.b3dm'):
|
|||
|
|
full_path = os.path.join(base_path, uri)
|
|||
|
|
if os.path.exists(full_path):
|
|||
|
|
b3dm_files.append((full_path, final_transform))
|
|||
|
|
elif uri.endswith('.json'):
|
|||
|
|
sub_tileset_path = os.path.join(base_path, uri)
|
|||
|
|
if os.path.exists(sub_tileset_path):
|
|||
|
|
sub_files = self.parse_tileset_json(sub_tileset_path, final_transform)
|
|||
|
|
b3dm_files.extend(sub_files)
|
|||
|
|
|
|||
|
|
if 'children' in node:
|
|||
|
|
for child in node['children']:
|
|||
|
|
process_node(child, base_path, final_transform)
|
|||
|
|
|
|||
|
|
base_path = os.path.dirname(tileset_path)
|
|||
|
|
if 'root' in tileset_data:
|
|||
|
|
process_node(tileset_data['root'], base_path, parent_transform)
|
|||
|
|
|
|||
|
|
return b3dm_files
|
|||
|
|
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"解析tileset.json时出错: {e}")
|
|||
|
|
return []
|
|||
|
|
|
|||
|
|
def parse_b3dm_file(self, file_path):
|
|||
|
|
"""解析B3DM文件"""
|
|||
|
|
try:
|
|||
|
|
with open(file_path, 'rb') as f:
|
|||
|
|
magic = f.read(4)
|
|||
|
|
if magic != b'b3dm':
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
version = struct.unpack('<I', f.read(4))[0]
|
|||
|
|
byte_length = struct.unpack('<I', f.read(4))[0]
|
|||
|
|
feature_table_json_byte_length = struct.unpack('<I', f.read(4))[0]
|
|||
|
|
feature_table_binary_byte_length = struct.unpack('<I', f.read(4))[0]
|
|||
|
|
batch_table_json_byte_length = struct.unpack('<I', f.read(4))[0]
|
|||
|
|
batch_table_binary_byte_length = struct.unpack('<I', f.read(4))[0]
|
|||
|
|
|
|||
|
|
f.seek(28 + feature_table_json_byte_length + feature_table_binary_byte_length +
|
|||
|
|
batch_table_json_byte_length + batch_table_binary_byte_length)
|
|||
|
|
|
|||
|
|
gltf_data = f.read()
|
|||
|
|
return self.parse_gltf_data(gltf_data)
|
|||
|
|
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"解析B3DM文件 {file_path} 失败: {e}")
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
def parse_gltf_data(self, gltf_data):
|
|||
|
|
"""解析glTF数据"""
|
|||
|
|
try:
|
|||
|
|
if gltf_data[:4] == b'glTF':
|
|||
|
|
return self.parse_glb_data(gltf_data)
|
|||
|
|
else:
|
|||
|
|
gltf_json = json.loads(gltf_data.decode('utf-8'))
|
|||
|
|
return self.extract_vertices_from_gltf(gltf_json, None)
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"解析glTF数据失败: {e}")
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
def parse_glb_data(self, glb_data):
|
|||
|
|
"""解析GLB格式的glTF数据"""
|
|||
|
|
try:
|
|||
|
|
magic = glb_data[:4]
|
|||
|
|
if magic != b'glTF':
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
version = struct.unpack('<I', glb_data[4:8])[0]
|
|||
|
|
total_length = struct.unpack('<I', glb_data[8:12])[0]
|
|||
|
|
|
|||
|
|
offset = 12
|
|||
|
|
json_data = None
|
|||
|
|
binary_data = None
|
|||
|
|
|
|||
|
|
while offset < len(glb_data):
|
|||
|
|
if offset + 8 > len(glb_data):
|
|||
|
|
break
|
|||
|
|
|
|||
|
|
chunk_length = struct.unpack('<I', glb_data[offset:offset+4])[0]
|
|||
|
|
chunk_type = glb_data[offset+4:offset+8]
|
|||
|
|
chunk_data = glb_data[offset+8:offset+8+chunk_length]
|
|||
|
|
|
|||
|
|
if chunk_type == b'JSON':
|
|||
|
|
json_data = json.loads(chunk_data.decode('utf-8'))
|
|||
|
|
elif chunk_type == b'BIN\x00':
|
|||
|
|
binary_data = chunk_data
|
|||
|
|
|
|||
|
|
offset += 8 + chunk_length
|
|||
|
|
offset = (offset + 3) & ~3
|
|||
|
|
|
|||
|
|
if json_data:
|
|||
|
|
return self.extract_vertices_from_gltf(json_data, binary_data)
|
|||
|
|
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"解析GLB数据失败: {e}")
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
def extract_vertices_from_gltf(self, gltf_json, binary_data):
|
|||
|
|
"""从glTF JSON中提取顶点数据"""
|
|||
|
|
vertices = []
|
|||
|
|
|
|||
|
|
try:
|
|||
|
|
if 'extensionsUsed' in gltf_json and 'KHR_draco_mesh_compression' in gltf_json['extensionsUsed']:
|
|||
|
|
if DracoPy is None:
|
|||
|
|
print("警告: 检测到Draco压缩但DracoPy未安装")
|
|||
|
|
return vertices
|
|||
|
|
return self.extract_draco_vertices(gltf_json, binary_data)
|
|||
|
|
|
|||
|
|
if 'meshes' not in gltf_json:
|
|||
|
|
return vertices
|
|||
|
|
|
|||
|
|
for mesh in gltf_json['meshes']:
|
|||
|
|
for primitive in mesh['primitives']:
|
|||
|
|
if 'attributes' in primitive and 'POSITION' in primitive['attributes']:
|
|||
|
|
position_accessor_index = primitive['attributes']['POSITION']
|
|||
|
|
|
|||
|
|
if 'accessors' in gltf_json and position_accessor_index < len(gltf_json['accessors']):
|
|||
|
|
accessor = gltf_json['accessors'][position_accessor_index]
|
|||
|
|
buffer_view_index = accessor['bufferView']
|
|||
|
|
|
|||
|
|
if 'bufferViews' in gltf_json and buffer_view_index < len(gltf_json['bufferViews']):
|
|||
|
|
buffer_view = gltf_json['bufferViews'][buffer_view_index]
|
|||
|
|
buffer_index = buffer_view['buffer']
|
|||
|
|
byte_offset = buffer_view.get('byteOffset', 0) + accessor.get('byteOffset', 0)
|
|||
|
|
|
|||
|
|
if binary_data and buffer_index == 0:
|
|||
|
|
component_type = accessor['componentType']
|
|||
|
|
count = accessor['count']
|
|||
|
|
|
|||
|
|
if component_type == 5126: # FLOAT
|
|||
|
|
vertex_data = struct.unpack(f'<{count*3}f',
|
|||
|
|
binary_data[byte_offset:byte_offset+count*12])
|
|||
|
|
for i in range(0, len(vertex_data), 3):
|
|||
|
|
vertices.append([vertex_data[i], vertex_data[i+1], vertex_data[i+2]])
|
|||
|
|
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"提取顶点数据失败: {e}")
|
|||
|
|
|
|||
|
|
return vertices
|
|||
|
|
|
|||
|
|
def extract_draco_vertices(self, gltf_json, binary_data):
|
|||
|
|
"""提取Draco压缩的顶点数据"""
|
|||
|
|
vertices = []
|
|||
|
|
|
|||
|
|
try:
|
|||
|
|
if 'meshes' not in gltf_json:
|
|||
|
|
return vertices
|
|||
|
|
|
|||
|
|
for mesh in gltf_json['meshes']:
|
|||
|
|
for primitive in mesh['primitives']:
|
|||
|
|
if 'extensions' in primitive and 'KHR_draco_mesh_compression' in primitive['extensions']:
|
|||
|
|
draco_ext = primitive['extensions']['KHR_draco_mesh_compression']
|
|||
|
|
buffer_view_index = draco_ext['bufferView']
|
|||
|
|
|
|||
|
|
if 'bufferViews' in gltf_json and buffer_view_index < len(gltf_json['bufferViews']):
|
|||
|
|
buffer_view = gltf_json['bufferViews'][buffer_view_index]
|
|||
|
|
byte_offset = buffer_view.get('byteOffset', 0)
|
|||
|
|
byte_length = buffer_view['byteLength']
|
|||
|
|
|
|||
|
|
if binary_data:
|
|||
|
|
draco_data = binary_data[byte_offset:byte_offset+byte_length]
|
|||
|
|
|
|||
|
|
mesh_data = DracoPy.decode(draco_data)
|
|||
|
|
if hasattr(mesh_data, 'points'):
|
|||
|
|
points = mesh_data.points
|
|||
|
|
for point in points:
|
|||
|
|
vertices.append([float(point[0]), float(point[1]), float(point[2])])
|
|||
|
|
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"解码Draco数据失败: {e}")
|
|||
|
|
|
|||
|
|
return vertices
|
|||
|
|
|
|||
|
|
def calculate_triangle_angles(self, p1, p2, p3):
|
|||
|
|
"""计算三角形的三个内角（度）"""
|
|||
|
|
# 计算三边长度
|
|||
|
|
a = np.linalg.norm(p2 - p3) # 边a对应角A(p1)
|
|||
|
|
b = np.linalg.norm(p1 - p3) # 边b对应角B(p2)
|
|||
|
|
c = np.linalg.norm(p1 - p2) # 边c对应角C(p3)
|
|||
|
|
|
|||
|
|
# 避免除零错误
|
|||
|
|
if a == 0 or b == 0 or c == 0:
|
|||
|
|
return [0, 0, 0]
|
|||
|
|
|
|||
|
|
# 使用余弦定理计算角度
|
|||
|
|
try:
|
|||
|
|
# 角A = arccos((b²+c²-a²)/(2bc))
|
|||
|
|
cos_A = (b*b + c*c - a*a) / (2*b*c)
|
|||
|
|
cos_B = (a*a + c*c - b*b) / (2*a*c)
|
|||
|
|
cos_C = (a*a + b*b - c*c) / (2*a*b)
|
|||
|
|
|
|||
|
|
# 限制余弦值范围，避免数值误差
|
|||
|
|
cos_A = np.clip(cos_A, -1.0, 1.0)
|
|||
|
|
cos_B = np.clip(cos_B, -1.0, 1.0)
|
|||
|
|
cos_C = np.clip(cos_C, -1.0, 1.0)
|
|||
|
|
|
|||
|
|
angle_A = np.arccos(cos_A) * 180 / np.pi
|
|||
|
|
angle_B = np.arccos(cos_B) * 180 / np.pi
|
|||
|
|
angle_C = np.arccos(cos_C) * 180 / np.pi
|
|||
|
|
|
|||
|
|
return [angle_A, angle_B, angle_C]
|
|||
|
|
except:
|
|||
|
|
return [0, 0, 0]
|
|||
|
|
|
|||
|
|
def is_valid_triangle(self, p1, p2, p3, min_angle=10.0, max_aspect_ratio=10.0):
|
|||
|
|
"""验证三角形质量，基于角度约束和长宽比"""
|
|||
|
|
angles = self.calculate_triangle_angles(p1, p2, p3)
|
|||
|
|
|
|||
|
|
# 检查最小角度约束（参考C#代码中的10度限制）
|
|||
|
|
if min(angles) < min_angle:
|
|||
|
|
return False
|
|||
|
|
|
|||
|
|
# 计算三边长度
|
|||
|
|
a = np.linalg.norm(p2 - p3)
|
|||
|
|
b = np.linalg.norm(p1 - p3)
|
|||
|
|
c = np.linalg.norm(p1 - p2)
|
|||
|
|
|
|||
|
|
# 检查长宽比约束
|
|||
|
|
max_side = max(a, b, c)
|
|||
|
|
min_side = min(a, b, c)
|
|||
|
|
if min_side > 0 and max_side / min_side > max_aspect_ratio:
|
|||
|
|
return False
|
|||
|
|
|
|||
|
|
return True
|
|||
|
|
|
|||
|
|
def calculate_circumcenter_and_radius(self, p1, p2, p3):
|
|||
|
|
"""计算三角形外接圆圆心和半径（高精度算法）"""
|
|||
|
|
try:
|
|||
|
|
x1, y1 = p1[0], p1[1]
|
|||
|
|
x2, y2 = p2[0], p2[1]
|
|||
|
|
x3, y3 = p3[0], p3[1]
|
|||
|
|
|
|||
|
|
# 使用C#代码中的高精度外接圆计算公式
|
|||
|
|
d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
|
|||
|
|
|
|||
|
|
if abs(d) < 1e-10: # 三点共线
|
|||
|
|
return None, float('inf')
|
|||
|
|
|
|||
|
|
ux = ((x1*x1 + y1*y1) * (y2 - y3) + (x2*x2 + y2*y2) * (y3 - y1) + (x3*x3 + y3*y3) * (y1 - y2)) / d
|
|||
|
|
uy = ((x1*x1 + y1*y1) * (x3 - x2) + (x2*x2 + y2*y2) * (x1 - y3) + (x3*x3 + y3*y3) * (x2 - x1)) / d
|
|||
|
|
|
|||
|
|
# 计算半径
|
|||
|
|
radius = np.sqrt((ux - x1)**2 + (uy - y1)**2)
|
|||
|
|
|
|||
|
|
return np.array([ux, uy]), radius
|
|||
|
|
except:
|
|||
|
|
return None, float('inf')
|
|||
|
|
|
|||
|
|
def calculate_volume_delaunay(self, vertices, base_elevation=None, min_angle=10.0, use_quality_filter=True):
|
|||
|
|
"""使用Delaunay三角剖分计算体积"""
|
|||
|
|
if len(vertices) < 4:
|
|||
|
|
print("顶点数量不足，无法进行三角剖分")
|
|||
|
|
return 0.0
|
|||
|
|
|
|||
|
|
try:
|
|||
|
|
points = np.array(vertices)
|
|||
|
|
|
|||
|
|
if base_elevation is None:
|
|||
|
|
base_elevation = np.min(points[:, 2])
|
|||
|
|
|
|||
|
|
print(f"Delaunay方法使用基准面高程: {base_elevation:.2f} 米")
|
|||
|
|
print(f"质量控制参数: 最小角度={min_angle}°, 质量过滤={'开启' if use_quality_filter else '关闭'}")
|
|||
|
|
|
|||
|
|
adjusted_points = points.copy()
|
|||
|
|
adjusted_points[:, 2] = np.maximum(points[:, 2] - base_elevation, 0)
|
|||
|
|
|
|||
|
|
print("正在进行Delaunay三角剖分...")
|
|||
|
|
tri = Delaunay(adjusted_points)
|
|||
|
|
|
|||
|
|
valid_simplices = []
|
|||
|
|
total_simplices = len(tri.simplices)
|
|||
|
|
|
|||
|
|
if use_quality_filter:
|
|||
|
|
print("正在进行三角形质量检查...")
|
|||
|
|
for simplex in tqdm(tri.simplices, desc="质量检查", unit="个"):
|
|||
|
|
p0 = adjusted_points[simplex[0]]
|
|||
|
|
p1 = adjusted_points[simplex[1]]
|
|||
|
|
p2 = adjusted_points[simplex[2]]
|
|||
|
|
p3 = adjusted_points[simplex[3]]
|
|||
|
|
|
|||
|
|
faces = [(p0, p1, p2), (p0, p1, p3), (p0, p2, p3), (p1, p2, p3)]
|
|||
|
|
valid_faces = 0
|
|||
|
|
|
|||
|
|
for face in faces:
|
|||
|
|
if self.is_valid_triangle(face[0], face[1], face[2], min_angle):
|
|||
|
|
valid_faces += 1
|
|||
|
|
|
|||
|
|
if valid_faces >= 3:
|
|||
|
|
valid_simplices.append(simplex)
|
|||
|
|
|
|||
|
|
print(f"质量过滤: {len(valid_simplices)}/{total_simplices} 个四面体通过质量检查")
|
|||
|
|
else:
|
|||
|
|
valid_simplices = tri.simplices
|
|||
|
|
|
|||
|
|
total_volume = 0.0
|
|||
|
|
valid_volume_count = 0
|
|||
|
|
|
|||
|
|
print(f"计算 {len(valid_simplices)} 个四面体的体积...")
|
|||
|
|
|
|||
|
|
for simplex in tqdm(valid_simplices, desc="计算四面体体积", unit="个"):
|
|||
|
|
p0 = adjusted_points[simplex[0]]
|
|||
|
|
p1 = adjusted_points[simplex[1]]
|
|||
|
|
p2 = adjusted_points[simplex[2]]
|
|||
|
|
p3 = adjusted_points[simplex[3]]
|
|||
|
|
|
|||
|
|
v1 = p1 - p0
|
|||
|
|
v2 = p2 - p0
|
|||
|
|
v3 = p3 - p0
|
|||
|
|
|
|||
|
|
det = np.linalg.det(np.array([v1, v2, v3]))
|
|||
|
|
volume = abs(det) / 6.0
|
|||
|
|
|
|||
|
|
if volume > 1e-12:
|
|||
|
|
total_volume += volume
|
|||
|
|
valid_volume_count += 1
|
|||
|
|
|
|||
|
|
print(f"有效体积计算: {valid_volume_count}/{len(valid_simplices)} 个四面体")
|
|||
|
|
|
|||
|
|
return total_volume
|
|||
|
|
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f"Delaunay三角剖分计算体积失败: {e}")
|
|||
|
|
return 0.0
|
|||
|
|
|
|||
|
|
|
|||
|
|
def load_and_filter_vertices(self, tileset_path):
|
|||
|
|
"""加载并过滤指定范围内的顶点"""
|
|||
|
|
print(f"开始处理tileset: {tileset_path}")
|
|||
|
|
|
|||
|
|
# 解析tileset获取所有B3DM文件
|
|||
|
|
b3dm_files = self.parse_tileset_json(tileset_path)
|
|||
|
|
print(f"找到 {len(b3dm_files)} 个B3DM文件")
|
|||
|
|
|
|||
|
|
if not b3dm_files:
|
|||
|
|
print("未找到任何B3DM文件")
|
|||
|
|
return False
|
|||
|
|
|
|||
|
|
# 处理每个B3DM文件
|
|||
|
|
processed_count = 0
|
|||
|
|
total_vertices = 0
|
|||
|
|
filtered_count = 0
|
|||
|
|
|
|||
|
|
for i, (b3dm_file, transform_matrix) in enumerate(tqdm(b3dm_files, desc="处理B3DM文件", unit="文件")):
|
|||
|
|
# print(f"处理文件 {i+1}/{len(b3dm_files)}: {os.path.basename(b3dm_file)}")
|
|||
|
|
|
|||
|
|
vertices = self.parse_b3dm_file(b3dm_file)
|
|||
|
|
if vertices:
|
|||
|
|
# 应用变换矩阵
|
|||
|
|
if transform_matrix:
|
|||
|
|
vertices = self.apply_transform_to_vertices(vertices, transform_matrix)
|
|||
|
|
|
|||
|
|
# 过滤范围内的顶点
|
|||
|
|
for vertex in vertices:
|
|||
|
|
total_vertices += 1
|
|||
|
|
if self.is_point_in_bounds(vertex):
|
|||
|
|
self.filtered_vertices.append(vertex)
|
|||
|
|
filtered_count += 1
|
|||
|
|
|
|||
|
|
self.all_vertices.extend(vertices)
|
|||
|
|
processed_count += 1
|
|||
|
|
# tqdm.write(f" {os.path.basename(b3dm_file)}: 提取到 {len(vertices)} 个顶点")
|
|||
|
|
|
|||
|
|
print(f"\n成功处理 {processed_count}/{len(b3dm_files)} 个文件")
|
|||
|
|
print(f"总计提取 {total_vertices} 个顶点")
|
|||
|
|
print(f"范围内顶点 {filtered_count} 个")
|
|||
|
|
|
|||
|
|
return len(self.filtered_vertices) > 0
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calculate_volume(self, tileset_path, base_elevation=None, min_angle=10.0, use_quality_filter=True):
|
|||
|
|
"""计算指定范围内三维模型的体积
|
|||
|
|
|
|||
|
|
Args:
|
|||
|
|
tileset_path: tileset.json文件路径
|
|||
|
|
base_elevation: 基准面高程
|
|||
|
|
min_angle: 最小角度约束（度）
|
|||
|
|
use_quality_filter: 是否启用质量过滤
|
|||
|
|
"""
|
|||
|
|
if not self.load_and_filter_vertices(tileset_path):
|
|||
|
|
print("未找到范围内的顶点数据")
|
|||
|
|
return 0.0
|
|||
|
|
|
|||
|
|
print(f"\n开始计算体积，使用Delaunay三角剖分方法")
|
|||
|
|
print(f"参与计算的顶点数: {len(self.filtered_vertices)}")
|
|||
|
|
|
|||
|
|
points = np.array(self.filtered_vertices)
|
|||
|
|
if base_elevation is None:
|
|||
|
|
base_elevation = np.min(points[:, 2])
|
|||
|
|
|
|||
|
|
print(f"统一基准面高程: {base_elevation:.2f} 米")
|
|||
|
|
|
|||
|
|
volume = self.calculate_volume_delaunay(self.filtered_vertices, base_elevation, min_angle, use_quality_filter)
|
|||
|
|
|
|||
|
|
print(f"\n计算结果:")
|
|||
|
|
print(f"体积: {volume:.6f} 立方米")
|
|||
|
|
print(f"体积: {volume/1000000:.6f} 立方千米")
|
|||
|
|
|
|||
|
|
return volume
|
|||
|
|
|
|||
|
|
def main():
|
|||
|
|
if len(sys.argv) < 3:
|
|||
|
|
print("用法: python volume_calculator.py <lboundary.json路径> <tileset.json路径> [基准面高程] [最小角度] [质量过滤]")
|
|||
|
|
print("基准面高程: 基准面高程(米)，默认使用最低点")
|
|||
|
|
print("最小角度: 三角形最小角度约束(度)，默认10.0")
|
|||
|
|
print("质量过滤: 是否启用质量过滤(true/false)，默认true")
|
|||
|
|
print("示例: python volume_calculator.py boundary.json tileset.json 100.0 15.0 true")
|
|||
|
|
return
|
|||
|
|
|
|||
|
|
location_file = sys.argv[1]
|
|||
|
|
tileset_path = sys.argv[2]
|
|||
|
|
|
|||
|
|
# 解析可选参数
|
|||
|
|
base_elevation = None
|
|||
|
|
min_angle = 10.0
|
|||
|
|
use_quality_filter = True
|
|||
|
|
|
|||
|
|
if len(sys.argv) > 3:
|
|||
|
|
try:
|
|||
|
|
base_elevation = float(sys.argv[3])
|
|||
|
|
except ValueError:
|
|||
|
|
print("错误：基准面高程必须是数字")
|
|||
|
|
return
|
|||
|
|
|
|||
|
|
if len(sys.argv) > 4:
|
|||
|
|
try:
|
|||
|
|
min_angle = float(sys.argv[4])
|
|||
|
|
except ValueError:
|
|||
|
|
print("错误：最小角度必须是数字")
|
|||
|
|
return
|
|||
|
|
|
|||
|
|
if len(sys.argv) > 5:
|
|||
|
|
use_quality_filter = sys.argv[5].lower() in ['true', '1', 'yes', 'on']
|
|||
|
|
|
|||
|
|
if not os.path.exists(location_file):
|
|||
|
|
print(f"错误: 位置文件不存在 {location_file}")
|
|||
|
|
return
|
|||
|
|
|
|||
|
|
if not os.path.exists(tileset_path):
|
|||
|
|
print(f"错误: tileset文件不存在 {tileset_path}")
|
|||
|
|
return
|
|||
|
|
|
|||
|
|
calculator = VolumeCalculator(location_file)
|
|||
|
|
volume = calculator.calculate_volume(tileset_path, base_elevation, min_angle, use_quality_filter)
|
|||
|
|
|
|||
|
|
if volume > 0:
|
|||
|
|
print("\n体积计算完成!")
|
|||
|
|
else:
|
|||
|
|
print("\n体积计算失败!")
|
|||
|
|
|
|||
|
|
if __name__ == "__main__":
|
|||
|
|
main()
|