Revision: 76840
http://sourceforge.net/p/brlcad/code/76840
Author: starseeker
Date: 2020-08-19 01:05:00 +0000 (Wed, 19 Aug 2020)
Log Message:
-----------
Remove a few more old, unused files
Modified Paths:
--------------
brlcad/trunk/doc/legal/embedded/CMakeLists.txt
brlcad/trunk/src/libbrep/CMakeLists.txt
Removed Paths:
-------------
brlcad/trunk/doc/legal/embedded/nurbs_fit.txt
brlcad/trunk/src/libbrep/opennurbs_fit.cpp
brlcad/trunk/src/libbrep/opennurbs_fit.h
brlcad/trunk/src/libbrep/surface_tree_queue_tests.patch
Modified: brlcad/trunk/doc/legal/embedded/CMakeLists.txt
===================================================================
--- brlcad/trunk/doc/legal/embedded/CMakeLists.txt 2020-08-19 00:59:34 UTC
(rev 76839)
+++ brlcad/trunk/doc/legal/embedded/CMakeLists.txt 2020-08-19 01:05:00 UTC
(rev 76840)
@@ -37,7 +37,6 @@
mt19937ar.txt
naca.txt
normalize.txt
- nurbs_fit.txt
obr.txt
OpenNURBS.txt
osg.txt
Deleted: brlcad/trunk/doc/legal/embedded/nurbs_fit.txt
===================================================================
--- brlcad/trunk/doc/legal/embedded/nurbs_fit.txt 2020-08-19 00:59:34 UTC
(rev 76839)
+++ brlcad/trunk/doc/legal/embedded/nurbs_fit.txt 2020-08-19 01:05:00 UTC
(rev 76840)
@@ -1,32 +0,0 @@
-Copyright (c) 2011, Thomas Mörwald, Jonathan Balzer, Inc.
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without
-modification, are permitted provided that the following conditions
-are met:
-
- * Redistributions of source code must retain the above copyright
- notice, this list of conditions and the following disclaimer.
- * Redistributions in binary form must reproduce the above
- copyright notice, this list of conditions and the following
- disclaimer in the documentation and/or other materials provided
- with the distribution.
- * Neither the name of Thomas Mörwald or Jonathan Balzer nor the names of its
- contributors may be used to endorse or promote products derived
- from this software without specific prior written permission.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
-FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
-COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
-INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
-BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
-LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
-ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-POSSIBILITY OF SUCH DAMAGE.
-
-file:/src/libbrep/opennurbs_fit.cpp
-file:/src/libbrep/opennurbs_fit.h
Modified: brlcad/trunk/src/libbrep/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/libbrep/CMakeLists.txt 2020-08-19 00:59:34 UTC (rev
76839)
+++ brlcad/trunk/src/libbrep/CMakeLists.txt 2020-08-19 01:05:00 UTC (rev
76840)
@@ -46,7 +46,6 @@
shape_recognition/util.cpp
ssx_event.cpp
tools/util.cpp
- #opennurbs_fit.cpp
)
set(libbrep_ignored_files
@@ -55,12 +54,9 @@
cdt/cdt.h
cdt/mesh.h
cdt/RTree.h
- opennurbs_fit.h
- opennurbs_fit.cpp
PullbackCurve.h
shape_recognition/shape_recognition.h
shape_recognition/torus.cpp
- surface_tree_queue_tests.patch
)
CMAKEFILES(${libbrep_ignored_files})
Deleted: brlcad/trunk/src/libbrep/opennurbs_fit.cpp
===================================================================
--- brlcad/trunk/src/libbrep/opennurbs_fit.cpp 2020-08-19 00:59:34 UTC (rev
76839)
+++ brlcad/trunk/src/libbrep/opennurbs_fit.cpp 2020-08-19 01:05:00 UTC (rev
76840)
@@ -1,1557 +0,0 @@
-/** @file opennurbs_fit.cpp
- *
- * Extensions to the openNURBS library, based off of Thomas Mörwald's
- * surface fitting code in the Point Cloud Library (pcl). His code is
- # BSD licensed:
- *
- * Copyright (c) 2011, Thomas Mörwald, Jonathan Balzer, Inc.
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- *
- * * Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * * Redistributions in binary form must reproduce the above
- * copyright notice, this list of conditions and the following
- * disclaimer in the documentation and/or other materials provided
- * with the distribution.
- * * Neither the name of Thomas Mörwald or Jonathan Balzer nor the names of
its
- * contributors may be used to endorse or promote products derived
- * from this software without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
- * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
- * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
- * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
- * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
- * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- *
- */
-
-#include "common.h"
-
-#include <iostream>
-#include <stdexcept>
-#include <map>
-#include <algorithm>
-
-#include "opennurbs_fit.h"
-#include "vmath.h"
-
-using namespace on_fit;
-
-// ************************
-// * from nurbs_tools.cpp *
-// ************************
-
-void
-NurbsTools::downsample_random (const vector_vec3d &data1, vector_vec3d &data2,
unsigned size)
-{
- if (data1.size () <= size && size > 0) {
- data2 = data1;
- return;
- }
-
- unsigned s = data1.size ();
- data2.clear ();
-
- for (unsigned i = 0; i < size; i++) {
- unsigned rnd = unsigned (s * (double (rand ()) / RAND_MAX));
- data2.push_back (data1[rnd]);
- }
-}
-
-
-void
-NurbsTools::downsample_random (vector_vec3d &data, unsigned size)
-{
- if (data.size () <= size && size > 0)
- return;
-
- unsigned s = data.size ();
-
- vector_vec3d data_tmp;
-
- for (unsigned i = 0; i < size; i++) {
- unsigned rnd = unsigned ((s - 1) * (double (rand ()) / RAND_MAX));
- data_tmp.push_back (data[rnd]);
- }
-
- data = data_tmp;
-}
-
-
-unsigned
-NurbsTools::getClosestPoint (const ON_2dPoint &p, const vector_vec2d &data)
-{
- if (data.empty ())
- throw std::runtime_error ("[NurbsTools::getClosestPoint(2d)] Data
empty.\n");
-
- unsigned idx (0);
- double dist2 (DBL_MAX);
- for (unsigned i = 0; i < data.size (); i++) {
- ON_2dVector v = (data[i] - p);
- double d2 = ON_DotProduct(v, v); // Was the squaredNorm in Eigen
- if (d2 < dist2) {
- idx = i;
- dist2 = d2;
- }
- }
- return idx;
-}
-
-
-unsigned
-NurbsTools::getClosestPoint (const ON_3dPoint &p, const vector_vec3d &data)
-{
- if (data.empty ())
- throw std::runtime_error ("[NurbsTools::getClosestPoint(2d)] Data
empty.\n");
-
- unsigned idx (0);
- double dist2 (DBL_MAX);
- for (unsigned i = 0; i < data.size (); i++) {
- ON_3dVector v = (data[i] - p);
- double d2 = ON_DotProduct(v, v); // Was the squaredNorm in Eigen
- if (d2 < dist2) {
- idx = i;
- dist2 = d2;
- }
- }
- return idx;
-}
-
-
-unsigned
-NurbsTools::getClosestPoint (const ON_2dPoint &p, const ON_2dVector &dir,
const vector_vec2d &data,
- unsigned &idxcp)
-{
- if (data.empty ())
- throw std::runtime_error ("[NurbsTools::getClosestPoint(2d)] Data
empty.\n");
-
- unsigned idx (0);
- idxcp = 0;
- double dist2 (0.0);
- double dist2cp (DBL_MAX);
- for (unsigned i = 0; i < data.size (); i++) {
- ON_2dVector v = (data[i] - p);
- double d2 = ON_DotProduct(v, v);
-
- if (d2 < dist2cp) {
- idxcp = i;
- dist2cp = d2;
- }
-
- if (NEAR_ZERO(d2, SMALL_FASTF))
- return i;
-
- v.Unitize();
-
- double d1 = ON_DotProduct(dir, v);
- if (d1 / d2 > dist2) {
- idx = i;
- dist2 = d1 / d2;
- }
- }
- return idx;
-}
-
-
-ON_3dVector
-NurbsTools::computeMean (const vector_vec3d &data)
-{
- ON_3dVector u(0.0, 0.0, 0.0);
-
- unsigned s = data.size ();
- double ds = 1.0 / s;
-
- for (unsigned i = 0; i < s; i++)
- u += (data[i] * ds);
-
- return u;
-}
-
-
-ON_2dVector
-NurbsTools::computeMean (const vector_vec2d &data)
-{
- ON_2dVector u (0.0, 0.0);
-
- unsigned s = data.size ();
- double ds = 1.0 / s;
-
- for (unsigned i = 0; i < s; i++)
- u += (data[i] * ds);
-
- return u;
-}
-
-
-void
-NurbsTools::pca (const vector_vec3d &data, ON_3dVector &mean, Eigen::Matrix3d
&eigenvectors,
- Eigen::Vector3d &eigenvalues)
-{
- if (data.empty ()) {
- printf ("[NurbsTools::pca] Error, data is empty\n");
- abort ();
- }
-
- mean = computeMean (data);
-
- unsigned s = data.size ();
-
- ON_Matrix Q(3, s);
-
- for (unsigned i = 0; i < s; i++) {
- Q[0][i] = data[i].x - mean.x;
- Q[1][i] = data[i].y - mean.y;
- Q[2][i] = data[i].z - mean.z;
- }
-
- ON_Matrix Qt = Q;
- Qt.Transpose();
-
- ON_Matrix oC;
- oC.Multiply(Q, Qt);
-
- Eigen::Matrix3d C(3, 3);
- for (unsigned i = 0; i < 3; i++) {
- for (unsigned j = 0; j < 3; j++) {
- C(i, j) = oC[i][j];
- }
- }
-
- Eigen::SelfAdjointEigenSolver < Eigen::Matrix3d > eigensolver (C);
- if (eigensolver.info () != Eigen::Success) {
- printf ("[NurbsTools::pca] Can not find eigenvalues.\n");
- abort ();
- }
-
- for (int i = 0; i < 3; ++i) {
- eigenvalues (i) = eigensolver.eigenvalues () (2 - i);
- if (i == 2)
- eigenvectors.col (2) = eigenvectors.col (0).cross (eigenvectors.col
(1));
- else
- eigenvectors.col (i) = eigensolver.eigenvectors ().col (2 - i);
- }
-}
-
-
-// ******************************
-// * from nurbs_solve_eigen.cpp *
-// ******************************
-
-void
-NurbsSolve::assign (unsigned rows, unsigned cols, unsigned dims)
-{
- m_Ksparse.resize(rows, cols);
- //std::cout << "m_Ksparse(rows=" << rows << ", cols=" << cols << "\n";
- m_Ksparse.setZero();
- m_xeig = new ON_Matrix(cols, dims);
- m_xeig->Zero();
- //std::cout << "m_xeig(rows=" << cols << ", cols=" << dims << "\n";
- m_feig = new ON_Matrix(rows, dims);
- m_feig->Zero();
- //std::cout << "m_feig(rows=" << rows << ", cols=" << dims << "\n";
-}
-
-
-void
-NurbsSolve::K (unsigned i, unsigned j, double v)
-{
- m_Ksparse.insert(i, j) = v;
- //std::cout << "m_Ksparse count: " << m_Ksparse.nonZeros() << "\n";
- //std::cout << "m_Ksparse item(" << i << ", " << j << "):" <<
m_Ksparse.coeffRef(i, j) << "\n";
-}
-void
-NurbsSolve::x (unsigned i, unsigned j, double v)
-{
- (*m_xeig)[i][j] = v;
-}
-void
-NurbsSolve::f (unsigned i, unsigned j, double v)
-{
- (*m_feig)[i][j] = v;
-}
-
-
-double
-NurbsSolve::K (unsigned i, unsigned j)
-{
- return m_Ksparse.coeffRef(i, j);
-}
-double
-NurbsSolve::x (unsigned i, unsigned j)
-{
- return (*m_xeig)[i][j];
-}
-double
-NurbsSolve::f (unsigned i, unsigned j)
-{
- return (*m_feig)[i][j];
-}
-
-
-void
-NurbsSolve::resizeF (unsigned rows)
-{
- ON_Matrix *new_feig = new ON_Matrix(rows, (*m_feig).ColCount());
- for (int i = 0; i < (*m_feig).RowCount(); i++) {
- for (int j = 0; j < (*m_feig).ColCount(); j++) {
- (*new_feig)[i][j] = (*m_feig)[i][j];
- }
- }
- delete m_feig;
- m_feig = new_feig;
-}
-
-
-bool
-solveSparseLinearSystemLQ (Eigen::SparseMatrix<double>* A, Eigen::MatrixXd* b,
Eigen::MatrixXd* x)
-{
- Eigen::SparseMatrix<double> At(A->cols(), A->rows());
- Eigen::MatrixXd Atb(A->cols(), b->cols());
-
- At = A->transpose();
- Eigen::SparseMatrix<double> AtA = At * (*A);
- Atb = At * (*b);
-
- Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
- solver.compute(AtA);
- if (solver.info()!=Eigen::Success) {
- // decomposition failed
- std::cout << "decomposition failed\n";
- }
-
- (*x) = solver.solve(Atb);
-
- if (solver.info()!=Eigen::Success) {
- std::cout << "solver failed: " << solver.info() << "\n";
- return false;
- } else {
- return true;
- }
-}
-
-
-bool
-NurbsSolve::solve ()
-{
-
- clock_t time_start, time_end;
- time_start = clock ();
-
- Eigen::MatrixXd e_m_xeig = Eigen::MatrixXd::Zero ((*m_xeig).RowCount(),
(*m_xeig).ColCount());
- Eigen::MatrixXd e_m_feig = Eigen::MatrixXd::Zero ((*m_feig).RowCount(),
(*m_feig).ColCount());
-
- for (int i = 0; i < (*m_xeig).RowCount(); i++) {
- for (int j = 0; j < (*m_xeig).ColCount(); j++) {
- e_m_xeig (i, j) = (*m_xeig)[i][j];
- }
- }
- for (int i = 0; i < (*m_feig).RowCount(); i++) {
- for (int j = 0; j < (*m_feig).ColCount(); j++) {
- e_m_feig (i, j) = (*m_feig)[i][j];
- }
- }
-
- bool success = solveSparseLinearSystemLQ (&m_Ksparse, &e_m_feig,
&e_m_xeig);
- if (!success) {
- std::cout << "solver failed!\n";
- return false;
- }
-
- for (int i = 0; i < (*m_xeig).RowCount(); i++) {
- for (int j = 0; j < (*m_xeig).ColCount(); j++) {
- (*m_xeig)[i][j] = e_m_xeig (i, j);
- }
- }
-
- time_end = clock ();
-
- if (!m_quiet) {
- double solve_time = (double)(time_end - time_start) /
(double)(CLOCKS_PER_SEC);
- printf ("[NurbsSolve[Eigen::SimplicialLDLT]::solve_()] solution found!
(%f sec)\n", solve_time);
- }
-
- return true;
-}
-
-
-// ********************************
-// * from fitting_surface_pdm.cpp *
-// ********************************
-
-FittingSurface::FittingSurface (int order, NurbsDataSurface *in_m_data,
ON_3dVector z)
-{
- if (order < 2)
- throw std::runtime_error ("[FittingSurface::FittingSurface] Error order
too low (order<2).");
-
- ON::Begin ();
-
- this->m_data = in_m_data;
- m_nurbs = initNurbsPCA (order, m_data, z);
-
- this->init ();
-}
-
-
-FittingSurface::FittingSurface (NurbsDataSurface *in_m_data, const
ON_NurbsSurface &ns)
-{
- ON::Begin ();
-
- this->m_nurbs = ON_NurbsSurface (ns);
- this->m_data = in_m_data;
-
- this->init ();
-}
-
-
-void
-FittingSurface::refine (int dim)
-{
- std::vector<double> xi;
- std::vector<double> elements = getElementVector (m_nurbs, dim);
-
- for (unsigned i = 0; i < elements.size () - 1; i++)
- xi.push_back (elements[i] + 0.5 * (elements[i + 1] - elements[i]));
-
- for (unsigned i = 0; i < xi.size (); i++)
- m_nurbs.InsertKnot (dim, xi[i], 1);
-
- m_elementsU = getElementVector (m_nurbs, 0);
- m_elementsV = getElementVector (m_nurbs, 1);
- m_minU = m_elementsU[0];
- m_minV = m_elementsV[0];
- m_maxU = m_elementsU[m_elementsU.size () - 1];
- m_maxV = m_elementsV[m_elementsV.size () - 1];
-}
-
-
-void
-FittingSurface::assemble (Parameter param)
-{
- clock_t time_start, time_end;
- time_start = clock ();
-
- int nBnd = m_data->boundary.size ();
- int nInt = m_data->interior.size ();
- int nCurInt = param.regularisation_resU * param.regularisation_resV;
- int nCurBnd = 2 * param.regularisation_resU + 2 *
param.regularisation_resV;
- int nCageReg = (m_nurbs.m_cv_count[0] - 2) * (m_nurbs.m_cv_count[1] - 2);
- int nCageRegBnd = 2 * (m_nurbs.m_cv_count[0] - 1) + 2 *
(m_nurbs.m_cv_count[1] - 1);
-
- if (param.boundary_weight <= 0.0)
- nBnd = 0;
- if (param.interior_weight <= 0.0)
- nInt = 0;
- if (param.boundary_regularisation <= 0.0)
- nCurBnd = 0;
- if (param.interior_regularisation <= 0.0)
- nCurInt = 0;
- if (param.interior_smoothness <= 0.0)
- nCageReg = 0;
- if (param.boundary_smoothness <= 0.0)
- nCageRegBnd = 0;
-
- int ncp = m_nurbs.m_cv_count[0] * m_nurbs.m_cv_count[1];
- int nrows = nBnd + nInt + nCurInt + nCurBnd + nCageReg + nCageRegBnd;
-
- m_solver.assign (nrows, ncp, 3);
-
- unsigned row = 0;
-
- // boundary points should lie on edges of surface
- if (nBnd > 0)
- assembleBoundary (param.boundary_weight, row);
-
- // interior points should lie on surface
- if (nInt > 0)
- assembleInterior (param.interior_weight, row);
-
- // minimal curvature on surface
- if (nCurInt > 0) {
- if (m_nurbs.Order (0) < 3 || m_nurbs.Order (1) < 3)
- printf ("[FittingSurface::assemble] Error insufficient NURBS order
to add curvature regularisation.\n");
- else
- addInteriorRegularisation (2, param.regularisation_resU,
param.regularisation_resV,
- param.interior_regularisation /
param.regularisation_resU, row);
- }
-
- // minimal curvature on boundary
- if (nCurBnd > 0) {
- if (m_nurbs.Order (0) < 3 || m_nurbs.Order (1) < 3)
- printf ("[FittingSurface::assemble] Error insufficient NURBS order
to add curvature regularisation.\n");
- else
- addBoundaryRegularisation (2, param.regularisation_resU,
param.regularisation_resV,
- param.boundary_regularisation /
param.regularisation_resU, row);
- }
-
- // cage regularisation
- if (nCageReg > 0)
- addCageInteriorRegularisation (param.interior_smoothness, row);
-
- if (nCageRegBnd > 0) {
- addCageBoundaryRegularisation (param.boundary_smoothness, NORTH, row);
- addCageBoundaryRegularisation (param.boundary_smoothness, SOUTH, row);
- addCageBoundaryRegularisation (param.boundary_smoothness, WEST, row);
- addCageBoundaryRegularisation (param.boundary_smoothness, EAST, row);
- addCageCornerRegularisation (param.boundary_smoothness * 2.0, row);
- }
-
- time_end = clock ();
- if (!m_quiet) {
- double solve_time = (double)(time_end - time_start) /
(double)(CLOCKS_PER_SEC);
- printf ("[FittingSurface::assemble()] (assemble (%d, %d): %f sec)\n",
nrows, ncp, solve_time);
- }
-}
-
-
-void
-FittingSurface::init ()
-{
- m_elementsU = getElementVector (m_nurbs, 0);
- m_elementsV = getElementVector (m_nurbs, 1);
- m_minU = m_elementsU[0];
- m_minV = m_elementsV[0];
- m_maxU = m_elementsU[m_elementsU.size () - 1];
- m_maxV = m_elementsV[m_elementsV.size () - 1];
-
- in_max_steps = 100;
- in_accuracy = 1e-4;
-
- m_quiet = true;
-}
-
-
-void
-FittingSurface::solve (double damp)
-{
- if (m_solver.solve ())
- updateSurf (damp);
-}
-
-
-void
-FittingSurface::updateSurf (double damp)
-{
- int ncp = m_nurbs.m_cv_count[0] * m_nurbs.m_cv_count[1];
-
- for (int A = 0; A < ncp; A++) {
-
- int I = gl2gr (A);
- int J = gl2gc (A);
-
- ON_3dPoint cp_prev;
- m_nurbs.GetCV (I, J, cp_prev);
-
- ON_3dPoint cp;
- cp.x = cp_prev.x + damp * (m_solver.x (A, 0) - cp_prev.x);
- cp.y = cp_prev.y + damp * (m_solver.x (A, 1) - cp_prev.y);
- cp.z = cp_prev.z + damp * (m_solver.x (A, 2) - cp_prev.z);
-
- m_nurbs.SetCV (I, J, cp);
-
- }
-
-}
-
-
-void
-FittingSurface::setInvMapParams (unsigned max_steps, double accuracy)
-{
- this->in_max_steps = max_steps;
- this->in_accuracy = accuracy;
-}
-
-
-std::vector<double>
-FittingSurface::getElementVector (const ON_NurbsSurface &nurbs, int dim) // !
-{
- std::vector<double> result;
-
- int idx_min = 0;
- int idx_max = nurbs.KnotCount (dim) - 1;
- if (nurbs.IsClosed (dim)) {
- idx_min = nurbs.Order (dim) - 2;
- idx_max = nurbs.KnotCount (dim) - nurbs.Order (dim) + 1;
- }
-
- const double* knots = nurbs.Knot (dim);
-
- result.push_back (knots[idx_min]);
-
- //for (int E=(m_nurbs.Order(0)-2);
E<(m_nurbs.KnotCount(0)-m_nurbs.Order(0)+2); E++) {
- for (int E = idx_min + 1; E <= idx_max; E++) {
-
- if (!NEAR_EQUAL(knots[E], knots[E - 1], SMALL_FASTF)) // do not count
double knots
- result.push_back (knots[E]);
-
- }
-
- return result;
-}
-
-
-void
-FittingSurface::assembleInterior (double wInt, unsigned &row)
-{
- m_data->interior_line_start.clear ();
- m_data->interior_line_end.clear ();
- m_data->interior_error.clear ();
- m_data->interior_normals.clear ();
- unsigned nInt = m_data->interior.size ();
- for (unsigned p = 0; p < nInt; p++) {
- ON_3dVector &pcp = m_data->interior[p];
-
- // inverse mapping
- ON_2dVector params;
- ON_3dPoint pt;
- ON_3dVector tu, tv, n;
- double error;
- if (p < m_data->interior_param.size ()) {
- params = inverseMapping (m_nurbs, pcp, m_data->interior_param[p],
error, pt, tu, tv, in_max_steps, in_accuracy);
- m_data->interior_param[p] = params;
- } else {
- params = findClosestElementMidPoint (m_nurbs, pcp);
- params = inverseMapping (m_nurbs, pcp, params, error, pt, tu, tv,
in_max_steps, in_accuracy);
- m_data->interior_param.push_back (params);
- }
- m_data->interior_error.push_back (error);
-
- n = ON_CrossProduct(tu, tv);
- n.Unitize();
-
- m_data->interior_normals.push_back (n);
- m_data->interior_line_start.push_back (pcp);
- m_data->interior_line_end.push_back (pt);
-
- double w (wInt);
- if (p < m_data->interior_weight.size ())
- w = m_data->interior_weight[p];
-
- addPointConstraint (m_data->interior_param[p], m_data->interior[p], w,
row);
- }
-}
-
-
-void
-FittingSurface::assembleBoundary (double wBnd, unsigned &row)
-{
- m_data->boundary_line_start.clear ();
- m_data->boundary_line_end.clear ();
- m_data->boundary_error.clear ();
- m_data->boundary_normals.clear ();
- unsigned nBnd = m_data->boundary.size ();
- for (unsigned p = 0; p < nBnd; p++) {
- ON_3dVector &pcp = m_data->boundary[p];
-
- double error;
- ON_3dPoint pt;
- ON_3dVector tu, tv, n;
- ON_2dVector params = inverseMappingBoundary (m_nurbs, pcp, error, pt,
tu, tv, in_max_steps, in_accuracy);
- m_data->boundary_error.push_back (error);
-
- if (p < m_data->boundary_param.size ()) {
- m_data->boundary_param[p] = params;
- } else {
- m_data->boundary_param.push_back (params);
- }
-
- n = ON_CrossProduct(tu, tv);
- n.Unitize();
-
- m_data->boundary_normals.push_back (n);
- m_data->boundary_line_start.push_back (pcp);
- m_data->boundary_line_end.push_back (pt);
-
- double w (wBnd);
- if (p < m_data->boundary_weight.size ())
- w = m_data->boundary_weight[p];
-
- addPointConstraint (m_data->boundary_param[p], m_data->boundary[p], w,
row);
- }
-}
-
-
-ON_NurbsSurface
-FittingSurface::initNurbs4Corners (int order, ON_3dPoint ll, ON_3dPoint lr,
ON_3dPoint ur, ON_3dPoint ul)
-{
- std::map<int, std::map<int, ON_3dPoint> > cv_map;
-
- double dc = 1.0 / (order - 1);
-
- for (int i = 0; i < order; i++) {
- double di = dc * i;
- cv_map[i][0] = ll + (lr - ll) * di;
- cv_map[0][i] = ll + (ul - ll) * di;
- cv_map[i][order - 1] = ul + (ur - ul) * di;
- cv_map[order - 1][i] = lr + (ur - lr) * di;
- }
-
- for (int i = 1; i < order - 1; i++) {
- for (int j = 1; j < order - 1; j++) {
- ON_3dPoint du = cv_map[0][j] + (cv_map[order - 1][j] -
cv_map[0][j]) * dc * i;
- ON_3dPoint dv = cv_map[i][0] + (cv_map[i][order - 1] -
cv_map[i][0]) * dc * j;
- cv_map[i][j] = du * 0.5 + dv * 0.5;
- }
- }
-
- ON_NurbsSurface nurbs (3, false, order, order, order, order);
- nurbs.MakeClampedUniformKnotVector (0, 1.0);
- nurbs.MakeClampedUniformKnotVector (1, 1.0);
-
- for (int i = 0; i < order; i++) {
- for (int j = 0; j < order; j++) {
- nurbs.SetCV (i, j, cv_map[i][j]);
- }
- }
- return nurbs;
-}
-
-
-ON_NurbsSurface
-FittingSurface::initNurbsPCA (int order, NurbsDataSurface *m_data, ON_3dVector
z)
-{
- ON_3dVector mean;
- Eigen::Matrix3d eigenvectors;
- Eigen::Vector3d eigenvalues;
-
- unsigned s = m_data->interior.size ();
-
- NurbsTools::pca (m_data->interior, mean, eigenvectors, eigenvalues);
-
- m_data->mean = mean;
- //m_data->eigenvectors = (*eigenvectors);
-
- bool flip (false);
- Eigen::Vector3d ez(z[0], z[1], z[2]);
- if (eigenvectors.col (2).dot (ez) < 0.0)
- flip = true;
-
- eigenvalues = eigenvalues / s; // seems that the eigenvalues are dependent
on the number of points (???)
-
- ON_3dVector sigma(sqrt(eigenvalues[0]), sqrt(eigenvalues[1]),
sqrt(eigenvalues[2]));
-
- ON_NurbsSurface nurbs (3, false, order, order, order, order);
- nurbs.MakeClampedUniformKnotVector (0, 1.0);
- nurbs.MakeClampedUniformKnotVector (1, 1.0);
-
- // +- 2 sigma -> 95, 45 % aller Messwerte
- double dcu = (4.0 * sigma[0]) / (nurbs.Order (0) - 1);
- double dcv = (4.0 * sigma[1]) / (nurbs.Order (1) - 1);
-
- ON_3dVector cv_t, cv;
- Eigen::Vector3d ecv_t, ecv;
- Eigen::Vector3d emean(mean[0], mean[1], mean[2]);
- for (int i = 0; i < nurbs.Order (0); i++) {
- for (int j = 0; j < nurbs.Order (1); j++) {
- cv[0] = -2.0 * sigma[0] + dcu * i;
- cv[1] = -2.0 * sigma[1] + dcv * j;
- cv[2] = 0.0;
- ecv (0) = -2.0 * sigma[0] + dcu * i;
- ecv (1) = -2.0 * sigma[1] + dcv * j;
- ecv (2) = 0.0;
- ecv_t = eigenvectors * ecv + emean;
- cv_t[0] = ecv_t (0);
- cv_t[1] = ecv_t (1);
- cv_t[2] = ecv_t (2);
- if (flip)
- nurbs.SetCV (nurbs.Order (0) - 1 - i, j, ON_3dPoint (cv_t[0],
cv_t[1], cv_t[2]));
- else
- nurbs.SetCV (i, j, ON_3dPoint (cv_t[0], cv_t[1], cv_t[2]));
- }
- }
- return nurbs;
-}
-
-
-ON_NurbsSurface
-FittingSurface::initNurbsPCABoundingBox (int order, NurbsDataSurface *m_data,
ON_3dVector z)
-{
- ON_3dVector mean;
- Eigen::Matrix3d eigenvectors;
- Eigen::Vector3d eigenvalues;
-
- unsigned s = m_data->interior.size ();
- m_data->interior_param.clear ();
-
- NurbsTools::pca (m_data->interior, mean, eigenvectors, eigenvalues);
-
- m_data->mean = mean;
- //m_data->eigenvectors = (*eigenvectors);
-
- bool flip (false);
- Eigen::Vector3d ez(z[0], z[1], z[2]);
- if (eigenvectors.col (2).dot (ez) < 0.0)
- flip = true;
-
- eigenvalues = eigenvalues / s; // seems that the eigenvalues are dependent
on the number of points (???)
- Eigen::Matrix3d eigenvectors_inv = eigenvectors.inverse ();
-
- ON_3dVector v_max(0.0, 0.0, 0.0);
- ON_3dVector v_min(DBL_MAX, DBL_MAX, DBL_MAX);
- Eigen::Vector3d emean(mean[0], mean[1], mean[2]);
- for (unsigned i = 0; i < s; i++) {
- Eigen::Vector3d eint(m_data->interior[i][0], m_data->interior[i][1],
m_data->interior[i][2]);
- Eigen::Vector3d ep = eigenvectors_inv * (eint - emean);
- ON_3dPoint p(ep (0), ep (1), ep(2));
- m_data->interior_param.push_back (ON_2dPoint(p[0], p[1]));
-
- V_MAX(v_max[0], p[0]);
- V_MAX(v_max[1], p[1]);
- V_MAX(v_max[2], p[2]);
-
- V_MIN(v_min[0], p[0]);
- V_MIN(v_min[1], p[1]);
- V_MIN(v_min[2], p[2]);
- }
-
- for (unsigned i = 0; i < s; i++) {
- ON_2dVector &p = m_data->interior_param[i];
- if (v_max[0] > v_min[0] && v_max[0] > v_min[0]) {
- p[0] = (p[0] - v_min[0]) / (v_max[0] - v_min[0]);
- p[1] = (p[1] - v_min[1]) / (v_max[1] - v_min[1]);
- } else {
- throw std::runtime_error ("[NurbsTools::initNurbsPCABoundingBox]
Error: v_max <= v_min");
- }
- }
-
- ON_NurbsSurface nurbs (3, false, order, order, order, order);
-
- nurbs.MakeClampedUniformKnotVector (0, 1.0);
- nurbs.MakeClampedUniformKnotVector (1, 1.0);
-
- double dcu = (v_max[0] - v_min[0]) / (nurbs.Order (0) - 1);
- double dcv = (v_max[1] - v_min[1]) / (nurbs.Order (1) - 1);
-
- ON_3dPoint cv_t, cv;
- Eigen::Vector3d ecv_t2, ecv2;
- Eigen::Vector3d emean2(mean[0], mean[1], mean[2]);
- for (int i = 0; i < nurbs.Order (0); i++) {
- for (int j = 0; j < nurbs.Order (1); j++) {
- cv[0] = v_min[0] + dcu * i;
- cv[1] = v_min[1] + dcv * j;
- cv[2] = 0.0;
- ecv2 (0) = cv[0];
- ecv2 (1) = cv[1];
- ecv2 (2) = cv[2];
- ecv_t2 = eigenvectors * ecv2 + emean2;
- if (flip)
- nurbs.SetCV (nurbs.Order (0) - 1 - i, j, cv_t);
- else
- nurbs.SetCV (i, j, cv_t);
- }
- }
- return nurbs;
-}
-
-
-void
-FittingSurface::addPointConstraint (const ON_2dVector ¶ms, const
ON_3dPoint &point, double weight,
- unsigned &row)
-{
- double *N0 = new double[m_nurbs.Order (0) * m_nurbs.Order (0)];
- double *N1 = new double[m_nurbs.Order (1) * m_nurbs.Order (1)];
-
- int E = ON_NurbsSpanIndex (m_nurbs.m_order[0], m_nurbs.m_cv_count[0],
m_nurbs.m_knot[0], params[0], 0, 0);
- int F = ON_NurbsSpanIndex (m_nurbs.m_order[1], m_nurbs.m_cv_count[1],
m_nurbs.m_knot[1], params[1], 0, 0);
-
- ON_EvaluateNurbsBasis (m_nurbs.Order (0), m_nurbs.m_knot[0] + E,
params[0], N0);
- ON_EvaluateNurbsBasis (m_nurbs.Order (1), m_nurbs.m_knot[1] + F,
params[1], N1);
-
- m_solver.f (row, 0, point[0] * weight);
- m_solver.f (row, 1, point[1] * weight);
- m_solver.f (row, 2, point[2] * weight);
-
- for (int i = 0; i < m_nurbs.Order (0); i++) {
-
- for (int j = 0; j < m_nurbs.Order (1); j++) {
-
- m_solver.K (row, lrc2gl (E, F, i, j), weight * N0[i] * N1[j]);
-
- } // j
-
- } // i
-
- row++;
-
- delete [] N0;
- delete [] N1;
-
-}
-
-
-void
-FittingSurface::addCageInteriorRegularisation (double weight, unsigned &row)
-{
- for (int i = 1; i < (m_nurbs.m_cv_count[0] - 1); i++) {
- for (int j = 1; j < (m_nurbs.m_cv_count[1] - 1); j++) {
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- m_solver.K (row, grc2gl (i + 0, j + 0), -4.0 * weight);
- m_solver.K (row, grc2gl (i + 0, j - 1), 1.0 * weight);
- m_solver.K (row, grc2gl (i + 0, j + 1), 1.0 * weight);
- m_solver.K (row, grc2gl (i - 1, j + 0), 1.0 * weight);
- m_solver.K (row, grc2gl (i + 1, j + 0), 1.0 * weight);
-
- row++;
- }
- }
-}
-
-
-void
-FittingSurface::addCageBoundaryRegularisation (double weight, int side,
unsigned &row)
-{
- int i = 0;
- int j = 0;
-
- switch (side) {
- case SOUTH:
- j = m_nurbs.m_cv_count[1] - 1;
- case NORTH:
- for (i = 1; i < (m_nurbs.m_cv_count[0] - 1); i++) {
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- m_solver.K (row, grc2gl (i + 0, j), -2.0 * weight);
- m_solver.K (row, grc2gl (i - 1, j), 1.0 * weight);
- m_solver.K (row, grc2gl (i + 1, j), 1.0 * weight);
-
- row++;
- }
- break;
-
- case EAST:
- i = m_nurbs.m_cv_count[0] - 1;
- case WEST:
- for (j = 1; j < (m_nurbs.m_cv_count[1] - 1); j++) {
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- m_solver.K (row, grc2gl (i, j + 0), -2.0 * weight);
- m_solver.K (row, grc2gl (i, j - 1), 1.0 * weight);
- m_solver.K (row, grc2gl (i, j + 1), 1.0 * weight);
-
- row++;
- }
- break;
- }
-}
-
-
-void
-FittingSurface::addCageCornerRegularisation (double weight, unsigned &row)
-{
- { // NORTH-WEST
- int i = 0;
- int j = 0;
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- m_solver.K (row, grc2gl (i + 0, j + 0), -2.0 * weight);
- m_solver.K (row, grc2gl (i + 1, j + 0), 1.0 * weight);
- m_solver.K (row, grc2gl (i + 0, j + 1), 1.0 * weight);
-
- row++;
- }
-
- { // NORTH-EAST
- int i = m_nurbs.m_cv_count[0] - 1;
- int j = 0;
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- m_solver.K (row, grc2gl (i + 0, j + 0), -2.0 * weight);
- m_solver.K (row, grc2gl (i - 1, j + 0), 1.0 * weight);
- m_solver.K (row, grc2gl (i + 0, j + 1), 1.0 * weight);
-
- row++;
- }
-
- { // SOUTH-EAST
- int i = m_nurbs.m_cv_count[0] - 1;
- int j = m_nurbs.m_cv_count[1] - 1;
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- m_solver.K (row, grc2gl (i + 0, j + 0), -2.0 * weight);
- m_solver.K (row, grc2gl (i - 1, j + 0), 1.0 * weight);
- m_solver.K (row, grc2gl (i + 0, j - 1), 1.0 * weight);
-
- row++;
- }
-
- { // SOUTH-WEST
- int i = 0;
- int j = m_nurbs.m_cv_count[1] - 1;
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- m_solver.K (row, grc2gl (i + 0, j + 0), -2.0 * weight);
- m_solver.K (row, grc2gl (i + 1, j + 0), 1.0 * weight);
- m_solver.K (row, grc2gl (i + 0, j - 1), 1.0 * weight);
-
- row++;
- }
-
-}
-
-
-void
-FittingSurface::addInteriorRegularisation (int order, int resU, int resV,
double UNUSED(weight), unsigned &row)
-{
- double *N0 = new double[m_nurbs.Order (0) * m_nurbs.Order (0)];
- double *N1 = new double[m_nurbs.Order (1) * m_nurbs.Order (1)];
-
- double dU = (m_maxU - m_minU) / resU;
- double dV = (m_maxV - m_minV) / resV;
-
- for (int u = 0; u < resU; u++) {
- for (int v = 0; v < resV; v++) {
-
- ON_2dPoint params;
- params[0] = m_minU + u * dU + 0.5 * dU;
- params[1] = m_minV + v * dV + 0.5 * dV;
-
- // printf("%f %f, %f %f\n", m_minU, dU, params(0),
params(1));
-
- int E = ON_NurbsSpanIndex (m_nurbs.m_order[0],
m_nurbs.m_cv_count[0], m_nurbs.m_knot[0], params[0], 0, 0);
- int F = ON_NurbsSpanIndex (m_nurbs.m_order[1],
m_nurbs.m_cv_count[1], m_nurbs.m_knot[1], params[1], 0, 0);
-
- ON_EvaluateNurbsBasis (m_nurbs.Order (0), m_nurbs.m_knot[0] + E,
params[0], N0);
- ON_EvaluateNurbsBasis (m_nurbs.Order (1), m_nurbs.m_knot[1] + F,
params[1], N1);
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (0),
m_nurbs.m_knot[0] + E, order, N0); // derivative order?
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (1),
m_nurbs.m_knot[1] + F, order, N1);
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- for (int i = 0; i < m_nurbs.Order (0); i++) {
-
- for (int j = 0; j < m_nurbs.Order (1); j++) {
-
- //m_solver.K (row, lrc2gl (E, F, i, j),
- // weight * (N0[order * m_nurbs.Order (0) + i] *
N1[j] + N0[i] * N1[order * m_nurbs.Order (1) + j]));
-
- } // i
-
- } // j
-
- row++;
-
- }
- }
- delete [] N0;
- delete [] N1;
-}
-
-
-void
-FittingSurface::addBoundaryRegularisation (int order, int resU, int resV,
double weight, unsigned &row)
-{
- double *N0 = new double[m_nurbs.Order (0) * m_nurbs.Order (0)];
- double *N1 = new double[m_nurbs.Order (1) * m_nurbs.Order (1)];
-
- double dU = (m_maxU - m_minU) / resU;
- double dV = (m_maxV - m_minV) / resV;
-
- for (int u = 0; u < resU; u++) {
-
- ON_2dPoint params;
- params[0] = m_minU + u * dU + 0.5 * dU;
- params[1] = m_minV;
-
- int E = ON_NurbsSpanIndex (m_nurbs.m_order[0], m_nurbs.m_cv_count[0],
m_nurbs.m_knot[0], params[0], 0, 0);
- int F = ON_NurbsSpanIndex (m_nurbs.m_order[1], m_nurbs.m_cv_count[1],
m_nurbs.m_knot[1], params[1], 0, 0);
-
- ON_EvaluateNurbsBasis (m_nurbs.Order (0), m_nurbs.m_knot[0] + E,
params[0], N0);
- ON_EvaluateNurbsBasis (m_nurbs.Order (1), m_nurbs.m_knot[1] + F,
params[1], N1);
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (0), m_nurbs.m_knot[0]
+ E, order, N0); // derivative order?
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (1), m_nurbs.m_knot[1]
+ F, order, N1);
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- for (int i = 0; i < m_nurbs.Order (0); i++) {
-
- for (int j = 0; j < m_nurbs.Order (1); j++) {
-
- m_solver.K (row, lrc2gl (E, F, i, j),
- weight * (N0[order * m_nurbs.Order (0) + i] * N1[j]
+ N0[i] * N1[order * m_nurbs.Order (1) + j]));
-
- } // i
-
- } // j
-
- row++;
-
- }
-
- for (int u = 0; u < resU; u++) {
-
- ON_2dPoint params;
- params[0] = m_minU + u * dU + 0.5 * dU;
- params[1] = m_maxV;
-
- int E = ON_NurbsSpanIndex (m_nurbs.m_order[0], m_nurbs.m_cv_count[0],
m_nurbs.m_knot[0], params[0], 0, 0);
- int F = ON_NurbsSpanIndex (m_nurbs.m_order[1], m_nurbs.m_cv_count[1],
m_nurbs.m_knot[1], params[1], 0, 0);
-
- ON_EvaluateNurbsBasis (m_nurbs.Order (0), m_nurbs.m_knot[0] + E,
params[0], N0);
- ON_EvaluateNurbsBasis (m_nurbs.Order (1), m_nurbs.m_knot[1] + F,
params[1], N1);
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (0), m_nurbs.m_knot[0]
+ E, order, N0); // derivative order?
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (1), m_nurbs.m_knot[1]
+ F, order, N1);
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- for (int i = 0; i < m_nurbs.Order (0); i++) {
-
- for (int j = 0; j < m_nurbs.Order (1); j++) {
-
- m_solver.K (row, lrc2gl (E, F, i, j),
- weight * (N0[order * m_nurbs.Order (0) + i] * N1[j]
+ N0[i] * N1[order * m_nurbs.Order (1) + j]));
-
- } // i
-
- } // j
-
- row++;
-
- }
-
- for (int v = 0; v < resV; v++) {
-
- ON_2dPoint params;
- params[0] = m_minU;
- params[1] = m_minV + v * dV + 0.5 * dV;
-
- int E = ON_NurbsSpanIndex (m_nurbs.m_order[0], m_nurbs.m_cv_count[0],
m_nurbs.m_knot[0], params[0], 0, 0);
- int F = ON_NurbsSpanIndex (m_nurbs.m_order[1], m_nurbs.m_cv_count[1],
m_nurbs.m_knot[1], params[1], 0, 0);
-
- ON_EvaluateNurbsBasis (m_nurbs.Order (0), m_nurbs.m_knot[0] + E,
params[0], N0);
- ON_EvaluateNurbsBasis (m_nurbs.Order (1), m_nurbs.m_knot[1] + F,
params[1], N1);
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (0), m_nurbs.m_knot[0]
+ E, order, N0); // derivative order?
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (1), m_nurbs.m_knot[1]
+ F, order, N1);
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- for (int i = 0; i < m_nurbs.Order (0); i++) {
-
- for (int j = 0; j < m_nurbs.Order (1); j++) {
-
- m_solver.K (row, lrc2gl (E, F, i, j),
- weight * (N0[order * m_nurbs.Order (0) + i] * N1[j]
+ N0[i] * N1[order * m_nurbs.Order (1) + j]));
-
- } // i
-
- } // j
-
- row++;
-
- }
-
- for (int v = 0; v < resV; v++) {
-
- ON_2dPoint params;
- params[0] = m_maxU;
- params[1] = m_minV + v * dV + 0.5 * dV;
-
- int E = ON_NurbsSpanIndex (m_nurbs.m_order[0], m_nurbs.m_cv_count[0],
m_nurbs.m_knot[0], params[0], 0, 0);
- int F = ON_NurbsSpanIndex (m_nurbs.m_order[1], m_nurbs.m_cv_count[1],
m_nurbs.m_knot[1], params[1], 0, 0);
-
- ON_EvaluateNurbsBasis (m_nurbs.Order (0), m_nurbs.m_knot[0] + E,
params[0], N0);
- ON_EvaluateNurbsBasis (m_nurbs.Order (1), m_nurbs.m_knot[1] + F,
params[1], N1);
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (0), m_nurbs.m_knot[0]
+ E, order, N0); // derivative order?
- ON_EvaluateNurbsBasisDerivatives (m_nurbs.Order (1), m_nurbs.m_knot[1]
+ F, order, N1);
-
- m_solver.f (row, 0, 0.0);
- m_solver.f (row, 1, 0.0);
- m_solver.f (row, 2, 0.0);
-
- for (int i = 0; i < m_nurbs.Order (0); i++) {
-
- for (int j = 0; j < m_nurbs.Order (1); j++) {
-
- m_solver.K (row, lrc2gl (E, F, i, j),
- weight * (N0[order * m_nurbs.Order (0) + i] * N1[j]
+ N0[i] * N1[order * m_nurbs.Order (1) + j]));
-
- } // i
-
- } // j
-
- row++;
-
- }
- delete [] N0;
- delete [] N1;
-}
-
-
-ON_2dPoint
-FittingSurface::inverseMapping (const ON_NurbsSurface &nurbs, const ON_3dPoint
&pt, const ON_2dPoint &hint, double &error,
- ON_3dPoint &p, ON_3dVector &tu, ON_3dVector
&tv, int maxSteps, double accuracy, bool quiet)
-{
-
- double pointAndTangents[9];
-
- ON_2dVector current, delta;
- ON_Matrix A(2, 2);
- ON_2dVector b;
- ON_3dVector r;
- std::vector<double> elementsU = getElementVector (nurbs, 0);
- std::vector<double> elementsV = getElementVector (nurbs, 1);
- double minU = elementsU[0];
- double minV = elementsV[0];
- double maxU = elementsU[elementsU.size () - 1];
- double maxV = elementsV[elementsV.size () - 1];
-
- current = hint;
-
- for (int k = 0; k < maxSteps; k++) {
-
- nurbs.Evaluate (current[0], current[1], 1, 3, pointAndTangents);
- p[0] = pointAndTangents[0];
- p[1] = pointAndTangents[1];
- p[2] = pointAndTangents[2];
- tu[0] = pointAndTangents[3];
- tu[1] = pointAndTangents[4];
- tu[2] = pointAndTangents[5];
- tv[0] = pointAndTangents[6];
- tv[1] = pointAndTangents[7];
- tv[2] = pointAndTangents[8];
-
- r = p - pt;
-
- b[0] = -ON_DotProduct(r, tu);
- b[1] = -ON_DotProduct(r, tv);
-
- A[0][0] = ON_DotProduct(tu, tu);
- A[0][1] = ON_DotProduct(tu, tv);
- A[1][0] = A[0][1];
- A[1][1] = ON_DotProduct(tv, tv);
-
- Eigen::Vector2d eb(b[0], b[1]);
- Eigen::Vector2d edelta;
- Eigen::Matrix2d eA;
-
- eA (0, 0) = A[0][0];
- eA (0, 1) = A[0][1];
- eA (1, 0) = A[1][0];
- eA (1, 1) = A[1][1];
-
- edelta = eA.ldlt ().solve (eb);
-
- delta[0] = edelta (0);
- delta[1] = edelta (1);
-
- if (sqrt(ON_DotProduct(delta, delta)) < accuracy) {
-
- error = sqrt(ON_DotProduct(r, r));
- return current;
-
- } else {
- current = current + delta;
-
- CLAMP(current[0], minU, maxU);
- CLAMP(current[1], minV, maxV);
- }
-
- }
-
- error = sqrt(ON_DotProduct(r, r));
- if (!quiet) {
- printf ("[FittingSurface::inverseMapping] Warning: Method did not
converge (%e %d)\n", accuracy, maxSteps);
- printf (" %f %f ... %f %f\n", hint[0], hint[1], current[0],
current[1]);
- }
-
- return current;
-
-}
-
-
-ON_2dPoint
-FittingSurface::findClosestElementMidPoint (const ON_NurbsSurface &nurbs,
const ON_3dPoint &pt)
-{
- ON_2dPoint hint;
- ON_3dVector r;
- std::vector<double> elementsU = getElementVector (nurbs, 0);
- std::vector<double> elementsV = getElementVector (nurbs, 1);
-
- double d_shortest (DBL_MAX);
- for (unsigned i = 0; i < elementsU.size () - 1; i++) {
- for (unsigned j = 0; j < elementsV.size () - 1; j++) {
- double points[3];
- double d;
-
- double xi = elementsU[i] + 0.5 * (elementsU[i + 1] - elementsU[i]);
- double eta = elementsV[j] + 0.5 * (elementsV[j + 1] - elementsV[j]);
-
- nurbs.Evaluate (xi, eta, 0, 3, points);
- r[0] = points[0] - pt[0];
- r[1] = points[1] - pt[1];
- r[2] = points[2] - pt[2];
-
- d = ON_DotProduct(r, r);
-
- if ((i == 0 && j == 0) || d < d_shortest) {
- d_shortest = d;
- hint[0] = xi;
- hint[1] = eta;
- }
- }
- }
-
- return hint;
-}
-
-
-ON_2dPoint
-FittingSurface::inverseMappingBoundary (const ON_NurbsSurface &nurbs, const
ON_3dPoint &pt, double &error, ON_3dPoint &p,
- ON_3dVector &tu, ON_3dVector &tv, int
maxSteps, double accuracy, bool quiet)
-{
-
- ON_2dPoint result;
- double min_err = 100.0;
- std::vector<myvec> ini_points;
- double err_tmp;
- ON_3dPoint p_tmp;
- ON_3dVector tu_tmp, tv_tmp;
-
- std::vector<double> elementsU = getElementVector (nurbs, 0);
- std::vector<double> elementsV = getElementVector (nurbs, 1);
-
- // NORTH - SOUTH
- for (unsigned i = 0; i < (elementsV.size () - 1); i++) {
- ini_points.push_back (myvec (WEST, elementsV[i] + 0.5 * (elementsV[i +
1] - elementsV[i])));
- ini_points.push_back (myvec (EAST, elementsV[i] + 0.5 * (elementsV[i +
1] - elementsV[i])));
- }
-
- // WEST - EAST
- for (unsigned i = 0; i < (elementsU.size () - 1); i++) {
- ini_points.push_back (myvec (NORTH, elementsU[i] + 0.5 * (elementsU[i +
1] - elementsU[i])));
- ini_points.push_back (myvec (SOUTH, elementsU[i] + 0.5 * (elementsU[i +
1] - elementsU[i])));
- }
-
- for (unsigned i = 0; i < ini_points.size (); i++) {
-
- ON_2dPoint params = inverseMappingBoundary (nurbs, pt,
ini_points[i].side, ini_points[i].hint, err_tmp, p_tmp,
- tu_tmp, tv_tmp, maxSteps,
accuracy, quiet);
-
- if (i == 0 || err_tmp < min_err) {
- min_err = err_tmp;
- result = params;
- p = p_tmp;
- tu = tu_tmp;
- tv = tv_tmp;
- }
- }
-
- error = min_err;
- return result;
-
-}
-
-
-ON_2dPoint
-FittingSurface::inverseMappingBoundary (const ON_NurbsSurface &nurbs, const
ON_3dPoint &pt, int side, double hint,
- double &error, ON_3dPoint &p,
ON_3dVector &tu, ON_3dVector &tv, int maxSteps,
- double accuracy, bool quiet)
-{
-
- double pointAndTangents[9];
- double current, delta;
- ON_3dVector r, t;
- ON_2dPoint params;
-
- current = hint;
-
- std::vector<double> elementsU = getElementVector (nurbs, 0);
- std::vector<double> elementsV = getElementVector (nurbs, 1);
- double minU = elementsU[0];
- double minV = elementsV[0];
- double maxU = elementsU[elementsU.size () - 1];
- double maxV = elementsV[elementsV.size () - 1];
-
- for (int k = 0; k < maxSteps; k++) {
-
- switch (side) {
-
- case WEST:
-
- params[0] = minU;
- params[1] = current;
- nurbs.Evaluate (minU, current, 1, 3, pointAndTangents);
- p[0] = pointAndTangents[0];
- p[1] = pointAndTangents[1];
- p[2] = pointAndTangents[2];
- tu[0] = pointAndTangents[3];
- tu[1] = pointAndTangents[4];
- tu[2] = pointAndTangents[5];
- tv[0] = pointAndTangents[6];
- tv[1] = pointAndTangents[7];
- tv[2] = pointAndTangents[8];
-
- t = tv; // use tv
-
- break;
- case SOUTH:
-
- params[0] = current;
- params[1] = maxV;
- nurbs.Evaluate (current, maxV, 1, 3, pointAndTangents);
- p[0] = pointAndTangents[0];
- p[1] = pointAndTangents[1];
- p[2] = pointAndTangents[2];
- tu[0] = pointAndTangents[3];
- tu[1] = pointAndTangents[4];
- tu[2] = pointAndTangents[5];
- tv[0] = pointAndTangents[6];
- tv[1] = pointAndTangents[7];
- tv[2] = pointAndTangents[8];
-
- t = tu; // use tu
-
- break;
- case EAST:
-
- params[0] = maxU;
- params[1] = current;
- nurbs.Evaluate (maxU, current, 1, 3, pointAndTangents);
- p[0] = pointAndTangents[0];
- p[1] = pointAndTangents[1];
- p[2] = pointAndTangents[2];
- tu[0] = pointAndTangents[3];
- tu[1] = pointAndTangents[4];
- tu[2] = pointAndTangents[5];
- tv[0] = pointAndTangents[6];
- tv[1] = pointAndTangents[7];
- tv[2] = pointAndTangents[8];
-
- t = tv; // use tv
-
- break;
- case NORTH:
-
- params[0] = current;
- params[1] = minV;
- nurbs.Evaluate (current, minV, 1, 3, pointAndTangents);
- p[0] = pointAndTangents[0];
- p[1] = pointAndTangents[1];
- p[2] = pointAndTangents[2];
- tu[0] = pointAndTangents[3];
- tu[1] = pointAndTangents[4];
- tu[2] = pointAndTangents[5];
- tv[0] = pointAndTangents[6];
- tv[1] = pointAndTangents[7];
- tv[2] = pointAndTangents[8];
-
- t = tu; // use tu
-
- break;
- default:
- throw std::runtime_error
("[FittingSurface::inverseMappingBoundary] ERROR: Specify a boundary!");
-
- }
-
- r[0] = pointAndTangents[0] - pt[0];
- r[1] = pointAndTangents[1] - pt[1];
- r[2] = pointAndTangents[2] - pt[2];
-
- delta = -0.5 * ON_DotProduct(r, t) / ON_DotProduct(t, t);
-
- if (fabs (delta) < accuracy) {
-
- error = sqrt(ON_DotProduct(r, r));
- return params;
-
- } else {
-
- current = current + delta;
-
- bool stop = false;
-
- switch (side) {
-
- case WEST:
- case EAST:
- if (current < minV) {
- params[1] = minV;
- stop = true;
- } else if (current > maxV) {
- params[1] = maxV;
- stop = true;
- }
-
- break;
-
- case NORTH:
- case SOUTH:
- if (current < minU) {
- params[0] = minU;
- stop = true;
- } else if (current > maxU) {
- params[0] = maxU;
- stop = true;
- }
-
- break;
- }
-
- if (stop) {
- error = sqrt(ON_DotProduct(r, r));
- return params;
- }
-
- }
-
- }
-
- error = sqrt(ON_DotProduct(r, r));
- if (!quiet)
- printf (
- "[FittingSurface::inverseMappingBoundary] Warning: Method did not
converge! (residual: %f, delta: %f, params: %f %f)\n",
- error, delta, params[0], params[1]);
-
- return params;
-}
-
-
-/*
- * Local Variables:
- * tab-width: 8
- * mode: C++
- * indent-tabs-mode: t
- * c-file-style: "stroustrup"
- * coding: utf-8
- * End:
- * ex: shiftwidth=4 tabstop=8
- */
Deleted: brlcad/trunk/src/libbrep/opennurbs_fit.h
===================================================================
--- brlcad/trunk/src/libbrep/opennurbs_fit.h 2020-08-19 00:59:34 UTC (rev
76839)
+++ brlcad/trunk/src/libbrep/opennurbs_fit.h 2020-08-19 01:05:00 UTC (rev
76840)
@@ -1,575 +0,0 @@
-/** @File opennurbs_fit.h
- *
- * Extensions to the openNURBS library, based off of Thomas Mörwald's
- * surface fitting code in the Point Cloud Library (pcl). His code is
- * BSD licensed:
- *
- * Copyright (c) 2011, Thomas Mörwald, Jonathan Balzer, Inc.
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- *
- * * Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * * Redistributions in binary form must reproduce the above
- * copyright notice, this list of conditions and the following
- * disclaimer in the documentation and/or other materials provided
- * with the distribution.
- * * Neither the name of Thomas Mörwald or Jonathan Balzer nor the names of
its
- * contributors may be used to endorse or promote products derived
- * from this software without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
- * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
- * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
- * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
- * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
- * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
- * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
- * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- *
- */
-
-#ifndef LIBBREP_OPENNURBS_FIT_H
-#define LIBBREP_OPENNURBS_FIT_H
-
-#include "common.h"
-
-/* system headers */
-#include <vector>
-#include <list>
-#include <map>
-
-/* library headers */
-#include "bio.h" /* needed to include windows.h with protections */
-#define ON_NO_WINDOWS 1 /* don't let opennurbs include windows.h */
-#include "opennurbs.h"
-
-
-#if defined(__GNUC__) && (__GNUC__ == 4 && __GNUC_MINOR__ < 6) &&
!defined(__clang__)
-# pragma message "Disabling GCC shadow warnings via pragma due to Eigen
headers..."
-# pragma message "Disabling GCC float equality comparison warnings via pragma
due to Eigen headers..."
-# pragma message "Disabling GCC inline failure warnings via pragma due to
Eigen headers..."
-#endif
-#if defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >=
6)) && !defined(__clang__)
-# pragma GCC diagnostic push
-#endif
-#if defined(__clang__)
-# pragma clang diagnostic push
-#endif
-#if defined(__GNUC__) && (__GNUC__ == 4 && __GNUC_MINOR__ >= 3) &&
!defined(__clang__)
-# pragma GCC diagnostic ignored "-Wshadow"
-# pragma GCC diagnostic ignored "-Wfloat-equal"
-# pragma GCC diagnostic ignored "-Winline"
-#endif
-#if defined(__clang__)
-# pragma clang diagnostic ignored "-Wshadow"
-# pragma clang diagnostic ignored "-Wfloat-equal"
-# pragma clang diagnostic ignored "-Winline"
-#endif
-#undef Success
-#include <Eigen/StdVector>
-#undef Success
-#include <Eigen/Dense>
-#undef Success
-#include <Eigen/Sparse>
-#undef Success
-#include <Eigen/SparseCholesky>
-#if defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >=
6)) && !defined(__clang__)
-# pragma GCC diagnostic pop
-#endif
-#if defined(__clang__)
-# pragma clang diagnostic pop
-#endif
-
-namespace on_fit
-{
-
-// **********************
-// * from nurbs_data.h *
-// **********************
-
-typedef std::vector<ON_2dVector > vector_vec2d;
-typedef std::vector<ON_3dVector > vector_vec3d;
-
-/** \brief Data structure for NURBS surface fitting
- * (FittingSurface, FittingSurfaceTDM, FittingCylinder, GlobalOptimization,
GlobalOptimizationTDM) */
-struct NurbsDataSurface
-{
- ON_Matrix *eigenvectors;
- ON_3dVector mean;
-
- vector_vec3d interior; ///<< input
- std::vector<double> interior_weight; ///<< input
- std::vector<double> interior_error; ///>> output
- vector_vec2d interior_param; ///>> output
- vector_vec3d interior_line_start; ///>> output
- vector_vec3d interior_line_end; ///>> output
- vector_vec3d interior_normals; ///>> output
-
- vector_vec3d boundary; ///<< input
- std::vector<double> boundary_weight; ///<< input
- std::vector<double> boundary_error; ///>> output
- vector_vec2d boundary_param; ///>> output
- vector_vec3d boundary_line_start; ///>> output
- vector_vec3d boundary_line_end; ///>> output
- vector_vec3d boundary_normals; ///>> output
-
- vector_vec3d common_boundary_point;
- std::vector<unsigned> common_boundary_idx;
- vector_vec2d common_boundary_param;
-
- /** \brief Clear all interior data */
- inline void
- clear_interior ()
- {
- interior.clear ();
- interior_weight.clear ();
- interior_error.clear ();
- interior_param.clear ();
- interior_line_start.clear ();
- interior_line_end.clear ();
- interior_normals.clear ();
- }
-
- /** \brief Clear all boundary data */
- inline void
- clear_boundary ()
- {
- boundary.clear ();
- boundary_weight.clear ();
- boundary_error.clear ();
- boundary_param.clear ();
- boundary_line_start.clear ();
- boundary_line_end.clear ();
- boundary_normals.clear ();
- }
-
- /** \brief Clear all common data */
- inline void
- clear_common_boundary ()
- {
- common_boundary_point.clear ();
- common_boundary_idx.clear ();
- common_boundary_param.clear ();
- }
-
-};
-
-
-// **********************
-// * from nurbs_tools.h *
-// **********************
-enum
-{
- NORTH = 1, NORTHEAST = 2, EAST = 3, SOUTHEAST = 4, SOUTH = 5, SOUTHWEST =
6, WEST = 7, NORTHWEST = 8
-};
-
-
-/** \brief Some useful tools for initialization, point search, ... */
-class NurbsTools
-{
-public:
-
- /** \brief Get the closest point with respect to 'point'
- * \param[in] point The point to which the closest point is searched for.
- * \param[in] data Vector containing the set of points for searching. */
- static unsigned
- getClosestPoint (const ON_2dPoint &point, const vector_vec2d &data);
-
- /** \brief Get the closest point with respect to 'point'
- * \param[in] point The point to which the closest point is searched for.
- * \param[in] data Vector containing the set of points for searching. */
- static unsigned
- getClosestPoint (const ON_3dPoint &point, const vector_vec3d &data);
-
- /** \brief Get the closest point with respect to 'point' in Non-Euclidean
metric
- * \brief Related paper: TODO
- * \param[in] point The point to which the closest point is searched for.
- * \param[in] dir The direction defining 'inside' and 'outside'
- * \param[in] data Vector containing the set of points for searching.
- * \param[out] idxcp Closest point with respect to Euclidean metric. */
- static unsigned
- getClosestPoint (const ON_2dPoint &point, const ON_2dVector &dir, const
vector_vec2d &data,
- unsigned &idxcp);
-
- /** \brief Compute the mean of a set of points
- * \param[in] data Set of points. */
- static ON_3dVector
- computeMean (const vector_vec3d &data);
- /** \brief Compute the mean of a set of points
- * \param[in] data Set of points. */
- static ON_2dVector
- computeMean (const vector_vec2d &data);
-
- /** \brief PCA - principal-component-analysis
- * \param[in] data Set of points.
- * \param[out] mean The mean of the set of points.
- * \param[out] eigenvectors Matrix containing column-wise the
eigenvectors of the set of points.
- * \param[out] eigenvalues The eigenvalues of the set of points with
respect to the eigenvectors. */
- static void
- pca (const vector_vec3d &data, ON_3dVector &mean, Eigen::Matrix3d
&eigenvectors,
- Eigen::Vector3d &eigenvalues);
-
- /** \brief Downsample data points to a certain size.
- * \param[in] data1 The original set of points.
- * \param[out] data2 The downsampled set of points of size 'size'.
- * \param[in] size The desired size of the resulting set of points.
*/
- static void
- downsample_random (const vector_vec3d &data1, vector_vec3d &data2,
unsigned size);
- /** \brief Downsample data points to a certain size.
- * \param[in/out] data1 The set of points for downsampling;
- * will be replaced by the resulting set of points of size 'size'.
- * \param[in] size The desired size of the resulting set of points.
*/
- static void
- downsample_random (vector_vec3d &data1, unsigned size);
-
-};
-
-
-// **********************
-// * from nurbs_solve.h *
-// **********************
-
-/** \brief Solving the linear system of equations using Eigen or UmfPack.
- * (can be defined in on_nurbs.cmake)*/
-class NurbsSolve
-{
-public:
- /** \brief Empty constructor */
- NurbsSolve () :
- m_quiet (true)
- {
- }
-
- /** \brief Assign size and dimension (2D, 3D) of system of equations. */
- void
- assign (unsigned rows, unsigned cols, unsigned dims);
-
- /** \brief Set value for system matrix K (stiffness matrix, basis
functions) */
- void
- K (unsigned i, unsigned j, double v);
- /** \brief Set value for state vector x (control points) */
- void
- x (unsigned i, unsigned j, double v);
- /** \brief Set value for target vector f (force vector) */
- void
- f (unsigned i, unsigned j, double v);
-
- /** \brief Get value for system matrix K (stiffness matrix, basis
functions) */
- double
- K (unsigned i, unsigned j);
- /** \brief Get value for state vector x (control points) */
- double
- x (unsigned i, unsigned j);
- /** \brief Get value for target vector f (force vector) */
- double
- f (unsigned i, unsigned j);
-
- /** \brief Resize target vector f (force vector) */
- void
- resizeF (unsigned rows);
-
- /** \brief Print system matrix K (stiffness matrix, basis functions) */
- void
- printK ();
- /** \brief Print state vector x (control points) */
- void
- printX ();
- /** \brief Print target vector f (force vector) */
- void
- printF ();
-
- /** \brief Solves the system of equations with respect to x.
- * - Using UmfPack incredibly speeds up this function. */
- bool
- solve ();
-
- /** \brief Compute the difference between solution K*x and target f */
- ON_Matrix *diff ();
-
- /** \brief Enable/Disable debug outputs in console. */
- inline void
- setQuiet (bool val)
- {
- m_quiet = val;
- }
-
-private:
- bool m_quiet;
- Eigen::SparseMatrix<double> m_Ksparse;
- ON_Matrix *m_xeig;
- ON_Matrix *m_feig;
-
-};
-
-
-// ******************************
-// * from fitting_surface_pdm.h *
-// ******************************
-
-/** \brief Fitting a B-Spline surface to 3D point-clouds using
point-distance-minimization
- * Based on paper: TODO
- * \author Thomas Moerwald
- * \ingroup surface */
-class FittingSurface
-{
-public:
- ON_NurbsSurface m_nurbs;
- NurbsDataSurface *m_data;
-
- class myvec
- {
- public:
- int side;
- double hint;
- myvec (int in_side, double in_hint)
- {
- this->side = in_side;
- this->hint = in_hint;
- }
- };
-
- /** \brief Parameters for fitting */
- struct Parameter
- {
- double interior_weight;
- double interior_smoothness;
- double interior_regularisation;
-
- double boundary_weight;
- double boundary_smoothness;
- double boundary_regularisation;
-
- unsigned regularisation_resU;
- unsigned regularisation_resV;
-
- Parameter (double intW = 1.0, double intS = 0.000001, double intR =
0.0, double bndW = 1.0,
- double bndS = 0.000001, double bndR = 0.0, unsigned regU =
0, unsigned regV = 0) :
- interior_weight (intW), interior_smoothness (intS),
interior_regularisation (intR), boundary_weight (bndW),
- boundary_smoothness (bndS), boundary_regularisation (bndR),
regularisation_resU (regU),
- regularisation_resV (regV)
- {
-
- }
-
- };
-
- /** \brief Constructor initializing with the B-Spline surface given in
argument 2.
- * \param[in] data pointer to the 3D point-cloud data to be fit.
- * \param[in] ns B-Spline surface used for fitting. */
- FittingSurface (NurbsDataSurface *data, const ON_NurbsSurface &ns);
-
- /** \brief Constructor initializing B-Spline surface using
initNurbsPCA(...).
- * \param[in] order the polynomial order of the B-Spline surface.
- * \param[in] data pointer to the 2D point-cloud data to be fit.
- * \param[in] z vector defining front face of surface. */
- FittingSurface (int order, NurbsDataSurface *data, ON_3dVector z =
ON_3dVector(0.0, 0.0, 1.0));
-
- /** \brief Refines surface by inserting a knot in the middle of each
element.
- * \param[in] dim dimension of refinement (0, 1) */
- void
- refine (int dim);
-
- /** \brief Assemble the system of equations for fitting
- * - for large point-clouds this is time consuming.
- * - should be done once before refinement to initialize the starting
points for point inversion. */
- virtual void
- assemble (Parameter param = Parameter ());
-
- /** \brief Solve system of equations using Eigen or UmfPack (can be
defined in on_nurbs.cmake),
- * and updates B-Spline surface if a solution can be obtained. */
- virtual void
- solve (double damp = 1.0);
-
- /** \brief Update surface according to the current system of equations.
- * \param[in] damp damping factor from one iteration to the other. */
- virtual void
- updateSurf (double damp);
-
- /** \brief Set parameters for inverse mapping
- * \param[in] in_max_steps maximum number of iterations.
- * \param[in] in_accuracy stops iteration if specified accuracy is
reached. */
- void
- setInvMapParams (unsigned in_max_steps, double in_accuracy);
-
- /** \brief Get the elements of a B-Spline surface.*/
- static std::vector<double>
- getElementVector (const ON_NurbsSurface &nurbs, int dim);
-
- /** \brief Inverse mapping / point inversion: Given a point pt, this
function finds the closest
- * point on the B-Spline surface using Newtons method and point-distance
(L2-, Euclidean norm).
- * \param[in] nurbs the B-Spline surface.
- * \param[in] pt the point to which the closest point on the surface will
be computed.
- * \param[in] hint the starting point in parametric domain (warning: may
lead to convergence at local minima).
- * \param[in] error the distance between the point pt and p after
convergence.
- * \param[in] p closest point on surface.
- * \param[in] tu the tangent vector at point p in u-direction.
- * \param[in] tv the tangent vector at point p in v-direction.
- * \param[in] maxSteps maximum number of iterations.
- * \param[in] accuracy convergence criteria: if error is lower than
accuracy the function returns
- * \return closest point on surface in parametric domain.*/
- static ON_2dPoint
- inverseMapping (const ON_NurbsSurface &nurbs, const ON_3dPoint &pt, const
ON_2dPoint &hint,
- double &error, ON_3dPoint &p, ON_3dVector &tu, ON_3dVector
&tv, int maxSteps = 100,
- double accuracy = 1e-6, bool quiet = true);
-
- /** \brief Given a point pt, the function finds the closest midpoint of
the elements of the surface.
- * \param[in] nurbs the B-Spline surface.
- * \param[in] pt the point to which the closest midpoint of the elements
will be computed.
- * return closest midpoint in parametric domain. */
- static ON_2dPoint
- findClosestElementMidPoint (const ON_NurbsSurface &nurbs, const ON_3dPoint
&pt);
-
- /** \brief Inverse mapping / point inversion: Given a point pt, this
function finds the closest
- * point on the boundary of the B-Spline surface using Newtons method and
point-distance (L2-, Euclidean norm).
- * \param[in] nurbs the B-Spline surface.
- * \param[in] pt the point to which the closest point on the surface will
be computed.
- * \param[in] error the distance between the point pt and p after
convergence.
- * \param[in] p closest boundary point on surface.
- * \param[in] tu the tangent vector at point p in u-direction.
- * \param[in] tv the tangent vector at point p in v-direction.
- * \param[in] maxSteps maximum number of iterations.
- * \param[in] accuracy convergence criteria: if error is lower than
accuracy the function returns
- * \return closest point on surface in parametric domain.*/
- static ON_2dPoint
- inverseMappingBoundary (const ON_NurbsSurface &nurbs, const ON_3dPoint
&pt, double &error,
- ON_3dPoint &p, ON_3dVector &tu, ON_3dVector &tv,
int maxSteps = 100,
- double accuracy = 1e-6, bool quiet = true);
-
- /** \brief Inverse mapping / point inversion: Given a point pt, this
function finds the closest
- * point on one side of the boundary of the B-Spline surface using Newtons
method and
- * point-distance (L2-, Euclidean norm).
- * \param[in] nurbs the B-Spline surface.
- * \param[in] pt the point to which the closest point on the surface will
be computed.
- * \param[in] side the side of the boundary (NORTH, SOUTH, EAST, WEST)
- * \param[in] hint the starting point in parametric domain (warning: may
lead to convergence at local minima).
- * \param[in] error the distance between the point pt and p after
convergence.
- * \param[in] p closest boundary point on surface.
- * \param[in] tu the tangent vector at point p in u-direction.
- * \param[in] tv the tangent vector at point p in v-direction.
- * \param[in] maxSteps maximum number of iterations.
- * \param[in] accuracy convergence criteria: if error is lower than
accuracy the function returns
- * \return closest point on surface in parametric domain.*/
- static ON_2dPoint
- inverseMappingBoundary (const ON_NurbsSurface &nurbs, const ON_3dPoint
&pt, int side, double hint,
- double &error, ON_3dPoint &p, ON_3dVector &tu,
ON_3dVector &tv,
- int maxSteps = 100, double accuracy = 1e-6, bool
quiet = true);
-
- /** \brief Initializing a B-Spline surface using 4 corners */
- static ON_NurbsSurface
- initNurbs4Corners (int order, ON_3dPoint ll, ON_3dPoint lr, ON_3dPoint ur,
ON_3dPoint ul);
-
- /** \brief Initializing a B-Spline surface using
principal-component-analysis and eigen values */
- static ON_NurbsSurface
- initNurbsPCA (int order, NurbsDataSurface *data, ON_3dVector z =
ON_3dVector(0.0, 0.0, 1.0));
-
- /** \brief Initializing a B-Spline surface using
principal-component-analysis and bounding box of points */
- static ON_NurbsSurface
- initNurbsPCABoundingBox (int order, NurbsDataSurface *data, ON_3dVector z
= ON_3dVector(0.0, 0.0, 1.0));
-
- /** \brief Enable/Disable debug outputs in console. */
- inline void
- setQuiet (bool val)
- {
- m_quiet = val;
- m_solver.setQuiet(val);
- }
-
-protected:
-
- /** \brief Initialisation of member variables */
- void
- init ();
-
- /** \brief Assemble point-to-surface constraints for interior points. */
- virtual void
- assembleInterior (double wInt, unsigned &row);
-
- /** \brief Assemble point-to-surface constraints for boundary points. */
- virtual void
- assembleBoundary (double wBnd, unsigned &row);
-
- /** \brief Add minimization constraint: point-to-surface distance
(point-distance-minimization). */
- virtual void
- addPointConstraint (const ON_2dVector ¶ms, const ON_3dPoint &point,
double weight, unsigned &row);
- // void addBoundaryPointConstraint(double paramU, double paramV, double
weight, unsigned &row);
-
- /** \brief Add minimization constraint: interior smoothness by control
point regularisation. */
- virtual void
- addCageInteriorRegularisation (double weight, unsigned &row);
-
- /** \brief Add minimization constraint: boundary smoothness by control
point regularisation. */
- virtual void
- addCageBoundaryRegularisation (double weight, int side, unsigned &row);
-
- /** \brief Add minimization constraint: corner smoothness by control point
regularisation. */
- virtual void
- addCageCornerRegularisation (double weight, unsigned &row);
-
- /** \brief Add minimization constraint: interior smoothness by derivatives
regularisation. */
- virtual void
- addInteriorRegularisation (int order, int resU, int resV, double weight,
unsigned &row);
-
- /** \brief Add minimization constraint: boundary smoothness by derivatives
regularisation. */
- virtual void
- addBoundaryRegularisation (int order, int resU, int resV, double weight,
unsigned &row);
-
- NurbsSolve m_solver;
-
- bool m_quiet;
-
- std::vector<double> m_elementsU;
- std::vector<double> m_elementsV;
-
- double m_minU;
- double m_minV;
- double m_maxU;
- double m_maxV;
-
- int in_max_steps;
- double in_accuracy;
-
- // index routines
- int
- grc2gl (int I, int J)
- {
- return m_nurbs.CVCount (1) * I + J;
- } // global row/col index to global lexicographic index
- int
- lrc2gl (int E, int F, int i, int j)
- {
- return grc2gl (E + i, F + j);
- } // local row/col index to global lexicographic index
- int
- gl2gr (int A)
- {
- return (int)(A / m_nurbs.CVCount (1));
- } // global lexicographic in global row index
- int
- gl2gc (int A)
- {
- return (int)(A % m_nurbs.CVCount (1));
- } // global lexicographic in global col index
-};
-
-
-}
-#endif
-
-/*
- * Local Variables:
- * tab-width: 8
- * mode: C++
- * indent-tabs-mode: t
- * c-file-style: "stroustrup"
- * coding: utf-8
- * End:
- * ex: shiftwidth=4 tabstop=8
- */
Deleted: brlcad/trunk/src/libbrep/surface_tree_queue_tests.patch
===================================================================
--- brlcad/trunk/src/libbrep/surface_tree_queue_tests.patch 2020-08-19
00:59:34 UTC (rev 76839)
+++ brlcad/trunk/src/libbrep/surface_tree_queue_tests.patch 2020-08-19
01:05:00 UTC (rev 76840)
@@ -1,2159 +0,0 @@
-Index: src/libbrep/libbrep_curvetree.cpp
-===================================================================
---- src/libbrep/libbrep_curvetree.cpp (revision 0)
-+++ src/libbrep/libbrep_curvetree.cpp (revision 57887)
-@@ -0,0 +1,622 @@
-+/* L I B B R E P _ C U R V E T R E E . C P P
-+ * BRL-CAD
-+ *
-+ * Copyright (c) 2013-2020 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 libbrep_curvetree.cpp
-+ *
-+ * Brief description
-+ *
-+ */
-+
-+#include "libbrep_curvetree.h"
-+
-+#include "vmath.h"
-+#include <iostream>
-+
-+#include <vector>
-+#include <map>
-+#include <set>
-+#include <queue>
-+#include <limits>
-+
-+// Return 0 if there are no tangent points in the interval, 1
-+// if there is a horizontal tangent, two if there is a vertical
-+// tangent, 3 if a normal split is in order.
-+int ON_Curve_Has_Tangent(const ON_Curve* curve, double min, double max) {
-+
-+ bool tanx1, tanx2, x_changed;
-+ bool tany1, tany2, y_changed;
-+ bool slopex, slopey;
-+ double xdelta, ydelta;
-+ ON_3dVector tangent1, tangent2;
-+ ON_3dPoint p1, p2;
-+ ON_Interval t(min, max);
-+
-+ 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);
-+
-+ x_changed =(tanx1 != tanx2);
-+ y_changed =(tany1 != tany2);
-+
-+ if (x_changed && y_changed) return 3; //horz & vert
-+ if (x_changed) return 1;//need to get vertical tangent
-+ if (y_changed) return 2;//need to find horizontal tangent
-+
-+ 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 we have no slope change
-+ // in x or y, we have a tangent line
-+ if (NEAR_ZERO(xdelta, TOL) || NEAR_ZERO(ydelta, TOL)) return 0;
-+
-+ if ((slopex != tanx1) || (slopey != tany1)) return 3;
-+
-+ return 0;
-+}
-+
-+int ON_Curve_HV_Split(const ON_Curve *curve, double min, double max, double
*mid) {
-+ // Don't split if min ~= max - that's infinite loop territory
-+ if (NEAR_ZERO(max - min, ON_ZERO_TOLERANCE)) return 0;
-+ // Go HV tangent hunting
-+ int tan_result = ON_Curve_Has_Tangent(curve, min, max);
-+ if (tan_result == 3) return 1;
-+ if (tan_result == 0) return 0;
-+ if (tan_result == 1 || tan_result == 2) {
-+ // If we have a tangent point at the min or max, that's OK
-+ ON_3dVector hv_tangent = curve->TangentAt(max);
-+ if (NEAR_ZERO(hv_tangent.x, TOL2) || NEAR_ZERO(hv_tangent.y, TOL2))
tan_result = 0;
-+ hv_tangent = curve->TangentAt(min);
-+ if (NEAR_ZERO(hv_tangent.x, TOL2) || NEAR_ZERO(hv_tangent.y, TOL2))
tan_result = 0;
-+ // The tangent point is in a bad spot - find it and subdivide on it
-+ if(tan_result) {
-+ bool tanmin;
-+ (tan_result == 1) ? (tanmin = (hv_tangent[X] < 0.0)) : (tanmin =
(hv_tangent[Y] < 0.0));
-+ double hv_mid;
-+ double hv_min = min;
-+ double hv_max = max;
-+ // find the tangent point using binary search
-+ while (fabs(hv_min - hv_max) > TOL2) {
-+ hv_mid = (hv_max + hv_min)/2.0;
-+ hv_tangent = curve->TangentAt(hv_mid);
-+ if ((tan_result == 1 && NEAR_ZERO(hv_tangent[X], TOL2))
-+ ||(tan_result == 2 && NEAR_ZERO(hv_tangent[Y], TOL2))) {
-+ (*mid) = hv_mid;
-+ return 1;
-+ }
-+ if ((tan_result == 1 && ((hv_tangent[X] < 0.0) == tanmin)) ||
-+ (tan_result == 2 && ((hv_tangent[Y] < 0.0) == tanmin)
)) {
-+ hv_min = hv_mid;
-+ } else {
-+ hv_max = hv_mid;
-+ }
-+ }
-+ return 0;
-+ } else {
-+ // The tangent was at the min and/or max, we don't need to split on
account
-+ // of this case
-+ return 0;
-+ }
-+ }
-+ // Shouldn't get here - something went drastically wrong if we did
-+ return -1;
-+}
-+
-+double ON_Curve_Solve_For_V(const ON_Curve *curve, double u, double min,
double max, double tol) {
-+
-+ ON_3dVector Tan_start, Tan_end;
-+ ON_3dPoint p;
-+ double guess, dT;
-+ double Ta = min;
-+ double Tb = max;
-+ ON_3dPoint A = curve->PointAt(Ta);
-+ ON_3dPoint B = curve->PointAt(Tb);
-+ double dU = fabs(A.x - B.x);
-+ if (dU <= tol) {
-+ if (A.y <= B.y) {return A.y;} else {return B.y;};
-+ }
-+ Tan_start = curve->TangentAt(Ta);
-+ Tan_end = curve->TangentAt(Tb);
-+ dT = Tb - Ta;
-+
-+ /* Use quick binary subdivision until derivatives at end points in 'u'
are within 5 percent */
-+ while (!(fabs(dU) <= tol) && !(fabs(dT) <= tol)) {
-+ guess = Ta + dT/2;
-+ p = curve->PointAt(guess);
-+
-+ /* if we hit 'u' exactly, done deal */
-+ if (fabs(p.x-u) <= SMALL_FASTF) return p.y;
-+
-+ if (p.x > u) {
-+ /* v is behind us, back up the end */
-+ Tb = guess;
-+ B = p;
-+ Tan_end = curve->TangentAt(Tb);
-+ } else {
-+ /* v is in front, move start forward */
-+ Ta = guess;
-+ A = p;
-+ Tan_start = curve->TangentAt(Ta);
-+ }
-+ dT = Tb - Ta;
-+ dU = B.x - A.x;
-+ }
-+
-+ dU = fabs(B.x - A.x);
-+ if (dU <= tol) { //vertical
-+ if (A.y <= B.y) {return A.y;} else {return B.y;};
-+ }
-+
-+ guess = Ta + (u - A.x) * dT/dU;
-+ p = curve->PointAt(guess);
-+
-+ int cnt=0;
-+ while ((cnt < 1000) && (!(fabs(p.x-u) <= tol))) {
-+ if (p.x < u) {
-+ Ta = guess;
-+ A = p;
-+ } else {
-+ Tb = guess;
-+ B = p;
-+ }
-+ dU = fabs(B.x - A.x);
-+ if (dU <= tol) { //vertical
-+ if (A.y <= B.y) {return A.y;} else {return B.y;};
-+ }
-+
-+ dT = Tb - Ta;
-+ guess =Ta + (u - A[X]) * dT/dU;
-+ p = curve->PointAt(guess);
-+ cnt++;
-+ }
-+ if (cnt > 999) {
-+ std::cout << "ON_Curve_Solve_For_V(): estimate of 'v' given a trim
curve and 'u' did not converge within iteration bound(" << cnt << ")\n";
-+ }
-+ return p.y;
-+}
-+
-+int ON_Curve_Relative_Size(const ON_Curve *curve, const ON_Surface *srf,
double tmin, double tmax, double factor) {
-+ int retval = 0;
-+ ON_Interval u = srf->Domain(0);
-+ ON_Interval v = srf->Domain(1);
-+ ON_3dPoint a(u[0], v[0], 0.0);
-+ ON_3dPoint b(u[1], v[1], 0.0);
-+ double ab = a.DistanceTo(b);
-+ double cd = curve->PointAt(tmin).DistanceTo(curve->PointAt(tmax));
-+ if (cd > factor*ab) retval = 1;
-+ return retval;
-+}
-+
-+long int ON_CurveTree::New_Trim_Node(int parent, int trim_index, int tmin,
int tmax) {
-+ // Parent node of new node
-+ nodes.push_back(parent);
-+ // Curve parameter information
-+ nodes.push_back(tmin);
-+ nodes.push_back(tmax);
-+ // As yet, the node has no children
-+ nodes.push_back(0);
-+ nodes.push_back(0);
-+ // knots/HV tangents info node
-+ nodes.push_back(0);
-+ // Populate uv_bbox_values
-+ nodes.push_back(uv_bbox_values.size());
-+ for (int i = 0; i < 2; ++i)
uv_bbox_values.push_back(std::numeric_limits<double>::infinity());
-+ for (int i = 2; i < 4; ++i) uv_bbox_values.push_back(-1 *
std::numeric_limits<double>::infinity());
-+ // Trim index in brep
-+ nodes.push_back(trim_index);
-+
-+ return nodes.size() - TRIM_NODE_STEP;
-+}
-+
-+unsigned int ON_CurveTree::Depth(long int node_id) {
-+ if (node_id < first_trim_node) return 0;
-+ int depth = 1;
-+ int parent = nodes[node_id];
-+ while (parent != 0 && parent >= first_trim_node) {
-+ depth++;
-+ parent = nodes[parent];
-+ }
-+ return depth;
-+}
-+
-+int ON_Curve_Flat(const ON_Curve *curve, double min, double max, double tol){
-+ int retval = 1;
-+ // If min ~= max, call it flat
-+ if (NEAR_ZERO(max - min, ON_ZERO_TOLERANCE)) return retval;
-+ ON_3dVector tangent_start = curve->TangentAt(min);
-+ ON_3dVector tangent_end = curve->TangentAt(max);
-+ if (fabs(tangent_start * tangent_end) < tol) retval = 0;
-+ return retval;
-+}
-+
-+void ON_CurveTree::Subdivide_Trim_Node(long int first_node_id, const
ON_BrepTrim *trim) {
-+
-+ // Initialize
-+ ON_3dPoint pmin, pmax;
-+ std::set<long int> leaf_nodes;
-+ const ON_Curve *curve = trim->TrimCurveOf();
-+ const ON_Surface *srf = trim->SurfaceOf();
-+ double c_min, c_max;
-+ int knotcnt = curve->SpanCount();
-+ double *knots = (double *)malloc(sizeof(double) * (knotcnt + 1));
-+ (void)curve->GetSpanVector(knots);
-+ (void)curve->GetDomain(&c_min, &c_max);
-+
-+ // Populate the root node t_coords
-+ t_coords.push_back(c_min);
-+ t_coords.push_back(c_max);
-+ nodes.at(first_node_id + 1) = t_coords.size() - 2;
-+ nodes.at(first_node_id + 2) = t_coords.size() - 1;
-+
-+ // Use a queue to process nodes
-+ std::queue<long int> trim_node_queue;
-+
-+ // Prime the queue with the root node
-+ trim_node_queue.push(first_node_id);
-+
-+ // Begin the core of the tree building process. Pop a node
-+ // off the queue, test it to see if it must be split further,
-+ // and proceed accordingly.
-+ while (!trim_node_queue.empty()) {
-+ int split = 0;
-+ long int node = trim_node_queue.front();
-+ trim_node_queue.pop();
-+ double node_tmin = t_coords.at(nodes.at(node+1));
-+ double node_tmax = t_coords.at(nodes.at(node+2));
-+ //std::cout << "tmin: " << node_tmin << ", tmax: " << node_tmax <<
"\n";
-+ ON_Interval node_ival(node_tmin, node_tmax);
-+
-+ // Determine the depth of this node
-+ int node_depth = this->Depth(node);
-+ //std::cout << "Node depth: " << node_depth << "\n";
-+
-+ // Use a variable for the midpoint, since it may be overridden
-+ // if we have knots or tangent points to split on.
-+ double t_mid = node_ival.Mid();
-+
-+ // If we have knots in the interval, we need to split
-+ split += ON_Interval_Find_Split_Knot(&node_ival, knotcnt, knots,
&t_mid);
-+ if (split) nodes.at(node+5) = 1;
-+ if (split) {
-+ //std::cout << "Knot split:" << node << "," << t_mid << "\n";
-+ //if(NEAR_ZERO(t_mid-node_tmin, TOL2) || NEAR_ZERO(t_mid-node_tmax,
TOL2)) std::cout << "Arrgh! Knots\n";
-+ }
-+ // If we aren't already splitting, check for horizontal and vertical
tangents
-+ if (!split) {
-+ split += ON_Curve_HV_Split(curve, node_tmin, node_tmax, &t_mid);
-+ if (split) nodes.at(node+5) = 1;
-+ if (split) {
-+ //std::cout << "HV split:" << node << "," << t_mid << "\n";
-+ //if(NEAR_ZERO(t_mid-node_tmin, TOL2) ||
NEAR_ZERO(t_mid-node_tmax, TOL2)) std::cout << "Arrgh! HV\n";
-+ }
-+ }
-+
-+ // If we aren't already splitting, check flatness
-+ if (!split) {
-+ split += !ON_Curve_Flat(curve, node_tmin, node_tmax,
BREP_CURVE_FLATNESS);
-+ if (split) {
-+ //std::cout << "Flatness:" << node << "\n";
-+ //if(NEAR_ZERO(t_mid-node_tmin, TOL2) ||
NEAR_ZERO(t_mid-node_tmax, TOL2)) std::cout << "Arrgh! Flat\n";
-+ }
-+ }
-+
-+ // We want subdivided curves to be small relative to their
-+ // parent surface.
-+ if (!split) {
-+ split += ON_Curve_Relative_Size(curve, srf, node_tmin, node_tmax,
BREP_TRIM_SUB_FACTOR);
-+ if (split) {
-+ // std::cout << "Relative:" << node << "\n";
-+ //if(NEAR_ZERO(t_mid-node_tmin, TOL2) ||
NEAR_ZERO(t_mid-node_tmax, TOL2)) std::cout << "Arrgh! Relative\n";
-+ }
-+ }
-+
-+ // If we DO need to split, proceed
-+ if (split) {
-+ long int left_child, right_child;
-+ // Push the new midpoint onto tcoords
-+ t_coords.push_back(t_mid);
-+ left_child = this->New_Trim_Node(node, nodes.at(node+7),
nodes.at(node+1), t_coords.size() - 1);
-+ nodes.at(node+3) = left_child;
-+ trim_node_queue.push(left_child);
-+ right_child = this->New_Trim_Node(node, nodes.at(node+7),
t_coords.size() - 1, nodes.at(node+2));
-+ nodes.at(node+4) = right_child;
-+ trim_node_queue.push(right_child);
-+ } else {
-+ double pmin_u, pmin_v, pmax_u, pmax_v;
-+ leaf_nodes.insert(node);
-+ pmin = curve->PointAt(node_tmin);
-+ pmax = curve->PointAt(node_tmax);
-+ (pmin.x > pmax.x) ? (pmin_u = pmax.x) : (pmin_u = pmin.x);
-+ (pmax.x > pmin.x) ? (pmax_u = pmax.x) : (pmax_u = pmin.x);
-+ (pmin.y > pmax.y) ? (pmin_v = pmax.y) : (pmin_v = pmin.y);
-+ (pmax.y > pmin.y) ? (pmax_v = pmax.y) : (pmax_v = pmin.y);
-+ uv_bbox_values.at(nodes.at(node+6)) = pmin_u;
-+ uv_bbox_values.at(nodes.at(node+6)+1) = pmin_v;
-+ uv_bbox_values.at(nodes.at(node+6)+2) = pmax_u;
-+ uv_bbox_values.at(nodes.at(node+6)+3) = pmax_v;
-+ //std::cout << "node " << node << " bbox:" << pmin_u << "," <<
pmin_v << "," << pmax_u << "," << pmax_v << "\n";
-+ }
-+
-+ }
-+ // Clean up
-+ free(knots);
-+
-+ // Sync bboxes up to the root trim node
-+ std::set<long int> node_set;
-+ std::set<long int>::iterator p_it;
-+ std::set<long int> *parent_nodes, *current_nodes, *temp;
-+ double *p_umin, *p_vmin, *p_umax, *p_vmax;
-+ double *c_umin, *c_vmin, *c_umax, *c_vmax;
-+ current_nodes = &leaf_nodes;
-+ parent_nodes = &node_set;
-+ while (!current_nodes->empty()) {
-+ for (p_it = current_nodes->begin(); p_it != current_nodes->end();
++p_it) {
-+ if ((*p_it) > first_node_id) {
-+ long int parent_id = nodes.at((*p_it));
-+ p_umin = &(uv_bbox_values[nodes.at(parent_id+6)]);
-+ p_vmin = &(uv_bbox_values[nodes.at(parent_id+6)+1]);
-+ p_umax = &(uv_bbox_values[nodes.at(parent_id+6)+2]);
-+ p_vmax = &(uv_bbox_values[nodes.at(parent_id+6)+3]);
-+ //std::cout << "pnode " << parent_id << " bbox:" << *p_umin <<
"," << *p_vmin << "," << *p_umax << "," << *p_vmax << "\n";
-+ c_umin = &(uv_bbox_values[nodes.at((*p_it)+6)]);
-+ c_vmin = &(uv_bbox_values[nodes.at((*p_it)+6)+1]);
-+ c_umax = &(uv_bbox_values[nodes.at((*p_it)+6)+2]);
-+ c_vmax = &(uv_bbox_values[nodes.at((*p_it)+6)+3]);
-+ //std::cout << "cnode " << (*p_it) << " bbox:" << *c_umin <<
"," << *c_vmin << "," << *c_umax << "," << *c_vmax << "\n";
-+ if (*c_umin < *p_umin) (*p_umin) = (*c_umin);
-+ if (*c_vmin < *p_vmin) (*p_vmin) = (*c_vmin);
-+ if (*c_umax > *p_umax) (*p_umax) = (*c_umax);
-+ if (*c_vmax > *p_vmax) (*p_vmax) = (*c_vmax);
-+ if (parent_id > first_node_id) parent_nodes->insert(parent_id);
-+ //std::cout << "pnode u" << parent_id << " bbox:" << *p_umin <<
"," << *p_vmin << "," << *p_umax << "," << *p_vmax << "\n";
-+ }
-+ }
-+ current_nodes->clear();
-+ temp = parent_nodes;
-+ parent_nodes = current_nodes;
-+ current_nodes = temp;
-+ }
-+}
-+
-+// When testing a curve for trimming contributions, the first stage is to
check
-+// the curve tree node to see whether it contains a segment from below the
knot
-+// and HV tangent breakdowns of the curve. If it does not, we need to walk
down
-+// the tree until we have the set of nodes that represent corresponding
sub-segments
-+// that are below the knot and HV tangent breakdowns.
-+//
-+// Once we have an appropriate segment set, do the pnpoly intersection checks.
-+// For anything where the point is not in the segment bounding box, the test
is
-+// simple - if the point is within the bounding box of the node, walk down the
-+// tree until either a determination can be made or the point is found to be
-+// within a leaf node.
-+//
-+// If the point is within a leaf node, the more expensive search against the
-+// actual curve to determine 2D ray intersection status is necessary.
-+bool ON_CurveTree::CurveTrim(double u, double v, long int node) {
-+ bool c = false;
-+ std::set<long int> test_nodes;
-+ std::set<long int>::iterator tn_it;
-+ std::queue<long int> q1, q2;
-+ std::queue<long int> *current_queue, *next_queue, *tmp;
-+ // Get the node or nodes we need into test_nodes
-+ q1.push(node);
-+ current_queue = &q1;
-+ next_queue= &q2;
-+ while (!current_queue->empty()) {
-+ while (!current_queue->empty()) {
-+ long int current_node = current_queue->front();
-+ current_queue->pop();
-+ if (nodes.at(current_node+5) != 1) {
-+ test_nodes.insert(current_node);
-+ } else {
-+ next_queue->push(nodes.at(current_node+3));
-+ next_queue->push(nodes.at(current_node+4));
-+ }
-+ }
-+ tmp = current_queue;
-+ current_queue = next_queue;
-+ next_queue = tmp;
-+ }
-+ // Check node bounding boxes vs uv point
-+ double *c_umin, *c_vmin, *c_umax, *c_vmax;
-+ for(tn_it = test_nodes.begin(); tn_it != test_nodes.end(); ++tn_it) {
-+ c_umin = &uv_bbox_values[nodes.at((*tn_it)+6)];
-+ c_vmin = &uv_bbox_values[nodes.at((*tn_it)+6)+1];
-+ c_umax = &uv_bbox_values[nodes.at((*tn_it)+6)+2];
-+ c_vmax = &uv_bbox_values[nodes.at((*tn_it)+6)+3];
-+ // Are we inside the u interval of the bbox?
-+ if ((u <= *c_umax) && (u > *c_umin)) {
-+ // Are we below the bbox? if so, a v+ ray is guaranteed to
intersect
-+ if (v <= *c_vmin) {
-+ c = !c;
-+ } else {
-+ // We're not below it - are we inside it? If not, no
intersection possible,
-+ // otherwise, need closer look
-+ if (v < *c_vmax) {
-+ // Descend the tree until either we're outside
-+ // the bboxes and can make a determination, or
-+ // we find the leaf node containing the uv point
-+ // and do the full curve test.
-+ current_queue->push((*tn_it));
-+ while (!current_queue->empty()) {
-+ while (!current_queue->empty()) {
-+ long int current_node = current_queue->front();
-+ current_queue->pop();
-+ c_umin = &uv_bbox_values[nodes.at(current_node+6)];
-+ c_vmin = &uv_bbox_values[nodes.at(current_node+6)+1];
-+ c_umax = &uv_bbox_values[nodes.at(current_node+6)+2];
-+ c_vmax = &uv_bbox_values[nodes.at(current_node+6)+3];
-+ if ((u <= *c_umax) && (u >= *c_umin)) {
-+ if (v <= *c_vmin) {
-+ c = !c;
-+ } else {
-+ if (v < *c_vmax) {
-+ // If we are at a leaf, it's time to
do a
-+ // close check. Otherwise, queue the
children
-+ if (nodes.at(current_node + 3)) {
-+
next_queue->push(nodes.at(current_node + 3));
-+
next_queue->push(nodes.at(current_node + 4));
-+ } else {
-+ // We're close to a leaf
-+ const ON_Curve *curve =
this->ctree_face->Brep()->m_C2[(int)nodes.at(current_node + 7)];
-+ double v_curve =
ON_Curve_Solve_For_V(curve, u, nodes.at(current_node + 1),
nodes.at(current_node + 2), TOL3);
-+ if (v <= v_curve) c = !c;
-+ }
-+ }
-+ }
-+ }
-+ }
-+ tmp = current_queue;
-+ current_queue = next_queue;
-+ next_queue = tmp;
-+ }
-+ }
-+ }
-+ }
-+ }
-+ return c;
-+}
-+
-+// Remember - unlike other
-+// trimming loops, points inside the outer loop are untrimmed and
-+// points outside it ARE trimmed. Flip test results accordingly.
-+bool ON_CurveTree::IsTrimmed(double u, double v, std::vector<long int>
*subset, bool pnpoly_state) {
-+ //std::cout << "u :" << u << ", v :" << v << "\n";
-+ int last_loop = 0;
-+ std::set<long int> loops_to_check;
-+ long int node_pos = 0;
-+ // First stage - identify which trimming loop(s) we care about, based on
loop UV bounding boxes
-+ while (last_loop != 1) {
-+ double *p_umin = &uv_bbox_values[nodes.at(node_pos+2)];
-+ double *p_vmin = &uv_bbox_values[nodes.at(node_pos+2)+1];
-+ double *p_umax = &uv_bbox_values[nodes.at(node_pos+2)+2];
-+ double *p_vmax = &uv_bbox_values[nodes.at(node_pos+2)+3];
-+ //std::cout << "loop " << node_pos << " bbox:" << (*p_umin) << "," <<
(*p_vmin) << "," << (*p_umax) << "," << (*p_vmax) << "\n";
-+ // Is uv point inside bbox?
-+ if (!( u < (*p_umin) || v < (*p_vmin) || u > (*p_umax) || v >
(*p_vmax) )) loops_to_check.insert(node_pos);
-+ // Jump the starting position to the next loop.
-+ node_pos = nodes.at(node_pos);
-+ if (node_pos == 0) last_loop = 1;
-+ }
-+ std::set<long int>::iterator l_it;
-+ if (loops_to_check.size() > 1) {
-+ std::cout << "Loops possibly containing UV pt (" << u << "," << v <<
"):";
-+ for(l_it = loops_to_check.begin(); l_it != loops_to_check.end();
++l_it) {
-+ std::cout << " " << (*l_it) << " ";
-+ }
-+ std::cout << "\n";
-+ for(l_it = loops_to_check.begin(); l_it != loops_to_check.end(); ++l_it) {
-+ long int current_node_pos = (*l_it) + 3;
-+ long int final_node_pos = (*l_it) + nodes.at((*l_it)+1) + 3;
-+ while (current_node_pos <= final_node_pos) {
-+ //std::cout << "loop " << (*l_it) << ", trim " <<
nodes.at(current_node_pos) << "\n";
-+ current_node_pos++;
-+ }
-+ }
-+ std::cout << "\n";
-+}
-+}
-+
-+ON_CurveTree::ON_CurveTree(const ON_BrepFace* face)
-+{
-+ // Save the face pointer associated with this tree, in case we
-+ // need to access data from the face
-+ ctree_face = face;
-+
-+ // We can find out in advance how many node elements we will have due to
-+ // trimming loops - populate those elements up front.
-+ long int elements = 0;
-+ for (int loop_index = 0; loop_index < face->LoopCount(); loop_index++) {
-+ elements += face->Loop(loop_index)->TrimCount();
-+ }
-+ elements += (face->LoopCount()) * NON_TRIM_LOOP_STEP;
-+ // Import to initialize to zero - used later as stopping criteria
-+ nodes.assign(elements, 0);
-+
-+ // Any node elements after those already in place will pertain to trim
-+ // nodes - to be able to easily distinguish when an index is referring
-+ // to a loop or a trim, record the index which marks the cross-over point
-+ // in the vector. The test is then a simple < or >= check.
-+ first_trim_node = nodes.size();
-+
-+ // Now that we have the setup, add actual information and handle the trims
-+ long int node_pos = 0;
-+ for (int loop_index = 0; loop_index < face->LoopCount(); loop_index++) {
-+ ON_BrepLoop *loop = face->Loop(loop_index);
-+ nodes.at(node_pos + 1) = loop->TrimCount();
-+ // Populate uv_bbox_values
-+ nodes.at(node_pos + 2) = uv_bbox_values.size();
-+ for (int i = 0; i < 2; ++i)
uv_bbox_values.push_back(std::numeric_limits<double>::infinity());
-+ for (int i = 2; i < 4; ++i) uv_bbox_values.push_back(-1 *
std::numeric_limits<double>::infinity());
-+
-+ for (int trim_index = 0; trim_index < loop->TrimCount(); trim_index++) {
-+ int curve_index = loop->Trim(trim_index)->TrimCurveIndexOf();
-+ long int current_trim_node =
node_pos+NON_TRIM_LOOP_STEP+trim_index;
-+ // Add "leaf" trim node from loop
-+ long int trim_root = this->New_Trim_Node(node_pos, curve_index,
0, 0);
-+ nodes.at(current_trim_node) = trim_root;
-+ // Build KDtree below "leaf" trim
-+ this->Subdivide_Trim_Node(trim_root,
face->Loop(loop_index)->Trim(trim_index));
-+ }
-+ // Now that we've handled the trims, build the loop bbox
-+ double *p_umin = &uv_bbox_values[nodes.at(node_pos+2)];
-+ double *p_vmin = &uv_bbox_values[nodes.at(node_pos+2)+1];
-+ double *p_umax = &uv_bbox_values[nodes.at(node_pos+2)+2];
-+ double *p_vmax = &uv_bbox_values[nodes.at(node_pos+2)+3];
-+ for (int trim_index = 0; trim_index < loop->TrimCount(); trim_index++) {
-+ long int current_trim_node =
nodes.at(node_pos+NON_TRIM_LOOP_STEP+trim_index);
-+ double *c_umin = &uv_bbox_values[nodes.at(current_trim_node+6)];
-+ double *c_vmin = &uv_bbox_values[nodes.at(current_trim_node+6)+1];
-+ double *c_umax = &uv_bbox_values[nodes.at(current_trim_node+6)+2];
-+ double *c_vmax = &uv_bbox_values[nodes.at(current_trim_node+6)+3];
-+ if (*c_umin < *p_umin) (*p_umin) = (*c_umin);
-+ if (*c_vmin < *p_vmin) (*p_vmin) = (*c_vmin);
-+ if (*c_umax > *p_umax) (*p_umax) = (*c_umax);
-+ if (*c_vmax > *p_vmax) (*p_vmax) = (*c_vmax);
-+ }
-+ std::cout << "loop " << node_pos << " bbox:" << (*p_umin) << "," <<
(*p_vmin) << "," << (*p_umax) << "," << (*p_vmax) << "\n";
-+
-+ // Identify the size of this loop node and store it in the first entry
-+ // to enable "walking" of the loops
-+ // Jump the starting position to the next loop.
-+ if (loop_index != face->LoopCount() - 1)
-+ nodes.at(node_pos) = node_pos + loop->TrimCount() +
NON_TRIM_LOOP_STEP + 1;
-+ node_pos = nodes.at(node_pos);
-+ }
-+}
-+
-+ON_CurveTree::~ON_CurveTree() {
-+}
-+
-+// Local Variables:
-+// tab-width: 8
-+// mode: C++
-+// c-basic-offset: 4
-+// indent-tabs-mode: t
-+// c-file-style: "stroustrup"
-+// End:
-+// ex: shiftwidth=4 tabstop=8
-
-Property changes on: src/libbrep/libbrep_curvetree.cpp
-___________________________________________________________________
-Added: svn:mime-type
- + text/plain
-Added: svn:eol-style
- + native
-
-Index: src/libbrep/libbrep_surfacetree.cpp
-===================================================================
---- src/libbrep/libbrep_surfacetree.cpp (revision 0)
-+++ src/libbrep/libbrep_surfacetree.cpp (revision 57887)
-@@ -0,0 +1,439 @@
-+/* L I B B R E P _ S U R F A C E T R E E . C P P
-+ * BRL-CAD
-+ *
-+ * Copyright (c) 2013-2020 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 libbrep_surfacetree.cpp
-+ *
-+ * Brief description
-+ *
-+ */
-+
-+#include <vector>
-+#include <map>
-+#include <set>
-+#include <queue>
-+#include <iostream>
-+
-+#include "libbrep_surfacetree.h"
-+
-+long int ON_SurfaceTree::New_Node(int parent, int umin, int umax, int vmin,
int vmax) {
-+ // Parent node of new node
-+ nodes.push_back(parent);
-+ // UV information
-+ nodes.push_back(umin);
-+ nodes.push_back(umax);
-+ nodes.push_back(vmin);
-+ nodes.push_back(vmax);
-+ // As yet, the node has no children.
-+ for (int i = 0; i < 4; ++i) nodes.push_back(0);
-+ ON_BoundingBox surf_bbox;
-+ bboxes.Append(surf_bbox);
-+ nodes.push_back(bboxes.Count() - 1);
-+ nodes.push_back(0); // We haven't tested for trimming curves yet.
-+ nodes.push_back(0); // No temporary nodes as yet.
-+ return nodes.size() - SURF_NODE_STEP;
-+}
-+
-+// Split a node to produce new nodes in both U and V directions
-+void ON_SurfaceTree::Node_Split_UV(double u_mid, double v_mid, std::pair<int,
std::vector<int> *> *node, std::queue<std::vector<int> *> *frame_arrays,
std::queue<std::pair<int, std::vector<int> *> > *node_queue) {
-+ std::vector<int> *parent_frame_index = node->second;
-+
-+ // Four new frame index arrays will be needed
-+ std::vector<int> *q1f = NULL;
-+ std::vector<int> *q2f = NULL;
-+ std::vector<int> *q3f = NULL;
-+ std::vector<int> *q4f = NULL;
-+
-+ // Add the midpoint U and V coordinates, which will be the
-+ // only UV values not already present in the double array
-+ uv_coords.push_back(u_mid);
-+ int umid_id = uv_coords.size() - 1;
-+ uv_coords.push_back(v_mid);
-+ int vmid_id = uv_coords.size() - 1;
-+
-+ // Create child nodes. Reuse frame arrays if any are available to reuse
-+
-+ // Create the SW node
-+ long int q1 = this->New_Node(node->first, nodes[node->first+1], umid_id,
nodes[node->first+3], vmid_id);
-+ nodes[node->first+5] = q1;
-+ if (!frame_arrays->empty()) {q1f = frame_arrays->front();
frame_arrays->pop();} else {q1f = new std::vector<int>;}
-+ q1f->assign(13, -1);
-+ q1f->at(0) = parent_frame_index->at(0);
-+ q1f->at(1) = parent_frame_index->at(9);
-+ q1f->at(2) = parent_frame_index->at(4);
-+ q1f->at(3) = parent_frame_index->at(10);
-+ q1f->at(4) = parent_frame_index->at(5);
-+ node_queue->push(std::make_pair(q1,q1f));
-+
-+ // Create the SE node
-+ long int q2 = this->New_Node(node->first, umid_id, nodes[node->first+2],
nodes[node->first+3], vmid_id);
-+ nodes[node->first+6] = q2;
-+ if (!frame_arrays->empty()) {q2f = frame_arrays->front();
frame_arrays->pop();} else {q2f = new std::vector<int>;}
-+ q2f->assign(13, -1);
-+ q2f->at(0) = parent_frame_index->at(9);
-+ q2f->at(1) = parent_frame_index->at(1);
-+ q2f->at(2) = parent_frame_index->at(12);
-+ q2f->at(3) = parent_frame_index->at(4);
-+ q2f->at(4) = parent_frame_index->at(7);
-+ node_queue->push(std::make_pair(q2,q2f));
-+
-+ // Create the NE node
-+ long int q3 = this->New_Node(node->first, umid_id, nodes[node->first+2],
vmid_id, nodes[node->first+4]);
-+ nodes[node->first+7] = q3;
-+ if (!frame_arrays->empty()) {q3f = frame_arrays->front();
frame_arrays->pop();} else {q3f = new std::vector<int>;}
-+ q3f->assign(13, -1);
-+ q3f->at(0) = parent_frame_index->at(4);
-+ q3f->at(1) = parent_frame_index->at(12);
-+ q3f->at(2) = parent_frame_index->at(2);
-+ q3f->at(3) = parent_frame_index->at(11);
-+ q3f->at(4) = parent_frame_index->at(8);
-+ node_queue->push(std::make_pair(q3,q3f));
-+
-+ // Create the NW node
-+ long int q4 = this->New_Node(node->first, nodes[node->first+1], umid_id,
vmid_id, nodes[node->first+4]);
-+ nodes[node->first+8] = q4;
-+ if (!frame_arrays->empty()) {q4f = frame_arrays->front();
frame_arrays->pop();} else {q4f = new std::vector<int>;}
-+ q4f->assign(13, -1);
-+ q4f->at(0) = parent_frame_index->at(10);
-+ q4f->at(1) = parent_frame_index->at(4);
-+ q4f->at(2) = parent_frame_index->at(11);
-+ q4f->at(3) = parent_frame_index->at(3);
-+ q4f->at(4) = parent_frame_index->at(6);
-+ node_queue->push(std::make_pair(q4,q4f));
-+}
-+
-+// Create two child nodes, splitting in the U direction
-+void ON_SurfaceTree::Node_Split_U(double u_mid, std::pair<int,
std::vector<int> *> *node, std::queue<std::vector<int> *> *frame_arrays,
std::queue<std::pair<int, std::vector<int> *> > *node_queue) {
-+ std::vector<int> *parent_frame_index = node->second;
-+
-+ // Two new frame index arrays will be needed
-+ std::vector<int> *q1f = NULL;
-+ std::vector<int> *q2f = NULL;
-+
-+ // Add the midpoint U coordinate, which will be the
-+ // only UV value not already present in the double array
-+ uv_coords.push_back(u_mid);
-+ int umid_id = uv_coords.size() - 1;
-+
-+ // Create child nodes. Reuse frame arrays if any are available to reuse
-+
-+ // Create the W node
-+ long int q1 = this->New_Node(node->first, nodes[node->first+1], umid_id,
nodes[node->first+3], nodes[node->first+4]);
-+ nodes[node->first+5] = q1;
-+ if (!frame_arrays->empty()) {q1f = frame_arrays->front();
frame_arrays->pop();} else {q1f = new std::vector<int>;}
-+ q1f->assign(13, -1);
-+ q1f->at(0) = parent_frame_index->at(0);
-+ q1f->at(1) = parent_frame_index->at(9);
-+ q1f->at(2) = parent_frame_index->at(11);
-+ q1f->at(3) = parent_frame_index->at(3);
-+ q1f->at(10) = parent_frame_index->at(10);
-+ q1f->at(12) = parent_frame_index->at(4);
-+ node_queue->push(std::make_pair(q1,q1f));
-+
-+ // Create the E node
-+ long int q2 = this->New_Node(node->first, umid_id, nodes[node->first+2],
nodes[node->first+3], nodes[node->first+4]);
-+ nodes[node->first+6] = q2;
-+ if (!frame_arrays->empty()) {q2f = frame_arrays->front();
frame_arrays->pop();} else {q2f = new std::vector<int>;}
-+ q2f->assign(13, -1);
-+ q2f->at(0) = parent_frame_index->at(9);
-+ q2f->at(1) = parent_frame_index->at(1);
-+ q2f->at(2) = parent_frame_index->at(2);
-+ q2f->at(3) = parent_frame_index->at(11);
-+ q2f->at(10) = parent_frame_index->at(4);
-+ q2f->at(12) = parent_frame_index->at(12);
-+ node_queue->push(std::make_pair(q2,q2f));
-+
-+}
-+
-+// Create two child nodes, splitting in the V direction
-+void ON_SurfaceTree::Node_Split_V(double v_mid, std::pair<int,
std::vector<int> *> *node, std::queue<std::vector<int> *> *frame_arrays,
std::queue<std::pair<int, std::vector<int> *> > *node_queue) {
-+ std::vector<int> *parent_frame_index = node->second;
-+
-+ // Two new frame index arrays will be needed
-+ std::vector<int> *q1f = NULL;
-+ std::vector<int> *q2f = NULL;
-+
-+ // Add the midpoint V coordinate, which will be the
-+ // only UV value not already present in the double array
-+ uv_coords.push_back(v_mid);
-+ int vmid_id = uv_coords.size() - 1;
-+
-+ // Create child nodes. Reuse frame arrays if any are available to reuse
-+
-+ // Create the S node
@@ Diff output truncated at 100000 characters. @@
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
_______________________________________________
BRL-CAD Source Commits mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/brlcad-commits
