11import math
22from typing import Optional
33
4+ import numpy as np
5+ from compas_cgal .meshing import pull_points_on_mesh
6+ from numpy import any
47from numpy import asarray
8+ from numpy import isnan
9+ from numpy import nan
510from scipy .interpolate import griddata
611
712from compas .datastructures import Mesh
13+ from compas .geometry import distance_point_point
814from compas_tna .diagrams import FormDiagram
915from compas_tna .envelope import Envelope
1016
1117
12- def griddata_project (xy : list [list [float ]], xyz_target : list [list [float ]], method = "linear" ) :
13- """Project a point cloud onto a target point cloud using griddata interpolation.
18+ def griddata_project (xy : list [list [float ]], xyz_target : list [list [float ]]) -> list [ float ] :
19+ """Project a set of XY points onto a target point cloud using griddata interpolation.
1420
1521 Parameters
1622 ----------
1723 xy : list[list[float]]
1824 The XY coordinates of the points to project.
1925 xyz_target : list[list[float]]
2026 The XYZ coordinates of the target points.
21- method : str, optional
22- The method to use for interpolation. Default is "linear".
2327
2428 Returns
2529 -------
2630 list[float]
27- The projected Z coordinates.
31+ The interpolated ( projected) Z coordinates.
2832 """
29- xy = asarray (xy )
30- xy_target = asarray (xyz_target )[:, :2 ]
31- z_target = asarray (xyz_target )[:, 2 ]
32- return griddata (xy_target , z_target , xy , method = method ).tolist ()
33+ xy = np .asarray (xy )
34+ xy_target = np .asarray (xyz_target )[:, :2 ]
35+ z_target = np .asarray (xyz_target )[:, 2 ]
36+
37+ # Initial interpolation
38+ z_proj = griddata (xy_target , z_target , xy , method = "linear" , fill_value = nan )
39+
40+ # Handle NaNs (outside convex hull)
41+ mask = isnan (z_proj )
42+ if any (mask ):
43+ z_proj [mask ] = griddata (xy_target , z_target , xy [mask ], method = "nearest" )
44+
45+ return z_proj .tolist ()
3346
3447
3548def interpolate_middle_mesh (intrados : Mesh , extrados : Mesh ) -> Mesh :
@@ -50,27 +63,24 @@ def interpolate_middle_mesh(intrados: Mesh, extrados: Mesh) -> Mesh:
5063 Mesh
5164 The interpolated middle mesh with proper normal-based thickness stored.
5265 """
53- # Use the intrados as base topology
54- middle = intrados .copy ()
55-
56- # Get Z coordinates from both surfaces based on the same XY coordinates
57- zi = middle .vertices_attribute ("z" )
58- ze = griddata_project (middle .vertices_attributes ("xy" ), extrados .vertices_attributes ("xyz" ), method = "linear" )
59-
60- # First loop: set middle Z as average of intrados and extrados
61- for i , key in enumerate (middle .vertices ()):
62- middle_z = (zi [i ] + ze [i ]) / 2.0
63- middle .vertex_attribute (key , "z" , middle_z )
64-
65- # Second loop: calculate and set thickness using correct normals
66- for i , key in enumerate (middle .vertices ()):
67- _ , _ , nz = middle .vertex_normal (key )
68- z_diff = ze [i ] - zi [i ]
69- if abs (nz ) > 0.1 :
70- thickness = abs (z_diff ) * abs (nz )
71- else :
72- thickness = abs (z_diff )
73- middle .vertex_attribute (key , "thickness" , thickness )
66+
67+ points , faces = intrados .to_vertices_and_faces ()
68+
69+ normals = [intrados .vertex_normal (key ) for key in intrados .vertices ()]
70+
71+ pulled_points = pull_points_on_mesh (points , normals , extrados )
72+
73+ midpoints = []
74+ thicknesses = []
75+
76+ for i , xyz in enumerate (pulled_points ):
77+ midpoints .append ([(xyz [0 ] + points [i ][0 ]) / 2 , (xyz [1 ] + points [i ][1 ]) / 2 , (xyz [2 ] + points [i ][2 ]) / 2 ])
78+ thicknesses .append (distance_point_point (xyz , points [i ]))
79+
80+ middle = Mesh .from_vertices_and_faces (midpoints , faces )
81+
82+ for i , vertex in enumerate (middle .vertices ()):
83+ middle .vertex_attribute (vertex , "thickness" , thicknesses [i ])
7484
7585 return middle
7686
@@ -134,46 +144,6 @@ def offset_from_middle(middle: Mesh, fixed_xy: bool = True) -> tuple[Mesh, Mesh]
134144 return intrados , extrados
135145
136146
137- def project_mesh_to_target_vertica_nearest (mesh : Mesh , target : Mesh ) -> None :
138- """Project a mesh vertically (in Z direction) onto a target mesh.
139-
140- Parameters
141- ----------
142- mesh : Mesh
143- The mesh to be projected.
144- target : Mesh
145- The target mesh to project onto.
146-
147- Returns
148- -------
149- None
150- The mesh is modified in place.
151- """
152- # Get target mesh vertices for simple vertical projection
153- target_vertices = list (target .vertices ())
154- target_points = [target .vertex_point (v ) for v in target_vertices ]
155-
156- for vertex in mesh .vertices ():
157- point = mesh .vertex_point (vertex )
158-
159- # Find the closest target vertex in XY plane
160- min_distance = float ("inf" )
161- closest_z = point .z
162-
163- for target_point in target_points :
164- # Calculate XY distance (ignore Z)
165- xy_distance = ((point .x - target_point .x ) ** 2 + (point .y - target_point .y ) ** 2 ) ** 0.5
166-
167- if xy_distance < min_distance :
168- min_distance = xy_distance
169- closest_z = target_point .z
170-
171- # Update vertex to closest Z value
172- new_point = point .copy ()
173- new_point .z = closest_z
174- mesh .vertex_attributes (vertex , "xyz" , new_point )
175-
176-
177147def project_mesh_to_target_vertical (mesh : Mesh , target : Mesh ) -> None :
178148 """Project a mesh vertically (in Z direction) onto a target mesh.
179149
@@ -189,7 +159,7 @@ def project_mesh_to_target_vertical(mesh: Mesh, target: Mesh) -> None:
189159 None
190160 The mesh is modified in place.
191161 """
192- z_target = griddata_project (mesh .vertices_attributes ("xy" ), target .vertices_attributes ("xyz" ), method = "linear" )
162+ z_target = griddata_project (mesh .vertices_attributes ("xy" ), target .vertices_attributes ("xyz" ))
193163
194164 for key , i in enumerate (mesh .vertices ()):
195165 mesh .vertex_attribute (key , "z" , z_target [i ])
@@ -586,7 +556,10 @@ def apply_fill_weight_to_formdiagram(self, formdiagram: FormDiagram) -> None:
586556 else :
587557 print (f"Warning: Vertex { vertex } not found in projected meshes" )
588558
589- print (f"Fill weight applied to form diagram. Total load: { fill_weight } " )
559+ sum_pz = sum (abs (asarray (formdiagram .vertices_attribute ("pz" ))))
560+
561+ print (f"Fill weight applied to form diagram. Fill load: { fill_weight :.1f} " )
562+ print (f"Total Load: { sum_pz :.1f} " )
590563
591564 # =============================================================================
592565 # Envelope and target projection operations
0 commit comments