22
22
"outersurf" ,
23
23
"surfvolume" ,
24
24
"insurface" ,
25
+ "remeshsurf" ,
25
26
]
26
27
27
28
##====================================================================================
31
32
import numpy as np
32
33
import os
33
34
import re
34
- import platform
35
+ import sys
35
36
import subprocess
36
37
from iso2mesh .trait import (
37
38
surfinterior ,
40
41
elemvolume ,
41
42
meshreorient ,
42
43
finddisconnsurf ,
44
+ maxsurf ,
45
+ volface ,
43
46
)
44
47
from iso2mesh .utils import *
45
48
from iso2mesh .io import saveoff , readoff , saveinr , readtetgen , savesurfpoly , readmedit
49
52
meshresample ,
50
53
removeisolatedsurf ,
51
54
removeisolatednode ,
55
+ qmeshcut ,
52
56
)
53
57
54
58
##====================================================================================
@@ -572,7 +576,9 @@ def surf2volz(node, face, xi, yi, zi):
572
576
573
577
for i in iz :
574
578
plane = np .array ([[0 , 100 , zi [i ]], [100 , 0 , zi [i ]], [0 , 0 , zi [i ]]])
575
- bcutpos , bcutvalue , bcutedges = qmeshcut (face [:, :3 ], node , node [:, 0 ], plane )
579
+ bcutpos , bcutvalue , bcutedges , _ , _ = qmeshcut (
580
+ face [:, :3 ], node , node [:, 0 ], plane
581
+ )
576
582
577
583
if bcutpos .size == 0 :
578
584
continue
@@ -622,6 +628,7 @@ def surf2vol(node, face, xi, yi, zi, **kwargs):
622
628
img: a volumetric binary image at the position of ndgrid(xi, yi, zi)
623
629
v2smap (optional): a 4x4 matrix denoting the Affine transformation to map voxel coordinates back to the mesh space.
624
630
"""
631
+ from scipy .ndimage import binary_fill_holes
625
632
626
633
opt = kwargs
627
634
label = opt .get ("label" , 0 )
@@ -651,7 +658,7 @@ def surf2vol(node, face, xi, yi, zi, **kwargs):
651
658
im |= np .moveaxis (surf2volz (node [:, [1 , 2 , 0 ]], fc [:, :3 ], yi , zi , xi ), 0 , 1 )
652
659
653
660
if opt .get ("fill" , 0 ) or label :
654
- im = imfill (im , "holes" )
661
+ im = binary_fill_holes (im )
655
662
if label :
656
663
im = im .astype (elabel .dtype ) * lbl
657
664
@@ -851,7 +858,7 @@ def cgalv2m(vol, opt, maxvol):
851
858
return node , elem , face
852
859
853
860
854
- def cgals2m (v , f , opt , maxvol , * args ):
861
+ def cgals2m (v , f , opt , maxvol , ** kwargs ):
855
862
"""
856
863
Convert a triangular surface to a tetrahedral mesh using CGAL mesher.
857
864
@@ -888,7 +895,6 @@ def cgals2m(v, f, opt, maxvol, *args):
888
895
ssize = 6
889
896
approx = 0.5
890
897
reratio = 3
891
- flags = args_to_dict (* args )
892
898
893
899
if not isinstance (opt , dict ):
894
900
ssize = opt
@@ -899,7 +905,7 @@ def cgals2m(v, f, opt, maxvol, *args):
899
905
approx = opt .get ("distbound" , approx )
900
906
reratio = opt .get ("reratio" , reratio )
901
907
902
- if flags .get ("DoRepair " , 0 ) == 1 :
908
+ if kwargs .get ("dorepair " , 0 ) == 1 :
903
909
v , f = meshcheckrepair (v , f )
904
910
905
911
saveoff (v , f , mwpath ("pre_cgalpoly.off" ))
@@ -1149,3 +1155,68 @@ def insurface(node, face, points):
1149
1155
tf = tf .astype (int )
1150
1156
1151
1157
return tf
1158
+
1159
+
1160
+ def remeshsurf (node , face , opt ):
1161
+ """
1162
+ remeshsurf(node, face, opt)
1163
+
1164
+ Remesh a triangular surface, output is guaranteed to be free of self-intersecting elements.
1165
+ This function can both downsample or upsample a mesh.
1166
+
1167
+ Parameters:
1168
+ node: list of nodes on the input surface mesh, 3 columns for x, y, z
1169
+ face: list of triangular elements on the surface, [n1, n2, n3, region_id]
1170
+ opt: function parameters
1171
+ opt.gridsize: resolution for the voxelization of the mesh
1172
+ opt.closesize: if there are openings, set the closing diameter
1173
+ opt.elemsize: the size of the element of the output surface
1174
+ If opt is a scalar, it defines the elemsize and gridsize = opt / 4
1175
+
1176
+ Returns:
1177
+ newno: list of nodes on the resulting surface mesh, 3 columns for x, y, z
1178
+ newfc: list of triangular elements on the surface, [n1, n2, n3, region_id]
1179
+ """
1180
+ from scipy .ndimage import binary_fill_holes
1181
+
1182
+ # Step 1: convert the old surface to a volumetric image
1183
+ p0 = np .min (node , axis = 0 )
1184
+ p1 = np .max (node , axis = 0 )
1185
+
1186
+ if isinstance (opt , dict ):
1187
+ dx = opt .get ("gridsize" , None )
1188
+ else :
1189
+ dx = opt / 4
1190
+
1191
+ x_range = np .arange (p0 [0 ] - dx , p1 [0 ] + dx , dx )
1192
+ y_range = np .arange (p0 [1 ] - dx , p1 [1 ] + dx )
1193
+ z_range = np .arange (p0 [2 ] - dx , p1 [2 ] + dx )
1194
+
1195
+ img = surf2vol (node , face , x_range , y_range , z_range )
1196
+
1197
+ # Compute surface edges
1198
+ eg = surfedge (face )
1199
+
1200
+ closesize = 0
1201
+ if eg .size > 0 and isinstance (opt , dict ):
1202
+ closesize = opt .get ("closesize" , 0 )
1203
+
1204
+ # Step 2: fill holes in the volumetric binary image
1205
+ img = binary_fill_holes (img )
1206
+
1207
+ # Step 3: convert the filled volume to a new surface
1208
+ if isinstance (opt , dict ):
1209
+ if "elemsize" in opt :
1210
+ opt ["radbound" ] = opt ["elemsize" ] / dx
1211
+ newno , newfc , _ , _ = v2s (img , 0.5 , opt , "cgalsurf" )
1212
+ else :
1213
+ opt = {"radbound" : opt / dx }
1214
+ newno , newfc , _ , _ = v2s (img , 0.5 , opt , "cgalsurf" )
1215
+
1216
+ # Adjust new nodes to match original coordinates
1217
+ newno [:, 0 :3 ] *= dx
1218
+ newno [:, 0 ] += p0 [0 ]
1219
+ newno [:, 1 ] += p0 [1 ]
1220
+ newno [:, 2 ] += p0 [2 ]
1221
+
1222
+ return newno , newfc
0 commit comments