370 lines
13 KiB
Python
370 lines
13 KiB
Python
"""
|
|
Georeferenzierungs-Modul
|
|
Transformation mit Passpunkten (mind. 3 Punkte)
|
|
NUR Rotation und Translation - KEINE Maßstabsänderung
|
|
"""
|
|
|
|
import math
|
|
import numpy as np
|
|
from typing import List, Tuple, Optional, Dict
|
|
from dataclasses import dataclass
|
|
from .cor_generator import CORPoint
|
|
|
|
|
|
@dataclass
|
|
class ControlPoint:
|
|
"""Passpunkt mit Soll- und Ist-Koordinaten"""
|
|
name: str
|
|
# Ist-Koordinaten (lokales System)
|
|
local_x: float
|
|
local_y: float
|
|
local_z: float
|
|
# Soll-Koordinaten (Zielsystem, z.B. UTM)
|
|
target_x: float
|
|
target_y: float
|
|
target_z: float
|
|
# Residuen nach Transformation
|
|
residual_x: float = 0.0
|
|
residual_y: float = 0.0
|
|
residual_z: float = 0.0
|
|
|
|
@property
|
|
def residual_2d(self) -> float:
|
|
"""2D-Residuum"""
|
|
return math.sqrt(self.residual_x**2 + self.residual_y**2)
|
|
|
|
@property
|
|
def residual_3d(self) -> float:
|
|
"""3D-Residuum"""
|
|
return math.sqrt(self.residual_x**2 + self.residual_y**2 + self.residual_z**2)
|
|
|
|
|
|
@dataclass
|
|
class GeoreferenceResult:
|
|
"""Ergebnis der Georeferenzierung"""
|
|
# Transformationsparameter
|
|
translation_x: float = 0.0
|
|
translation_y: float = 0.0
|
|
translation_z: float = 0.0
|
|
rotation_rad: float = 0.0
|
|
|
|
# Qualitätsparameter
|
|
rmse_x: float = 0.0
|
|
rmse_y: float = 0.0
|
|
rmse_z: float = 0.0
|
|
rmse_2d: float = 0.0
|
|
rmse_3d: float = 0.0
|
|
max_residual_2d: float = 0.0
|
|
max_residual_3d: float = 0.0
|
|
|
|
# Passpunkte mit Residuen
|
|
control_points: List[ControlPoint] = None
|
|
|
|
def __post_init__(self):
|
|
if self.control_points is None:
|
|
self.control_points = []
|
|
|
|
@property
|
|
def rotation_gon(self) -> float:
|
|
"""Rotation in Gon"""
|
|
return self.rotation_rad * 200.0 / math.pi
|
|
|
|
@property
|
|
def rotation_deg(self) -> float:
|
|
"""Rotation in Grad"""
|
|
return math.degrees(self.rotation_rad)
|
|
|
|
|
|
class Georeferencer:
|
|
"""
|
|
Georeferenzierung mit mindestens 3 Passpunkten
|
|
Verwendet Methode der kleinsten Quadrate für Überbestimmung
|
|
NUR Rotation und Translation (4 Parameter)
|
|
"""
|
|
|
|
def __init__(self):
|
|
self.control_points: List[ControlPoint] = []
|
|
self.result: Optional[GeoreferenceResult] = None
|
|
self.points_to_transform: List[CORPoint] = []
|
|
self.transformed_points: List[CORPoint] = []
|
|
|
|
def add_control_point(self, name: str,
|
|
local_x: float, local_y: float, local_z: float,
|
|
target_x: float, target_y: float, target_z: float):
|
|
"""Fügt einen Passpunkt hinzu"""
|
|
self.control_points.append(ControlPoint(
|
|
name=name,
|
|
local_x=local_x, local_y=local_y, local_z=local_z,
|
|
target_x=target_x, target_y=target_y, target_z=target_z
|
|
))
|
|
|
|
def clear_control_points(self):
|
|
"""Löscht alle Passpunkte"""
|
|
self.control_points = []
|
|
self.result = None
|
|
|
|
def set_points_to_transform(self, points: List[CORPoint]):
|
|
"""Setzt die zu transformierenden Punkte"""
|
|
self.points_to_transform = points.copy()
|
|
|
|
def compute_transformation(self) -> GeoreferenceResult:
|
|
"""
|
|
Berechnet die Transformationsparameter
|
|
4-Parameter-Transformation: Rotation + Translation (kein Maßstab)
|
|
Methode der kleinsten Quadrate
|
|
"""
|
|
if len(self.control_points) < 3:
|
|
raise ValueError("Mindestens 3 Passpunkte erforderlich!")
|
|
|
|
n = len(self.control_points)
|
|
|
|
# Koordinaten extrahieren
|
|
local_coords = np.array([[cp.local_x, cp.local_y] for cp in self.control_points])
|
|
target_coords = np.array([[cp.target_x, cp.target_y] for cp in self.control_points])
|
|
|
|
# Schwerpunkte berechnen
|
|
local_centroid = np.mean(local_coords, axis=0)
|
|
target_centroid = np.mean(target_coords, axis=0)
|
|
|
|
# Zum Schwerpunkt verschieben
|
|
local_centered = local_coords - local_centroid
|
|
target_centered = target_coords - target_centroid
|
|
|
|
# Rotation berechnen mit SVD (Procrustes ohne Skalierung)
|
|
H = local_centered.T @ target_centered
|
|
U, S, Vt = np.linalg.svd(H)
|
|
|
|
# Rotationsmatrix
|
|
R = Vt.T @ U.T
|
|
|
|
# Reflexion korrigieren falls nötig
|
|
if np.linalg.det(R) < 0:
|
|
Vt[-1, :] *= -1
|
|
R = Vt.T @ U.T
|
|
|
|
# Rotationswinkel aus Matrix
|
|
rotation_rad = math.atan2(R[1, 0], R[0, 0])
|
|
|
|
# Translation berechnen
|
|
translation = target_centroid - R @ local_centroid
|
|
|
|
# Ergebnis speichern
|
|
self.result = GeoreferenceResult(
|
|
translation_x=translation[0],
|
|
translation_y=translation[1],
|
|
rotation_rad=rotation_rad
|
|
)
|
|
|
|
# Z-Translation (Mittelwert der Höhendifferenzen)
|
|
z_diffs = [cp.target_z - cp.local_z for cp in self.control_points]
|
|
self.result.translation_z = np.mean(z_diffs)
|
|
|
|
# Residuen berechnen
|
|
self._compute_residuals()
|
|
|
|
return self.result
|
|
|
|
def _compute_residuals(self):
|
|
"""Berechnet Residuen für alle Passpunkte"""
|
|
if self.result is None:
|
|
return
|
|
|
|
cos_r = math.cos(self.result.rotation_rad)
|
|
sin_r = math.sin(self.result.rotation_rad)
|
|
|
|
sum_res_x2 = 0.0
|
|
sum_res_y2 = 0.0
|
|
sum_res_z2 = 0.0
|
|
max_res_2d = 0.0
|
|
max_res_3d = 0.0
|
|
|
|
for cp in self.control_points:
|
|
# Transformation anwenden
|
|
x_trans = cos_r * cp.local_x - sin_r * cp.local_y + self.result.translation_x
|
|
y_trans = sin_r * cp.local_x + cos_r * cp.local_y + self.result.translation_y
|
|
z_trans = cp.local_z + self.result.translation_z
|
|
|
|
# Residuen
|
|
cp.residual_x = cp.target_x - x_trans
|
|
cp.residual_y = cp.target_y - y_trans
|
|
cp.residual_z = cp.target_z - z_trans
|
|
|
|
sum_res_x2 += cp.residual_x**2
|
|
sum_res_y2 += cp.residual_y**2
|
|
sum_res_z2 += cp.residual_z**2
|
|
|
|
max_res_2d = max(max_res_2d, cp.residual_2d)
|
|
max_res_3d = max(max_res_3d, cp.residual_3d)
|
|
|
|
n = len(self.control_points)
|
|
self.result.rmse_x = math.sqrt(sum_res_x2 / n)
|
|
self.result.rmse_y = math.sqrt(sum_res_y2 / n)
|
|
self.result.rmse_z = math.sqrt(sum_res_z2 / n)
|
|
self.result.rmse_2d = math.sqrt((sum_res_x2 + sum_res_y2) / n)
|
|
self.result.rmse_3d = math.sqrt((sum_res_x2 + sum_res_y2 + sum_res_z2) / n)
|
|
self.result.max_residual_2d = max_res_2d
|
|
self.result.max_residual_3d = max_res_3d
|
|
self.result.control_points = self.control_points
|
|
|
|
def transform_points(self) -> List[CORPoint]:
|
|
"""Transformiert alle Punkte"""
|
|
if self.result is None:
|
|
raise ValueError("Transformation noch nicht berechnet!")
|
|
|
|
self.transformed_points = []
|
|
|
|
cos_r = math.cos(self.result.rotation_rad)
|
|
sin_r = math.sin(self.result.rotation_rad)
|
|
|
|
for p in self.points_to_transform:
|
|
x_trans = cos_r * p.x - sin_r * p.y + self.result.translation_x
|
|
y_trans = sin_r * p.x + cos_r * p.y + self.result.translation_y
|
|
z_trans = p.z + self.result.translation_z
|
|
|
|
self.transformed_points.append(CORPoint(
|
|
name=p.name,
|
|
x=x_trans,
|
|
y=y_trans,
|
|
z=z_trans
|
|
))
|
|
|
|
return self.transformed_points
|
|
|
|
def transform_single_point(self, x: float, y: float, z: float) -> Tuple[float, float, float]:
|
|
"""Transformiert einen einzelnen Punkt"""
|
|
if self.result is None:
|
|
raise ValueError("Transformation noch nicht berechnet!")
|
|
|
|
cos_r = math.cos(self.result.rotation_rad)
|
|
sin_r = math.sin(self.result.rotation_rad)
|
|
|
|
x_trans = cos_r * x - sin_r * y + self.result.translation_x
|
|
y_trans = sin_r * x + cos_r * y + self.result.translation_y
|
|
z_trans = z + self.result.translation_z
|
|
|
|
return x_trans, y_trans, z_trans
|
|
|
|
def get_transformation_report(self) -> str:
|
|
"""Erstellt einen detaillierten Transformationsbericht"""
|
|
if self.result is None:
|
|
return "Keine Transformation berechnet."
|
|
|
|
lines = []
|
|
lines.append("=" * 70)
|
|
lines.append("GEOREFERENZIERUNG - TRANSFORMATIONSBERICHT")
|
|
lines.append("=" * 70)
|
|
lines.append("")
|
|
lines.append("TRANSFORMATIONSPARAMETER (4-Parameter, ohne Maßstab)")
|
|
lines.append("-" * 70)
|
|
lines.append(f" Translation X: {self.result.translation_x:>15.4f} m")
|
|
lines.append(f" Translation Y: {self.result.translation_y:>15.4f} m")
|
|
lines.append(f" Translation Z: {self.result.translation_z:>15.4f} m")
|
|
lines.append(f" Rotation: {self.result.rotation_gon:>15.6f} gon")
|
|
lines.append(f" {self.result.rotation_deg:>15.6f}°")
|
|
lines.append(f" {self.result.rotation_rad:>15.8f} rad")
|
|
lines.append(f" Maßstab: {1.0:>15.6f} (fest)")
|
|
lines.append("")
|
|
lines.append("QUALITÄTSPARAMETER")
|
|
lines.append("-" * 70)
|
|
lines.append(f" RMSE X: {self.result.rmse_x*1000:>15.2f} mm")
|
|
lines.append(f" RMSE Y: {self.result.rmse_y*1000:>15.2f} mm")
|
|
lines.append(f" RMSE Z: {self.result.rmse_z*1000:>15.2f} mm")
|
|
lines.append(f" RMSE 2D: {self.result.rmse_2d*1000:>15.2f} mm")
|
|
lines.append(f" RMSE 3D: {self.result.rmse_3d*1000:>15.2f} mm")
|
|
lines.append(f" Max. Residuum 2D: {self.result.max_residual_2d*1000:>15.2f} mm")
|
|
lines.append(f" Max. Residuum 3D: {self.result.max_residual_3d*1000:>15.2f} mm")
|
|
lines.append("")
|
|
lines.append("PASSPUNKTE UND RESIDUEN")
|
|
lines.append("-" * 70)
|
|
lines.append(f"{'Punkt':<10} {'vX [mm]':>10} {'vY [mm]':>10} {'vZ [mm]':>10} {'v2D [mm]':>10} {'v3D [mm]':>10}")
|
|
lines.append("-" * 70)
|
|
|
|
for cp in self.result.control_points:
|
|
lines.append(f"{cp.name:<10} {cp.residual_x*1000:>10.2f} {cp.residual_y*1000:>10.2f} "
|
|
f"{cp.residual_z*1000:>10.2f} {cp.residual_2d*1000:>10.2f} {cp.residual_3d*1000:>10.2f}")
|
|
|
|
lines.append("-" * 70)
|
|
lines.append(f"Anzahl Passpunkte: {len(self.control_points)}")
|
|
lines.append(f"Redundanz (2D): {2 * len(self.control_points) - 4}")
|
|
lines.append("")
|
|
lines.append("=" * 70)
|
|
|
|
return "\n".join(lines)
|
|
|
|
def get_control_points_comparison(self) -> str:
|
|
"""Erstellt eine Vergleichstabelle für Passpunkte"""
|
|
if self.result is None:
|
|
return "Keine Transformation berechnet."
|
|
|
|
lines = []
|
|
lines.append("PASSPUNKT-KOORDINATEN: LOKAL → ZIEL")
|
|
lines.append("=" * 100)
|
|
lines.append(f"{'Punkt':<8} {'X_lokal':>12} {'Y_lokal':>12} {'Z_lokal':>10} | "
|
|
f"{'X_Ziel':>12} {'Y_Ziel':>12} {'Z_Ziel':>10}")
|
|
lines.append("-" * 100)
|
|
|
|
for cp in self.control_points:
|
|
lines.append(f"{cp.name:<8} {cp.local_x:>12.3f} {cp.local_y:>12.3f} {cp.local_z:>10.3f} | "
|
|
f"{cp.target_x:>12.3f} {cp.target_y:>12.3f} {cp.target_z:>10.3f}")
|
|
|
|
lines.append("=" * 100)
|
|
return "\n".join(lines)
|
|
|
|
def export_report(self, filepath: str):
|
|
"""Exportiert den vollständigen Bericht"""
|
|
report = self.get_transformation_report()
|
|
with open(filepath, 'w', encoding='utf-8') as f:
|
|
f.write(report)
|
|
|
|
|
|
class RigidBodyTransformation:
|
|
"""
|
|
Alternative Implementierung: Rigide Körpertransformation
|
|
6 Parameter: 3 Translationen + 3 Rotationen
|
|
Für 3D-Daten mit mehreren Rotationsachsen
|
|
"""
|
|
|
|
def __init__(self):
|
|
self.control_points: List[ControlPoint] = []
|
|
self.translation = np.zeros(3)
|
|
self.rotation_matrix = np.eye(3)
|
|
|
|
def compute(self, local_points: np.ndarray, target_points: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
|
|
"""
|
|
Berechnet rigide Transformation mit Kabsch-Algorithmus
|
|
Keine Maßstabsänderung
|
|
"""
|
|
# Schwerpunkte
|
|
centroid_local = np.mean(local_points, axis=0)
|
|
centroid_target = np.mean(target_points, axis=0)
|
|
|
|
# Zentrieren
|
|
local_centered = local_points - centroid_local
|
|
target_centered = target_points - centroid_target
|
|
|
|
# Kovarianzmatrix
|
|
H = local_centered.T @ target_centered
|
|
|
|
# SVD
|
|
U, S, Vt = np.linalg.svd(H)
|
|
|
|
# Rotation
|
|
R = Vt.T @ U.T
|
|
|
|
# Reflexionskorrektur
|
|
if np.linalg.det(R) < 0:
|
|
Vt[-1, :] *= -1
|
|
R = Vt.T @ U.T
|
|
|
|
# Translation
|
|
t = centroid_target - R @ centroid_local
|
|
|
|
self.rotation_matrix = R
|
|
self.translation = t
|
|
|
|
return R, t
|
|
|
|
def transform(self, points: np.ndarray) -> np.ndarray:
|
|
"""Transformiert Punkte"""
|
|
return (self.rotation_matrix @ points.T).T + self.translation
|