2222 "outersurf" ,
2323 "surfvolume" ,
2424 "insurface" ,
25+ "remeshsurf" ,
2526]
2627
2728##====================================================================================
3132import numpy as np
3233import os
3334import re
34- import platform
35+ import sys
3536import subprocess
3637from iso2mesh .trait import (
3738 surfinterior ,
4041 elemvolume ,
4142 meshreorient ,
4243 finddisconnsurf ,
44+ maxsurf ,
45+ volface ,
4346)
4447from iso2mesh .utils import *
4548from iso2mesh .io import saveoff , readoff , saveinr , readtetgen , savesurfpoly , readmedit
4952 meshresample ,
5053 removeisolatedsurf ,
5154 removeisolatednode ,
55+ qmeshcut ,
5256)
5357
5458##====================================================================================
@@ -572,7 +576,9 @@ def surf2volz(node, face, xi, yi, zi):
572576
573577 for i in iz :
574578 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+ )
576582
577583 if bcutpos .size == 0 :
578584 continue
@@ -622,6 +628,7 @@ def surf2vol(node, face, xi, yi, zi, **kwargs):
622628 img: a volumetric binary image at the position of ndgrid(xi, yi, zi)
623629 v2smap (optional): a 4x4 matrix denoting the Affine transformation to map voxel coordinates back to the mesh space.
624630 """
631+ from scipy .ndimage import binary_fill_holes
625632
626633 opt = kwargs
627634 label = opt .get ("label" , 0 )
@@ -651,7 +658,7 @@ def surf2vol(node, face, xi, yi, zi, **kwargs):
651658 im |= np .moveaxis (surf2volz (node [:, [1 , 2 , 0 ]], fc [:, :3 ], yi , zi , xi ), 0 , 1 )
652659
653660 if opt .get ("fill" , 0 ) or label :
654- im = imfill (im , "holes" )
661+ im = binary_fill_holes (im )
655662 if label :
656663 im = im .astype (elabel .dtype ) * lbl
657664
@@ -851,7 +858,7 @@ def cgalv2m(vol, opt, maxvol):
851858 return node , elem , face
852859
853860
854- def cgals2m (v , f , opt , maxvol , * args ):
861+ def cgals2m (v , f , opt , maxvol , ** kwargs ):
855862 """
856863 Convert a triangular surface to a tetrahedral mesh using CGAL mesher.
857864
@@ -888,7 +895,6 @@ def cgals2m(v, f, opt, maxvol, *args):
888895 ssize = 6
889896 approx = 0.5
890897 reratio = 3
891- flags = args_to_dict (* args )
892898
893899 if not isinstance (opt , dict ):
894900 ssize = opt
@@ -899,7 +905,7 @@ def cgals2m(v, f, opt, maxvol, *args):
899905 approx = opt .get ("distbound" , approx )
900906 reratio = opt .get ("reratio" , reratio )
901907
902- if flags .get ("DoRepair " , 0 ) == 1 :
908+ if kwargs .get ("dorepair " , 0 ) == 1 :
903909 v , f = meshcheckrepair (v , f )
904910
905911 saveoff (v , f , mwpath ("pre_cgalpoly.off" ))
@@ -1149,3 +1155,68 @@ def insurface(node, face, points):
11491155 tf = tf .astype (int )
11501156
11511157 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