[brlcad-commits] SF.net SVN: brlcad:[51752] brlcad/trunk/src
Open Source Solid Modeling CAD
Brought to you by:
brlcad
From: <sta...@us...> - 2012-08-02 17:21:44
|
Revision: 51752 http://brlcad.svn.sourceforge.net/brlcad/?rev=51752&view=rev Author: starseeker Date: 2012-08-02 17:21:35 +0000 (Thu, 02 Aug 2012) Log Message: ----------- Start moving files into libnurbs. Modified Paths: -------------- brlcad/trunk/src/CMakeLists.txt brlcad/trunk/src/conv/step/CMakeLists.txt brlcad/trunk/src/conv/step/OpenNurbsInterfaces.cpp brlcad/trunk/src/conv/step/PullbackCurve.cpp brlcad/trunk/src/librt/CMakeLists.txt brlcad/trunk/src/librt/primitives/brep/brep_debug.h brlcad/trunk/src/librt/primitives/sketch/sketch_tess.cpp Added Paths: ----------- brlcad/trunk/src/libnurbs/CMakeLists.txt brlcad/trunk/src/libnurbs/opennurbs_ext.cpp brlcad/trunk/src/libnurbs/opennurbs_ext.h brlcad/trunk/src/libnurbs/opennurbs_fit.cpp brlcad/trunk/src/libnurbs/opennurbs_fit.h Removed Paths: ------------- brlcad/trunk/src/librt/opennurbs_ext.cpp brlcad/trunk/src/librt/opennurbs_ext.h brlcad/trunk/src/librt/opennurbs_fit.cpp brlcad/trunk/src/librt/opennurbs_fit.h Modified: brlcad/trunk/src/CMakeLists.txt =================================================================== --- brlcad/trunk/src/CMakeLists.txt 2012-08-02 16:58:07 UTC (rev 51751) +++ brlcad/trunk/src/CMakeLists.txt 2012-08-02 17:21:35 UTC (rev 51752) @@ -32,6 +32,7 @@ libbu libbn libsysv + libnurbs librt ) @@ -156,7 +157,7 @@ endforeach(subdir ${level_${current_level}_dirs}) endwhile(${current_level} LESS ${HIGHEST_TARGET_LEVEL}) -CMAKEFILES(README external libnurbs) +CMAKEFILES(README external) CMAKEFILES(Makefile.am) # Local Variables: Modified: brlcad/trunk/src/conv/step/CMakeLists.txt =================================================================== --- brlcad/trunk/src/conv/step/CMakeLists.txt 2012-08-02 16:58:07 UTC (rev 51751) +++ brlcad/trunk/src/conv/step/CMakeLists.txt 2012-08-02 17:21:35 UTC (rev 51752) @@ -461,6 +461,7 @@ set(stepg_LIBS libbu + libnurbs librt libwdb stepcore Modified: brlcad/trunk/src/conv/step/OpenNurbsInterfaces.cpp =================================================================== --- brlcad/trunk/src/conv/step/OpenNurbsInterfaces.cpp 2012-08-02 16:58:07 UTC (rev 51751) +++ brlcad/trunk/src/conv/step/OpenNurbsInterfaces.cpp 2012-08-02 17:21:35 UTC (rev 51752) @@ -98,7 +98,7 @@ /* FIXME: should not be peeking into a private header and cannot be a * public header (of librt). */ -#include "../../librt/opennurbs_ext.h" +#include "../../libnurbs/opennurbs_ext.h" ON_Brep * Modified: brlcad/trunk/src/conv/step/PullbackCurve.cpp =================================================================== --- brlcad/trunk/src/conv/step/PullbackCurve.cpp 2012-08-02 16:58:07 UTC (rev 51751) +++ brlcad/trunk/src/conv/step/PullbackCurve.cpp 2012-08-02 17:21:35 UTC (rev 51752) @@ -40,7 +40,7 @@ /* FIXME: should not be peeking into a private header and cannot be a * public header (of librt). */ -#include "../../librt/opennurbs_ext.h" +#include "../../libnurbs/opennurbs_ext.h" /* interface header */ #include "PullbackCurve.h" Added: brlcad/trunk/src/libnurbs/CMakeLists.txt =================================================================== --- brlcad/trunk/src/libnurbs/CMakeLists.txt (rev 0) +++ brlcad/trunk/src/libnurbs/CMakeLists.txt 2012-08-02 17:21:35 UTC (rev 51752) @@ -0,0 +1,63 @@ +# Include directories needed by libnurbs users +set(NURBS_INCLUDE_DIRS + ${BRLCAD_BINARY_DIR}/include + ${BRLCAD_SOURCE_DIR}/include + ${BU_INCLUDE_DIRS} + ${BN_INCLUDE_DIRS} + ${OPENNURBS_INCLUDE_DIR} + ${ZLIB_INCLUDE_DIR} + ) + +# locally used but not needed by users of the library +set(NURBS_LOCAL_INCLUDE_DIRS + ${CMAKE_CURRENT_SOURCE_DIR} + ${REGEX_INCLUDE_DIR} + ${TNT_INCLUDE_DIR} + #${BRLCAD_SOURCE_DIR}/src/other/Eigen + ) + +BRLCAD_LIB_INCLUDE_DIRS(nurbs NURBS_INCLUDE_DIRS NURBS_LOCAL_INCLUDE_DIRS) + +set(LIBNURBS_SOURCES + opennurbs_ext.cpp + #opennurbs_fit.cpp + ) + +set(librt_ignored_files + opennurbs_ext.h + opennurbs_fit.h + opennurbs_fit.cpp + ) +CMAKEFILES(${librt_ignored_files}) + +# Initialize librt_DEFINES in case of reconfiguration +set(libnurbs_DEFINES "") + +# Set libnurbs compile definitions +get_property(libnurbs_DEFINES GLOBAL PROPERTY libnurbs_DEFINES) +list(APPEND libnurbs_DEFINES "OBJ_BREP=1") +#list(APPEND libnurbs_DEFINES "EIGEN_NO_DEBUG") +#list(APPEND libnurbs_DEFINES "EIGEN_USE_NEW_STDVECTOR") +if(CPP_DLL_DEFINES) + list(APPEND libnurbs_DEFINES "ON_DLL_IMPORTS") +endif(CPP_DLL_DEFINES) +set_property(GLOBAL PROPERTY libnurbs_DEFINES "${libnurbs_DEFINES}") + +BRLCAD_ADDLIB(libnurbs "${LIBNURBS_SOURCES}" "libbn;libbu;${OPENNURBS_LIBRARY};${ZLIB_LIBRARY};${WINSOCK_LIB};${RPCRT_LIB};${STDCXX_LIBRARIES}" NOSTRICTCXX) +if(BUILD_STATIC_LIBS) + set_property(TARGET libnurbs-static APPEND PROPERTY COMPILE_DEFINITIONS "OBJ_BREP=1") +endif(BUILD_STATIC_LIBS) + +set_target_properties(libnurbs PROPERTIES VERSION 20.0.1 SOVERSION 20) +if(CPP_DLL_DEFINES) + if(BRLCAD_BUILD_STATIC_LIBS AND BRLCAD_ENABLE_BRLCAD_LIBRARY) + SET_PROPERTY(TARGET librt-static APPEND PROPERTY COMPILE_DEFINITIONS "ON_DLL_IMPORTS") + endif(BRLCAD_BUILD_STATIC_LIBS AND BRLCAD_ENABLE_BRLCAD_LIBRARY) +endif(CPP_DLL_DEFINES) + +# Local Variables: +# tab-width: 8 +# mode: cmake +# indent-tabs-mode: t +# End: +# ex: shiftwidth=2 tabstop=8 Property changes on: brlcad/trunk/src/libnurbs/CMakeLists.txt ___________________________________________________________________ Added: svn:mime-type + text/plain Added: svn:eol-style + native Copied: brlcad/trunk/src/libnurbs/opennurbs_ext.cpp (from rev 51750, brlcad/trunk/src/librt/opennurbs_ext.cpp) =================================================================== --- brlcad/trunk/src/libnurbs/opennurbs_ext.cpp (rev 0) +++ brlcad/trunk/src/libnurbs/opennurbs_ext.cpp 2012-08-02 17:21:35 UTC (rev 51752) @@ -0,0 +1,3400 @@ +/* O P E N N U R B S _ E X T . C P P + * BRL-CAD + * + * Copyright (c) 2007-2012 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. + */ +/** @file opennurbs_ext.cpp + * + * Implementation of routines openNURBS left out. + * + */ + +#include "common.h" + +#include "bio.h" +#include <assert.h> +#include <vector> + +#include "tnt.h" +#include "jama_lu.h" + +#include "vmath.h" +#include "bu.h" + +#include "brep.h" +#include "dvec.h" + +/* our interface header */ +#include "./opennurbs_ext.h" + + +#define RANGE_HI 0.55 +#define RANGE_LO 0.45 +#define UNIVERSAL_SAMPLE_COUNT 1001 + +#define BBOX_GROW 0.0 + +/// grows 3D BBox along each axis by this factor +#define BBOX_GROW_3D 0.1 + +/// arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bounds) +#define TOL 0.000001 + +/// another arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bounds) +#define TOL2 0.00001 + + +extern int brep_getSurfacePoint(const ON_3dPoint&, ON_2dPoint&, brlcad::BBNode*); + +bool +ON_NearZero(double x, double tolerance) { + return (x > -tolerance) && (x < tolerance); +} + + +namespace brlcad { + +inline void +distribute(const int count, const ON_3dVector* v, double x[], double y[], double z[]) +{ + for (int i = 0; i < count; i++) { + x[i] = v[i].x; + y[i] = v[i].y; + z[i] = v[i].z; + } +} + + +//-------------------------------------------------------------------------------- +// CurveTree +CurveTree::CurveTree(ON_BrepFace* face) : + m_face(face), m_adj_face_index(-99) +{ + m_root = initialLoopBBox(); + + for (int li = 0; li < face->LoopCount(); li++) { + bool innerLoop = (li > 0) ? true : false; + ON_BrepLoop* loop = face->Loop(li); + // for each trim + for (int ti = 0; ti < loop->m_ti.Count(); ti++) { + int adj_face_index = -99; + ON_BrepTrim& trim = face->Brep()->m_T[loop->m_ti[ti]]; + + if (trim.m_ei != -1) { // does not lie on a portion of a singular surface side + ON_BrepEdge& edge = face->Brep()->m_E[trim.m_ei]; + switch (trim.m_type) { + case ON_BrepTrim::unknown: + bu_log("ON_BrepTrim::unknown on Face:%d\n", face->m_face_index); + break; + case ON_BrepTrim::boundary: + //bu_log("ON_BrepTrim::boundary on Face:%d\n", face->m_face_index); + break; + case ON_BrepTrim::mated: + if (edge.m_ti.Count() == 2) { + if (face->m_face_index == face->Brep()->m_T[edge.m_ti[0]].FaceIndexOf()) { + adj_face_index = face->Brep()->m_T[edge.m_ti[1]].FaceIndexOf(); + } else { + adj_face_index = face->Brep()->m_T[edge.m_ti[0]].FaceIndexOf(); + } + } else { + bu_log("Mated Edge should have 2 adjacent faces, right? Face(%d) has %d trim indexes\n", face->m_face_index, edge.m_ti.Count()); + } + break; + case ON_BrepTrim::seam: + if (edge.m_ti.Count() == 2) { + if ((face->m_face_index == face->Brep()->m_T[edge.m_ti[0]].FaceIndexOf()) && (face->m_face_index == face->Brep()->m_T[edge.m_ti[1]].FaceIndexOf())) { + adj_face_index = face->m_face_index; + } else { + bu_log("Seamed Edge should have 1 faces sharing the trim so trim index should be one, right? Face(%d) has %d trim indexes\n", face->m_face_index, edge.m_ti.Count()); + bu_log("Face(%d) has %d, %d trim indexes\n", face->m_face_index, face->Brep()->m_T[edge.m_ti[0]].FaceIndexOf(), face->Brep()->m_T[edge.m_ti[1]].FaceIndexOf()); + } + } else if (edge.m_ti.Count() == 1) { + adj_face_index = face->m_face_index; + } else { + bu_log("Seamed Edge should have 1 faces sharing the trim so trim index should be one, right? Face(%d) has %d trim indexes\n", face->m_face_index, edge.m_ti.Count()); + } + break; + case ON_BrepTrim::singular: + bu_log("ON_BrepTrim::singular on Face:%d\n", face->m_face_index); + break; + case ON_BrepTrim::crvonsrf: + bu_log("ON_BrepTrim::crvonsrf on Face:%d\n", face->m_face_index); + break; + case ON_BrepTrim::ptonsrf: + bu_log("ON_BrepTrim::ptonsrf on Face:%d\n", face->m_face_index); + break; + case ON_BrepTrim::slit: + bu_log("ON_BrepTrim::slit on Face:%d\n", face->m_face_index); + break; + default: + bu_log("ON_BrepTrim::default on Face:%d\n", face->m_face_index); + } + } + const ON_Curve* trimCurve = trim.TrimCurveOf(); + double min, max; + (void) trimCurve->GetDomain(&min, &max); + ON_Interval t(min, max); + + TRACE("need to subdivide"); + // divide on param interval + + if (!trimCurve->IsLinear()) { + int knotcnt = trimCurve->SpanCount(); + double *knots = new double[knotcnt + 1]; + + trimCurve->GetSpanVector(knots); + std::list<double> splitlist; + for (int knot_index = 1; knot_index <= knotcnt; knot_index++) { + ON_Interval range(knots[knot_index - 1], knots[knot_index]); + + getHVTangents(trimCurve, range, splitlist); + } + for (std::list<double>::iterator l = splitlist.begin(); l != splitlist.end(); l++) { + double xmax = *l; + if (!NEAR_EQUAL(xmax, min, TOL)) { + m_root->addChild(subdivideCurve(trimCurve, adj_face_index, min, xmax, innerLoop, 0)); + } + min = xmax; + } + delete [] knots; + } else { + int knotcnt = trimCurve->SpanCount(); + double *knots = new double[knotcnt + 1]; + + trimCurve->GetSpanVector(knots); + for (int knot_index = 1; knot_index <= knotcnt; knot_index++) { + double xmax = knots[knot_index]; + if (!NEAR_EQUAL(xmax, min, TOL)) { + m_root->addChild(subdivideCurve(trimCurve, adj_face_index, min, xmax, innerLoop, 0)); + } + min = xmax; + } + delete [] knots; + } + + if (!NEAR_EQUAL(max, min, TOL)) { + m_root->addChild(subdivideCurve(trimCurve, adj_face_index, min, max, innerLoop, 0)); + } + } + } + getLeaves(m_sortedX); + m_sortedX.sort(sortX); + getLeaves(m_sortedY); + m_sortedY.sort(sortY); + + return; +} + + +CurveTree::~CurveTree() +{ + delete m_root; +} + + +BRNode* +CurveTree::getRootNode() const +{ + return m_root; +} + + +int +CurveTree::depth() +{ + return m_root->depth(); +} + + +ON_2dPoint +CurveTree::getClosestPointEstimate(const ON_3dPoint& pt) +{ + return m_root->getClosestPointEstimate(pt); +} + + +ON_2dPoint +CurveTree::getClosestPointEstimate(const ON_3dPoint& pt, ON_Interval& u, ON_Interval& v) +{ + return m_root->getClosestPointEstimate(pt, u, v); +} + + +void +CurveTree::getLeaves(std::list<BRNode*>& out_leaves) +{ + m_root->getLeaves(out_leaves); +} + + +void +CurveTree::getLeavesAbove(std::list<BRNode*>& out_leaves, const ON_Interval& u, const ON_Interval& v) +{ + point_t bmin, bmax; + double dist; + for (std::list<BRNode*>::iterator i = m_sortedX.begin(); i != m_sortedX.end(); i++) { + BRNode* br = dynamic_cast<BRNode*>(*i); + br->GetBBox(bmin, bmax); + + dist = TOL;//0.03*DIST_PT_PT(bmin, bmax); + if (bmax[X]+dist < u[0]) + continue; + if (bmin[X]-dist < u[1]) { + if (bmax[Y]+dist > v[0]) { + out_leaves.push_back(br); + } + } + } +} + + +void +CurveTree::getLeavesAbove(std::list<BRNode*>& out_leaves, const ON_2dPoint& pt, fastf_t tol) +{ + point_t bmin, bmax; + for (std::list<BRNode*>::iterator i = m_sortedX.begin(); i != m_sortedX.end(); i++) { + BRNode* br = dynamic_cast<BRNode*>(*i); + br->GetBBox(bmin, bmax); + + if (bmax[X]+tol < pt.x) + continue; + if (bmin[X]-tol < pt.x) { + if (bmax[Y]+tol > pt.y) { + out_leaves.push_back(br); + } + } + } +} + + +void +CurveTree::getLeavesRight(std::list<BRNode*>& out_leaves, const ON_Interval& u, const ON_Interval& v) +{ + point_t bmin, bmax; + double dist; + for (std::list<BRNode*>::iterator i = m_sortedX.begin(); i != m_sortedX.end(); i++) { + BRNode* br = dynamic_cast<BRNode*>(*i); + br->GetBBox(bmin, bmax); + + dist = TOL;//0.03*DIST_PT_PT(bmin, bmax); + if (bmax[Y]+dist < v[0]) + continue; + if (bmin[Y]-dist < v[1]) { + if (bmax[X]+dist > u[0]) { + out_leaves.push_back(br); + } + } + } +} + + +void +CurveTree::getLeavesRight(std::list<BRNode*>& out_leaves, const ON_2dPoint& pt, fastf_t tol) +{ + point_t bmin, bmax; + for (std::list<BRNode*>::iterator i = m_sortedX.begin(); i != m_sortedX.end(); i++) { + BRNode* br = dynamic_cast<BRNode*>(*i); + br->GetBBox(bmin, bmax); + + if (bmax[Y]+tol < pt.y) + continue; + if (bmin[Y]-tol < pt.y) { + if (bmax[X]+tol > pt.x) { + out_leaves.push_back(br); + } + } + } +} + + +fastf_t +CurveTree::getVerticalTangent(const ON_Curve *curve, fastf_t min, fastf_t max) +{ + fastf_t mid; + ON_3dVector tangent; + bool tanmin; + + // first lets check end points + tangent = curve->TangentAt(max); + if (NEAR_ZERO(tangent.x, TOL2)) + return max; + tangent = curve->TangentAt(min); + if (NEAR_ZERO(tangent.x, TOL2)) + return min; + + tanmin = (tangent[X] < 0.0); + while ((max-min) > TOL2) { + mid = (max + min)/2.0; + tangent = curve->TangentAt(mid); + if (NEAR_ZERO(tangent[X], TOL2)) { + return mid; + } + if ((tangent[X] < 0.0) == tanmin) { + min = mid; + } else { + max = mid; + } + } + return min; +} + + +fastf_t +CurveTree::getHorizontalTangent(const ON_Curve *curve, fastf_t min, fastf_t max) +{ + fastf_t mid; + ON_3dVector tangent; + bool tanmin; + + // first lets check end points + tangent = curve->TangentAt(max); + if (NEAR_ZERO(tangent.y, TOL2)) + return max; + tangent = curve->TangentAt(min); + if (NEAR_ZERO(tangent.y, TOL2)) + return min; + + tanmin = (tangent[Y] < 0.0); + while ((max-min) > TOL2) { + mid = (max + min)/2.0; + tangent = curve->TangentAt(mid); + if (NEAR_ZERO(tangent[Y], TOL2)) { + return mid; + } + if ((tangent[Y] < 0.0) == tanmin) { + min = mid; + } else { + max = mid; + } + } + return min; +} + + +bool +CurveTree::getHVTangents(const ON_Curve* curve, ON_Interval& t, std::list<fastf_t>& list) +{ + bool tanx1, tanx2, tanx_changed; + bool tany1, tany2, tany_changed; + bool tan_changed; + ON_3dVector tangent1, tangent2; + ON_3dPoint p1, p2; + + tangent1 = curve->TangentAt(t[0]); + tangent2 = curve->TangentAt(t[1]); + + tanx1 = (tangent1[X] < 0.0); + tanx2 = (tangent2[X] < 0.0); + tany1 = (tangent1[Y] < 0.0); + tany2 = (tangent2[Y] < 0.0); + + tanx_changed =(tanx1 != tanx2); + tany_changed =(tany1 != tany2); + + tan_changed = tanx_changed || tany_changed; + + if (tan_changed) { + if (tanx_changed && tany_changed) {//horz & vert simply split + double midpoint = (t[1]+t[0])/2.0; + ON_Interval left(t[0], midpoint); + ON_Interval right(midpoint, t[1]); + getHVTangents(curve, left, list); + getHVTangents(curve, right, list); + return true; + } else if (tanx_changed) {//find horz + double x = getVerticalTangent(curve, t[0], t[1]); + list.push_back(x); + } else { //find vert + double x = getHorizontalTangent(curve, t[0], t[1]); + list.push_back(x); + } + } else { // check point slope for change + bool slopex, slopex_changed; + bool slopey, slopey_changed; + bool slope_changed; + fastf_t xdelta, ydelta; + + p1 = curve->PointAt(t[0]); + p2 = curve->PointAt(t[1]); + + xdelta = (p2[X] - p1[X]); + slopex = (xdelta < 0.0); + ydelta = (p2[Y] - p1[Y]); + slopey = (ydelta < 0.0); + + if (NEAR_ZERO(xdelta, TOL) || + NEAR_ZERO(ydelta, TOL)) { + return true; + } + + slopex_changed = (slopex != tanx1); + slopey_changed = (slopey != tany1); + + slope_changed = slopex_changed || slopey_changed; + + if (slope_changed) { //2 horz or 2 vert changes simply split + double midpoint = (t[1]+t[0])/2.0; + ON_Interval left(t[0], midpoint); + ON_Interval right(midpoint, t[1]); + getHVTangents(curve, left, list); + getHVTangents(curve, right, list); + return true; + } + } + return true; +} + + +BRNode* +CurveTree::curveBBox(const ON_Curve* curve, int adj_face_index, ON_Interval& t, bool isLeaf, bool innerTrim, const ON_BoundingBox& bb) +{ + BRNode* node; + int vdot = 1; + + if (isLeaf) { + TRACE("creating leaf: u(" << u.Min() << ", " << u.Max() << ") v(" << v.Min() << ", " << v.Max() << ")"); + node = new BRNode(curve, adj_face_index, bb, m_face, t, vdot, innerTrim); + } else { + node = new BRNode(bb); + } + + return node; + +} + + +BRNode* +CurveTree::initialLoopBBox() +{ + ON_BoundingBox bb; + for (int i = 0; i < m_face->LoopCount(); i++) { + ON_BrepLoop* loop = m_face->Loop(i); + if (loop->m_type == ON_BrepLoop::outer) { + if (loop->GetBBox(bb[0], bb[1], 0)) { + TRACE("BBox for Loop min<" << bb[0][0] << ", " << bb[0][1] ", " << bb[0][2] << ">"); + TRACE("BBox for Loop max<" << bb[1][0] << ", " << bb[1][1] ", " << bb[1][2] << ">"); + } + break; + } + } + BRNode* node = new BRNode(bb); + return node; +} + + +BRNode* +CurveTree::subdivideCurve(const ON_Curve* curve, int adj_face_index, double min, double max, bool innerTrim, int divDepth) +{ + ON_Interval dom = curve->Domain(); + ON_3dPoint points[2]; + points[0] = curve->PointAt(min); + points[1] = curve->PointAt(max); + point_t minpt, maxpt; + VSETALL(minpt, MAX_FASTF); + VSETALL(maxpt, -MAX_FASTF); + for (int i = 0; i < 2; i++) + VMINMAX(minpt, maxpt, ((double*)points[i])); + points[0]=ON_3dPoint(minpt); + points[1]=ON_3dPoint(maxpt); + ON_BoundingBox bb(points[0], points[1]); + + ON_Interval t(min, max); + if (isLinear(curve, min, max) || divDepth >= BREP_MAX_LN_DEPTH) { + double delta = (max - min)/(BREP_BB_CRV_PNT_CNT-1); + point_t pnts[BREP_BB_CRV_PNT_CNT]; + ON_3dPoint pnt; + for (int i=0;i<BREP_BB_CRV_PNT_CNT-1;i++) { + pnt = curve->PointAt(min + delta*i); + VSET(pnts[i], pnt[0], pnt[1], pnt[2]); + } + pnt = curve->PointAt(max); + VSET(pnts[BREP_BB_CRV_PNT_CNT-1], pnt[0], pnt[1], pnt[2]); + + VSETALL(minpt, MAX_FASTF); + VSETALL(maxpt, -MAX_FASTF); + for (int i = 0; i < BREP_BB_CRV_PNT_CNT; i++) + VMINMAX(minpt, maxpt, ((double*)pnts[i])); + + VMOVE(pnt, minpt); + bb.Set(pnt, false); + VMOVE(pnt, maxpt); + bb.Set(pnt, true); + return curveBBox(curve, adj_face_index, t, true, innerTrim, bb); + } + + // else subdivide + BRNode* parent = curveBBox(curve, adj_face_index, t, false, innerTrim, bb); + double mid = (max+min)/2.0; + BRNode* l = subdivideCurve(curve, adj_face_index, min, mid, innerTrim, divDepth+1); + BRNode* r = subdivideCurve(curve, adj_face_index, mid, max, innerTrim, divDepth+1); + parent->addChild(l); + parent->addChild(r); + return parent; +} + + +/** + * Determine whether a given curve segment is linear + */ +bool +CurveTree::isLinear(const ON_Curve* curve, double min, double max) +{ + ON_3dVector tangent_start = curve->TangentAt(min); + ON_3dVector tangent_end = curve->TangentAt(max); + double vdot = tangent_start * tangent_end; + if (vdot < BREP_CURVE_FLATNESS) + return false; + + ON_3dPoint pmin = curve->PointAt(min); + ON_3dPoint pmax = curve->PointAt(max); + + const ON_Surface* surf = m_face->SurfaceOf(); + ON_Interval u = surf->Domain(0); + ON_Interval v = surf->Domain(1); + point_t a, b; + VSET(a, u[0], v[0], 0.0); + VSET(b, u[1], v[1], 0.0); + double dd = DIST_PT_PT(a, b); + double cd = DIST_PT_PT(pmin, pmax); + + if (cd > BREP_TRIM_SUB_FACTOR*dd) + return false; + + double delta = (max - min)/(BREP_BB_CRV_PNT_CNT-1); + ON_3dPoint points[BREP_BB_CRV_PNT_CNT]; + for (int i=0;i<BREP_BB_CRV_PNT_CNT-1;i++) { + points[i] = curve->PointAt(min + delta*i); + } + points[BREP_BB_CRV_PNT_CNT-1] = curve->PointAt(max); + + ON_3dVector A; + ON_3dVector B; + vdot = 1.0; + A = points[BREP_BB_CRV_PNT_CNT-1] - points[0]; + A.Unitize(); + for (int i=1;i<BREP_BB_CRV_PNT_CNT-1;i++) { + B = points[i] - points[0]; + B.Unitize(); + vdot = vdot * (A * B); + if (vdot < BREP_CURVE_FLATNESS) + return false; //already failed + } + + return vdot >= BREP_CURVE_FLATNESS; +} + + +//-------------------------------------------------------------------------------- +// SurfaceTree +SurfaceTree::SurfaceTree(ON_BrepFace* face, bool removeTrimmed, int depthLimit) + : m_removeTrimmed(removeTrimmed), + m_face(face) +{ + // first, build the Curve Tree + if (removeTrimmed) + ctree = new CurveTree(m_face); + else + ctree = NULL; + + // build the surface bounding volume hierarchy + const ON_Surface* surf = face->SurfaceOf(); + TRACE("Creating surface tree for: " << face->m_face_index); + ON_Interval u = surf->Domain(0); + ON_Interval v = surf->Domain(1); + // Populate initial corner and normal arrays for use in + // tree build + ON_3dPoint corners[9]; + ON_3dVector normals[9]; + double uq = u.Length()*0.25; + double vq = v.Length()*0.25; + + /////////////////////////////////////////////////////////////////////// + ON_Plane frames[9]; + surf->FrameAt(u.Min(), v.Min(), frames[0]); + surf->FrameAt(u.Max(), v.Min(), frames[1]); + surf->FrameAt(u.Max(), v.Max(), frames[2]); + surf->FrameAt(u.Min(), v.Max(), frames[3]); + surf->FrameAt(u.Mid(), v.Mid(), frames[4]); + surf->FrameAt(u.Mid() - uq, v.Mid() - vq, frames[5]); + surf->FrameAt(u.Mid() - uq, v.Mid() + vq, frames[6]); + surf->FrameAt(u.Mid() + uq, v.Mid() - vq, frames[7]); + surf->FrameAt(u.Mid() + uq, v.Mid() + vq, frames[8]); + + corners[0] = frames[0].origin; + normals[0] = frames[0].zaxis; + corners[1] = frames[1].origin; + normals[1] = frames[1].zaxis; + corners[2] = frames[2].origin; + normals[2] = frames[2].zaxis; + corners[3] = frames[3].origin; + normals[3] = frames[3].zaxis; + corners[4] = frames[4].origin; + normals[4] = frames[4].zaxis; + corners[5] = frames[5].origin; + normals[5] = frames[5].zaxis; + corners[6] = frames[6].origin; + normals[6] = frames[6].zaxis; + corners[7] = frames[7].origin; + normals[7] = frames[7].zaxis; + corners[8] = frames[8].origin; + normals[8] = frames[8].zaxis; + + m_root = subdivideSurfaceByKnots(surf, u, v, frames, corners, normals, 0, depthLimit); + if (m_root) { + m_root->BuildBBox(); + } + TRACE("u: [" << u[0] << ", " << u[1] << "]"); + TRACE("v: [" << v[0] << ", " << v[1] << "]"); + TRACE("m_root: " << m_root); +} + + +SurfaceTree::~SurfaceTree() +{ + delete ctree; + delete m_root; +} + + +BBNode* +SurfaceTree::getRootNode() const +{ + return m_root; +} + + +int +SurfaceTree::depth() +{ + return m_root->depth(); +} + + +ON_2dPoint +SurfaceTree::getClosestPointEstimate(const ON_3dPoint& pt) +{ + return m_root->getClosestPointEstimate(pt); +} + + +ON_2dPoint +SurfaceTree::getClosestPointEstimate(const ON_3dPoint& pt, ON_Interval& u, ON_Interval& v) +{ + return m_root->getClosestPointEstimate(pt, u, v); +} + + +void +SurfaceTree::getLeaves(std::list<BBNode*>& out_leaves) +{ + m_root->getLeaves(out_leaves); +} + + +const ON_Surface * +SurfaceTree::getSurface() +{ + return m_face->SurfaceOf(); +} + + +int +SurfaceTree::getSurfacePoint(const ON_3dPoint& pt, ON_2dPoint& uv, const ON_3dPoint& from, double tolerance) const +{ + std::list<BBNode*> nodes; + (void)m_root->getLeavesBoundingPoint(from, nodes); + + double min_dist = MAX_FASTF; + ON_2dPoint curr_uv(0.0, 0.0); + bool found = false; + + std::list<BBNode*>::iterator i; + for (i = nodes.begin(); i != nodes.end(); i++) { + BBNode* node = (*i); + if (::brep_getSurfacePoint(pt, curr_uv, node)) { + ON_3dPoint fp = m_face->SurfaceOf()->PointAt(curr_uv.x, curr_uv.y); + double dist = fp.DistanceTo(pt); + if (NEAR_ZERO(dist, BREP_SAME_POINT_TOLERANCE)) { + uv = curr_uv; + found = true; + return 1; //close enough to same point so no sense in looking for one closer + } else if (NEAR_ZERO(dist, tolerance)) { + if (dist < min_dist) { + uv = curr_uv; + min_dist = dist; + found = true; //within tolerance but may be a point closer so keep looking + } + } + } + } + if (found) { + return 1; + } + return -1; +} + + +//static int bb_cnt=0; +BBNode* +SurfaceTree::surfaceBBox(const ON_Surface *localsurf, bool isLeaf, ON_3dPoint *m_corners, ON_3dVector *m_normals, const ON_Interval& u, const ON_Interval& v) +{ + point_t min, max, buffer; + ON_BoundingBox bbox = localsurf->BoundingBox(); + + VSETALL(buffer, BREP_EDGE_MISS_TOLERANCE); + //VSETALL(min, MAX_FASTF); + //VSETALL(max, -MAX_FASTF); + //for (int i = 0; i < 9; i++) { + //VMINMAX(min, max, ((double*)m_corners[i])); + //bu_log("in c%d_%d sph %f %f %f %f\n", bb_cnt, i, m_corners[i].x, m_corners[i].y, m_corners[i].z, 1.0); + //} + + //bu_log("in bb%d rpp %f %f %f %f %f %f\n", bb_cnt, min[0], max[0], min[1], max[1], min[2], max[2]); + VMOVE(min, bbox.Min()); + VMOVE(max, bbox.Max()); + + //bu_log("in b%d rpp %f %f %f %f %f %f\n", bb_cnt++, bbox.m_min.x, bbox.m_max.x, bbox.m_min.y, bbox.m_max.y, bbox.m_min.z, bbox.m_max.z); + //bu_log("in bc%d rpp %f %f %f %f %f %f\n", bb_cnt++, min[0], max[0], min[1], max[1], min[2], max[2]); + + // calculate the estimate point on the surface: i.e. use the point + // on the surface defined by (u.Mid(), v.Mid()) as a heuristic for + // finding the uv domain bounding a portion of the surface close + // to an arbitrary model-space point (see + // getClosestPointEstimate()) + ON_3dPoint estimate; + ON_3dVector normal; + estimate = m_corners[4]; + normal = m_normals[4]; + + BBNode* node; + if (isLeaf) { +/* vect_t delta; */ + + VSUB2(min, min, buffer); + VADD2(max, max, buffer); +/* VSUB2(delta, max, min); + VSCALE(delta, delta, BBOX_GROW_3D); + VSUB2(min, min, delta); + VADD2(max, max, delta); +*/ + TRACE("creating leaf: u(" << u.Min() << ", " << u.Max() << + ") v(" << v.Min() << ", " << v.Max() << ")"); + node = new BBNode(ctree, ON_BoundingBox(ON_3dPoint(min), + ON_3dPoint(max)), + m_face, + u, v); + node->prepTrims(); + + } else { + node = new BBNode(ctree, ON_BoundingBox(ON_3dPoint(min), + ON_3dPoint(max))); + } + + node->m_estimate = estimate; + node->m_normal = normal; + node->m_face = m_face; + node->m_u = u; + node->m_v = v; + return node; +} + + +BBNode* +initialBBox(CurveTree* ctree, const ON_Surface* surf, const ON_BrepFace* face, const ON_Interval& u, const ON_Interval& v) +{ + ON_BoundingBox bb = surf->BoundingBox(); + BBNode* node = new BBNode(ctree, bb, face, u, v, false, false); + ON_3dPoint estimate; + ON_3dVector normal; + if (!surf->EvNormal(surf->Domain(0).Mid(), surf->Domain(1).Mid(), estimate, normal)) { + bu_bomb("Could not evaluate estimate point on surface"); + } + node->m_estimate = estimate; + node->m_normal = normal; + node->m_face = face; + node->m_u = u; + node->m_v = v; + return node; +} + + +BBNode* +SurfaceTree::subdivideSurfaceByKnots(const ON_Surface *localsurf, + const ON_Interval& u, + const ON_Interval& v, + ON_Plane frames[], + ON_3dPoint corners[], + ON_3dVector normals[], + int divDepth, + int depthLimit) +{ + const ON_Surface* surf = m_face->SurfaceOf(); + ON_Interval usurf = surf->Domain(0); + ON_Interval vsurf = surf->Domain(1); + BBNode* parent = initialBBox(ctree, localsurf, m_face, u, v); + BBNode* quads[4]; + int spanu_cnt = localsurf->SpanCount(0); + int spanv_cnt = localsurf->SpanCount(1); + double *spanu = NULL; + double *spanv = NULL; + spanu = new double[spanu_cnt+1]; + spanv = new double[spanv_cnt+1]; + localsurf->GetSpanVector(0, spanu); + localsurf->GetSpanVector(1, spanv); + + /////////////////////////////////// + double uq = u.Length()*0.25; + double vq = v.Length()*0.25; + surf->FrameAt(u.Mid() - uq, v.Mid() - vq, frames[5]); + surf->FrameAt(u.Mid() - uq, v.Mid() + vq, frames[6]); + surf->FrameAt(u.Mid() + uq, v.Mid() - vq, frames[7]); + surf->FrameAt(u.Mid() + uq, v.Mid() + vq, frames[8]); + + corners[5] = frames[5].origin; + normals[5] = frames[5].zaxis; + corners[6] = frames[6].origin; + normals[6] = frames[6].zaxis; + corners[7] = frames[7].origin; + normals[7] = frames[7].zaxis; + corners[8] = frames[8].origin; + normals[8] = frames[8].zaxis; + + if ((spanu_cnt > 1) && (spanv_cnt > 1)) { + double usplit = spanu[(spanu_cnt+1)/2]; + double vsplit = spanv[(spanv_cnt+1)/2]; + ON_Interval firstu(u.Min(), usplit); + ON_Interval secondu(usplit, u.Max()); + ON_Interval firstv(v.Min(), vsplit); + ON_Interval secondv(vsplit, v.Max()); + + ON_Surface *north = NULL; + ON_Surface *south = NULL; + ON_Surface *q2surf = NULL; + ON_Surface *q3surf = NULL; + ON_Surface *q1surf = NULL; + ON_Surface *q0surf = NULL; + + ON_BoundingBox box = localsurf->BoundingBox(); + + int dir = 1; + bool split = localsurf->Split(dir, localsurf->Domain(dir).Mid(), south, north); + + /* FIXME: this needs to be handled more gracefully */ + if (!split || !south || !north) { + bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)south, (void *)north); + delete parent; + return NULL; + } + + split = localsurf->Split(dir, vsplit, south, north); + + /* FIXME: this needs to be handled more gracefully */ + if (!split || !south || !north) { + bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)south, (void *)north); + delete parent; + return NULL; + } + + south->ClearBoundingBox(); + north->ClearBoundingBox(); + + dir = 0; + split = south->Split(dir, usplit, q0surf, q1surf); + + /* FIXME: this needs to be handled more gracefully */ + if (!split || !q0surf || !q1surf) { + bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)q0surf, (void *)q1surf); + delete parent; + return NULL; + } + + delete south; + q0surf->ClearBoundingBox(); + q1surf->ClearBoundingBox(); + split = north->Split(dir, usplit, q3surf, q2surf); + + /* FIXME: this needs to be handled more gracefully */ + if (!split || !q3surf || !q2surf) { + bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)q3surf, (void *)q2surf); + delete parent; + return NULL; + } + + delete north; + q3surf->ClearBoundingBox(); + q2surf->ClearBoundingBox(); + + /********************************************************************* + * In order to avoid fairly expensive re-calculation of 3d points at + * uv coordinates, all values that are shared between children at + * the same depth of a surface subdivision are pre-computed and + * passed as paramters. + * + * The majority of these points are already evaluated in the process + * of testing whether a subdivision has produced a leaf node. These + * values are in the normals and corners arrays and have index values + * corresponding to the values of the figure on the left below. There + * are four other shared values that are precomputed in the sharedcorners + * and sharednormals arrays; their index values in those arrays are + * illustrated in the figure on the right: + * + * + * 3-------------------2 +---------2---------+ + * | | | | + * | 6 8 | | | + * | | | | + * V| 4 | 1 3 + * | | | | + * | 5 7 | | | + * | | | | + * 0-------------------1 +---------0---------+ + * U U + * + * Values inherited from Values pre-prepared in + * parent subdivision shared arrays + * + * + * When the four subdivisions are made, the parent parameters + * are passed to the children as follows (values from the + * shared arrays are prefaced with an s): + * + * 3--------------S2 S2--------------2 + * | | | | + * | | | | + * V | 6 | | 8 | + * | | | | + * | | | | + * S1--------------4 4--------------S3 + * U U + * + * Quadrant 3 Quadrant 2 + * + * S1--------------4 4--------------S3 + * | | | | + * | | | | + * V | 5 | | 7 | + * | | | | + * | | | | + * 0--------------S0 S0--------------1 + * U U + * + * Quadrant 0 Quadrant 1 + * + * + **********************************************************************/ + + ON_Plane sharedframes[4]; + ON_3dPoint sharedcorners[4]; + ON_3dVector sharednormals[4]; + localsurf->FrameAt(usplit, v.Min(), sharedframes[0]); + localsurf->FrameAt(u.Min(), vsplit, sharedframes[1]); + localsurf->FrameAt(usplit, v.Max(), sharedframes[2]); + localsurf->FrameAt(u.Max(), vsplit, sharedframes[3]); + localsurf->FrameAt(usplit, vsplit, frames[4]); + + sharedcorners[0] = sharedframes[0].origin; + sharednormals[0] = sharedframes[0].zaxis; + sharedcorners[1] = sharedframes[1].origin; + sharednormals[1] = sharedframes[1].zaxis; + sharedcorners[2] = sharedframes[2].origin; + sharednormals[2] = sharedframes[2].zaxis; + sharedcorners[3] = sharedframes[3].origin; + sharednormals[3] = sharedframes[3].zaxis; + corners[4] = frames[4].origin; + normals[4] = frames[4].zaxis; + + ON_Plane *newframes; + ON_3dPoint *newcorners; + ON_3dVector *newnormals; + newframes = (ON_Plane *)bu_malloc(9*sizeof(ON_Plane), "new frames"); + newcorners = (ON_3dPoint *)bu_malloc(9*sizeof(ON_3dPoint), "new corners"); + newnormals = (ON_3dVector *)bu_malloc(9*sizeof(ON_3dVector), "new normals"); + newframes[0] = frames[0]; + newcorners[0] = corners[0]; + newnormals[0] = normals[0]; + newframes[1] = sharedframes[0]; + newcorners[1] = sharedcorners[0]; + newnormals[1] = sharednormals[0]; + newframes[2] = frames[4]; + newcorners[2] = corners[4]; + newnormals[2] = normals[4]; + newframes[3] = sharedframes[1]; + newcorners[3] = sharedcorners[1]; + newnormals[3] = sharednormals[1]; + newframes[4] = frames[5]; + newcorners[4] = corners[5]; + newnormals[4] = normals[5]; + quads[0] = subdivideSurfaceByKnots(q0surf, firstu, firstv, newframes, newcorners, newnormals, divDepth+1, depthLimit); + delete q0surf; + newframes[0] = sharedframes[0]; + newcorners[0] = sharedcorners[0]; + newnormals[0] = sharednormals[0]; + newframes[1] = frames[1]; + newcorners[1] = corners[1]; + newnormals[1] = normals[1]; + newframes[2] = sharedframes[3]; + newcorners[2] = sharedcorners[3]; + newnormals[2] = sharednormals[3]; + newframes[3] = frames[4]; + newcorners[3] = corners[4]; + newnormals[3] = normals[4]; + newframes[4] = frames[7]; + newcorners[4] = corners[7]; + newnormals[4] = normals[7]; + quads[1] = subdivideSurfaceByKnots(q1surf, secondu, firstv, newframes, newcorners, newnormals, divDepth+1, depthLimit); + delete q1surf; + newframes[0] = frames[4]; + newcorners[0] = corners[4]; + newnormals[0] = normals[4]; + newframes[1] = sharedframes[3]; + newcorners[1] = sharedcorners[3]; + newnormals[1] = sharednormals[3]; + newframes[2] = frames[2]; + newcorners[2] = corners[2]; + newnormals[2] = normals[2]; + newframes[3] = sharedframes[2]; + newcorners[3] = sharedcorners[2]; + newnormals[3] = sharednormals[2]; + newframes[4] = frames[8]; + newcorners[4] = corners[8]; + newnormals[4] = normals[8]; + quads[2] = subdivideSurfaceByKnots(q2surf, secondu, secondv, newframes, newcorners, newnormals, divDepth+1, depthLimit); + delete q2surf; + newframes[0] = sharedframes[1]; + newcorners[0] = sharedcorners[1]; + newnormals[0] = sharednormals[1]; + newframes[1] = frames[4]; + newcorners[1] = corners[4]; + newnormals[1] = normals[4]; + newframes[2] = sharedframes[2]; + newcorners[2] = sharedcorners[2]; + newnormals[2] = sharednormals[2]; + newframes[3] = frames[3]; + newcorners[3] = corners[3]; + newnormals[3] = normals[3]; + newframes[4] = frames[6]; + newcorners[4] = corners[6]; + newnormals[4] = normals[6]; + quads[3] = subdivideSurfaceByKnots(q3surf, firstu, secondv, newframes, newcorners, newnormals, divDepth+1, depthLimit); + delete q3surf; + bu_free(newframes, "free subsurface frames array"); + bu_free(newcorners, "free subsurface corners array"); + bu_free(newnormals, "free subsurface normals array"); + + parent->m_trimmed = true; + parent->m_checkTrim = false; + + for (int i = 0; i < 4; i++) { + if (!quads[i]) { + continue; + } + if (!(quads[i]->m_trimmed)) { + parent->m_trimmed = false; + } + if (quads[i]->m_checkTrim) { + parent->m_checkTrim = true; + } + } + if (m_removeTrimmed) { + for (int i = 0; i < 4; i++) { + if (!quads[i]) { + continue; + } + if (!(quads[i]->m_trimmed)) { + parent->addChild(quads[i]); + } + } + } else { + for (int i = 0; i < 4; i++) { + parent->addChild(quads[i]); + } + } + } else if (spanu_cnt > 1) { + double usplit = spanu[(spanu_cnt+1)/2]; + ON_Interval firstu(u.Min(), usplit); + ON_Interval secondu(usplit, u.Max()); + ////////////////////////////////////// + ON_Surface *east = NULL; + ON_Surface *west = NULL; + + ON_BoundingBox box = localsurf->BoundingBox(); + + int dir = 0; + bool split = localsurf->Split(dir, usplit, east, west); + + /* FIXME: this needs to be handled more gracefully */ + if (!split || !east || !west) { + bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)east, (void *)west); + delete parent; + return NULL; + } + + east->ClearBoundingBox(); + west->ClearBoundingBox(); + + ////////////////////////////////// + /********************************************************************* + * In order to avoid fairly expensive re-calculation of 3d points at + * uv coordinates, all values that are shared between children at + * the same depth of a surface subdivision are pre-computed and + * passed as paramters. + * + * The majority of these points are already evaluated in the process + * of testing whether a subdivision has produced a leaf node. These + * values are in the normals and corners arrays and have index values + * corresponding to the values of the figure on the left below. There + * are four other shared values that are precomputed in the sharedcorners + * and sharednormals arrays; their index values in those arrays are + * illustrated in the figure on the right: + * + * + * 3-------------------2 +---------2---------+ + * | | | | + * | 6 8 | | | + * | | | | + * V| 4 | 1 3 + * | | | | + * | 5 7 | | | + * | | | | + * 0-------------------1 +---------0---------+ + * U U + * + * Values inherited from Values pre-prepared in + * parent subdivision shared arrays + * + * + * When the four subdivisions are made, the parent parameters + * are passed to the children as follows (values from the + * shared arrays are prefaced with an s): + * + * 3--------------S1 S1--------------2 + * | | | | + * | | | | + * V | 4 | | 5 | + * | | | | + * | | | | + * 0--------------S0 S0--------------1 + * U U + * + * East West + * + * + **********************************************************************/ + /* When the four subdivisions are made, the parent parameters + * are passed to the children as follows (values from the + * shared arrays are prefaced with an s): + * + * 3--------------S2 S2--------------2 + * | | | | + * | | | | + * V | 6 | | 8 | + * | | | | + * | | | | + * S1--------------4 4--------------S3 + * U U + * + * Quadrant 3 Quadrant 2 + * + * S1--------------4 4--------------S3 + * | | | | + * | | | | + * V | 5 | | 7 | + * | | | | + * | | | | + * 0--------------S0 S0--------------1 + * U U + * + * Quadrant 0 Quadrant 1 + * + * + **********************************************************************/ + + ON_Plane sharedframes[4]; + ON_3dPoint sharedcorners[4]; + ON_3dVector sharednormals[4]; + localsurf->FrameAt(usplit, v.Min(), sharedframes[0]); + localsurf->FrameAt(usplit, v.Max(), sharedframes[1]); + + sharedcorners[0] = sharedframes[0].origin; + sharednormals[0] = sharedframes[0].zaxis; + sharedcorners[1] = sharedframes[1].origin; + sharednormals[1] = sharedframes[1].zaxis; + + ON_Plane *newframes; + ON_3dPoint *newcorners; + ON_3dVector *newnormals; + newframes = (ON_Plane *) bu_malloc(9 * sizeof(ON_Plane), + "new frames"); + newcorners = (ON_3dPoint *) bu_malloc(9 * sizeof(ON_3dPoint), + "new corners"); + newnormals = (ON_3dVector *) bu_malloc(9 * sizeof(ON_3dVector), + "new normals"); + newframes[0] = frames[0]; + newcorners[0] = corners[0]; + newnormals[0] = normals[0]; + newframes[1] = sharedframes[0]; + newcorners[1] = sharedcorners[0]; + newnormals[1] = sharednormals[0]; + newframes[2] = sharedframes[1]; + newcorners[2] = sharedcorners[1]; + newnormals[2] = sharednormals[1]; + newframes[3] = frames[3]; + newcorners[3] = corners[3]; + newnormals[3] = normals[3]; + localsurf->FrameAt(firstu.Mid(), v.Mid(), newframes[4]); + + newcorners[4] = newframes[4].origin; + newnormals[4] = newframes[4].zaxis; + + //ON_BoundingBox bbox = q0surf->BoundingBox(); + //bu_log("%d - in bbq0 rpp %f %f %f %f %f %f\n", divDepth, bbox.m_min.x, bbox.m_max.x, bbox.m_min.y, bbox.m_max.y, bbox.m_min.z, bbox.m_max.z); + quads[0] = subdivideSurfaceByKnots(east, firstu, v, newframes, newcorners, newnormals, divDepth+1, depthLimit); + delete east; + + newframes[0] = sharedframes[0]; + newcorners[0] = sharedcorners[0]; + newnormals[0] = sharednormals[0]; + newframes[1] = frames[1]; + newcorners[1] = corners[1]; + newnormals[1] = normals[1]; + newframes[2] = frames[2]; + newcorners[2] = corners[2]; + newnormals[2] = normals[2]; + newframes[3] = sharedframes[1]; + newcorners[3] = sharedcorners[1]; + newnormals[3] = sharednormals[1]; + localsurf->FrameAt(secondu.Mid(), v.Mid(), newframes[4]); + + newcorners[4] = newframes[4].origin; + newnormals[4] = newframes[4].zaxis; + + //bbox = q1surf->BoundingBox(); + //bu_log("%d - in bbq1 rpp %f %f %f %f %f %f\n", divDepth, bbox.m_min.x, bbox.m_max.x, bbox.m_min.y, bbox.m_max.y, bbox.m_min.z, bbox.m_max.z); + quads[1] = subdivideSurfaceByKnots(west, secondu, v, newframes, newcorners, newnormals, divDepth+1, depthLimit); + delete west; + + bu_free(newframes, "free subsurface frames array"); + bu_free(newcorners, "free subsurface corners array"); + bu_free(newnormals, "free subsurface normals array"); + + parent->m_trimmed = true; + parent->m_checkTrim = false; + + for (int i = 0; i < 2; i++) { + if (!quads[i]) { + continue; + } + if (!(quads[i]->m_trimmed)) { + parent->m_trimmed = false; + } + if (quads[i]->m_checkTrim) { + parent->m_checkTrim = true; + } + } + if (m_removeTrimmed) { + for (int i = 0; i < 2; i++) { + if (!quads[i]) { + continue; + } + if (!(quads[i]->m_trimmed)) { + parent->addChild(quads[i]); + } + } + } else { + for (int i = 0; i < 2; i++) { + parent->addChild(quads[i]); + } + } + } else if (spanv_cnt > 1) { + double vsplit = spanv[(spanv_cnt+1)/2]; + ON_Interval firstv(v.Min(), vsplit); + ON_Interval secondv(vsplit, v.Max()); + + ////////////////////////////////////// + ON_Surface *north = NULL; + ON_Surface *south = NULL; + + ON_BoundingBox box = localsurf->BoundingBox(); + + int dir = 1; + bool split = localsurf->Split(dir, vsplit, south, north); + + /* FIXME: this needs to be handled more gracefully */ + if (!split || !south || !north) { + bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)south, (void *)north); + delete parent; + return NULL; + } + + south->ClearBoundingBox(); + north->ClearBoundingBox(); + + ////////////////////////////////// + /********************************************************************* + * In order to avoid fairly expensive re-calculation of 3d points at + * uv coordinates, all values that are shared between children at + * the same depth of a surface subdivision are pre-computed and + * passed as paramters. + * + * The majority of these points are already evaluated in the process + * of testing whether a subdivision has produced a leaf node. These + * values are in the normals and corners arrays and have index values + * corresponding to the values of the figure on the left below. There + * are four other shared values that are precomputed in the sharedcorners + * and sharednormals arrays; their index values in those arrays are + * illustrated in the figure on the right: + * + * + * 3-------------------2 +---------2---------+ + * | | | | + * | 6 8 | | | + * | | | | + * V| 4 | 1 3 + * | | | | + * | 5 7 | | | + * | | | | + * 0-------------------1 +---------0---------+ + * U U + * + * Values inherited from Values pre-prepared in + * parent subdivision shared arrays + * + * + * When the four subdivisions are made, the parent parameters + * are passed to the children as follows (values from the + * shared arrays are prefaced with an s): + * + * 3--------------------2 + * | | + * | | + * V | 5 | + * | | + * | | + * S0-------------------S1 + * U + * + * North + * + * S0-------------------S1 + * | | + * | | + * V | 4 | + * | | + * | | + * 0--------------------1 + * U + * + * South + * + * + **********************************************************************/ + + ON_Plane sharedframes[2]; + ON_3dPoint sharedcorners[2]; + ON_3dVector sharednormals[2]; + localsurf->FrameAt(u.Min(), vsplit, sharedframes[0]); + localsurf->FrameAt(u.Max(), vsplit, sharedframes[1]); + + sharedcorners[0] = sharedframes[0].origin; + sharednormals[0] = sharedframes[0].zaxis; + sharedcorners[1] = sharedframes[1].origin; + sharednormals[1] = sharedframes[1].zaxis; + + ON_Plane *newframes; + ON_3dPoint *newcorners; + ON_3dVector *newnormals; + newframes = (ON_Plane *) bu_malloc(9 * sizeof(ON_Plane), + "new frames"); + newcorners = (ON_3dPoint *) bu_malloc(9 * sizeof(ON_3dPoint), + "new corners"); + newnormals = (ON_3dVector *) bu_malloc(9 * sizeof(ON_3dVector), + "new normals"); + newframes[0] = frames[0]; + newcorners[0] = corners[0]; + newnormals[0] = normals[0]; + newframes[1] = frames[1]; + newcorners[1] = corners[1]; + newnormals[1] = normals[1]; + newframes[2] = sharedframes[1]; + newcorners[2] = sharedcorners[1]; + newnormals[2] = sharednormals[1]; + newframes[3] = sharedframes[0]; + newcorners[3] = sharedcorners[0]; + newnormals[3] = sharednormals[0]; + localsurf->FrameAt(u.Mid(), firstv.Mid(), newframes[4]); + + newcorners[4] = newframes[4].origin; + newnormals[4] = newframes[4].zaxis; + //ON_BoundingBox bbox = q0surf->BoundingBox(); + //bu_log("%d - in bbq0 rpp %f %f %f %f %f %f\n", divDepth, bbox.m_min.x, bbox.m_max.x, bbox.m_min.y, bbox.m_max.y, bbox.m_min.z, bbox.m_max.z); + quads[0] = subdivideSurfaceByKnots(south, u, firstv, newframes, newcorners, newnormals, divDepth+1, depthLimit); + delete south; + + newframes[0] = sharedframes[0]; + newcorners[0] = sharedcorners[0]; + newnormals[0] = sharednormals[0]; + newframes[1] = sharedframes[1]; + newcorners[1] = sharedcorners[1]; + newnormals[1] = sharednormals[1]; + newframes[2] = frames[2]; + newcorners[2] = corners[2]; + newnormals[2] = normals[2]; + newframes[3] = frames[3]; + newcorners[3] = corners[3]; + newnormals[3] = normals[3]; + localsurf->FrameAt(u.Mid(), secondv.Mid(), newframes[4]); + + newcorners[4] = newframes[4].origin; + newnormals[4] = newframes[4].zaxis; + //bbox = q1surf->BoundingBox(); + //bu_log("%d - in bbq1 rpp %f %f %f %f %f %f\n", divDepth, bbox.m_min.x, bbox.m_max.x, bbox.m_min.y, bbox.m_max.y, bbox.m_min.z, bbox.m_max.z); + quads[1] = subdivideSurfaceByKnots(north, u, secondv, newframes, newcorners, newnormals, divDepth+1, depthLimit); + delete north; + + bu_free(newframes, "free subsurface frames array"); + bu_free(newcorners, "free subsurface corners array"); + bu_free(newnormals, "free subsurface normals array"); + + parent->m_trimmed = true; + parent->m_checkTrim = false; + + for (int i = 0; i < 2; i++) { + if (!quads[i]) { + continue; + } + if (!(quads[i]->m_trimmed)) { + parent->m_trimmed = false; + } + if (quads[i]->m_checkTrim) { + parent->m_checkTrim = true; + } + } + if (m_removeTrimmed) { + for (int i = 0; i < 2; i++) { + if (!quads[i]) { + continue; + } + if (!(quads[i]->m_trimmed)) { + parent->addChild(quads[i]); + } + } + } else { + for (int i = 0; i < 2; i++) { + parent->addChild(quads[i]); + } + } + } else { + //return surfaceBBox(localsurf, true, corners, normals, u, v); + //parent->addChild(subdivideSurface(localsurf, u, v, frames, corners, normals, 0)); + ((ON_Surface *)localsurf)->ClearBoundingBox(); + delete parent; + return subdivideSurface(localsurf, u, v, frames, corners, normals, 0, depthLimit); + } + delete [] spanu; + delete [] spanv; + + parent->BuildBBox(); + return parent; + +} + + +BBNode* +SurfaceTree::subdivideSurface(const ON_Surface *localsurf, + const ON_Interval& u, + const ON_Interval& v, + ON_Plane frames[], + ON_3dPoint corners[], + ON_3dVector normals[], + int divDepth, + int depthLimit) +{ + const ON_Surface* surf = m_face->SurfaceOf(); + ON_Interval usurf = surf->Domain(0); + ON_Interval vsurf = surf->Domain(1); + + double uq = u.Length()*0.25; + double vq = v.Length()*0.25; + localsurf->FrameAt(u.Mid() - uq, v.Mid() - vq, frames[5]); + localsurf->FrameAt(u.Mid() - uq, v.Mid() + vq, frames[6]); + localsurf->FrameAt(u.Mid() + uq, v.Mid() - vq, frames[7]); + localsurf->FrameAt(u.Mid() + uq, v.Mid() + vq, frames[8]); + + corners[5] = frames[5].origin; + normals[5] = frames[5].zaxis; + corners[6] = frames[6].origin; + normals[6] = frames[6].zaxis; + corners[7] = frames[7].origin; + normals[7] = frames[7].zaxis; + corners[8] = frames[8].origin; + normals[8] = frames[8].zaxis; + + double width, height; + double ratio = 5.0; + localsurf->GetSurfaceSize(&width, &height); + if (((width/height < ratio) && (width/height > 1.0/ratio) && isFlat(frames) && isStraight(frames)) + || (divDepth >= depthLimit)) { //BREP_MAX_FT_DEPTH))) { + return surfaceBBox(localsurf, true, corners, normals, u, v); + } else { + bool isUFlat = isFlatU(frames); + bool isVFlat = isFlatV(frames); + + BBNode* parent = (divDepth == 0) ? initialBBox(ctree, localsurf, m_face, u, v) : + surfaceBBox(localsurf, false, corners, normals, u, v); + BBNode* quads[4]; + ON_Interval first(0, 0.5); + ON_Interval second(0.5, 1.0); + + if ((!isVFlat || (width/height > ratio)) && (!isUFlat || (height/width > ratio))) { + ON_Surface *north = NULL; + ON_Surface *south = NULL; + ON_Surface *q2surf = NULL; + ON_Surface *q3surf = NULL; + ON_Surface *q1surf = NULL; + ON_Surface *q0surf = NULL; + + ON_BoundingBox box = localsurf->BoundingBox(); + + int dir = 1; + bool split = localsurf->Split(dir, localsurf->Domain(dir).Mid(), south, north); + + /* FIXME: this needs to be handled more gracefully */ + if (!split || !south || !north) { + bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)south, (void *)north); + delete parent; + return NULL; + } + + south->ClearBoundingBox(); + north->ClearBoundingBox(); + + dir = 0; + split = south->Split(dir, south->Domain(dir).Mid(), q0surf, q1surf); + + /* FIXME: this needs to be handled more gracefully */ + if (!split || !q0surf || !q1surf) { + bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)q0surf, (void *)q1surf); + delete parent; + return NULL; + } + + delete south; + q0surf->ClearBoundingBox(); + q1surf->ClearBoundingBox(); + split = north->Split(dir, north->Domain(dir).Mid(), q3surf, q2surf); + + /* FIXME: this needs to be handled more gracefully */ + if (!split || !q3surf || !q2surf) { + bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)q3surf, (void *)q2surf); + delete parent; + return NULL; + } + + delete north; + q3surf->ClearBoundingBox(); + q2surf->ClearBoundingBox(); + /********************************************************************* + * In order to avoid fairly expensive re-calculation of 3d points at + * uv coordinates, all values that are shared between children at + * the same depth of a surface subdivision are pre-computed and + * passed as paramters. + * + * The majority of these points are already evaluated in the process + * of testing whether a subdivision has produced a leaf node. These + * values are in the normals and corners arrays and have index values + * corresponding to the values of the figure on the left below. There + * are four other shared values that are precomputed in the sharedcorners + * and sharednormals arrays; their index values in those arrays are + * illustrated in the figure on the right: + * + * + * 3-------------------2 +---------2---------+ + * | | | | + * | 6 8 | | | + * | | | | + * V| 4 | 1 3 + * | | | | + * | 5 7 | | | + * | | | | + * 0-------------------1 +---------0---------+ + * U U + * + * Values inherited from Values pre-prepared in + * parent subdivision shared arrays + * + * + * When the four subdivisions are made, the parent parameters + * are passed to the children as follows (values from the + * shared arrays are prefaced with an s): + * + * 3--------------S2 S2--------------2 + * | | | | + * | | | | + * V | 6 | | 8 | + * | | | | + * | | | | + * S1--------------4 4--------------S3 + * U U + * + * Quadrant 3 Quadrant 2 + * + * S1--------------4 4--------------S3 + * | | | | + * | | | | + * V | 5 | | 7 | + * | | | | + * | | | | + * 0--------------S0 S0--------------1 + * U U + * + * Quadrant 0 Quadrant 1 + * + * + **********************************************************************/ + + ON_Plane sharedframes[4]; + ON_3dPoint sharedcorners[4]; + ON_3dVector sharednormals[4]; + localsurf->FrameAt(u.Mid(), v.Min(), sharedframes[0]); + localsurf->FrameAt(u.Min(), v.Mid(), sharedframes[1]); + localsurf->FrameAt(u.Mid(), v.Max(), sharedframes[2]); + localsurf->FrameAt(u.Max(), v.Mid(), sharedframes[3]); + + sharedcorners[0] = sharedframes[0].origin; + sharednormals[0] = sharedframes[0].zaxis; + sharedcorners[1] = sharedframes[1].origin; + sharednormals[1] = sharedframes[1].zaxis; + sharedcorners[2] = sharedframes[2].origin; + sharednormals[2] = sharedframes[2].zaxis; + sharedcorners[3] = sharedframes[3].origin; + sharednormals[3] = sharedframes[3].zaxis; + + ON_Plane *newframes; + ON_3dPoint *newcorners; + ON_3dVector *newnormals; + newframes = (ON_Plane *)bu_malloc(9*sizeof(ON_Plane), "new frames"); + newcorners = (ON_3dPoint *)bu_malloc(9*sizeof(ON_3dPoint), "new corners"); + newnormals = (ON_3dVector *)bu_malloc(9*sizeof(ON_3dVector), "new normals"); + newframes[0] = frames[0]; + newcorners[0] = corners[0]; + newnormals[0] = normals[0]; + newframes[1] = sharedframes[0]; + newcorners[... [truncated message content] |