diff --git a/src/librt/CMakeLists.txt b/src/librt/CMakeLists.txt index b4b948298a1..79746d29768 100644 --- a/src/librt/CMakeLists.txt +++ b/src/librt/CMakeLists.txt @@ -213,6 +213,8 @@ set(LIBRT_SOURCES primitives/tgc/tgc.c primitives/tgc/tgc_brep.cpp primitives/tgc/tgc_mirror.c + primitives/tgc/tgc_tess.c + primitives/tgc/tgc_tnurb.c primitives/tor/tor.c primitives/tor/tor_brep.cpp primitives/tor/tor_mirror.c diff --git a/src/librt/primitives/tgc/tgc.c b/src/librt/primitives/tgc/tgc.c index 4b86dc13851..e358fa252d7 100644 --- a/src/librt/primitives/tgc/tgc.c +++ b/src/librt/primitives/tgc/tgc.c @@ -69,15 +69,11 @@ extern int rt_rec_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_ static void rt_tgc_rotate(fastf_t *A, fastf_t *B, fastf_t *Hv, fastf_t *Rot, fastf_t *Inv, struct tgc_specific *tgc); static void rt_tgc_shear(const fastf_t *vect, int axis, fastf_t *Shr, fastf_t *Trn, fastf_t *Inv); static void rt_tgc_scale(fastf_t a, fastf_t b, fastf_t h, fastf_t *Scl, fastf_t *Inv); -static void nmg_tgc_disk(struct faceuse *fu, fastf_t *rmat, fastf_t height, int flip); -static void nmg_tgc_nurb_cyl(struct faceuse *fu, fastf_t *top_mat, fastf_t *bot_mat); void rt_pnt_sort(register fastf_t *t, int npts); #define VLARGE 1000000.0 -#define MAX_RATIO 10.0 /* maximum allowed height-to-width ration for triangles */ -#define RAT M_SQRT1_2 #define T_OUT 0 #define T_IN 1 @@ -120,33 +116,6 @@ const struct bu_structparse rt_tgc_parse[] = { { {'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL } }; - -static const fastf_t nmg_tgc_unitcircle[36] = { - 1.0, 0.0, 0.0, 1.0, - RAT, -RAT, 0.0, RAT, - 0.0, -1.0, 0.0, 1.0, - -RAT, -RAT, 0.0, RAT, - -1.0, 0.0, 0.0, 1.0, - -RAT, RAT, 0.0, RAT, - 0.0, 1.0, 0.0, 1.0, - RAT, RAT, 0.0, RAT, - 1.0, 0.0, 0.0, 1.0 -}; - - -static const fastf_t nmg_uv_unitcircle[27] = { - 1.0, .5, 1.0, - RAT, RAT, RAT, - .5, 1.0, 1.0, - 0.0, RAT, RAT, - 0.0, .5, 1.0, - 0.0, 0.0, RAT, - .5, 0.0, 1.0, - RAT, 0.0, RAT, - 1.0, .5, 1.0 -}; - - #ifdef USE_OPENCL /* largest data members first */ struct clt_tgc_specific { @@ -2207,1104 +2176,6 @@ rt_tgc_curve(register struct curvature *cvp, register struct hit *hitp, struct s } -/** - * Tessellation of the TGC. - * - * Returns - - * -1 failure - * 0 OK. *r points to nmgregion that holds this tessellation. - */ - -struct tgc_pts -{ - point_t pt; - vect_t tan_axb; - struct vertex *v; - char dont_use; -}; - - -/* version using tolerances */ -int -rt_tgc_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bg_tess_tol *ttol, const struct bn_tol *tol) -{ - struct shell *s; /* shell to hold facetted TGC */ - struct faceuse *fu, *fu_top, *fu_base; - struct rt_tgc_internal *tip; - fastf_t radius; /* bounding sphere radius */ - fastf_t max_radius; /* max of a, b, c, d */ - fastf_t h, a, b, c, d; /* lengths of TGC vectors */ - fastf_t inv_length; /* 1.0/length of a vector */ - vect_t unit_a, unit_b, unit_c, unit_d; /* units vectors in a, b, c, d directions */ - fastf_t rel, absolute, norm; /* interpreted tolerances */ - fastf_t alpha_tol; /* final tolerance for ellipse parameter */ - fastf_t abs_tol; /* handle invalid ttol->abs */ - size_t nells; /* total number of ellipses */ - size_t nsegs; /* number of vertices/ellipse */ - vect_t *A; /* array of A vectors for ellipses */ - vect_t *B; /* array of B vectors for ellipses */ - fastf_t *factors; /* array of ellipse locations along height vector */ - vect_t vtmp; - vect_t normal; /* normal vector */ - vect_t rev_norm; /* reverse normal */ - struct tgc_pts **pts; /* array of points (pts[ellipse#][seg#]) */ - struct bu_ptbl verts; /* table of vertices used for top and bottom faces */ - struct bu_ptbl faces; /* table of faceuses for nmg_gluefaces */ - struct vertex **v[3]; /* array for making triangular faces */ - - size_t i; - - VSETALL(unit_a, 0); - VSETALL(unit_b, 0); - VSETALL(unit_c, 0); - VSETALL(unit_d, 0); - - RT_CK_DB_INTERNAL(ip); - tip = (struct rt_tgc_internal *)ip->idb_ptr; - RT_TGC_CK_MAGIC(tip); - - if (ttol->abs > 0.0 && ttol->abs < tol->dist) { - bu_log("WARNING: tessellation tolerance is %fmm while calculational tolerance is %fmm\n", - ttol->abs, tol->dist); - bu_log("Cannot tessellate a TGC to finer tolerance than the calculational tolerance\n"); - abs_tol = tol->dist; - } else { - abs_tol = ttol->abs; - } - - h = MAGNITUDE(tip->h); - a = MAGNITUDE(tip->a); - if (2.0*a <= tol->dist) - a = 0.0; - b = MAGNITUDE(tip->b); - if (2.0*b <= tol->dist) - b = 0.0; - c = MAGNITUDE(tip->c); - if (2.0*c <= tol->dist) - c = 0.0; - d = MAGNITUDE(tip->d); - if (2.0*d <= tol->dist) - d = 0.0; - - if (ZERO(a) && ZERO(b) && (ZERO(c) || ZERO(d))) { - bu_log("Illegal TGC a, b, and c or d less than tolerance\n"); - return -1; - } else if (ZERO(c) && ZERO(d) && (ZERO(a) || ZERO(b))) { - bu_log("Illegal TGC c, d, and a or b less than tolerance\n"); - return -1; - } - - if (a > 0.0) { - inv_length = 1.0/a; - VSCALE(unit_a, tip->a, inv_length); - } - if (b > 0.0) { - inv_length = 1.0/b; - VSCALE(unit_b, tip->b, inv_length); - } - if (c > 0.0) { - inv_length = 1.0/c; - VSCALE(unit_c, tip->c, inv_length); - } - if (d > 0.0) { - inv_length = 1.0/d; - VSCALE(unit_d, tip->d, inv_length); - } - - /* get bounding sphere radius for relative tolerance */ - radius = h/2.0; - max_radius = 0.0; - if (a > max_radius) - max_radius = a; - if (b > max_radius) - max_radius = b; - if (c > max_radius) - max_radius = c; - if (d > max_radius) - max_radius = d; - - if (max_radius > radius) - radius = max_radius; - - if (abs_tol <= 0.0 && ttol->rel <= 0.0 && ttol->norm <= 0.0) { - /* no tolerances specified, use 10% relative tolerance */ - if ((radius * 0.2) < max_radius) - alpha_tol = 2.0 * acos(1.0 - 2.0 * radius * 0.1 / max_radius); - else - alpha_tol = M_PI_2; - } else { - if (abs_tol > 0.0) - absolute = 2.0 * acos(1.0 - abs_tol/max_radius); - else - absolute = M_PI_2; - - if (ttol->rel > 0.0) { - if (ttol->rel * 2.0 * radius < max_radius) - rel = 2.0 * acos(1.0 - ttol->rel * 2.0 * radius/max_radius); - else - rel = M_PI_2; - } else - rel = M_PI_2; - - if (ttol->norm > 0.0) { - fastf_t norm_top, norm_bot; - - if (anorm) * (a/b)); - else - norm_bot = 2.0 * atan(tan(ttol->norm) * (b/a)); - - if (cnorm) * (c/d)); - else - norm_top = 2.0 * atan(tan(ttol->norm) * (d/c)); - - if (norm_bot < norm_top) - norm = norm_bot; - else - norm = norm_top; - } else - norm = M_PI_2; - - if (absolute < rel) - alpha_tol = absolute; - else - alpha_tol = rel; - if (norm < alpha_tol) - alpha_tol = norm; - } - - /* get number of segments per quadrant */ - nsegs = (int)(M_PI_2 / alpha_tol + 0.9999); - if (nsegs < 2) - nsegs = 2; - - /* and for complete ellipse */ - nsegs *= 4; - - /* get number and placement of intermediate ellipses */ - { - fastf_t ratios[4], max_ratio; - fastf_t new_ratio = 0; - int which_ratio; - fastf_t len_ha, len_hb; - vect_t ha, hb; - fastf_t ang; - fastf_t sin_ang, cos_ang, cos_m_1_sq, sin_sq; - fastf_t len_A, len_B, len_C, len_D; - size_t bot_ell=0; - size_t top_ell=1; - int reversed=0; - - nells = 2; - - max_ratio = MAX_RATIO + 1.0; - - factors = (fastf_t *)bu_malloc(nells*sizeof(fastf_t), "factors"); - A = (vect_t *)bu_malloc(nells*sizeof(vect_t), "A vectors"); - B = (vect_t *)bu_malloc(nells*sizeof(vect_t), "B vectors"); - - factors[bot_ell] = 0.0; - factors[top_ell] = 1.0; - VMOVE(A[bot_ell], tip->a); - VMOVE(A[top_ell], tip->c); - VMOVE(B[bot_ell], tip->b); - VMOVE(B[top_ell], tip->d); - - /* make sure that AxB points in the general direction of H */ - VCROSS(vtmp, A[0], B[0]); - if (VDOT(vtmp, tip->h) < 0.0) { - VMOVE(A[bot_ell], tip->b); - VMOVE(A[top_ell], tip->d); - VMOVE(B[bot_ell], tip->a); - VMOVE(B[top_ell], tip->c); - reversed = 1; - } - ang = M_2PI/((double)nsegs); - sin_ang = sin(ang); - cos_ang = cos(ang); - cos_m_1_sq = (cos_ang - 1.0)*(cos_ang - 1.0); - sin_sq = sin_ang*sin_ang; - - VJOIN2(ha, tip->h, 1.0, tip->c, -1.0, tip->a); - VJOIN2(hb, tip->h, 1.0, tip->d, -1.0, tip->b); - len_ha = MAGNITUDE(ha); - len_hb = MAGNITUDE(hb); - - while (max_ratio > MAX_RATIO) { - fastf_t tri_width; - - len_A = MAGNITUDE(A[bot_ell]); - if (2.0*len_A <= tol->dist) - len_A = 0.0; - len_B = MAGNITUDE(B[bot_ell]); - if (2.0*len_B <= tol->dist) - len_B = 0.0; - len_C = MAGNITUDE(A[top_ell]); - if (2.0*len_C <= tol->dist) - len_C = 0.0; - len_D = MAGNITUDE(B[top_ell]); - if (2.0*len_D <= tol->dist) - len_D = 0.0; - - if ((len_B > 0.0 && len_D > 0.0) || - (len_B > 0.0 && (ZERO(len_D) && ZERO(len_C)))) - { - tri_width = sqrt(cos_m_1_sq*len_A*len_A + sin_sq*len_B*len_B); - ratios[0] = (factors[top_ell] - factors[bot_ell])*len_ha - /tri_width; - } else - ratios[0] = 0.0; - - if ((len_A > 0.0 && len_C > 0.0) || - (len_A > 0.0 && (ZERO(len_C) && ZERO(len_D)))) - { - tri_width = sqrt(sin_sq*len_A*len_A + cos_m_1_sq*len_B*len_B); - ratios[1] = (factors[top_ell] - factors[bot_ell])*len_hb - /tri_width; - } else - ratios[1] = 0.0; - - if ((len_D > 0.0 && len_B > 0.0) || - (len_D > 0.0 && (ZERO(len_A) && ZERO(len_B)))) - { - tri_width = sqrt(cos_m_1_sq*len_C*len_C + sin_sq*len_D*len_D); - ratios[2] = (factors[top_ell] - factors[bot_ell])*len_ha - /tri_width; - } else - ratios[2] = 0.0; - - if ((len_C > 0.0 && len_A > 0.0) || - (len_C > 0.0 && (ZERO(len_A) && ZERO(len_B)))) - { - tri_width = sqrt(sin_sq*len_C*len_C + cos_m_1_sq*len_D*len_D); - ratios[3] = (factors[top_ell] - factors[bot_ell])*len_hb - /tri_width; - } else - ratios[3] = 0.0; - - which_ratio = -1; - max_ratio = 0.0; - - for (i=0; i<4; i++) { - if (ratios[i] > max_ratio) { - max_ratio = ratios[i]; - which_ratio = i; - } - } - - if (ZERO(len_A) && ZERO(len_B) && ZERO(len_C) && ZERO(len_D)) { - if (top_ell == nells - 1) { - VMOVE(A[top_ell-1], A[top_ell]); - VMOVE(B[top_ell-1], A[top_ell]); - factors[top_ell-1] = factors[top_ell]; - } else if (bot_ell == 0) { - for (i=0; i 0.5) - new_ratio = 0.5; - } else if (which_ratio == 2 || which_ratio == 3) { - new_ratio = 1.0 - MAX_RATIO/max_ratio; - if (top_ell == nells - 1 && new_ratio < 0.5) - new_ratio = 0.5; - } else { - /* no MAX? */ - bu_bomb("rt_tgc_tess: Should never get here!!\n"); - } - - nells++; - factors = (fastf_t *)bu_realloc(factors, nells*sizeof(fastf_t), "factors"); - A = (vect_t *)bu_realloc(A, nells*sizeof(vect_t), "A vectors"); - B = (vect_t *)bu_realloc(B, nells*sizeof(vect_t), "B vectors"); - - for (i=nells-1; i>top_ell; i--) { - factors[i] = factors[i-1]; - VMOVE(A[i], A[i-1]); - VMOVE(B[i], B[i-1]); - } - - factors[top_ell] = factors[bot_ell] + - new_ratio*(factors[top_ell+1] - factors[bot_ell]); - - if (reversed) { - VBLEND2(A[top_ell], (1.0-factors[top_ell]), tip->b, factors[top_ell], tip->d); - VBLEND2(B[top_ell], (1.0-factors[top_ell]), tip->a, factors[top_ell], tip->c); - } else { - VBLEND2(A[top_ell], (1.0-factors[top_ell]), tip->a, factors[top_ell], tip->c); - VBLEND2(B[top_ell], (1.0-factors[top_ell]), tip->b, factors[top_ell], tip->d); - } - - if (which_ratio == 0 || which_ratio == 1) { - top_ell++; - bot_ell++; - } - - } - - } - - /* get memory for points */ - pts = (struct tgc_pts **)bu_calloc(nells, sizeof(struct tgc_pts *), "rt_tgc_tess: pts"); - for (i=0; iv); - } else if (i == nells-1 && ZERO(c) && ZERO(d)) { - VADD2(pts[i][j].pt, tip->v, tip->h); - } else { - VJOIN3(pts[i][j].pt, tip->v, h_factor, tip->h, cos_alpha, A[i], sin_alpha, B[i]); - } - - /* Storing the tangent here while sines and cosines are available */ - if (i == 0 && ZERO(a) && ZERO(b)) { - VCOMB2(pts[0][j].tan_axb, -sin_alpha, unit_c, cos_alpha, unit_d); - } else if (i == nells-1 && ZERO(c) && ZERO(d)) { - VCOMB2(pts[i][j].tan_axb, -sin_alpha, unit_a, cos_alpha, unit_b); - } else { - VCOMB2(pts[i][j].tan_axb, -sin_alpha, A[i], cos_alpha, B[i]); - } - } - } - - /* make sure no edges will be created with length < tol->dist */ - for (i=0; i tol->dist_sq) { - VMOVE(curr_pt, pts[i][j].pt); - } else { - /* don't use this point, it will create a too short edge */ - pts[i][j].dont_use = 'n'; - } - } - } - - /* make region, shell, vertex */ - *r = nmg_mrsv(m); - s = BU_LIST_FIRST(shell, &(*r)->s_hd); - - bu_ptbl_init(&verts, 64, " &verts "); - bu_ptbl_init(&faces, 64, " &faces "); - /* Make bottom face */ - if (a > 0.0 && b > 0.0) { - for (i=nsegs; i>0; i--) { - /* reverse order to get outward normal */ - if (!pts[0][i-1].dont_use) - bu_ptbl_ins(&verts, (long *)&pts[0][i-1].v); - } - - if (BU_PTBL_LEN(&verts) > 2) { - fu_base = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&verts), BU_PTBL_LEN(&verts)); - bu_ptbl_ins(&faces, (long *)fu_base); - } else - fu_base = (struct faceuse *)NULL; - } else - fu_base = (struct faceuse *)NULL; - - - /* Make top face */ - if (c > 0.0 && d > 0.0) { - bu_ptbl_reset(&verts); - for (i=0; i 2) { - fu_top = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&verts), BU_PTBL_LEN(&verts)); - bu_ptbl_ins(&faces, (long *)fu_top); - } else - fu_top = (struct faceuse *)NULL; - } else - fu_top = (struct faceuse *)NULL; - - /* Free table of vertices */ - bu_ptbl_free(&verts); - - /* Make triangular faces */ - for (i=0; i 0.0 || b > 0.0) { - if (!pts[i][k].dont_use) { - v[0] = curr_bot; - v[1] = &pts[i][k].v; - if (i+1 == nells-1 && ZERO(c) && ZERO(d)) - v[2] = &pts[i+1][0].v; - else - v[2] = curr_top; - fu = nmg_cmface(s, v, 3); - bu_ptbl_ins(&faces, (long *)fu); - curr_bot = &pts[i][k].v; - } - } - - if (i != nells-2 || c > 0.0 || d > 0.0) { - if (!pts[i+1][k].dont_use) { - v[0] = &pts[i+1][k].v; - v[1] = curr_top; - if (i == 0 && ZERO(a) && ZERO(b)) - v[2] = &pts[i][0].v; - else - v[2] = curr_bot; - fu = nmg_cmface(s, v, 3); - bu_ptbl_ins(&faces, (long *)fu); - curr_top = &pts[i+1][k].v; - } - } - } - } - - /* Assign geometry */ - for (i=0; iv); - } else if (i == nells-1 && ZERO(c) && ZERO(d)) { - if (j == 0) { - VADD2(pt_geom, tip->v, tip->h); - nmg_vertex_gv(pts[i][0].v, pt_geom); - } - } else if (pts[i][j].v) - nmg_vertex_gv(pts[i][j].v, pts[i][j].pt); - - /* Storing the tangent here while sines and cosines are available */ - if (i == 0 && ZERO(a) && ZERO(b)) { - VCOMB2(pts[0][j].tan_axb, -sin_alpha, unit_c, cos_alpha, unit_d); - } else if (i == nells-1 && ZERO(c) && ZERO(d)) { - VCOMB2(pts[i][j].tan_axb, -sin_alpha, unit_a, cos_alpha, unit_b); - } else { - VCOMB2(pts[i][j].tan_axb, -sin_alpha, A[i], cos_alpha, B[i]); - } - } - } - - /* Associate face plane equations */ - for (i=0; ivu_hd)) { - NMG_CK_VERTEXUSE(vu); - - fu = nmg_find_fu_of_vu(vu); - NMG_CK_FACEUSE(fu); - - /* don't need vertexuse normals for faces that are really flat */ - if (fu == fu_base || fu->fumate_p == fu_base || - fu == fu_top || fu->fumate_p == fu_top) - continue; - - if (fu->orientation == OT_SAME) - nmg_vertexuse_nv(vu, normal); - else if (fu->orientation == OT_OPPOSITE) - nmg_vertexuse_nv(vu, rev_norm); - } - } - } - } - - /* Finished with storage, so free it */ - bu_free((char *)factors, "rt_tgc_tess: factors"); - bu_free((char *)A, "rt_tgc_tess: A"); - bu_free((char *)B, "rt_tgc_tess: B"); - for (i=0; ilu_hd)) { - struct edgeuse *eu; - - NMG_CK_LOOPUSE(lu); - - if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC) - continue; - - for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) { - struct edge *e; - - NMG_CK_EDGEUSE(eu); - e = eu->e_p; - NMG_CK_EDGE(e); - e->is_real = 1; - } - } - } - - nmg_region_a(*r, tol); - - /* glue faces together */ - nmg_gluefaces((struct faceuse **)BU_PTBL_BASEADDR(&faces), BU_PTBL_LEN(&faces), &RTG.rtg_vlfree, tol); - bu_ptbl_free(&faces); - - return 0; -} - - -/** - * "Tessellate an TGC into a trimmed-NURB-NMG data structure. - * Computing NURB surfaces and trimming curves to interpolate - * the parameters of the TGC - * - * The process is to create the nmg topology of the TGC fill it - * in with a unit cylinder geometry (i.e. unitcircle at the top (0, 0, 1) - * unit cylinder of radius 1, and unitcircle at the bottom), and then - * scale it with a perspective matrix derived from the parameters of the - * tgc. The result is three trimmed nurb surfaces which interpolate the - * parameters of the original TGC. - * - * Returns - - * -1 failure - * 0 OK. *r points to nmgregion that holds this tessellation - */ - -int -rt_tgc_tnurb(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bn_tol *tol) -{ - struct rt_tgc_internal *tip; - - struct shell *s; - struct vertex *verts[2]; - struct vertex **vertp[4]; - struct faceuse * top_fu; - struct faceuse * cyl_fu; - struct faceuse * bot_fu; - vect_t uvw; - struct loopuse *lu; - struct edgeuse *eu; - struct edgeuse *top_eu; - struct edgeuse *bot_eu; - - mat_t mat; - mat_t imat, omat, top_mat, bot_mat; - vect_t anorm; - vect_t bnorm; - vect_t cnorm; - - - /* Get the internal representation of the tgc */ - - RT_CK_DB_INTERNAL(ip); - tip = (struct rt_tgc_internal *) ip->idb_ptr; - RT_TGC_CK_MAGIC(tip); - - /* Create the NMG Topology */ - - *r = nmg_mrsv(m); /* Make region, empty shell, vertex */ - s = BU_LIST_FIRST(shell, &(*r)->s_hd); - - - /* Create transformation matrix for the top cap surface*/ - - MAT_IDN(omat); - MAT_IDN(mat); - - omat[0] = MAGNITUDE(tip->c); - omat[5] = MAGNITUDE(tip->d); - omat[3] = tip->v[0] + tip->h[0]; - omat[7] = tip->v[1] + tip->h[1]; - omat[11] = tip->v[2] + tip->h[2]; - - bn_mat_mul(imat, mat, omat); - - VMOVE(anorm, tip->c); - VMOVE(bnorm, tip->d); - VCROSS(cnorm, tip->c, tip->d); - VUNITIZE(anorm); - VUNITIZE(bnorm); - VUNITIZE(cnorm); - - MAT_IDN(omat); - - VMOVE(&omat[0], anorm); - VMOVE(&omat[4], bnorm); - VMOVE(&omat[8], cnorm); - - - bn_mat_mul(top_mat, omat, imat); - - /* Create topology for top cap surface */ - - verts[0] = verts[1] = NULL; - vertp[0] = &verts[0]; - top_fu = nmg_cmface(s, vertp, 1); - - lu = BU_LIST_FIRST(loopuse, &top_fu->lu_hd); - NMG_CK_LOOPUSE(lu); - eu= BU_LIST_FIRST(edgeuse, &lu->down_hd); - NMG_CK_EDGEUSE(eu); - top_eu = eu; - - VSET(uvw, 0, 0, 0); - nmg_vertexuse_a_cnurb(eu->vu_p, uvw); - VSET(uvw, 1, 0, 0); - nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); - - /* Top cap surface */ - nmg_tgc_disk(top_fu, top_mat, 0.0, 0); - - /* Create transformation matrix for the bottom cap surface*/ - - MAT_IDN(omat); - MAT_IDN(mat); - - omat[0] = MAGNITUDE(tip->a); - omat[5] = MAGNITUDE(tip->b); - omat[3] = tip->v[0]; - omat[7] = tip->v[1]; - omat[11] = tip->v[2]; - - bn_mat_mul(imat, mat, omat); - - VMOVE(anorm, tip->a); - VMOVE(bnorm, tip->b); - VCROSS(cnorm, tip->a, tip->b); - VUNITIZE(anorm); - VUNITIZE(bnorm); - VUNITIZE(cnorm); - - MAT_IDN(omat); - - VMOVE(&omat[0], anorm); - VMOVE(&omat[4], bnorm); - VMOVE(&omat[8], cnorm); - - bn_mat_mul(bot_mat, omat, imat); - - /* Create topology for bottom cap surface */ - - vertp[0] = &verts[1]; - bot_fu = nmg_cmface(s, vertp, 1); - - lu = BU_LIST_FIRST(loopuse, &bot_fu->lu_hd); - NMG_CK_LOOPUSE(lu); - eu= BU_LIST_FIRST(edgeuse, &lu->down_hd); - NMG_CK_EDGEUSE(eu); - bot_eu = eu; - - VSET(uvw, 0, 0, 0); - nmg_vertexuse_a_cnurb(eu->vu_p, uvw); - VSET(uvw, 1, 0, 0); - nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); - - - nmg_tgc_disk(bot_fu, bot_mat, 0.0, 1); - - /* Create topology for cylinder surface */ - - vertp[0] = &verts[0]; - vertp[1] = &verts[0]; - vertp[2] = &verts[1]; - vertp[3] = &verts[1]; - cyl_fu = nmg_cmface(s, vertp, 4); - - nmg_tgc_nurb_cyl(cyl_fu, top_mat, bot_mat); - - /* fuse top cylinder edge to matching edge on body of cylinder */ - lu = BU_LIST_FIRST(loopuse, &cyl_fu->lu_hd); - - eu= BU_LIST_FIRST(edgeuse, &lu->down_hd); - - nmg_je(top_eu, eu); - - eu= BU_LIST_LAST(edgeuse, &lu->down_hd); - eu= BU_LIST_LAST(edgeuse, &eu->l); - - nmg_je(bot_eu, eu); - nmg_region_a(*r, tol); - - return 0; -} - - -static void -nmg_tgc_disk(struct faceuse *fu, fastf_t *rmat, fastf_t height, int flip) -{ - struct face_g_snurb * fg; - struct loopuse * lu; - struct edgeuse * eu; - struct edge_g_cnurb * eg; - fastf_t *mptr; - int i; - vect_t vect; - point_t point; - - nmg_face_g_snurb(fu, - 2, 2, /* u, v order */ - 4, 4, /* number of knots */ - NULL, NULL, /* initial knot vectors */ - 2, 2, /* n_rows, n_cols */ - RT_NURB_MAKE_PT_TYPE(3, 2, 0), - NULL); /* Initial mesh */ - - fg = fu->f_p->g.snurb_p; - - - fg->u.knots[0] = 0.0; - fg->u.knots[1] = 0.0; - fg->u.knots[2] = 1.0; - fg->u.knots[3] = 1.0; - - fg->v.knots[0] = 0.0; - fg->v.knots[1] = 0.0; - fg->v.knots[2] = 1.0; - fg->v.knots[3] = 1.0; - - if (flip) { - fg->ctl_points[0] = 1.; - fg->ctl_points[1] = -1.; - fg->ctl_points[2] = height; - - fg->ctl_points[3] = -1; - fg->ctl_points[4] = -1.; - fg->ctl_points[5] = height; - - fg->ctl_points[6] = 1.; - fg->ctl_points[7] = 1.; - fg->ctl_points[8] = height; - - fg->ctl_points[9] = -1.; - fg->ctl_points[10] = 1.; - fg->ctl_points[11] = height; - } else { - - fg->ctl_points[0] = -1.; - fg->ctl_points[1] = -1.; - fg->ctl_points[2] = height; - - fg->ctl_points[3] = 1; - fg->ctl_points[4] = -1.; - fg->ctl_points[5] = height; - - fg->ctl_points[6] = -1.; - fg->ctl_points[7] = 1.; - fg->ctl_points[8] = height; - - fg->ctl_points[9] = 1.; - fg->ctl_points[10] = 1.; - fg->ctl_points[11] = height; - } - - /* multiple the matrix to get orientation and scaling that we want */ - mptr = fg->ctl_points; - - i = fg->s_size[0] * fg->s_size[1]; - - for (; i> 0; i--) { - MAT4X3PNT(vect, rmat, mptr); - mptr[0] = vect[0]; - mptr[1] = vect[1]; - mptr[2] = vect[2]; - mptr += 3; - } - - lu = BU_LIST_FIRST(loopuse, &fu->lu_hd); - NMG_CK_LOOPUSE(lu); - eu= BU_LIST_FIRST(edgeuse, &lu->down_hd); - NMG_CK_EDGEUSE(eu); - - - if (!flip) { - nmg_nurb_s_eval(fu->f_p->g.snurb_p, - nmg_uv_unitcircle[0], nmg_uv_unitcircle[1], point); - nmg_vertex_gv(eu->vu_p->v_p, point); - } else { - nmg_nurb_s_eval(fu->f_p->g.snurb_p, - nmg_uv_unitcircle[12], nmg_uv_unitcircle[13], point); - nmg_vertex_gv(eu->vu_p->v_p, point); - } - - nmg_edge_g_cnurb(eu, 3, 12, NULL, 9, RT_NURB_MAKE_PT_TYPE(3, 3, 1), - NULL); - - eg = eu->g.cnurb_p; - eg->order = 3; - - eg->k.knots[0] = 0.0; - eg->k.knots[1] = 0.0; - eg->k.knots[2] = 0.0; - eg->k.knots[3] = .25; - eg->k.knots[4] = .25; - eg->k.knots[5] = .5; - eg->k.knots[6] = .5; - eg->k.knots[7] = .75; - eg->k.knots[8] = .75; - eg->k.knots[9] = 1.0; - eg->k.knots[10] = 1.0; - eg->k.knots[11] = 1.0; - - if (!flip) { - for (i = 0; i < 27; i++) - eg->ctl_points[i] = nmg_uv_unitcircle[i]; - } else { - - VSET(&eg->ctl_points[0], 0.0, .5, 1.0); - VSET(&eg->ctl_points[3], 0.0, 0.0, RAT); - VSET(&eg->ctl_points[6], 0.5, 0.0, 1.0); - VSET(&eg->ctl_points[9], RAT, 0.0, RAT); - VSET(&eg->ctl_points[12], 1.0, .5, 1.0); - VSET(&eg->ctl_points[15], RAT, RAT, RAT); - VSET(&eg->ctl_points[18], .5, 1.0, 1.0); - VSET(&eg->ctl_points[21], 0.0, RAT, RAT); - VSET(&eg->ctl_points[24], 0.0, .5, 1.0); - } -} - - -/* Create a cylinder with a top surface and a bottom surface - * defined by the ellipsoids at the top and bottom of the - * cylinder, the top_mat, and bot_mat are applied to a unit circle - * for the top row of the surface and the bot row of the surface - * respectively. - */ -static void -nmg_tgc_nurb_cyl(struct faceuse *fu, fastf_t *top_mat, fastf_t *bot_mat) -{ - struct face_g_snurb * fg; - struct loopuse * lu; - struct edgeuse * eu; - fastf_t * mptr; - int i; - hvect_t vect; - point_t uvw; - point_t point; - hvect_t hvect; - - nmg_face_g_snurb(fu, - 3, 2, - 12, 4, - NULL, NULL, - 2, 9, - RT_NURB_MAKE_PT_TYPE(4, 3, 1), - NULL); - - fg = fu->f_p->g.snurb_p; - - fg->v.knots[0] = 0.0; - fg->v.knots[1] = 0.0; - fg->v.knots[2] = 1.0; - fg->v.knots[3] = 1.0; - - fg->u.knots[0] = 0.0; - fg->u.knots[1] = 0.0; - fg->u.knots[2] = 0.0; - fg->u.knots[3] = .25; - fg->u.knots[4] = .25; - fg->u.knots[5] = .5; - fg->u.knots[6] = .5; - fg->u.knots[7] = .75; - fg->u.knots[8] = .75; - fg->u.knots[9] = 1.0; - fg->u.knots[10] = 1.0; - fg->u.knots[11] = 1.0; - - mptr = fg->ctl_points; - - for (i = 0; i < 9; i++) { - MAT4X4PNT(vect, top_mat, &nmg_tgc_unitcircle[i*4]); - mptr[0] = vect[0]; - mptr[1] = vect[1]; - mptr[2] = vect[2]; - mptr[3] = vect[3]; - mptr += 4; - } - - for (i = 0; i < 9; i++) { - MAT4X4PNT(vect, bot_mat, &nmg_tgc_unitcircle[i*4]); - mptr[0] = vect[0]; - mptr[1] = vect[1]; - mptr[2] = vect[2]; - mptr[3] = vect[3]; - mptr += 4; - } - - /* Assign edgeuse parameters & vertex geometry */ - - lu = BU_LIST_FIRST(loopuse, &fu->lu_hd); - NMG_CK_LOOPUSE(lu); - eu = BU_LIST_FIRST(edgeuse, &lu->down_hd); - NMG_CK_EDGEUSE(eu); - - /* March around the fu's loop assigning uv parameter values */ - - nmg_nurb_s_eval(fg, 0.0, 0.0, hvect); - HDIVIDE(point, hvect); - nmg_vertex_gv(eu->vu_p->v_p, point); /* 0, 0 vertex */ - - VSET(uvw, 0, 0, 0); - nmg_vertexuse_a_cnurb(eu->vu_p, uvw); - VSET(uvw, 1, 0, 0); - nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); - eu = BU_LIST_NEXT(edgeuse, &eu->l); - - VSET(uvw, 1, 0, 0); - nmg_vertexuse_a_cnurb(eu->vu_p, uvw); - VSET(uvw, 1, 1, 0); - nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); - eu = BU_LIST_NEXT(edgeuse, &eu->l); - - VSET(uvw, 1, 1, 0); - nmg_vertexuse_a_cnurb(eu->vu_p, uvw); - VSET(uvw, 0, 1, 0); - nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); - - nmg_nurb_s_eval(fg, 1., 1., hvect); - HDIVIDE(point, hvect); - nmg_vertex_gv(eu->vu_p->v_p, point); /* 4, 1 vertex */ - - eu = BU_LIST_NEXT(edgeuse, &eu->l); - - VSET(uvw, 0, 1, 0); - nmg_vertexuse_a_cnurb(eu->vu_p, uvw); - VSET(uvw, 0, 0, 0); - nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); - - /* Create the edge loop geometry */ - - lu = BU_LIST_FIRST(loopuse, &fu->lu_hd); - NMG_CK_LOOPUSE(lu); - eu = BU_LIST_FIRST(edgeuse, &lu->down_hd); - NMG_CK_EDGEUSE(eu); - - for (i = 0; i < 4; i++) { - nmg_edge_g_cnurb_plinear(eu); - eu = BU_LIST_NEXT(edgeuse, &eu->l); - } -} - - int rt_tgc_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip) { diff --git a/src/librt/primitives/tgc/tgc_tess.c b/src/librt/primitives/tgc/tgc_tess.c new file mode 100644 index 00000000000..93263c64f57 --- /dev/null +++ b/src/librt/primitives/tgc/tgc_tess.c @@ -0,0 +1,720 @@ +/* T G C _ T E S S . C + * BRL-CAD + * + * Copyright (c) 1985-2023 United States Government as represented by + * the U.S. Army Research Laboratory. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * version 2.1 as published by the Free Software Foundation. + * + * This library is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this file; see the file named COPYING for more + * information. + */ +/** @addtogroup primitives */ +/** @{ */ +/** @file primitives/tgc/tgc_tess.c + * + * Truncated General Cone Tessellation. + * + */ +/** @} */ + +#include "common.h" + +#include +#include +#include +#include "bio.h" + +#include "bu/cv.h" +#include "vmath.h" +#include "rt/db4.h" +#include "nmg.h" +#include "rt/geom.h" +#include "raytrace.h" + +#include "../../librt_private.h" + +#define MAX_RATIO 10.0 /* maximum allowed height-to-width ration for triangles */ + +/** + * Tessellation of the TGC. + * + * Returns - + * -1 failure + * 0 OK. *r points to nmgregion that holds this tessellation. + */ + +struct tgc_pts +{ + point_t pt; + vect_t tan_axb; + struct vertex *v; + char dont_use; +}; + + +/* version using tolerances */ +int +rt_tgc_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bg_tess_tol *ttol, const struct bn_tol *tol) +{ + struct shell *s; /* shell to hold facetted TGC */ + struct faceuse *fu, *fu_top, *fu_base; + struct rt_tgc_internal *tip; + fastf_t radius; /* bounding sphere radius */ + fastf_t max_radius; /* max of a, b, c, d */ + fastf_t h, a, b, c, d; /* lengths of TGC vectors */ + fastf_t inv_length; /* 1.0/length of a vector */ + vect_t unit_a, unit_b, unit_c, unit_d; /* units vectors in a, b, c, d directions */ + fastf_t rel, absolute, norm; /* interpreted tolerances */ + fastf_t alpha_tol; /* final tolerance for ellipse parameter */ + fastf_t abs_tol; /* handle invalid ttol->abs */ + size_t nells; /* total number of ellipses */ + size_t nsegs; /* number of vertices/ellipse */ + vect_t *A; /* array of A vectors for ellipses */ + vect_t *B; /* array of B vectors for ellipses */ + fastf_t *factors; /* array of ellipse locations along height vector */ + vect_t vtmp; + vect_t normal; /* normal vector */ + vect_t rev_norm; /* reverse normal */ + struct tgc_pts **pts; /* array of points (pts[ellipse#][seg#]) */ + struct bu_ptbl verts; /* table of vertices used for top and bottom faces */ + struct bu_ptbl faces; /* table of faceuses for nmg_gluefaces */ + struct vertex **v[3]; /* array for making triangular faces */ + + size_t i; + + VSETALL(unit_a, 0); + VSETALL(unit_b, 0); + VSETALL(unit_c, 0); + VSETALL(unit_d, 0); + + RT_CK_DB_INTERNAL(ip); + tip = (struct rt_tgc_internal *)ip->idb_ptr; + RT_TGC_CK_MAGIC(tip); + + if (ttol->abs > 0.0 && ttol->abs < tol->dist) { + bu_log("WARNING: tessellation tolerance is %fmm while calculational tolerance is %fmm\n", + ttol->abs, tol->dist); + bu_log("Cannot tessellate a TGC to finer tolerance than the calculational tolerance\n"); + abs_tol = tol->dist; + } else { + abs_tol = ttol->abs; + } + + h = MAGNITUDE(tip->h); + a = MAGNITUDE(tip->a); + if (2.0*a <= tol->dist) + a = 0.0; + b = MAGNITUDE(tip->b); + if (2.0*b <= tol->dist) + b = 0.0; + c = MAGNITUDE(tip->c); + if (2.0*c <= tol->dist) + c = 0.0; + d = MAGNITUDE(tip->d); + if (2.0*d <= tol->dist) + d = 0.0; + + if (ZERO(a) && ZERO(b) && (ZERO(c) || ZERO(d))) { + bu_log("Illegal TGC a, b, and c or d less than tolerance\n"); + return -1; + } else if (ZERO(c) && ZERO(d) && (ZERO(a) || ZERO(b))) { + bu_log("Illegal TGC c, d, and a or b less than tolerance\n"); + return -1; + } + + if (a > 0.0) { + inv_length = 1.0/a; + VSCALE(unit_a, tip->a, inv_length); + } + if (b > 0.0) { + inv_length = 1.0/b; + VSCALE(unit_b, tip->b, inv_length); + } + if (c > 0.0) { + inv_length = 1.0/c; + VSCALE(unit_c, tip->c, inv_length); + } + if (d > 0.0) { + inv_length = 1.0/d; + VSCALE(unit_d, tip->d, inv_length); + } + + /* get bounding sphere radius for relative tolerance */ + radius = h/2.0; + max_radius = 0.0; + if (a > max_radius) + max_radius = a; + if (b > max_radius) + max_radius = b; + if (c > max_radius) + max_radius = c; + if (d > max_radius) + max_radius = d; + + if (max_radius > radius) + radius = max_radius; + + if (abs_tol <= 0.0 && ttol->rel <= 0.0 && ttol->norm <= 0.0) { + /* no tolerances specified, use 10% relative tolerance */ + if ((radius * 0.2) < max_radius) + alpha_tol = 2.0 * acos(1.0 - 2.0 * radius * 0.1 / max_radius); + else + alpha_tol = M_PI_2; + } else { + if (abs_tol > 0.0) + absolute = 2.0 * acos(1.0 - abs_tol/max_radius); + else + absolute = M_PI_2; + + if (ttol->rel > 0.0) { + if (ttol->rel * 2.0 * radius < max_radius) + rel = 2.0 * acos(1.0 - ttol->rel * 2.0 * radius/max_radius); + else + rel = M_PI_2; + } else + rel = M_PI_2; + + if (ttol->norm > 0.0) { + fastf_t norm_top, norm_bot; + + if (anorm) * (a/b)); + else + norm_bot = 2.0 * atan(tan(ttol->norm) * (b/a)); + + if (cnorm) * (c/d)); + else + norm_top = 2.0 * atan(tan(ttol->norm) * (d/c)); + + if (norm_bot < norm_top) + norm = norm_bot; + else + norm = norm_top; + } else + norm = M_PI_2; + + if (absolute < rel) + alpha_tol = absolute; + else + alpha_tol = rel; + if (norm < alpha_tol) + alpha_tol = norm; + } + + /* get number of segments per quadrant */ + nsegs = (int)(M_PI_2 / alpha_tol + 0.9999); + if (nsegs < 2) + nsegs = 2; + + /* and for complete ellipse */ + nsegs *= 4; + + /* get number and placement of intermediate ellipses */ + { + fastf_t ratios[4], max_ratio; + fastf_t new_ratio = 0; + int which_ratio; + fastf_t len_ha, len_hb; + vect_t ha, hb; + fastf_t ang; + fastf_t sin_ang, cos_ang, cos_m_1_sq, sin_sq; + fastf_t len_A, len_B, len_C, len_D; + size_t bot_ell=0; + size_t top_ell=1; + int reversed=0; + + nells = 2; + + max_ratio = MAX_RATIO + 1.0; + + factors = (fastf_t *)bu_malloc(nells*sizeof(fastf_t), "factors"); + A = (vect_t *)bu_malloc(nells*sizeof(vect_t), "A vectors"); + B = (vect_t *)bu_malloc(nells*sizeof(vect_t), "B vectors"); + + factors[bot_ell] = 0.0; + factors[top_ell] = 1.0; + VMOVE(A[bot_ell], tip->a); + VMOVE(A[top_ell], tip->c); + VMOVE(B[bot_ell], tip->b); + VMOVE(B[top_ell], tip->d); + + /* make sure that AxB points in the general direction of H */ + VCROSS(vtmp, A[0], B[0]); + if (VDOT(vtmp, tip->h) < 0.0) { + VMOVE(A[bot_ell], tip->b); + VMOVE(A[top_ell], tip->d); + VMOVE(B[bot_ell], tip->a); + VMOVE(B[top_ell], tip->c); + reversed = 1; + } + ang = M_2PI/((double)nsegs); + sin_ang = sin(ang); + cos_ang = cos(ang); + cos_m_1_sq = (cos_ang - 1.0)*(cos_ang - 1.0); + sin_sq = sin_ang*sin_ang; + + VJOIN2(ha, tip->h, 1.0, tip->c, -1.0, tip->a); + VJOIN2(hb, tip->h, 1.0, tip->d, -1.0, tip->b); + len_ha = MAGNITUDE(ha); + len_hb = MAGNITUDE(hb); + + while (max_ratio > MAX_RATIO) { + fastf_t tri_width; + + len_A = MAGNITUDE(A[bot_ell]); + if (2.0*len_A <= tol->dist) + len_A = 0.0; + len_B = MAGNITUDE(B[bot_ell]); + if (2.0*len_B <= tol->dist) + len_B = 0.0; + len_C = MAGNITUDE(A[top_ell]); + if (2.0*len_C <= tol->dist) + len_C = 0.0; + len_D = MAGNITUDE(B[top_ell]); + if (2.0*len_D <= tol->dist) + len_D = 0.0; + + if ((len_B > 0.0 && len_D > 0.0) || + (len_B > 0.0 && (ZERO(len_D) && ZERO(len_C)))) + { + tri_width = sqrt(cos_m_1_sq*len_A*len_A + sin_sq*len_B*len_B); + ratios[0] = (factors[top_ell] - factors[bot_ell])*len_ha + /tri_width; + } else + ratios[0] = 0.0; + + if ((len_A > 0.0 && len_C > 0.0) || + (len_A > 0.0 && (ZERO(len_C) && ZERO(len_D)))) + { + tri_width = sqrt(sin_sq*len_A*len_A + cos_m_1_sq*len_B*len_B); + ratios[1] = (factors[top_ell] - factors[bot_ell])*len_hb + /tri_width; + } else + ratios[1] = 0.0; + + if ((len_D > 0.0 && len_B > 0.0) || + (len_D > 0.0 && (ZERO(len_A) && ZERO(len_B)))) + { + tri_width = sqrt(cos_m_1_sq*len_C*len_C + sin_sq*len_D*len_D); + ratios[2] = (factors[top_ell] - factors[bot_ell])*len_ha + /tri_width; + } else + ratios[2] = 0.0; + + if ((len_C > 0.0 && len_A > 0.0) || + (len_C > 0.0 && (ZERO(len_A) && ZERO(len_B)))) + { + tri_width = sqrt(sin_sq*len_C*len_C + cos_m_1_sq*len_D*len_D); + ratios[3] = (factors[top_ell] - factors[bot_ell])*len_hb + /tri_width; + } else + ratios[3] = 0.0; + + which_ratio = -1; + max_ratio = 0.0; + + for (i=0; i<4; i++) { + if (ratios[i] > max_ratio) { + max_ratio = ratios[i]; + which_ratio = i; + } + } + + if (ZERO(len_A) && ZERO(len_B) && ZERO(len_C) && ZERO(len_D)) { + if (top_ell == nells - 1) { + VMOVE(A[top_ell-1], A[top_ell]); + VMOVE(B[top_ell-1], A[top_ell]); + factors[top_ell-1] = factors[top_ell]; + } else if (bot_ell == 0) { + for (i=0; i 0.5) + new_ratio = 0.5; + } else if (which_ratio == 2 || which_ratio == 3) { + new_ratio = 1.0 - MAX_RATIO/max_ratio; + if (top_ell == nells - 1 && new_ratio < 0.5) + new_ratio = 0.5; + } else { + /* no MAX? */ + bu_bomb("rt_tgc_tess: Should never get here!!\n"); + } + + nells++; + factors = (fastf_t *)bu_realloc(factors, nells*sizeof(fastf_t), "factors"); + A = (vect_t *)bu_realloc(A, nells*sizeof(vect_t), "A vectors"); + B = (vect_t *)bu_realloc(B, nells*sizeof(vect_t), "B vectors"); + + for (i=nells-1; i>top_ell; i--) { + factors[i] = factors[i-1]; + VMOVE(A[i], A[i-1]); + VMOVE(B[i], B[i-1]); + } + + factors[top_ell] = factors[bot_ell] + + new_ratio*(factors[top_ell+1] - factors[bot_ell]); + + if (reversed) { + VBLEND2(A[top_ell], (1.0-factors[top_ell]), tip->b, factors[top_ell], tip->d); + VBLEND2(B[top_ell], (1.0-factors[top_ell]), tip->a, factors[top_ell], tip->c); + } else { + VBLEND2(A[top_ell], (1.0-factors[top_ell]), tip->a, factors[top_ell], tip->c); + VBLEND2(B[top_ell], (1.0-factors[top_ell]), tip->b, factors[top_ell], tip->d); + } + + if (which_ratio == 0 || which_ratio == 1) { + top_ell++; + bot_ell++; + } + + } + + } + + /* get memory for points */ + pts = (struct tgc_pts **)bu_calloc(nells, sizeof(struct tgc_pts *), "rt_tgc_tess: pts"); + for (i=0; iv); + } else if (i == nells-1 && ZERO(c) && ZERO(d)) { + VADD2(pts[i][j].pt, tip->v, tip->h); + } else { + VJOIN3(pts[i][j].pt, tip->v, h_factor, tip->h, cos_alpha, A[i], sin_alpha, B[i]); + } + + /* Storing the tangent here while sines and cosines are available */ + if (i == 0 && ZERO(a) && ZERO(b)) { + VCOMB2(pts[0][j].tan_axb, -sin_alpha, unit_c, cos_alpha, unit_d); + } else if (i == nells-1 && ZERO(c) && ZERO(d)) { + VCOMB2(pts[i][j].tan_axb, -sin_alpha, unit_a, cos_alpha, unit_b); + } else { + VCOMB2(pts[i][j].tan_axb, -sin_alpha, A[i], cos_alpha, B[i]); + } + } + } + + /* make sure no edges will be created with length < tol->dist */ + for (i=0; i tol->dist_sq) { + VMOVE(curr_pt, pts[i][j].pt); + } else { + /* don't use this point, it will create a too short edge */ + pts[i][j].dont_use = 'n'; + } + } + } + + /* make region, shell, vertex */ + *r = nmg_mrsv(m); + s = BU_LIST_FIRST(shell, &(*r)->s_hd); + + bu_ptbl_init(&verts, 64, " &verts "); + bu_ptbl_init(&faces, 64, " &faces "); + /* Make bottom face */ + if (a > 0.0 && b > 0.0) { + for (i=nsegs; i>0; i--) { + /* reverse order to get outward normal */ + if (!pts[0][i-1].dont_use) + bu_ptbl_ins(&verts, (long *)&pts[0][i-1].v); + } + + if (BU_PTBL_LEN(&verts) > 2) { + fu_base = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&verts), BU_PTBL_LEN(&verts)); + bu_ptbl_ins(&faces, (long *)fu_base); + } else + fu_base = (struct faceuse *)NULL; + } else + fu_base = (struct faceuse *)NULL; + + + /* Make top face */ + if (c > 0.0 && d > 0.0) { + bu_ptbl_reset(&verts); + for (i=0; i 2) { + fu_top = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&verts), BU_PTBL_LEN(&verts)); + bu_ptbl_ins(&faces, (long *)fu_top); + } else + fu_top = (struct faceuse *)NULL; + } else + fu_top = (struct faceuse *)NULL; + + /* Free table of vertices */ + bu_ptbl_free(&verts); + + /* Make triangular faces */ + for (i=0; i 0.0 || b > 0.0) { + if (!pts[i][k].dont_use) { + v[0] = curr_bot; + v[1] = &pts[i][k].v; + if (i+1 == nells-1 && ZERO(c) && ZERO(d)) + v[2] = &pts[i+1][0].v; + else + v[2] = curr_top; + fu = nmg_cmface(s, v, 3); + bu_ptbl_ins(&faces, (long *)fu); + curr_bot = &pts[i][k].v; + } + } + + if (i != nells-2 || c > 0.0 || d > 0.0) { + if (!pts[i+1][k].dont_use) { + v[0] = &pts[i+1][k].v; + v[1] = curr_top; + if (i == 0 && ZERO(a) && ZERO(b)) + v[2] = &pts[i][0].v; + else + v[2] = curr_bot; + fu = nmg_cmface(s, v, 3); + bu_ptbl_ins(&faces, (long *)fu); + curr_top = &pts[i+1][k].v; + } + } + } + } + + /* Assign geometry */ + for (i=0; iv); + } else if (i == nells-1 && ZERO(c) && ZERO(d)) { + if (j == 0) { + VADD2(pt_geom, tip->v, tip->h); + nmg_vertex_gv(pts[i][0].v, pt_geom); + } + } else if (pts[i][j].v) + nmg_vertex_gv(pts[i][j].v, pts[i][j].pt); + + /* Storing the tangent here while sines and cosines are available */ + if (i == 0 && ZERO(a) && ZERO(b)) { + VCOMB2(pts[0][j].tan_axb, -sin_alpha, unit_c, cos_alpha, unit_d); + } else if (i == nells-1 && ZERO(c) && ZERO(d)) { + VCOMB2(pts[i][j].tan_axb, -sin_alpha, unit_a, cos_alpha, unit_b); + } else { + VCOMB2(pts[i][j].tan_axb, -sin_alpha, A[i], cos_alpha, B[i]); + } + } + } + + /* Associate face plane equations */ + for (i=0; ivu_hd)) { + NMG_CK_VERTEXUSE(vu); + + fu = nmg_find_fu_of_vu(vu); + NMG_CK_FACEUSE(fu); + + /* don't need vertexuse normals for faces that are really flat */ + if (fu == fu_base || fu->fumate_p == fu_base || + fu == fu_top || fu->fumate_p == fu_top) + continue; + + if (fu->orientation == OT_SAME) + nmg_vertexuse_nv(vu, normal); + else if (fu->orientation == OT_OPPOSITE) + nmg_vertexuse_nv(vu, rev_norm); + } + } + } + } + + /* Finished with storage, so free it */ + bu_free((char *)factors, "rt_tgc_tess: factors"); + bu_free((char *)A, "rt_tgc_tess: A"); + bu_free((char *)B, "rt_tgc_tess: B"); + for (i=0; ilu_hd)) { + struct edgeuse *eu; + + NMG_CK_LOOPUSE(lu); + + if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC) + continue; + + for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) { + struct edge *e; + + NMG_CK_EDGEUSE(eu); + e = eu->e_p; + NMG_CK_EDGE(e); + e->is_real = 1; + } + } + } + + nmg_region_a(*r, tol); + + /* glue faces together */ + nmg_gluefaces((struct faceuse **)BU_PTBL_BASEADDR(&faces), BU_PTBL_LEN(&faces), &RTG.rtg_vlfree, tol); + bu_ptbl_free(&faces); + + return 0; +} + +/* + * Local Variables: + * mode: C + * tab-width: 8 + * indent-tabs-mode: t + * c-file-style: "stroustrup" + * End: + * ex: shiftwidth=4 tabstop=8 + */ diff --git a/src/librt/primitives/tgc/tgc_tnurb.c b/src/librt/primitives/tgc/tgc_tnurb.c new file mode 100644 index 00000000000..68d0db822db --- /dev/null +++ b/src/librt/primitives/tgc/tgc_tnurb.c @@ -0,0 +1,514 @@ +/* T G C _ T N U R B . C + * BRL-CAD + * + * Copyright (c) 1985-2023 United States Government as represented by + * the U.S. Army Research Laboratory. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * version 2.1 as published by the Free Software Foundation. + * + * This library is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this file; see the file named COPYING for more + * information. + */ +/** @addtogroup primitives */ +/** @{ */ +/** @file primitives/tgc/tgc_tnurb.c + * + * Truncated General Cone NMG nurb. + * + */ +/** @} */ + +#include "common.h" + +#include +#include +#include +#include "bio.h" + +#include "bu/cv.h" +#include "vmath.h" +#include "rt/db4.h" +#include "nmg.h" +#include "rt/geom.h" +#include "raytrace.h" + +#include "../../librt_private.h" + +#define RAT M_SQRT1_2 + +static const fastf_t nmg_tgc_unitcircle[36] = { + 1.0, 0.0, 0.0, 1.0, + RAT, -RAT, 0.0, RAT, + 0.0, -1.0, 0.0, 1.0, + -RAT, -RAT, 0.0, RAT, + -1.0, 0.0, 0.0, 1.0, + -RAT, RAT, 0.0, RAT, + 0.0, 1.0, 0.0, 1.0, + RAT, RAT, 0.0, RAT, + 1.0, 0.0, 0.0, 1.0 +}; + +static const fastf_t nmg_uv_unitcircle[27] = { + 1.0, .5, 1.0, + RAT, RAT, RAT, + .5, 1.0, 1.0, + 0.0, RAT, RAT, + 0.0, .5, 1.0, + 0.0, 0.0, RAT, + .5, 0.0, 1.0, + RAT, 0.0, RAT, + 1.0, .5, 1.0 +}; + +static void +nmg_tgc_disk(struct faceuse *fu, fastf_t *rmat, fastf_t height, int flip) +{ + struct face_g_snurb * fg; + struct loopuse * lu; + struct edgeuse * eu; + struct edge_g_cnurb * eg; + fastf_t *mptr; + int i; + vect_t vect; + point_t point; + + nmg_face_g_snurb(fu, + 2, 2, /* u, v order */ + 4, 4, /* number of knots */ + NULL, NULL, /* initial knot vectors */ + 2, 2, /* n_rows, n_cols */ + RT_NURB_MAKE_PT_TYPE(3, 2, 0), + NULL); /* Initial mesh */ + + fg = fu->f_p->g.snurb_p; + + + fg->u.knots[0] = 0.0; + fg->u.knots[1] = 0.0; + fg->u.knots[2] = 1.0; + fg->u.knots[3] = 1.0; + + fg->v.knots[0] = 0.0; + fg->v.knots[1] = 0.0; + fg->v.knots[2] = 1.0; + fg->v.knots[3] = 1.0; + + if (flip) { + fg->ctl_points[0] = 1.; + fg->ctl_points[1] = -1.; + fg->ctl_points[2] = height; + + fg->ctl_points[3] = -1; + fg->ctl_points[4] = -1.; + fg->ctl_points[5] = height; + + fg->ctl_points[6] = 1.; + fg->ctl_points[7] = 1.; + fg->ctl_points[8] = height; + + fg->ctl_points[9] = -1.; + fg->ctl_points[10] = 1.; + fg->ctl_points[11] = height; + } else { + + fg->ctl_points[0] = -1.; + fg->ctl_points[1] = -1.; + fg->ctl_points[2] = height; + + fg->ctl_points[3] = 1; + fg->ctl_points[4] = -1.; + fg->ctl_points[5] = height; + + fg->ctl_points[6] = -1.; + fg->ctl_points[7] = 1.; + fg->ctl_points[8] = height; + + fg->ctl_points[9] = 1.; + fg->ctl_points[10] = 1.; + fg->ctl_points[11] = height; + } + + /* multiple the matrix to get orientation and scaling that we want */ + mptr = fg->ctl_points; + + i = fg->s_size[0] * fg->s_size[1]; + + for (; i> 0; i--) { + MAT4X3PNT(vect, rmat, mptr); + mptr[0] = vect[0]; + mptr[1] = vect[1]; + mptr[2] = vect[2]; + mptr += 3; + } + + lu = BU_LIST_FIRST(loopuse, &fu->lu_hd); + NMG_CK_LOOPUSE(lu); + eu= BU_LIST_FIRST(edgeuse, &lu->down_hd); + NMG_CK_EDGEUSE(eu); + + + if (!flip) { + nmg_nurb_s_eval(fu->f_p->g.snurb_p, + nmg_uv_unitcircle[0], nmg_uv_unitcircle[1], point); + nmg_vertex_gv(eu->vu_p->v_p, point); + } else { + nmg_nurb_s_eval(fu->f_p->g.snurb_p, + nmg_uv_unitcircle[12], nmg_uv_unitcircle[13], point); + nmg_vertex_gv(eu->vu_p->v_p, point); + } + + nmg_edge_g_cnurb(eu, 3, 12, NULL, 9, RT_NURB_MAKE_PT_TYPE(3, 3, 1), + NULL); + + eg = eu->g.cnurb_p; + eg->order = 3; + + eg->k.knots[0] = 0.0; + eg->k.knots[1] = 0.0; + eg->k.knots[2] = 0.0; + eg->k.knots[3] = .25; + eg->k.knots[4] = .25; + eg->k.knots[5] = .5; + eg->k.knots[6] = .5; + eg->k.knots[7] = .75; + eg->k.knots[8] = .75; + eg->k.knots[9] = 1.0; + eg->k.knots[10] = 1.0; + eg->k.knots[11] = 1.0; + + if (!flip) { + for (i = 0; i < 27; i++) + eg->ctl_points[i] = nmg_uv_unitcircle[i]; + } else { + + VSET(&eg->ctl_points[0], 0.0, .5, 1.0); + VSET(&eg->ctl_points[3], 0.0, 0.0, RAT); + VSET(&eg->ctl_points[6], 0.5, 0.0, 1.0); + VSET(&eg->ctl_points[9], RAT, 0.0, RAT); + VSET(&eg->ctl_points[12], 1.0, .5, 1.0); + VSET(&eg->ctl_points[15], RAT, RAT, RAT); + VSET(&eg->ctl_points[18], .5, 1.0, 1.0); + VSET(&eg->ctl_points[21], 0.0, RAT, RAT); + VSET(&eg->ctl_points[24], 0.0, .5, 1.0); + } +} + + + + +/* Create a cylinder with a top surface and a bottom surface + * defined by the ellipsoids at the top and bottom of the + * cylinder, the top_mat, and bot_mat are applied to a unit circle + * for the top row of the surface and the bot row of the surface + * respectively. + */ +static void +nmg_tgc_nurb_cyl(struct faceuse *fu, fastf_t *top_mat, fastf_t *bot_mat) +{ + struct face_g_snurb * fg; + struct loopuse * lu; + struct edgeuse * eu; + fastf_t * mptr; + int i; + hvect_t vect; + point_t uvw; + point_t point; + hvect_t hvect; + + nmg_face_g_snurb(fu, + 3, 2, + 12, 4, + NULL, NULL, + 2, 9, + RT_NURB_MAKE_PT_TYPE(4, 3, 1), + NULL); + + fg = fu->f_p->g.snurb_p; + + fg->v.knots[0] = 0.0; + fg->v.knots[1] = 0.0; + fg->v.knots[2] = 1.0; + fg->v.knots[3] = 1.0; + + fg->u.knots[0] = 0.0; + fg->u.knots[1] = 0.0; + fg->u.knots[2] = 0.0; + fg->u.knots[3] = .25; + fg->u.knots[4] = .25; + fg->u.knots[5] = .5; + fg->u.knots[6] = .5; + fg->u.knots[7] = .75; + fg->u.knots[8] = .75; + fg->u.knots[9] = 1.0; + fg->u.knots[10] = 1.0; + fg->u.knots[11] = 1.0; + + mptr = fg->ctl_points; + + for (i = 0; i < 9; i++) { + MAT4X4PNT(vect, top_mat, &nmg_tgc_unitcircle[i*4]); + mptr[0] = vect[0]; + mptr[1] = vect[1]; + mptr[2] = vect[2]; + mptr[3] = vect[3]; + mptr += 4; + } + + for (i = 0; i < 9; i++) { + MAT4X4PNT(vect, bot_mat, &nmg_tgc_unitcircle[i*4]); + mptr[0] = vect[0]; + mptr[1] = vect[1]; + mptr[2] = vect[2]; + mptr[3] = vect[3]; + mptr += 4; + } + + /* Assign edgeuse parameters & vertex geometry */ + + lu = BU_LIST_FIRST(loopuse, &fu->lu_hd); + NMG_CK_LOOPUSE(lu); + eu = BU_LIST_FIRST(edgeuse, &lu->down_hd); + NMG_CK_EDGEUSE(eu); + + /* March around the fu's loop assigning uv parameter values */ + + nmg_nurb_s_eval(fg, 0.0, 0.0, hvect); + HDIVIDE(point, hvect); + nmg_vertex_gv(eu->vu_p->v_p, point); /* 0, 0 vertex */ + + VSET(uvw, 0, 0, 0); + nmg_vertexuse_a_cnurb(eu->vu_p, uvw); + VSET(uvw, 1, 0, 0); + nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); + eu = BU_LIST_NEXT(edgeuse, &eu->l); + + VSET(uvw, 1, 0, 0); + nmg_vertexuse_a_cnurb(eu->vu_p, uvw); + VSET(uvw, 1, 1, 0); + nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); + eu = BU_LIST_NEXT(edgeuse, &eu->l); + + VSET(uvw, 1, 1, 0); + nmg_vertexuse_a_cnurb(eu->vu_p, uvw); + VSET(uvw, 0, 1, 0); + nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); + + nmg_nurb_s_eval(fg, 1., 1., hvect); + HDIVIDE(point, hvect); + nmg_vertex_gv(eu->vu_p->v_p, point); /* 4, 1 vertex */ + + eu = BU_LIST_NEXT(edgeuse, &eu->l); + + VSET(uvw, 0, 1, 0); + nmg_vertexuse_a_cnurb(eu->vu_p, uvw); + VSET(uvw, 0, 0, 0); + nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); + + /* Create the edge loop geometry */ + + lu = BU_LIST_FIRST(loopuse, &fu->lu_hd); + NMG_CK_LOOPUSE(lu); + eu = BU_LIST_FIRST(edgeuse, &lu->down_hd); + NMG_CK_EDGEUSE(eu); + + for (i = 0; i < 4; i++) { + nmg_edge_g_cnurb_plinear(eu); + eu = BU_LIST_NEXT(edgeuse, &eu->l); + } +} + + + + +/** + * "Tessellate an TGC into a trimmed-NURB-NMG data structure. + * Computing NURB surfaces and trimming curves to interpolate + * the parameters of the TGC + * + * The process is to create the nmg topology of the TGC fill it + * in with a unit cylinder geometry (i.e. unitcircle at the top (0, 0, 1) + * unit cylinder of radius 1, and unitcircle at the bottom), and then + * scale it with a perspective matrix derived from the parameters of the + * tgc. The result is three trimmed nurb surfaces which interpolate the + * parameters of the original TGC. + * + * Returns - + * -1 failure + * 0 OK. *r points to nmgregion that holds this tessellation + */ + +int +rt_tgc_tnurb(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bn_tol *tol) +{ + struct rt_tgc_internal *tip; + + struct shell *s; + struct vertex *verts[2]; + struct vertex **vertp[4]; + struct faceuse * top_fu; + struct faceuse * cyl_fu; + struct faceuse * bot_fu; + vect_t uvw; + struct loopuse *lu; + struct edgeuse *eu; + struct edgeuse *top_eu; + struct edgeuse *bot_eu; + + mat_t mat; + mat_t imat, omat, top_mat, bot_mat; + vect_t anorm; + vect_t bnorm; + vect_t cnorm; + + + /* Get the internal representation of the tgc */ + + RT_CK_DB_INTERNAL(ip); + tip = (struct rt_tgc_internal *) ip->idb_ptr; + RT_TGC_CK_MAGIC(tip); + + /* Create the NMG Topology */ + + *r = nmg_mrsv(m); /* Make region, empty shell, vertex */ + s = BU_LIST_FIRST(shell, &(*r)->s_hd); + + + /* Create transformation matrix for the top cap surface*/ + + MAT_IDN(omat); + MAT_IDN(mat); + + omat[0] = MAGNITUDE(tip->c); + omat[5] = MAGNITUDE(tip->d); + omat[3] = tip->v[0] + tip->h[0]; + omat[7] = tip->v[1] + tip->h[1]; + omat[11] = tip->v[2] + tip->h[2]; + + bn_mat_mul(imat, mat, omat); + + VMOVE(anorm, tip->c); + VMOVE(bnorm, tip->d); + VCROSS(cnorm, tip->c, tip->d); + VUNITIZE(anorm); + VUNITIZE(bnorm); + VUNITIZE(cnorm); + + MAT_IDN(omat); + + VMOVE(&omat[0], anorm); + VMOVE(&omat[4], bnorm); + VMOVE(&omat[8], cnorm); + + + bn_mat_mul(top_mat, omat, imat); + + /* Create topology for top cap surface */ + + verts[0] = verts[1] = NULL; + vertp[0] = &verts[0]; + top_fu = nmg_cmface(s, vertp, 1); + + lu = BU_LIST_FIRST(loopuse, &top_fu->lu_hd); + NMG_CK_LOOPUSE(lu); + eu= BU_LIST_FIRST(edgeuse, &lu->down_hd); + NMG_CK_EDGEUSE(eu); + top_eu = eu; + + VSET(uvw, 0, 0, 0); + nmg_vertexuse_a_cnurb(eu->vu_p, uvw); + VSET(uvw, 1, 0, 0); + nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); + + /* Top cap surface */ + nmg_tgc_disk(top_fu, top_mat, 0.0, 0); + + /* Create transformation matrix for the bottom cap surface*/ + + MAT_IDN(omat); + MAT_IDN(mat); + + omat[0] = MAGNITUDE(tip->a); + omat[5] = MAGNITUDE(tip->b); + omat[3] = tip->v[0]; + omat[7] = tip->v[1]; + omat[11] = tip->v[2]; + + bn_mat_mul(imat, mat, omat); + + VMOVE(anorm, tip->a); + VMOVE(bnorm, tip->b); + VCROSS(cnorm, tip->a, tip->b); + VUNITIZE(anorm); + VUNITIZE(bnorm); + VUNITIZE(cnorm); + + MAT_IDN(omat); + + VMOVE(&omat[0], anorm); + VMOVE(&omat[4], bnorm); + VMOVE(&omat[8], cnorm); + + bn_mat_mul(bot_mat, omat, imat); + + /* Create topology for bottom cap surface */ + + vertp[0] = &verts[1]; + bot_fu = nmg_cmface(s, vertp, 1); + + lu = BU_LIST_FIRST(loopuse, &bot_fu->lu_hd); + NMG_CK_LOOPUSE(lu); + eu= BU_LIST_FIRST(edgeuse, &lu->down_hd); + NMG_CK_EDGEUSE(eu); + bot_eu = eu; + + VSET(uvw, 0, 0, 0); + nmg_vertexuse_a_cnurb(eu->vu_p, uvw); + VSET(uvw, 1, 0, 0); + nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw); + + + nmg_tgc_disk(bot_fu, bot_mat, 0.0, 1); + + /* Create topology for cylinder surface */ + + vertp[0] = &verts[0]; + vertp[1] = &verts[0]; + vertp[2] = &verts[1]; + vertp[3] = &verts[1]; + cyl_fu = nmg_cmface(s, vertp, 4); + + nmg_tgc_nurb_cyl(cyl_fu, top_mat, bot_mat); + + /* fuse top cylinder edge to matching edge on body of cylinder */ + lu = BU_LIST_FIRST(loopuse, &cyl_fu->lu_hd); + + eu= BU_LIST_FIRST(edgeuse, &lu->down_hd); + + nmg_je(top_eu, eu); + + eu= BU_LIST_LAST(edgeuse, &lu->down_hd); + eu= BU_LIST_LAST(edgeuse, &eu->l); + + nmg_je(bot_eu, eu); + nmg_region_a(*r, tol); + + return 0; +} + +/* + * Local Variables: + * mode: C + * tab-width: 8 + * indent-tabs-mode: t + * c-file-style: "stroustrup" + * End: + * ex: shiftwidth=4 tabstop=8 + */