diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index d3f98d5..ae868bc 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -5,6 +5,9 @@ on: branches: - main - fix-ci + pull_request: + branches: + - main env: diff --git a/python/cadconvert/bezier.py b/python/cadconvert/bezier.py new file mode 100644 index 0000000..eec446a --- /dev/null +++ b/python/cadconvert/bezier.py @@ -0,0 +1,39 @@ +import fem_tabulator as feta +import numpy as np +import quad_curve_utils as qr +import tqdm + +def eval_bc(verts, faces, bc_i, denom:int): + vf = (verts[faces[bc_i[:,0]]]) + return np.einsum('sed,se->sd', vf, bc_i[:,1:])/denom + +def bezier_fit_matrix(order : int, level : int) -> np.ndarray: + std_x, std_y = qr.quad_tuple_gen(level).T + bsv = qr.tp_sample_value(std_x/level, std_y/level, order=order) + bsv = bsv.reshape(len(bsv),-1) + return bsv + +def bezier_check_validity(mB,mT,F, quads, q2t, trim_types, quad_cp, order, valid_check, progress_bar = True): + all_b, all_t = qr.top_and_bottom_sample(mB,mT, F, quads, q2t, trim_types, level=1) + v4, f4 = qr.split_square(1) + + tup = feta.tuple_gen(order = order + 1, var_n=2) # elevated one order for quartic tetrahedra + grids = np.einsum('fed,Ee->fEd', v4[f4], np.asarray(tup)) + + grid_ids = np.ravel_multi_index(grids.reshape(-1,2).T, dims = (order + 2, + order + 2)).reshape(len(f4), -1) + valid_quad = np.ones(len(quad_cp), dtype=bool) + + A13 = bezier_fit_matrix(order, order+1) + if progress_bar: + pbar = tqdm.tqdm(quad_cp, desc='Bezier Quads Checking validity') + else: + pbar = quad_cp + for q,qcp in enumerate(pbar): + lagr = A13@qcp + for t,g in zip(f4, grid_ids): + if not (valid_check(all_b[q][t], all_t[q][t], + lagr[g], True)): + valid_quad[q] = False + break + return valid_quad \ No newline at end of file diff --git a/python/cadconvert/fem_tabulator.py b/python/cadconvert/fem_tabulator.py new file mode 100755 index 0000000..9b60c46 --- /dev/null +++ b/python/cadconvert/fem_tabulator.py @@ -0,0 +1,27 @@ +import sympy +import math +import numpy as np + + +def tuple_gen(order, var_n): + if var_n == 0: + return [[order]] + l = [] + for i in range(order + 1): + r = tuple_gen(order-i, var_n - 1) + l += [[i]+t for t in r] + l = sorted(l, key=lambda x: (-sum(i**2 for i in x), x[::-1])) + return l + +def bernstein_evaluator(x, y, z, codecs): + m = len(codecs[0]) # dim + 1 + n = codecs[0][0] # order + mc_dict = sympy.multinomial_coefficients(m, n) + mc = np.array([mc_dict[tuple(c)] for c in codecs]) + + w = 1-x-y-z + computed_powers = np.array([(w**i, x**i, y**i, z**i) + for i in range(n + 1)]) # make use of 0**0 == 1 + return mc[:,None]*np.array( + [np.prod([computed_powers[c, i] for i, c in enumerate(cod)], axis=0) for cod in codecs]) + diff --git a/python/cadconvert/occ_step.py b/python/cadconvert/occ_step.py new file mode 100644 index 0000000..c510252 --- /dev/null +++ b/python/cadconvert/occ_step.py @@ -0,0 +1,99 @@ +from OCC.Core.TColgp import TColgp_HArray1OfPnt2d, TColgp_Array1OfPnt2d,TColgp_Array2OfPnt +from OCC.Core.gp import gp_Pnt2d,gp_Vec, gp_Pnt +from OCC.Core.Geom import Geom_BezierSurface, Geom_BSplineSurface +from OCC.Core.TColGeom import TColGeom_Array2OfBezierSurface +from OCC.Core.GeomConvert import GeomConvert_CompBezierSurfacesToBSplineSurface +from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeFace +from OCC.Display.SimpleGui import init_display +from OCC.Core.STEPControl import STEPControl_Writer, STEPControl_AsIs +from OCC.Core.Interface import Interface_Static_SetCVal +import numpy as np + +def example_vis(surf): + display, start_display, add_menu, add_function_to_menu = init_display() + display.EraseAll() + display.DisplayShape(surf, update=True) + start_display() + +def cp_to_bz(cp): + """Tensor Product (Bezier) Control Points to OCC-BezierSurface object. + + Args: + cp (np.array): input control point, 4x4x3 + + Returns: + OCC-Geom_BezierSurface: OCC internal type for a Bezier Patch + """ + array1 = TColgp_Array2OfPnt(1, len(cp), 1, len(cp)) + for i in range(1,len(cp)+1): + for j in range(1,len(cp)+1): + array1.SetValue(i,j, gp_Pnt(*cp[i-1,j-1])) + BZ1 = Geom_BezierSurface(array1) + return BZ1 + +def cp_write_to_step(output_file, quad_cp): + step_writer = STEPControl_Writer() + Interface_Static_SetCVal("write.step.schema", "AP203") + + for cp in quad_cp: + assert len(cp) == 16 + b1 = cp_to_bz(cp.reshape(4,4,3)) + build = BRepBuilderAPI_MakeFace(b1, 1e-6) + step_writer.Transfer(build.Shape(),STEPControl_AsIs) + status = step_writer.Write(output_file) + + +def compose_bezier(bz_list): + bezierarray = TColGeom_Array2OfBezierSurface(1, len(bz_list), 1,1) + for i,b in enumerate(bz_list): + bezierarray.SetValue(i+1, 1, b) + BB = GeomConvert_CompBezierSurfacesToBSplineSurface(bezierarray) + if BB.IsDone(): + poles = BB.Poles().Array2() + uknots = BB.UKnots().Array1() + vknots = BB.VKnots().Array1() + umult = BB.UMultiplicities().Array1() + vmult = BB.VMultiplicities().Array1() + udeg = BB.UDegree() + vdeg = BB.VDegree() + BSPLSURF = Geom_BSplineSurface( poles, uknots, vknots, umult, vmult, udeg, vdeg, False, False) + + return BSPLSURF + else: + return None + +def test_compose(): + cp = np.zeros((4,4,3)) + cp2 = np.zeros((4,4,3)) + for i in range(4): + for j in range(4): + cp[i,j] = (i-3, j,(i-3)**2) + cp2[i,j] = (i,j,i**2) + b1 = cp_to_bz(cp) + b2 = cp_to_bz(cp2) + bsp = compose_bezier([b1,b2]) + + step_writer = STEPControl_Writer() + Interface_Static_SetCVal("write.step.schema", "AP203") + + build = BRepBuilderAPI_MakeFace(bsp, 1e-6) + step_writer.Transfer(build.Shape(),STEPControl_AsIs) + status = step_writer.Write('test.stp') + +def stripe_writer(out_file, all_stripes, quad_cp): + def rotate(cp, e): + return np.rot90(cp.reshape(4,4,3), k=-e) + step_writer = STEPControl_Writer() + Interface_Static_SetCVal("write.step.schema", "AP203") + + for stripe in all_stripes: + s0cp = np.array([rotate(quad_cp[f],e) for f,e in stripe]) + bzlist = [cp_to_bz(s) for s in s0cp] + if len(stripe) > 1: + bb = compose_bezier(bzlist) + else: + bb = bzlist[0] + step_writer.Transfer(BRepBuilderAPI_MakeFace(bb, 1e-6).Shape(), + STEPControl_AsIs) + status = step_writer.Write(out_file) + return status \ No newline at end of file diff --git a/python/cadconvert/quad_curve_utils.py b/python/cadconvert/quad_curve_utils.py new file mode 100644 index 0000000..eb7fe24 --- /dev/null +++ b/python/cadconvert/quad_curve_utils.py @@ -0,0 +1,496 @@ +import fem_tabulator as feta +import numpy as np +import igl +import scipy +import tqdm +from sksparse.cholmod import cholesky_AAt + +ref_quad = np.array([[0, 0], [1, 0], [1, 1], [0, 1.]]) + +def quadratic_minimize(A, b, known=None): + if known is None: + if scipy.sparse.issparse(A): + return cholesky_AAt(A.T)(A.T@b) + else: + return scipy.linalg.solve(A.T@A, A.T@b) + kid, kval = known + mask = np.zeros(A.shape[1]) + mask[kid] = 1 + L = A.T@A + rhs = A.T@b + Lii = L[mask == 0][:, mask == 0] + Lik = L[mask == 0][:, mask == 1] + x = np.zeros((A.shape[1], 3)) + x[mask == 1] = kval + x[mask == 0] = scipy.sparse.linalg.spsolve(Lii, + rhs[mask == 0] - Lik@kval) + return x + +def combine_tris(face_i, face_s): + fi, fs = set(face_i), set(face_s) + common = fi & fs + assert len(common) == 2 + id_list = list(face_i) + list(fs-common) + tid = list(face_i).index(list(fi-common)[0]) + permute = ([[0, 1, 3, 2], [0, 1, 2, 3], [0, 3, 1, 2]])[tid] + q = np.array(id_list)[permute] + + id_map = {k: m for m, k in enumerate(q)} + trim = [id_map[k] for k in face_i], [id_map[k] for k in face_s] + assert np.sum(q[trim[0]] - face_i) == 0 + return q, trim + +def quad_trim_assign(siblings, faces): + quad_assign = - np.ones(len(faces), dtype=int) + cnt = 0 + trimmers = np.zeros_like(faces, dtype=int) + quads = [] + for i, s in enumerate(siblings): + if s < 0: + quad_assign[i] = -1 + continue + if quad_assign[i] >= 0: + assert quad_assign[s] >= 0 + continue + q, trims = combine_tris(faces[i], faces[s]) + + trimmers[[i,s]] = trims + quads.append(q) + quad_assign[i] = quad_assign[s] = cnt + cnt += 1 + assert trimmers.max() == 3 + assert quad_assign.max() * 2 < len(faces) + + q2t = - np.ones((len(quads), 2), dtype=int) + for t, q in enumerate(quad_assign): + if q < 0: + continue + if q2t[q][0] < 0: + q2t[q][0] = t + else: + q2t[q][1] = t + return quad_assign, q2t, trimmers, np.asarray(quads) + + +def tp_sample_value(x, y, order=3): + tg1d = np.array([order - np.arange(order + 1), np.arange(order + 1)]).T + bas_x = feta.bernstein_evaluator(x, x*0, x*0, tg1d).T + bas_y = feta.bernstein_evaluator(y, y*0, y*0, tg1d).T + qt = quad_tuple_gen(order) + bas_xy = np.array([bas_x[:, i]*bas_y[:, j] for i, j in qt]).T + return bas_xy + + +def quad_tuple_gen(order): + x, y = np.unravel_index(range((order+1)**2), + shape=(order+1, order+1)) + return np.vstack([x, y]).T + + +def sample_for_quad_trim(tr0, tr1, order: int): + cache_key = (tuple(tr0), tuple(tr1), order) + if cache_key in sample_for_quad_trim.cache: + return sample_for_quad_trim.cache[cache_key] + TG = np.asarray(feta.tuple_gen(order=order, var_n=2)) + V0 = TG@ref_quad[tr0].astype(int) + V1 = TG@ref_quad[tr1].astype(int) + bc = TG@ref_quad[[0, 1, 3]].astype(int) + + dim = (order+1, order+1) + bc_collect = np.zeros(((order+1)**2, + 4), dtype=int) + # this implicitly correlate to the quad_tuple_gen which unravels. + index0 = np.ravel_multi_index(V0.T, dims=dim) + index1 = np.ravel_multi_index(V1.T, dims=dim) + + bc_collect[index0, 0] = 0 + bc_collect[index0, 2:] = bc + bc_collect[index1, 0] = 1 + bc_collect[index1, 2:] = bc + bc_collect[:, 1] = order - bc_collect[:, 2:].sum(axis=1) + + sample_for_quad_trim.cache[cache_key] = bc_collect + return bc_collect + + +sample_for_quad_trim.cache = dict() + + +def local_codecs_on_edge(level: int): + """List of codecs (partial) to resolve duplicates on edges for quads. + + Args: + level (int): number of points on the side. (order + 1) + """ + std_x, std_y = quad_tuple_gen(level).T + codes = [] + for (x, y) in zip(std_x, std_y): + if x % level == 0 or y % level == 0: + bc = [level - x, x, level - y, y] + zid = bc.index(0) + side = [[1, 2], [0, 3], [3, 2], [0, 1]][zid] + bc = bc[2:] if zid <= 1 else bc[:2] + code = tuple(side[:1]*bc[0] + side[1:]*bc[1]) + codes.append(code) + else: + codes.append(tuple([-1]*level)) + return codes + + +def quad_ho_F(quads, level: int): + std_cod = local_codecs_on_edge(level) + + global_ids = np.zeros((len(quads), len(std_cod)), dtype=int) + global_id_store = dict() + cnt = 0 + for q in range(len(quads)): + for i, c in enumerate(std_cod): + assert len(c) == level + if c[0] < 0: # internal. + node_id = cnt + cnt += 1 + else: + key = tuple(sorted(list(quads[q][cc] for cc in c))) + if key not in global_id_store: + global_id_store[key] = cnt + cnt += 1 + node_id = global_id_store[key] + global_ids[q, i] = node_id + assert global_ids[0].max() == len(global_ids[0]) - 1 + assert global_ids.shape[1] == (level+1)**2 + return global_ids + + +def highorder_tp_sv(cp, level=3, order=3): + def local_upsample(level: int): + usV, usF = igl.upsample(np.array([[0, 0], [1, 0], [1, 1], [0, 1.]]), + np.array([[0, 1, 2], [0, 2, 3]]), level) + bnd0 = igl.boundary_loop(usF) + usE = np.vstack([bnd0, np.roll(bnd0, -1)]).T + return usV, usF, usE + usV, usF, usE = local_upsample(level=level) + bas_uv = tp_sample_value(usV[:, 0], usV[:, 1], order) + sv = bas_uv@cp.reshape(len(cp), -1, 3) + return (sv, + np.vstack([usF+i*len(usV) for i in range(len(sv))]), + np.vstack([usE+i*len(usV) for i in range(len(sv))])) + + +def bezier_subd(cp): + # de Casteljau + l = [cp] + for i in range(len(cp)-1): + l.append((l[-1][1:] + l[-1][:-1])/2) + return np.asarray([i[0] for i in l]), np.asarray([i[-1] for i in l[::-1]]) + + +def bezier_eval(x, cp, order=3): + tg1d = np.array([order - np.arange(order + 1), np.arange(order + 1)]).T + bas_x = feta.bernstein_evaluator(x, x*0, x*0, tg1d).T + return bas_x@cp + +def quad_fit(V, F, quads, q2t, trim_types, level, order, bsv, query, regularizer=None): + """Least Square Fitting for QuadMesh, with Tensor Product Bezier Basis + + Args: + V ([type]): [description] + F ([type]): [description] + quads ([type]): [description] + q2t ([type]): Qx2, two triangles it points to + trim_types (npt.Array): used for converting quad parameter values to triangle and barycentric coordinates. + level ([type]): levels for upsample + order ([type]): [description] + bsv ([type]): [description] + query (Callable): [description] + + Returns: + [type]: [description] + """ + F_qh_sample = quad_ho_F(quads, level=level) + F_qh_order = quad_ho_F(quads, level=order) + assert F_qh_sample.shape[1] == bsv.shape[0] + + all_samples = np.zeros((F_qh_sample.max() + 1, 3)) + rows, cols, vals = [], [], [] + bsv = bsv.tocoo() + br, bc, bv = bsv.row, bsv.col, bsv.data + for q, (t0, t1) in enumerate(tqdm.tqdm(q2t, desc='Quad Fit')): + tbc0 = np.array(sample_for_quad_trim(trim_types[t0], trim_types[t1], level), + dtype=int) + tbc0[:, 0] = np.asarray([t0, t1])[tbc0[:, 0]] + sample_vals = query(tbc0, denom=level) + + all_samples[F_qh_sample[q]] = sample_vals + sample_q, order_q = F_qh_sample[q], F_qh_order[q] + rows = np.concatenate([rows, sample_q[br]]) + cols = np.concatenate([cols, order_q[bc]]) + vals = np.concatenate([vals, bv]) + if regularizer is not None: + offset = 0 + rnum = regularizer.shape[0] + regularizer = regularizer.tocoo() + r_rows, r_cols, r_vals = regularizer.row, regularizer.col, regularizer.data + for q, (t0, t1) in enumerate(tqdm.tqdm(q2t)): + order_q = F_qh_order[q] + rows = np.concatenate([rows, r_rows + q*rnum]) + cols = np.concatenate([cols, order_q[r_cols]]) + vals = np.concatenate([vals, r_vals]) + print(all_samples.shape) + all_samples = np.concatenate([all_samples, np.zeros((rnum*len(q2t),3))]) + + ijk, idx = np.unique(np.vstack([rows, cols]).astype(int), axis=1, return_index=True) + A = scipy.sparse.csr_matrix(arg1=(np.array(vals)[idx], + ijk)) + assert A.shape[0] == all_samples.shape[0] + sol = quadratic_minimize(A, all_samples) + res = np.linalg.norm(A@sol - all_samples, axis=1) + print('Residual', res.mean(), res.max()) + all_cp = sol[F_qh_order] + return all_cp, all_samples + + +def solo_cc_split(V, F, siblings, t2q, quads, q_cp, order: int, subd = None): + """Catmul-Clark style split for the solo triangles. + + Default, use de Casteljau (or de Boor) subdivision algorithm to assign for the subsequent constraints. + Returns information for constrained nodes and more. + Args: + V (np.array): TriMesh Vertices + F (np.array): TriMesh Faces + siblings (np.array): tris to sibling tris indices, only used -1 indicator here. + t2q ([type]): tris to quad indices. + quads ([type]): quad F matrix + q_cp ([type]): control points for the existing quads. + order (int): Polynomial order for the quad before the split. + + Returns: + known_cp: a map from tuples to control points. Used for subsequent fitting. + """ + assert len(q_cp) == len(quads) + solos = np.where(siblings == -1)[0] + if subd is None: + subd = bezier_subd + TT, TTi = igl.triangle_triangle_adjacency(F) + + edge_node_map = {tuple(sorted(k)): i + for i, k in enumerate(local_codecs_on_edge(order))} + edge_id = - np.ones_like(TT) + avail_id = F.max() + 1 + known_cod2cp = dict() + + ## CC split to quads, while recording known_cod2cp + newverts = [] + for f in solos: + for e in range(3): + v0, v1 = F[f, e], F[f, (e+1) % 3] + fo, eo = TT[f, e], TTi[f, e] + if fo < 0: # boundary, future more can go here + edge_id[f, e] = avail_id + newverts.append((V[v0]+V[v1])/2) + avail_id += 1 + continue + qid = t2q[fo] + if qid < 0: # opposite is not an existing quad + if edge_id[f, e] < 0: + edge_id[f, e] = edge_id[fo, eo] = avail_id + newverts.append((V[v0]+V[v1])/2) + avail_id += 1 + continue # simply register the available vertex without considering known_cod2cp assignments + assert edge_id[f, e] < 0, "Could not have visited twice." + + quad = list(quads[qid]) + qv0, qv1 = [quad.index(v0), quad.index(v1)] + # collect on this edge + + tup_list = [tuple(sorted([qv0]*(order - t1) + [qv1]*t1)) for t1 in range(order+1)] + cp_list = np.asarray([q_cp[qid][edge_node_map[t]] + for t in tup_list]) # Bezier control points of the original edge + cp0, cp1 = subd(cp_list) + for t1 in range(order + 1): + t0 = order - t1 + known_cod2cp[tuple(sorted([v0]*t0 + [avail_id]*t1))] = cp0[t1] + known_cod2cp[tuple(sorted([avail_id]*t0 + [v1]*t1))] = cp1[t1] + + edge_id[f, e] = avail_id + newverts.append((V[v0]+V[v1])/2) + avail_id += 1 + newquads = [] + for pid, fi in enumerate(solos, avail_id): # new internal id + v0, v1, v2 = F[fi] + for ei in range(3): + v0, e0, e2 = F[fi, ei], edge_id[fi, ei], edge_id[fi, (ei+2) % 3] + newquads += [[v0, e0, pid, e2]] + + assert np.asarray(newquads).max() + 1 == avail_id + len(solos) + return (np.concatenate([newverts, igl.barycenter(V, F[solos])]), + known_cod2cp, newquads) + + +def constrained_cc_fit(V, F, siblings, newquads, known_cp, level:int, order:int, bsv, query): + """Constrained Bezier fitting for the c-c split quads. + + Args: + V ([type]): [description] + F ([type]): [description] + siblings ([type]): [description] + newquads ([type]): [description] + known_cp (Map): codec -> control coeffs. + level (int): Sampling level for least squares. + order (int): Curved Order + bsv ([type]): [description] + query ([type]): lambda funuction, takes (tri_bc, a denominator scalar) returns sampled points. + + Returns: + [type]: [description] + """ + solos = np.where(siblings == -1)[0] + newF = F[siblings == -1] + assert len(newF)*3 == len(newquads) + F_sa = quad_ho_F(newquads, level=level) + F_or = quad_ho_F(newquads, level=order) + + std_codec = local_codecs_on_edge(order) + known_dict = dict() # relates xid in the system to the known control value. + for q, quad in enumerate(tqdm.tqdm(newquads, desc='Constrained CC')): + for e, code in enumerate(std_codec): + if code[0] < 0: + continue + key = tuple(sorted(list(quad[cc] for cc in code))) + if key in known_cp: + if F_or[q, e] in known_dict: + assert (known_dict[F_or[q, e]] == known_cp[key]).all() + else: + known_dict[F_or[q, e]] = known_cp[key] + + known_ids = np.ones(len(known_dict), dtype=int) + known_vals = np.ones((len(known_dict), 3)) + for i, (k, v) in enumerate(sorted(known_dict.items())): + known_ids[i], known_vals[i] = k, v + + all_samples = np.zeros((F_sa.max() + 1, 3)) + ijk, vals = [], [] + + # hardcoded pattern for catmull-clark + # used to warp the query barycentric coordinates. + refpts = np.array([[0, 0], + [3, 0], + [6, 0], + [3, 3], + [0, 6], + [0, 3], + [2, 2]]) # denom=6 + cc_q = np.array([[0, 1, 6, 5], [2, 3, 6, 1], [4, 5, 6, 3]]) + type2bc = [] + for t in range(3): + x, y = np.array(quad_tuple_gen(level)).T + x = x.reshape(-1, 1) + y = y.reshape(-1, 1) + assert x.max() <= level + qc = refpts[cc_q[t]].reshape(-1, 1, 2) + type2bc.append(((level-x)*qc[0] + x*qc[1])*(level-y) + + ((level-x)*qc[3] + x*qc[2])*y) + rows, cols, vals = [], [], [] + bsv = bsv.tocoo() + br, bc, bv = bsv.row, bsv.col, bsv.data + for qi, q in enumerate(tqdm.tqdm(newquads)): + tbc = np.zeros(((level+1)**2, 4), + dtype=int) + tbc[:, 0] = solos[qi//3] + tbc[:, 2:] = type2bc[qi % 3] + tbc[:, 1] = 6*(level)**2 - tbc[:, 2:].sum(axis=1) + assert tbc[:, 1].min() >= 0 + sample_vals = query(tbc, denom=6*(level)**2) + all_samples[F_sa[qi]] = sample_vals + sample_q, order_q = F_sa[qi], F_or[qi] + rows = np.concatenate([rows, sample_q[br]]) + cols = np.concatenate([cols, order_q[bc]]) + vals = np.concatenate([vals, bv]) + + ijk, idx = np.unique(np.vstack([rows, cols]).astype(int), axis=1, return_index=True) + A = scipy.sparse.csr_matrix(arg1=(np.array(vals)[idx], + ijk)) + + sol = quadratic_minimize(A, all_samples, (known_ids, known_vals)) + res = np.linalg.norm(A@sol - all_samples, axis=1) + print('Residual', res.mean(), res.max()) + all_cp = sol[F_or] + return all_cp + +def query(t_bc_i, denom): + # requires additional input of (mB,mT,F,aabb(refV,refF)) + all_f = query.F[t_bc_i[:,0]] + base_pts = np.einsum('sed,se->sd', query.mB[all_f], t_bc_i[:,1:])/denom + top_pts = np.einsum('sed,se->sd', query.mT[all_f], t_bc_i[:,1:])/denom + list_fid = np.zeros(len(base_pts),dtype=int) + list_bc = np.zeros((len(base_pts),3)) + for i, (b,t) in enumerate(zip(base_pts, top_pts)): + list_fid[i], u,v, _ = query.aabb.segment_hit(b,t, False) + list_bc[i] = (1-u-v,u,v) + return np.einsum('sed,se->sd', query.inpV[query.refF[list_fid]], list_bc) + + +def split_square(width): + x,y = quad_tuple_gen(width).T#np.meshgrid(np.arange(width+1), np.arange(width+1)) + tr = lambda i,j: np.ravel_multi_index((i,j), dims=(width+1,width+1)) + tris= [] + for i in range(width): + for j in range(width): + tris += [[tr(i,j), tr(i+1,j), tr(i,j+1)], [tr(i+1,j),tr(i+1,j+1), tr(i,j+1)]] + return np.vstack([x.flatten(), y.flatten()]).T, np.array(tris) + +def top_and_bottom_sample(mB, mT, F, quads, q2t, trim_types, level): + all_base = np.zeros((len(quads), (level+1)**2, 3)) + all_top = np.zeros((len(quads), (level+1)**2, 3)) + + for q, (t0, t1) in enumerate(q2t): + tbc0 = np.array(sample_for_quad_trim(trim_types[t0], trim_types[t1], level), + dtype=int) + tbc0[:, 0] = np.asarray([t0, t1])[tbc0[:, 0]] + all_base[q]= (np.einsum('fed,fe->fd',mB[F[tbc0[:, 0]]], tbc0[:,1:])/level) + all_top[q] = (np.einsum('fed,fe->fd',mT[F[tbc0[:, 0]]], tbc0[:,1:])/level) + + return np.array(all_base), np.array(all_top) + +import collections +def control_points_for_quad_edges(quads, valid, q_cp) -> dict: + """ + Based on valid, get the control points on the border and store to a dictionary + """ + assert len(q_cp) == len(quads) + assert len(quads) == len(valid) + known_cp = dict() + order = np.round(np.sqrt(q_cp.shape[1])).astype(int) - 1 + + edge_node_map = {tuple(sorted(k)): i + for i, k in enumerate(local_codecs_on_edge(order))} + # Collect some edges. + quad_connect = collections.defaultdict(lambda:[None,None]) + for q, _ in enumerate(quads): + # print(q) + if not valid[q]: + continue + for e in range(4): + v0, v1 = quads[q,e], quads[q,(e+1)%4] + ind = 0 + if v0 > v1: + v0,v1 = v1,v0 + ind = 1 + quad_connect[v0,v1][ind] = q + for (v0,v1), (q0,q1) in quad_connect.items(): + if (q0 is not None) and (q1 is not None): + continue + if q0 is None: + v0, v1 = v1, v0 + q0, q1 = q1, q0 + assert q0 is not None # readout v0,v1 + quad = list(quads[q0]) + qv0, qv1 = [quad.index(v0), quad.index(v1)] + tup1d = np.array([order - np.arange(order + 1), + np.arange(order + 1)]).T + + tup_list = [tuple(sorted([qv0]*t0 + [qv1]*t1)) for t0, t1 in tup1d] + cp_list = np.asarray([q_cp[q0][edge_node_map[t]] + for t in tup_list]) + known_cp[(v1,v0)] = cp_list[::-1] + return known_cp + diff --git a/python/cadconvert/quad_utils.py b/python/cadconvert/quad_utils.py new file mode 100644 index 0000000..3cbe4d5 --- /dev/null +++ b/python/cadconvert/quad_utils.py @@ -0,0 +1,267 @@ +import numpy as np +from collections import defaultdict +import scipy +import tqdm +import igl + + +def catmul_clark_split(v,f): + verts = list(v) + tt, tti = igl.triangle_triangle_adjacency(f) + edge_id = - np.ones_like(tt) + avail_id = len(v) + for fi, _ in enumerate(tt): + for ei in range(3): + if edge_id[fi,ei] == -1: + edge_id[fi,ei] = avail_id + edge_id[tt[fi,ei],tti[fi,ei]] = avail_id + v0, v1 = f[fi,ei], f[fi,(ei+1)%3] + verts.append((v[v0]+v[v1])/2) + avail_id += 1 + + bc = igl.barycenter(v,f) + quads = [] + for fi, [v0,v1,v2] in enumerate(f): + pid = fi + avail_id + for ei in range(3): + v0, e0, e2 = f[fi,ei], edge_id[fi,ei], edge_id[fi,(ei+2)%3] + quads += [[v0,e0,pid,e2]] + return np.asarray(verts + list(bc)), np.asarray(quads) + + +def q2e(f): + return np.array([f[:,i] for i in [[0, 1],[1,2],[2,3],[3,0]]]).reshape(-1,2) + + +def greedy_pairing(f, score, sharp_markers): + tt, tti = igl.triangle_triangle_adjacency(f) + occupied = -np.ones(len(f),dtype=int) + qid = 0 + pairs = [] + queue = [] + for fi in range(len(f)): + for e in range(3): + if sharp_markers[fi,e]: + continue + queue.append((score(fi,e), fi,e)) + queue = sorted(queue) + for _, fi,e in queue: + if occupied[fi] >= 0: + continue + fo = tt[fi,e] + if fo < 0: + continue + if occupied[fo] >= 0: + continue + occupied[fi] = fo + occupied[fo] = fi + q = list(f[fi]) + q.insert(e+1, f[fo][tti[fi,e]-1]) + pairs.append(q) + qid += 1 + return occupied, pairs + + +def crawling(hybrid): + def set_conn(v0,v1,x): + if v0 0: + if (len(bfs_queue)) > len(hybrid): + print('Wrong!InfLoop Warn') + assert False + + p = bfs_queue[0] + bfs_queue.pop(0) + t = p[-1][0] + ne = len(hybrid[t]) + + for i in range(ne): + e = hybrid[t][(i+1)%ne], hybrid[t][i] + fo = get_conn(*e) # oppo + if fo is None or fo in visited: + continue + visited.add(fo) + bfs_queue.append(p + [(fo,e)]) +# print(fo) +# print(fo, hybrid[fo]) + if len(hybrid[fo]) == 3: + return bfs_queue[-1] + return None + + def merge_tri_quad(tri, edge, quad): # this is a backward path + q, e0 = quad + v0,v1 = edge + qv = list(hybrid[q].copy()) +# print('In', tri, edge, qv) + x = list(set(tri)-set(edge))[0] # opposite vertex + qv.insert(qv.index(v0), x) +# print('Pent',qv, e0) + + if len(qv) == 4: + assert e0 is None + return None, None, qv + newtri = [qv[qv.index(e0[0])-1], e0[0], e0[1]] + qv.remove(e0[0]) + return newtri, e0, qv + + connect = defaultdict(lambda:[None,None]) + cnt_tri = 0 + for fi, f in enumerate(hybrid): + ne = len(f) + if ne == 3: + cnt_tri += 1 + for i in range(ne): + set_conn(f[i], f[(i+1)%ne],fi) + + for i in range(cnt_tri//2 + 1): + bachelor = np.where(np.array([len(h) for h in hybrid]) == 3)[0] + + if len(bachelor) == 0: + break + path = bfs_find_tri(bachelor[0]) + + tid, edge = path[-1] + tri = hybrid[tid] + new_quads = [] + pid = len(path) - 2 + while tri is not None: + tri, edge, new_q = merge_tri_quad(tri, edge, path[pid]) + pid -= 1 + new_quads.append(new_q) + + for p,_ in path: + f = hybrid[p] + ne = len(f) + for i in range(ne): + set_conn(f[i], f[(i+1)%ne], None) + hybrid[p] = [] + for f in new_quads: + ne = len(f) + for i in range(ne): + set_conn(f[i], f[(i+1)%ne], len(hybrid)) + hybrid.append(f) + return hybrid + +def split_edges_between_odd_components(f, cp): + tt,tti = igl.triangle_triangle_adjacency(f) + assert tt.min() >=0 + arr = defaultdict(lambda:set()) + for i in range(len(f)): + r0 = cp[i] + for j in range(3): + r1 = cp[tt[i,j]] + if r0 != r1: + v0, v1 = f[i,j], f[i,j-2] + arr[(min(r0,r1), max(r0,r1))].add((min(v0,v1), max(v0,v1))) + ijk = np.array([(k[0],k[1],1) for k, l in arr.items()],dtype=int) + n_comp = cp.max()+1 + comp_adj = scipy.sparse.coo_matrix(arg1=(ijk[:,2], ijk[:,:2].T), shape=(n_comp,n_comp)) + + odd_comps = np.where([np.count_nonzero(cp==i)%2 for i in range(max(cp) + 1)])[0] + assert len(odd_comps)%2 == 0 + + length, predec = scipy.sparse.csgraph.shortest_path(comp_adj, directed=False, indices=odd_comps, return_predecessors=True) + + def trace_path(prec, i,j): + assert prec[i] <0 + path = [j] + while True: + p = prec[path[-1]] + if p < 0: + break + path.append(p) + return path + + odd_pairing = [] + + + visited = np.zeros(len(odd_comps)) + for i,c in enumerate(odd_comps): + if visited[i] == 1: + continue + length[i][c] = np.inf + length[i][odd_comps[visited==1]] = np.inf + j = length[i][odd_comps].argmin() + visited[i] = visited[j] = 1 + path = trace_path(predec[i], c, odd_comps[j]) + odd_pairing.append((c, odd_comps[j], path)) + + split_pairs = [] + for i,j, path in odd_pairing: + split_pairs += [tuple(sorted(i)) for i in zip(path[:-1], path[1:])] + return [next(iter(arr[s])) for s in split_pairs] + + +def edge_dots(V, F): + FN = igl.per_face_normals(V,F, np.ones(3)) + tt, tti = igl.triangle_triangle_adjacency(F) + return np.einsum('fad,fd->fa',FN[tt],FN) + +# BFS traversal +def group_quads(quads): + def quad_connectivity(quads): + connect = defaultdict(lambda:[None,None]) + cnt_tri = 0 + def set_conn(v0,v1,x): + if v0= 0: + break + vertset = vertset.union(list(quads[f_oppo])) + if len(vertset) != 2*(len(stripe)+2): + break # not topology stripe. + stripe_paint[f_oppo] = next_color + v0id = list(quads[f_oppo]).index(v0) + e_next = (v0id + 1) %4 + f,e = f_oppo, e_next + stripe.append((f_oppo, v0id)) + next_color += 1 + all_stripes.append(stripe) + return stripe_paint, all_stripes \ No newline at end of file diff --git a/python/cadconvert/tensor-product-fit.py b/python/cadconvert/tensor-product-fit.py new file mode 100644 index 0000000..fa32e67 --- /dev/null +++ b/python/cadconvert/tensor-product-fit.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python +# coding: utf-8 +import fem_tabulator as feta +import h5py +import quad_curve_utils as qr +import quad_utils +import numpy as np +import igl +import sys +sys.path.append('/home/zhongshi/Workspace/bichon/python/debug') +import prism +import os +import scipy +import bezier as qb +import tqdm +from occ_step import cp_write_to_step +import occ_step + +def valid_pairing(faces, score, sharp_markers, valid_combine_func): + tt, tti = igl.triangle_triangle_adjacency(faces) + occupied = -np.ones(len(faces),dtype=int) + qid = 0 + pairs = [] + queue = [] + for fi in range(len(faces)): + for e in range(3): + if sharp_markers[fi,e]: + continue + queue.append((score(fi,e), fi,e)) + queue = sorted(queue) + for _, fi,e in queue: + if occupied[fi] >= 0: + continue + fo = tt[fi,e] + + if fo < 0: + continue + if occupied[fo] >= 0: + continue + + if not valid_combine_func(fi, fo): # combine fi with fo. + continue + + occupied[fi] = fo + occupied[fo] = fi + # q = list(faces[fi]) + # q.insert(e+1, faces[fo][tti[fi,e]-1]) + # pairs.append(q) + qid += 1 + return occupied, pairs + + +def main(input_file, output_file = None, order =3, level=6, post_check=False): + with h5py.File(input_file, 'r') as f: + V,F,refV,refF,inpV,mB,mT = map(lambda x:f[x][()], ('mV','mF','ref.V','ref.F','inpV','mbase','mtop')) + + ## Bezier fitting + A = scipy.sparse.coo_matrix(qb.bezier_fit_matrix(order, level)).tocsr() + query = qr.query + query.aabb, query.F, query.mB, query.mT, query.inpV, query.refF = prism.AABB(refV, refF), F, mB, mT, inpV, refF + + def valid_combine_func(fi, fo) -> bool: + # First fit + quad, trims = qr.combine_tris(F[fi], F[fo]) + tbc0 = np.array(qr.sample_for_quad_trim(trims[0], trims[1], level), + dtype=int) + tbc0[:, 0] = np.asarray([fi, fo])[tbc0[:, 0]] + sample_vals = query(tbc0, denom=level) + local_cp = qr.quadratic_minimize(A, sample_vals) + # Second, check + v = qb.bezier_check_validity(mB, mT, F[[fi,fo]], quad.reshape(-1,4), np.array([[0,1]]), + trims, np.array([local_cp]), order, + valid_check = prism.elevated_positive_check, + progress_bar=False) + + return v + siblings, _ = valid_pairing(F, score = lambda f,e: -np.linalg.norm(V[F[f,e]] - V[F[f,(e+1)%3]]), + sharp_markers=quad_utils.edge_dots(V,F) < np.cos(np.pi/4), + valid_combine_func=valid_combine_func) + print('empty siblings', np.count_nonzero(siblings == -1), '/', len(siblings)) + t2q, q2t,trim_types, quads = qr.quad_trim_assign(siblings, F) + + + quad_cp, samples = qr.quad_fit(V, F, quads, q2t, trim_types, level, order, A, query, None) + + if post_check: + valids = bezier_check_validity(mB, mT, F,quads, q2t, trim_types, quad_cp, 3) + + print('valids', np.count_nonzero(valids), '/', len(valids)) + # adjust quads, quads_cp, q2t, t2q based on validity + siblings[q2t[valids==False].flatten()] = -1 + quads = quads[valids] + quad_cp = quad_cp[valids] + q2t = q2t[valids] + t2q = np.ones_like(t2q) * -1 + for q, (t0,t1) in enumerate(q2t): + t2q[t0] = q + t2q[t1] = q + + new_v, known_cp, newquads = qr.solo_cc_split(V, F, siblings, t2q, quads, quad_cp, order, subd=None) + cc_cp = qr.constrained_cc_fit(V, F, siblings, newquads, known_cp, level, order, A, query) + if output_file is None: + output_file = f'/home/zhongshi/ntopo/ntopmodels/fit/{os.path.basename(input_file)}' + _, stripe0 = quad_utils.group_quads(quads) + _, stripe1 = quad_utils.group_quads(np.array(newquads)) + print(f'Quad Counts: {len(quads)} + {len(newquads)}') + print(f'Stripe Counts: {len(stripe0)} + {len(stripe1)}') + offset = lambda x: (x[0] + len(quads), x[1]) + occ_step.stripe_writer(output_file + '.stp', + stripe0 + [list(map(offset, stripe)) for stripe in stripe1], + np.vstack([quad_cp,cc_cp])) + + np.savez(output_file + '.npz', quad_cp = quad_cp, cc_cp = cc_cp, quads=quads, newquads=newquads, + stripe0=stripe0, stripe1=stripe1) + +def test_stripe(): + with np.load('temp.npz') as npl: + quads, quad_cp = npl['quads'], npl['quad_cp'] + + # quads -- quad_cp + stripe_paint, all_stripes = quad_utils.group_quads(quads) + print('num of stripes', len(all_stripes)) + stripe_writer(all_stripes, 'quad_cp, stripe.stp') + +if __name__ == '__main__': + import fire + fire.Fire(main) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 256463e..d9d3d86 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -72,3 +72,11 @@ if (ENABLE_ASAN) target_compile_options(cumin_bin PUBLIC "-fsanitize=address") target_link_options(cumin_bin PUBLIC "-fsanitize=address") endif() + +option(PYBICHON "Enable Python Binding" OFF) +if (PYBICHON) + find_package(PythonInterp QUIET) + prism_download_pybind11() + add_subdirectory(${PRISM_EXTERNAL}/pybind11/ pybind11) + add_subdirectory(python) +endif() \ No newline at end of file diff --git a/src/cumin/curve_utils.cpp b/src/cumin/curve_utils.cpp index f457bdb..fe25319 100755 --- a/src/cumin/curve_utils.cpp +++ b/src/cumin/curve_utils.cpp @@ -289,11 +289,7 @@ bool elevated_positive(const std::vector &base, if (recurse_check) { for (auto &t : tens) { // this is codec_bc. - if (!prism::curve::tetrahedron_inversion_check( - t, helper.volume_data.vol_codec, - helper.volume_data.vol_jac_codec, - helper.volume_data.vol_bern_from_lagr, - helper.volume_data.vol_jac_bern_from_lagr)) { + if (!prism::curve::tetrahedron_inversion_check(t)) { spdlog::debug("blocked by recursive {}", mf); return false; } diff --git a/src/cumin/curve_utils.hpp b/src/cumin/curve_utils.hpp index c243016..7c5ba4d 100644 --- a/src/cumin/curve_utils.hpp +++ b/src/cumin/curve_utils.hpp @@ -48,7 +48,7 @@ struct HelperTensors { inline HelperTensors &magic_matrices(int tri_order = -1, int level = -1) { if (!HelperTensors::tensors_) { if (tri_order == -1) { - throw 1; + throw std::runtime_error("HelperTensors::init(int tri_order, int level) must be called first"); } HelperTensors::init(tri_order, level); } @@ -56,7 +56,7 @@ inline HelperTensors &magic_matrices(int tri_order = -1, int level = -1) { if (HelperTensors::tensors_->tri_order != tri_order || HelperTensors::tensors_->level != level) { // unmatched. - throw 1; + throw std::runtime_error("HelperTensors::init(int tri_order, int level) called with mismatched parameters"); } } return *HelperTensors::tensors_; diff --git a/src/cumin/curve_validity.cpp b/src/cumin/curve_validity.cpp index 080af79..c511cb2 100755 --- a/src/cumin/curve_validity.cpp +++ b/src/cumin/curve_validity.cpp @@ -347,7 +347,7 @@ constexpr auto valid_curving = return false; }; -auto mean_curv_normals = [](const RowMatd &V, const RowMati &F) { +auto mean_curv_normals = [](const RowMatd &V, const RowMati &F) { // this is an inactive function for fitting normals. using namespace Eigen; RowMatd HN; SparseMatrix L, M, Minv; diff --git a/src/cumin/high_order_optimization.cpp b/src/cumin/high_order_optimization.cpp index 0bf5ffd..ba27c2a 100644 --- a/src/cumin/high_order_optimization.cpp +++ b/src/cumin/high_order_optimization.cpp @@ -630,7 +630,7 @@ int cutet_collapse(RowMatd &lagr, RowMati &p4T, double stop_energy) { auto n_node = codecs_o4.rows(); if (n_node != p4T.cols()) { spdlog::critical("tet order mismatch"); - throw 1; + throw std::runtime_error("tet order mismatch"); } std::vector vertices(lagr.rows()); @@ -842,11 +842,11 @@ int cutet_collapse(RowMatd &lagr, RowMati &p4T, double stop_energy) { auto [val, ign] = mips_energy(nodes35, vec_dxyz, false); if (std::isnan(val)) { spdlog::critical("PostNAN!!!"); - throw 1; + throw std::runtime_error("PostNAN!!!"); } if (val > 1e8) { spdlog::critical("large energy"); - throw 1; + throw std::runtime_error("large energy"); } } return true; @@ -1090,7 +1090,7 @@ int edge_collapsing(RowMatd &lagr, RowMati &p4T, double stop_energy) { auto n_node = codecs_o4.rows(); if (n_node != p4T.cols()) { spdlog::critical("tet order mismatch"); - throw 1; + throw std::runtime_error("tet order mismatch"); } std::vector vertices(lagr.rows()); diff --git a/src/cumin/stitch_surface_to_volume.cpp b/src/cumin/stitch_surface_to_volume.cpp index 96f4a7e..a8a385a 100644 --- a/src/cumin/stitch_surface_to_volume.cpp +++ b/src/cumin/stitch_surface_to_volume.cpp @@ -14,7 +14,7 @@ bool prism::curve::stitch_surface_to_volume( for (auto i = 0; i < base.rows(); i++) { if (base.row(i) != Vmsh.row(i)) { spdlog::critical("Stitch Precondition."); - throw 1; + throw std::runtime_error("Stitch Precondition."); } } spdlog::info("Setting up stitch: base {}, Fsh {}, Vmsh {} Tmsh {}", diff --git a/src/curve_in_shell.cpp b/src/curve_in_shell.cpp index cc6db69..184bccc 100755 --- a/src/curve_in_shell.cpp +++ b/src/curve_in_shell.cpp @@ -82,6 +82,7 @@ int main(int argc, char **argv) { {"skip_collapse", false}, {"skip_split", true}, {"skip_volume", false}, + {"danger_relax_precondition", false}, // this is a experiment switch: bypass thresholds in precondition, the result may or may not encounter floating point failures. }; config["tetfill"] = {{"tetwild", true}}; config["cutet"] = { diff --git a/src/pipeline_schedules.cpp b/src/pipeline_schedules.cpp index 826dff1..080664f 100644 --- a/src/pipeline_schedules.cpp +++ b/src/pipeline_schedules.cpp @@ -122,7 +122,8 @@ constexpr auto consistent_orientation = [](const RowMati &F) { }; bool preconditions(const RowMatd &V, const RowMati &F, - const std::string &filename) { + const std::string &filename, bool shortcircuit=true) { + spdlog::info("Preconditions: Checking consistency of input data, with short circuit {}", shortcircuit); if (F.rows() == 0) { spdlog::error("Precondition: Empty Mesh"); return false; @@ -149,26 +150,29 @@ bool preconditions(const RowMatd &V, const RowMati &F, double minarea = M.minCoeff() / 2; if (minarea < 1e-8) { spdlog::error("Precondition: Input {} area small {}", filename, minarea); - return false; + if (shortcircuit) + return false; } double minangle = min_dihedral_angles(V, F); spdlog::info("Minimum Angle {:.18f}", minangle); if (minangle < -1 + 0.0006091729809042379) { // smaller than 2 degree - spdlog::error("Precondition: Input {} flat angle", filename); - return false; + spdlog::error("Precondition: Input {} flat angle {}", filename, minangle); + if (shortcircuit) + return false; } // numerical self intersection prism::geogram::AABB tree(V, F, true); if (tree.numerical_self_intersection(1e-6)) { spdlog::error("Precondition: Input {} N-self intersects", filename); - return false; + if (shortcircuit) + return false; } std::vector vecV; std::vector vecF; eigen2vec(V, vecV); eigen2vec(F, vecF); if (!prism::spatial_hash::self_intersections(vecV, vecF).empty()) { - spdlog::error("Precondition: Input {} self intersects", filename); + spdlog::critical("Precondition: Input {} self intersects", filename); return false; } return true; @@ -333,7 +337,7 @@ void cutet_optim(RowMatd &lagr, RowMati &p4T, nlohmann::json config) { if (!(prism::curve::InversionCheck(lagr, p4T, codecs_o4_bc, codecs_o9_bc, bern_from_lagr_o4, bern_from_lagr_o9))) { spdlog::error("Inversion check failure: " + str + " "); - throw 1; + throw std::runtime_error("Inversion check failure"); } }; @@ -511,7 +515,7 @@ void feature_and_curve(std::string filename, std::string fgname, spdlog::info("V={}, F={}", V.rows(), F.rows()); put_in_unit_box(V); - if (!preconditions(V, F, filename)) return; + if (preconditions(V, F, filename, !control_cfg["danger_relax_precondition"]) ==false) return; } RowMati feature_edges; Eigen::VectorXi feature_corners; diff --git a/src/prism/PrismCage.cpp b/src/prism/PrismCage.cpp index 57570ca..52889ef 100755 --- a/src/prism/PrismCage.cpp +++ b/src/prism/PrismCage.cpp @@ -162,7 +162,7 @@ PrismCage::PrismCage(const RowMatd &vert, const RowMati &face, spdlog::error("Something wrong with reorder"); } } - + constraints_per_face.resize(ref.F.rows()); // singular split may increase. TODO: re-assign some of the constraints. for (int i = 0; i < ref.F.rows(); i++) { auto [s, mt, shift] = tetra_split_AorB({ref.F(i, 0), ref.F(i, 1), ref.F(i, 2)}); @@ -192,7 +192,7 @@ PrismCage::PrismCage(const RowMatd &vert, const RowMati &face, for (int i = 0; i < dsF.rows(); i++) { auto [s, mt, shift] = tetra_split_AorB({dsF(i, 0), dsF(i, 1), dsF(i, 2)}); for (int j = 0; j < 3; j++) dsF(i, j) = mt[j]; - if (mt[1] < num_cons || mt[2] < num_cons) throw 10; + if (mt[1] < num_cons || mt[2] < num_cons) throw std::runtime_error("only one vertex in each triangle can be singular."); } ref.aabb.reset( diff --git a/src/prism/bevel_utils.cpp b/src/prism/bevel_utils.cpp index f8d9d22..f554b19 100755 --- a/src/prism/bevel_utils.cpp +++ b/src/prism/bevel_utils.cpp @@ -317,7 +317,7 @@ void prism::bevel_utils::edge_based_bevel(const RowMatd& V, const RowMati& F, auto [e0, ve] = [&mark]() { for (auto j = 0; j < 3; j++) if (mark(j) >= 0) return std::pair(j, mark(j)); - throw 1; + assert(false); return std::pair(-1, -1); }(); auto v0 = Ff(e0); @@ -336,7 +336,7 @@ void prism::bevel_utils::edge_based_bevel(const RowMatd& V, const RowMati& F, auto e1 = [&mark]() { // the dull one for (auto j = 0; j < 3; j++) if (mark(j) < 0) return j; - throw 1; + assert(false && "should not be here"); return -1; }(); auto v0 = Ff((e1 + 2) % 3); diff --git a/src/prism/cage_check.cpp b/src/prism/cage_check.cpp index 7e0996f..b542a1f 100755 --- a/src/prism/cage_check.cpp +++ b/src/prism/cage_check.cpp @@ -54,7 +54,7 @@ bool prism::cage_check::verify_edge_based_track( prism::local_validity::identify_zig(meta, f); if (rej_id == -10) { spdlog::error("Not allow two features in same triangle"); - throw 1; + throw std::runtime_error("Not allow two features in same triangle"); return false; } if (rej_id >= 0) { diff --git a/src/prism/cage_utils.cpp b/src/prism/cage_utils.cpp index 7a0e340..317dc22 100755 --- a/src/prism/cage_utils.cpp +++ b/src/prism/cage_utils.cpp @@ -730,7 +730,7 @@ std::map, int> prism::cage_utils::split_singular_edges( for (auto fi = 0; fi < F.size(); fi++) { if (is_singular(F[fi][0]) && is_singular(F[fi][1]) && is_singular(F[fi][2])) { spdlog::error("Fully Singular Face"); - throw 1; + throw std::runtime_error("Fully Singular Face"); } for (auto ei = 0; ei < 3; ei++) { auto v0 = F[fi][(ei)], v1 = F[fi][(ei + 1) % 3]; diff --git a/src/prism/local_operations/validity_checks.cpp b/src/prism/local_operations/validity_checks.cpp index d98e6d9..5b55b22 100755 --- a/src/prism/local_operations/validity_checks.cpp +++ b/src/prism/local_operations/validity_checks.cpp @@ -371,6 +371,10 @@ distort_check(const std::vector &base, const std::set &combined_trackee, // indices to ref.F tracked const RowMatd &refV, const RowMati &refF, double distortion_bound, int num_freeze, bool bundled_intersection) { + // ANCHOR: 80% bottleneck for no-curve pipeline. + // reduce the predicates would go a long way + // Possible improvements: 1. check walls, instead of cells and fill in topologically + // 2. Parallel // NormalCheck spdlog::trace("In NC ct#{}, tris{}", combined_trackee.size(), tris.size()); igl::Timer timer; @@ -378,7 +382,7 @@ distort_check(const std::vector &base, assert(base.size() == top.size()); assert(base.size() == mid.size()); std::vector> distributed_refs(tris.size()); - for (int i = 0; i < tris.size(); i++) { + for (int i = 0; i < tris.size(); i++) { // over all the prism cells. auto [v0, v1, v2] = tris[i]; std::array base_vert{base[v0], base[v1], base[v2]}; std::array mid_vert{mid[v0], mid[v1], mid[v2]}; diff --git a/src/prism/predicates/inside_octahedron.cpp b/src/prism/predicates/inside_octahedron.cpp index dd49fe4..b0b0d6e 100755 --- a/src/prism/predicates/inside_octahedron.cpp +++ b/src/prism/predicates/inside_octahedron.cpp @@ -234,11 +234,7 @@ bool prism::pointless_triangle_intersect_octahedron( bool point_inside = false; for (auto j = singular? 4: 0; j<12; j++) { auto&t = TWELVE_TETRAS[j]; - // if (GEO::PCK::orient_3d(vecprism[t[0]].data(), vecprism[t[1]].data(), vecprism[t[2]].data(), - // vecprism[t[3]].data()) > 0) - // spdlog::info(" positive"); - // else - // spdlog::info("not positive"); + if (prism::predicates::point_in_tetrahedron( vectriangle[i], vecprism[t[0]], vecprism[t[1]], vecprism[t[2]], vecprism[t[3]])) { diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt new file mode 100644 index 0000000..542e347 --- /dev/null +++ b/src/python/CMakeLists.txt @@ -0,0 +1,65 @@ +################################################################################ +# Find Python +################################################################################ + +if(NOT ${PYTHONINTERP_FOUND}) + execute_process( + COMMAND + python -c "import sys; sys.stdout.write(sys.version)" + OUTPUT_VARIABLE PYTHON_VERSION) + + if(NOT PYTHON_VERSION) + message(FATAL_ERROR "Unable to run python") + endif() + + string(REGEX MATCH "^3\.*" IS_PYTHON3 ${PYTHON_VERSION}) + + if(NOT IS_PYTHON3) + message(FATAL_ERROR "Unable to find python 3") + else() + set(PYTHON_EXECUTABLE "python") + endif() +endif() + +################################################################################ +# Pybind11 +################################################################################ + +if(UNIX) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fPIC") + if(NOT ${U_CMAKE_BUILD_TYPE} MATCHES DEBUG) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fvisibility=hidden -flto") + endif() +endif() + + +pybind11_add_module(prism_python + module.cpp + prism.cpp + spatial.cpp + curve.cpp +) + +add_library(prism::python ALIAS prism_python) + +target_link_libraries(prism_python + PUBLIC igl::core prism_library +) +target_link_libraries(prism_python + PUBLIC cumin_library +) + +# Generate position independent code +set_target_properties(prism_python PROPERTIES POSITION_INDEPENDENT_CODE ON) + +# Output location +string(TOUPPER "${CMAKE_BUILD_TYPE}" U_CMAKE_BUILD_TYPE) +if(${U_CMAKE_BUILD_TYPE} MATCHES RELEASE) + set_target_properties(prism_python PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/python) + set_target_properties(prism_python PROPERTIES OUTPUT_NAME "prism") +else() + set_target_properties(prism_python PROPERTIES OUTPUT_NAME "prism") + set_target_properties(prism_python PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/python/debug) +endif() +set_target_properties(prism_python PROPERTIES CXX_VISIBILITY_PRESET hidden) diff --git a/src/python/curve.cpp b/src/python/curve.cpp new file mode 100644 index 0000000..3d757a7 --- /dev/null +++ b/src/python/curve.cpp @@ -0,0 +1,100 @@ +#include +#include +#include + +#include +#include +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include + +#include +#include +#include +#include +#include +#include +namespace py = pybind11; + +using namespace pybind11::literals; +auto find_order_tet = [](int rows) { + int order = 0; + while ((order + 1) * (order + 2) * (order + 3) < 6 * rows) { + order++; + } + if ((order + 1) * (order + 2) * (order + 3) != 6 * rows) { + throw std::runtime_error( + fmt::format("cp is not a tetrahedron control points, order={}, cp={}", + order, rows)); + } + return order; +}; +auto find_order_tri = [](int rows) { + int order = 0; + while ((order + 1) * (order + 2) < 2 * rows) { + order++; + } + if ((order + 1) * (order + 2) != 2 * rows) { + throw std::runtime_error( + fmt::format("cp is not a tetrahedron control points, order={}, cp={}", + order, rows)); + } + return order; +}; +void python_export_curve(py::module &m) { + m.def( + "tetrahedron_inversion_check", + [](const RowMatd &cp) { + prism::curve::magic_matrices(find_order_tet(cp.rows()) - 1); + return prism::curve::tetrahedron_inversion_check(cp); + }, + "", "cp"_a); + + m.def("clear_curve_cache", []() { + spdlog::info("clear curve cache. For repeated experiment of different orders."); + prism::curve::HelperTensors::tensors_.reset(); + }); + m.def( + "elevated_positive_check", + [](const RowMatd &f_base, const RowMatd &f_top, const RowMatd &lagcp, + bool recurse_check) -> bool { + if (f_base.rows() != 3 || f_top.rows() != 3) { + throw std::runtime_error( + "elevated_positive_check: f_base and f_top must be 3x3 matrices"); + } + auto order = find_order_tri(lagcp.rows()) - 1; + spdlog::trace("triangle order {}", order); + auto helper = + prism::curve::magic_matrices(order, 3); + auto &tri15lag_from_tri10bern = helper.elev_lag_from_bern; + auto &dxyz = helper.volume_data.vec_dxyz; + auto tri4_cod = codecs_gen_id(helper.tri_order + 1, 2); + auto tet4_cod = codecs_gen_id(helper.tri_order + 1, 3); + // 6T {35 x 3D} + auto tens = prism::curve::surface_to_decomposed_tetra( + f_base, lagcp, f_top, f_base.row(0) == f_top.row(0), true, tri4_cod, + tet4_cod); + for (auto &d : dxyz) { + for (auto &t : tens) { + Eigen::Matrix3d j = d * t; + auto det = j.determinant(); + if (det <= 0) { + spdlog::debug("negative {}", det); + return false; + } + } + } + if (recurse_check) { + for (auto &t : tens) { + // this is codec_bc. + if (!prism::curve::tetrahedron_inversion_check(t)) { + spdlog::debug("blocked by recursive"); + return false; + } + } + } + return true; + }, + "f_base"_a, "f_top"_a, "lagrcp"_a, "recurse_check"_a = false, "comment"); +} diff --git a/src/python/module.cpp b/src/python/module.cpp new file mode 100644 index 0000000..d7c4b0f --- /dev/null +++ b/src/python/module.cpp @@ -0,0 +1,39 @@ +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace py = pybind11; + +using namespace pybind11::literals; +// using namespace prism; + +extern void python_export_curve(py::module&); +extern void python_export_spatial(py::module&); +extern void python_export_prism(py::module&); + +PYBIND11_MODULE(prism, m) { + m.doc() = R"docstring( + ================ + pyPrism: Python Binding for Bijective Projection Shell and Curved Meshing + Techniques used in "Bijective Projection in a Shell" and "Bijective and Coarse High Order Tetrahedral Meshes" + Zhongshi Jiang, 2021 + ================ + )docstring"; + + m.def("foo", + [](const Eigen::MatrixXd& V, const Eigen::MatrixXi& F) { return V; }, + "foo", "vertices"_a, "facets"_a); + m.def("loglevel", [](int level) { + spdlog::set_level(static_cast(level)); + }, "set log level", "level"_a); + + python_export_spatial(m); + python_export_prism(m); + python_export_curve(m); +} diff --git a/src/python/prism.cpp b/src/python/prism.cpp new file mode 100644 index 0000000..d224c25 --- /dev/null +++ b/src/python/prism.cpp @@ -0,0 +1,51 @@ +#include +#include +#include + +#include +#include +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include + +#include +#include + +namespace py = pybind11; + +using namespace pybind11::literals; + + +void python_export_prism(py::module &m) +{ + py::class_ cage(m, "PrismCage"); + cage.def(py::init()) + .def_property_readonly("base", [](const PrismCage &cage) { + RowMatd mb; + vec2eigen(cage.base, mb); + return mb; + }) + .def_property_readonly("mid", [](const PrismCage &cage) { + RowMatd mb; + vec2eigen(cage.mid, mb); + return mb; + }) + .def_property_readonly("top", [](const PrismCage &cage) { + RowMatd mb; + vec2eigen(cage.top, mb); + return mb; + }) + .def_property_readonly("F", [](const PrismCage &cage) { + RowMati mb; + vec2eigen(cage.F, mb); + return mb; }) + .def_property_readonly("refV", [](const PrismCage &cage) { return cage.ref.V; }) + .def_property_readonly("refF", [](const PrismCage &cage) { return cage.ref.F; }) + .def("transfer", [](const PrismCage &cage, const RowMatd &pxV, const RowMati &pxF, const RowMatd &queryP) { + Eigen::VectorXi queryF; + RowMatd queryUV; + prism::correspond_bc(cage, pxV, pxF, queryP, queryF, queryUV); + return std::tuple(queryF, queryUV); + }); +} diff --git a/src/python/spatial.cpp b/src/python/spatial.cpp new file mode 100644 index 0000000..0ca3171 --- /dev/null +++ b/src/python/spatial.cpp @@ -0,0 +1,110 @@ +#include +#include +#include + +#include +#include +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include +#include + +namespace py = pybind11; + +using namespace pybind11::literals; + +void python_export_spatial(py::module& m) { + // AABB + py::class_ > AABB2(m, "AABB2", "using igl::AABB for 2D case, not very robust."); + AABB2.def(py::init<>()) + .def(py::init&>()) + .def("init", + [](igl::AABB& tree, const Eigen::MatrixXd& V, + const Eigen::MatrixXi& Ele) { + return tree.init(V, Ele, + Eigen::Matrix(), + Eigen::Matrix(), + Eigen::VectorXi(), 0); + }) + .def("squared_distance", + [](const igl::AABB& tree, + const Eigen::MatrixXd& V, const Eigen::MatrixXi& Ele, + const Eigen::MatrixXd& P, Eigen::MatrixXd& sqrD, + Eigen::MatrixXi& I, Eigen::MatrixXd& C) { + return tree.squared_distance(V, Ele, P, sqrD, I, C); + }) + .def("find", + [](const igl::AABB& tree, + const Eigen::MatrixXd& V, const Eigen::MatrixXi& Ele, + const Eigen::RowVectorXd& q, bool first = true) { + auto f = tree.find(V, Ele, q, first); + if (f.size() == 0) + return -1; + else + return f[0]; + }) + .def("find_all", [](const igl::AABB& tree, + const Eigen::MatrixXd& V, const Eigen::MatrixXi& Ele, + const Eigen::MatrixXd& q) { + Eigen::VectorXi Fid(q.rows()); + Fid.setConstant(-1); + Eigen::MatrixXd BC = Eigen::MatrixXd::Zero(q.rows(), 3); + for (int i = 0; i < q.rows(); i++) { + auto f = tree.find(V, Ele, q.row(i), true); + if (f.size() != 0) { + Fid(i) = f[0]; + } + } + int n = q.rows(); + Eigen::MatrixXd A(n, 2), B(n, 2), C(n, 2); + auto F = Ele; + for (int i = 0; i < n; i++) { + auto f = Fid(i); + if (f != -1) { + A.row(i) = V.row(F(f, 0)); + B.row(i) = V.row(F(f, 1)); + C.row(i) = V.row(F(f, 2)); + } else { + A.row(i) = V.row(0); + B.row(i) = V.row(0); + C.row(i) = V.row(0); + } + } + igl::barycentric_coordinates(q, A, B, C, BC); + return std::make_tuple(Fid, BC); + }); + + py::class_ AABB(m, "AABB"); + AABB.def(py::init()) + .def("intersects_triangle", + [](const prism::geogram::AABB& self, const Vec3d& P0, const Vec3d& P1, + const Vec3d& P2) { + return self.intersects_triangle({P0, P1, P2}); + }) + .def("segment_query", + [](const prism::geogram::AABB& self, const Vec3d& P0, const Vec3d& P1) { + return self.segment_query(P0, P1); + }) + .def("segment_hit", [](const prism::geogram::AABB& self, const Vec3d& P0, + const Vec3d& P1, bool ray = false) { + prism::Hit hit; + bool result = self.segment_hit(P0, P1, hit); + if (result) + return std::tuple(hit.id, hit.u, hit.v, hit.t); + else + return std::tuple(-1, 0., 0., -1.); + }); + + m.def("self_intersect", [] + ( const Eigen::MatrixXd & V, + const Eigen::MatrixXi& F)->bool { + std::vector vecV; + std::vector vecF; + eigen2vec(V, vecV); + eigen2vec(V, vecV); + auto pairs = prism::spatial_hash::self_intersections(vecV,vecF); + return pairs.size() > 0; + }, "use spatial hash for self intersection check", "V"_a, "F"_a); +}