diff --git a/compositemodel/include/GoTools/compositemodel/ParameterizeUtils.h b/compositemodel/include/GoTools/compositemodel/ParameterizeUtils.h new file mode 100644 index 000000000..e6c5065bd --- /dev/null +++ b/compositemodel/include/GoTools/compositemodel/ParameterizeUtils.h @@ -0,0 +1,74 @@ +/* + * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT, + * Applied Mathematics, Norway. + * + * Contact information: E-mail: tor.dokken@sintef.no + * SINTEF ICT, Department of Applied Mathematics, + * P.O. Box 124 Blindern, + * 0314 Oslo, Norway. + * + * This file is part of GoTools. + * + * GoTools is free software: you can redistribute it and/or modify + * it under the terms of the GNU Affero General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * GoTools 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 Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public + * License along with GoTools. If not, see + * . + * + * In accordance with Section 7(b) of the GNU Affero General Public + * License, a covered work must retain the producer line in every data + * file that is created or manipulated using GoTools. + * + * Other Usage + * You can be released from the requirements of the license by purchasing + * a commercial license. Buying such a license is mandatory as soon as you + * develop commercial activities involving the GoTools library without + * disclosing the source code of your own applications. + * + * This file may be used in accordance with the terms contained in a + * written agreement between you and SINTEF ICT. + */ + +#ifndef __PARAMETERIZEUTILS_H +#define __PARAMETERIZEUTILS_H + +#include "GoTools/parametrization/PrOrganizedPoints.h" +#include "GoTools/utils/Array.h" +#include + +namespace Go +{ + class ftPointSet; + + namespace ParameterizeUtils + { + void parameterizeGraph(shared_ptr graph, + std::vector& uv_pars); + + void parameterizeTriang(const double *xyz_points, int nmbp, + const int *triangles, int nmbt, + std::vector& uv_pars); + + void doParameterize(shared_ptr& op, + std::vector& uv_pars); + + bool recognizeCornerNodes(std::vector& bd_nodes, + std::vector& bd_pnts, + double bd_len, + std::vector& cc); + + void getPointsAndNormal(const double *xyz_points, int nmbp, + const int *triangles, int nmbt, + std::vector& point_and_normal); + } +} + +#endif diff --git a/compositemodel/include/GoTools/compositemodel/SurfaceModelUtils.h b/compositemodel/include/GoTools/compositemodel/SurfaceModelUtils.h index 13308b974..129170426 100644 --- a/compositemodel/include/GoTools/compositemodel/SurfaceModelUtils.h +++ b/compositemodel/include/GoTools/compositemodel/SurfaceModelUtils.h @@ -148,6 +148,24 @@ namespace Go void triangulateFaces(std::vector >& faces, shared_ptr& triang, double tol); + void triangulateModel(shared_ptr& model, double density, + shared_ptr& triang); + + void samplePointsModel(shared_ptr& model, + double density, + shared_ptr& triang); + + void getFaceInnerSamplePoints(shared_ptr& face, + double density, + shared_ptr& points); + + void getBoundarySamplePoints(shared_ptr& model, + double density, + shared_ptr& points); + + void performModelTriangulate(shared_ptr& model, + shared_ptr& triang, + double density); void reduceUnderlyingSurface(shared_ptr& bd_sf, std::vector >& cvs); diff --git a/compositemodel/include/GoTools/compositemodel/ftPointSet.h b/compositemodel/include/GoTools/compositemodel/ftPointSet.h index fb9b3e2ed..608239b4b 100644 --- a/compositemodel/include/GoTools/compositemodel/ftPointSet.h +++ b/compositemodel/include/GoTools/compositemodel/ftPointSet.h @@ -118,6 +118,15 @@ class ftSamplePoint /// Am I on sub surface boundary? bool isOnSubSurfaceBoundary() const { return ((at_boundary_ == 2) || (at_boundary_ == 1)); } + /// Am I on an inner boundary? + bool isOnInnerBoundary() const + { return at_boundary_ == 2; } + + virtual bool isCorner() const + { + return false; + } + /// Get number of neighbours. int getNmbNeighbour() const { return (int)next_.size();} diff --git a/compositemodel/src/ParameterizeUtils.C b/compositemodel/src/ParameterizeUtils.C new file mode 100644 index 000000000..52cd4ff1b --- /dev/null +++ b/compositemodel/src/ParameterizeUtils.C @@ -0,0 +1,372 @@ +/* + * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT, + * Applied Mathematics, Norway. + * + * Contact information: E-mail: tor.dokken@sintef.no + * SINTEF ICT, Department of Applied Mathematics, + * P.O. Box 124 Blindern, + * 0314 Oslo, Norway. + * + * This file is part of GoTools. + * + * GoTools is free software: you can redistribute it and/or modify + * it under the terms of the GNU Affero General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * GoTools 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 Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public + * License along with GoTools. If not, see + * . + * + * In accordance with Section 7(b) of the GNU Affero General Public + * License, a covered work must retain the producer line in every data + * file that is created or manipulated using GoTools. + * + * Other Usage + * You can be released from the requirements of the license by purchasing + * a commercial license. Buying such a license is mandatory as soon as you + * develop commercial activities involving the GoTools library without + * disclosing the source code of your own applications. + * + * This file may be used in accordance with the terms contained in a + * written agreement between you and SINTEF ICT. + */ + +#include "GoTools/compositemodel/ParameterizeUtils.h" +#include "GoTools/compositemodel/ftPointSet.h" +#include "GoTools/parametrization/PrPrmUniform.h" +#include "GoTools/parametrization/PrPrmExperimental.h" +#include "GoTools/parametrization/PrPrmLeastSquare.h" +#include "GoTools/parametrization/PrPrmMeanValue.h" +#include "GoTools/parametrization/PrPrmShpPres.h" +#include "GoTools/parametrization/PrTriangulation_OP.h" +#include "GoTools/parametrization/PrParametrizeBdy.h" +#include "GoTools/parametrization/PrParametrizeInt.h" +#include "GoTools/utils/Point.h" +#include +#include +#include +#include +//#define DEBUG + +using std::vector; +using std::make_shared; +using namespace Go; + +//=========================================================================== +void ParameterizeUtils::parameterizeGraph(shared_ptr graph, + vector& uv_pars) +//=========================================================================== +{ + int num = graph->size(); + ftSamplePoint *first = NULL; + ftSamplePoint *second = NULL; + for (int kj=0; kjisOnBoundary()) + { + first = curr; + vector adj = curr->getNeighbours(); + for (size_t kr=0; krisOnBoundary()) + { + second = adj[kr]; + break; + } + } + if (second) + break; + } + } + + // first and second may have the wrong order + if (first) + graph->setSecond(first); + if (second) + graph->setFirst(second); + + // We must make sure that the ftPointSet has the neighbour + // structure the way the parametrization code expects it. + graph->orderNeighbours(); + + // Parameterize + PrPrmUniform par; + PrParametrizeBdy bdy; + shared_ptr op = shared_ptr(graph); + + doParameterize(op, uv_pars); +} + +//=========================================================================== +void ParameterizeUtils::parameterizeTriang(const double *xyz_points, int nmbp, + const int *triangles, int nmbt, + vector& uv_pars) +//=========================================================================== +{ + shared_ptr prt = + make_shared(xyz_points, nmbp, triangles, nmbt); + PrParametrizeBdy bdy; + shared_ptr op = shared_ptr(prt); + + doParameterize(op, uv_pars); +} + +//=========================================================================== +void ParameterizeUtils::doParameterize(shared_ptr& op, + vector& uv_pars) +//=========================================================================== +{ + PrParametrizeBdy bdy; + bdy.attach(op); + bdy.parametrize(); + + vector bd_nodes; + vector bd_pnts; + int bnode = bdy.findBdyNode(); + bd_nodes.push_back(bnode); + bd_pnts.push_back(op->get3dNode(bnode)); + while (true) + { + int bnode2 = bdy.getNextBdyNode(bnode); + if (bnode2 == bd_nodes[0]) + break; + bnode = bnode2; + bd_nodes.push_back(bnode); + bd_pnts.push_back(op->get3dNode(bnode)); + } +#ifdef DEBUG + std::ofstream ofb("bd.g2"); + for (size_t kj=0; kj cc(4); + double bd_len = bdy.boundaryLength(bd_nodes[0], bd_nodes[0]); + bool found = recognizeCornerNodes(bd_nodes, bd_pnts, bd_len, cc); + if (!found) + { + MESSAGE("WARNING: Corners of parameter domain not found"); + //return; + + + cc[0] = bd_nodes[0]; + cc[1] = bd_nodes[bd_nodes.size()/4]; + cc[2] = bd_nodes[bd_nodes.size()/2]; + cc[3] = bd_nodes[3*bd_nodes.size()/4]; + } + + bdy.parametrize(cc[0],cc[1],cc[2],cc[3]); + +#ifdef DEBUG + vector bd; + vector inner; + int num = op->getNumNodes(); + for (int kj=0; kjisBoundary(kj)) + bd.push_back(op->get3dNode(kj)); + else + inner.push_back(op->get3dNode(kj)); + } + std::ofstream of2("triangvx.g2"); + (void)of2.precision(15); + int k2; + of2 << "400 1 0 4 0 255 0 255"<< std::endl; + of2 << inner.size() << std::endl; + for (k2=0; k2<(int)inner.size(); ++k2) + of2 << inner[k2][0] << " " << inner[k2][1] << " " << inner[k2][2] << std::endl; + of2 << std::endl; + of2 << "400 1 0 4 255 0 0 255"<< std::endl; + of2 << bd.size() << std::endl; + for (k2=0; k2<(int)bd.size(); ++k2) + of2 << bd[k2][0] << " " << bd[k2][1] << " " << bd[k2][2] << std::endl; + of2 << "400 1 0 4 255 0 255 0 "<< std::endl; + of2 << "4" << std::endl; + for (k2=0; k2<4; ++k2) + of2 << op->get3dNode(cc[k2]) << std::endl; +#endif + + //PrPrmUniform intr; + PrPrmMeanValue intr; + //PrPrmLeastSquare intr; + //PrPrmShpPres intr; + //PrPrmExperimental intr; + intr.attach(op); + intr.parametrize(); + + int nmbp = op->getNumNodes(); + uv_pars.reserve(2*nmbp); + for (int ki=0; kigetU(ki)); + uv_pars.push_back(op->getV(ki)); + } +} + +//=========================================================================== +bool ParameterizeUtils::recognizeCornerNodes(vector& bd_nodes, + vector& bd_pnts, + double bd_len, + vector& cc) +//=========================================================================== +{ + if (bd_nodes.size() < 4) + return false; // Do not make a suggestion for a degenerate surface +#ifdef DEBUG + std::ofstream ofbd("boundary_nodes.g2"); + ofbd << "400 1 0 4 255 0 0 255" << std::endl; + ofbd << bd_pnts.size() << std::endl; + for (size_t kh=0; kh 40) ? 2 : 1; + int ki, kj, kr; + vector bd_ang; + //vector bd_ang0; + //for (ki=nmb; ki perm(bd_nodes.size()); + for (ki=0; ki<(int)bd_nodes.size(); ++ki) + perm[ki] = ki; + + for (ki=0; ki<(int)bd_nodes.size(); ++ki) + { + double ang2 = bd_ang[perm[ki]]; // + bd_ang0[perm[ki]]; + for (kj=ki+1; kj<(int)bd_nodes.size(); ++kj) + { + double ang3 = bd_ang[perm[kj]]; // + bd_ang0[perm[kj]]; + if (ang3 > ang2) + { + std::swap(perm[ki], perm[kj]); + ang2 = ang3; + } + } + } + + // Dismiss very close corners + double threshold = 0.05*bd_len; + int nmb_perm = (int)perm.size(); + for (ki=0; ki<4; ki++) + { + for (kj=ki+1; kj& point_and_normal) +//=========================================================================== +{ + shared_ptr prt = + make_shared(xyz_points, nmbp, triangles, nmbt); + + int num = prt->getNumNodes(); + for (int ki=0; kiget3dNode(ki); + Point pos(node[0], node[1], node[2]); + vector neighbours; + prt->getNeighbours(ki, neighbours); + Vector3D node2 = prt->get3dNode(neighbours[neighbours.size()-1]); + Point pos2(node2[0], node2[1], node2[2]); + Point vec1 = pos2 - pos; + Point norm; + for (size_t kj=0; kjget3dNode(neighbours[kj]); + Point pos2(node2[0], node2[1], node2[2]); + Point vec2 = pos2 - pos; + Point norm2 = vec1.cross(vec2); + norm2.normalize(); + if (kj == 0) + norm = norm2; + else + norm += norm2; + vec1 = vec2; + } + norm.normalize(); + point_and_normal.insert(point_and_normal.end(), pos.begin(), pos.end()); + point_and_normal.insert(point_and_normal.end(), norm.begin(), norm.end()); + } +} diff --git a/compositemodel/src/SurfaceModelUtils.C b/compositemodel/src/SurfaceModelUtils.C index b12f7645f..6d6b3b68d 100644 --- a/compositemodel/src/SurfaceModelUtils.C +++ b/compositemodel/src/SurfaceModelUtils.C @@ -40,7 +40,9 @@ #include "GoTools/compositemodel/ftSurface.h" #include "GoTools/compositemodel/CompositeCurve.h" #include "GoTools/compositemodel/ftPointSet.h" +#include "GoTools/compositemodel/ftSurfaceSetPoint.h" #include "GoTools/compositemodel/AdaptSurface.h" +#include "GoTools/compositemodel/cmUtils.h" #include "GoTools/intersections/Identity.h" #include "GoTools/geometry/BoundedSurface.h" #include "GoTools/geometry/BoundedUtils.h" @@ -63,6 +65,8 @@ #include "GoTools/tesselator/TesselatorUtils.h" #include "GoTools/creators/CurveCreators.h" #include "GoTools/geometry/SISLconversion.h" +#include "GoTools/compositemodel/ttlPoint.h" +#include "GoTools/compositemodel/ttlTriang.h" #include "sislP.h" #include @@ -2477,6 +2481,695 @@ void SurfaceModelUtils::triangulateFaces(vector >& faces, } +//=========================================================================== +void SurfaceModelUtils::triangulateModel(shared_ptr& model, double density, + shared_ptr& triang) +//=========================================================================== +{ + samplePointsModel(model, density, triang); + performModelTriangulate(model, triang, density); +} + +//=========================================================================== +void SurfaceModelUtils::samplePointsModel(shared_ptr& model, + double density, + shared_ptr& triang) +//=========================================================================== +{ + // Compute sampling points at edges + getBoundarySamplePoints(model, density, triang); +#ifdef DEBUG + std::ofstream of("bd_nodes.g2"); + of << "400 1 0 4 0 255 0 255" << std::endl; + of << triang->size() << std::endl; + for (int ka=0; kasize(); ++ka) + of << (*triang)[ka]->getPoint() << std::endl; +#endif + + int bd_num = triang->size(); + + // Compute sampling points in the inner of the faces + int num_faces = model->nmbEntities(); + for (int ki=0; ki face = model->getFace(ki); + getFaceInnerSamplePoints(face, density, triang); + } + +#ifdef DEBUG + std::ofstream of2("inner_nodes.g2"); + of2 << "400 1 0 4 255 0 0 255" << std::endl; + of2 << triang->size()-bd_num << std::endl; + for (int ka=bd_num; kasize(); ++ka) + of2<< (*triang)[ka]->getPoint() << std::endl; + + std::ofstream of0("pretrimesh.g2"); + + vector edgs; + for (int ka=0; kasize(); ++ka) + { + Vector3D pos = (*triang)[ka]->getPoint(); + vector next = (*triang)[ka]->getNeighbours(); + for (size_t ki=0; kigetPoint()); + } + } + + of0 << "410 1 0 4 0 0 255 255" << std::endl; + of0 << edgs.size()/2 << std::endl; + for (size_t ki=0; ki& face, + double density, + shared_ptr& points) +//=========================================================================== +{ + shared_ptr surf = face->surface(); + shared_ptr faceb = face; + + // Get resolution + int min_num = 1; + int max_num = 100; + double lenfac = 0.9; + double dfac = 0.4; + + RectDomain dom = surf->containingDomain(); + int cv_dir = 0; + int pt_dir = 1; + double len_u1, len_u2, len_v1, len_v2; + GeometryTools::estimateIsoCurveLength(*surf, false, dom.umin(), len_u1); + GeometryTools::estimateIsoCurveLength(*surf, false, dom.umax(), len_u2); + GeometryTools::estimateIsoCurveLength(*surf, true, dom.vmin(), len_v1); + GeometryTools::estimateIsoCurveLength(*surf, true, dom.vmax(), len_v2); + double frac1 = std::min(len_u1, len_u2)/std::max(len_u1, len_u2); + double frac2 = std::min(len_v1, len_v2)/std::max(len_v1, len_v2); + int bd = 0; // Indicates inner point + double u_size, v_size; + surf->estimateSfSize(u_size, v_size); + + if (lenfac*v_size > u_size || + (v_size > u_size && lenfac*v_size <= u_size && frac1 > frac2)) + { + cv_dir = 1; + pt_dir = 0; + } + + // Fetch constant parameter curves in the selected curve direction + double u1 = (cv_dir == 0) ? dom.umin() : dom.vmin(); + double u2 = (cv_dir == 0) ? dom.umax() : dom.vmax(); + double len = (cv_dir == 0) ? u_size : v_size; + int num_cv = (int)(len/density) - 1; + num_cv = std::min(max_num, std::max(num_cv, min_num)); + double udel = (u2 - u1)/(double)(num_cv+2); + double tol1 = std::max(1.0e-7, 1.0e-5*(u2-u1)); + double par[2]; + + size_t kr; + par[cv_dir] = u1+udel; + double u_stop = u2 - 0.1*udel; // To avoid boundary points due to numerics + while (par[cv_dir] < u_stop) + { + vector > crvs = + surf->constParamCurves(par[cv_dir], cv_dir /*false*/); + + if (crvs.size() == 0) + { + par[cv_dir] += udel; + continue; // Outside domain of surface + } + + // Evaluate sampling points + for (kr=0; krestimatedCurveLength(); + if (len2 < dfac*density) + continue; + int nmb = (int)(len2/density) - 1; + nmb = std::min(max_num, std::max(nmb, min_num)); + double v1 = crvs[kr]->startparam(); + double v2 = crvs[kr]->endparam(); + double vdel = (v2 - v1)/(double)(nmb+1); + double tol2 = std::max(1.0e-7, 1.0e-5*(v2-v1)); + par[pt_dir] = v1 + vdel; + + double v_stop = v2 - 0.1*vdel; + while (par[pt_dir] < v_stop) + { + Point pos = crvs[kr]->point(par[pt_dir]); + Vector3D pnt3D(pos[0], pos[1], pos[2]); + Vector2D pntpar(par[0], par[1]); + shared_ptr ftpnt(new ftSurfaceSetPoint(pnt3D, + bd, + faceb, + pntpar)); + ftpnt->setPar(pntpar); + points->addEntry(ftpnt); + par[pt_dir] += vdel; + } + } + par[cv_dir] += udel; + } +} + +int getVertexIndex(vector >& vertices, + shared_ptr vx) +{ + int ix = -1; + for (size_t ki=0; ki& model, + double density, + shared_ptr& points) +//=========================================================================== +{ + int min_num = 1; + int max_num = 100; + int inner = 2; + int outer = 1; + + // Fetch all vertices in the model and define sample points corresponding + // to the vertices + vector > vertices, bd_vertices; + model->getAllVertices(vertices); + model->getBoundaryVertices(bd_vertices); + + for (size_t ki=0; kigetVertexPoint(); + Vector3D pos3d(vx_pos[0], vx_pos[1], vx_pos[2]); + shared_ptr sfsetpnt(new ftSurfaceSetPoint(pos3d, bd)); + + // Fetch associated faces and corresponding parameter values + vector > face_par = vertices[ki]->getFaces(); + + // Add face and parameter information to surface set point + for (size_t kj=0; kj face2 = model->fetchAsSharedPtr(face_par[kj].first); + shared_ptr face2b = face2; + sfsetpnt->addPair(face2b, par); + } + points->addEntry(sfsetpnt); + } + + vector > bd_edges = model->getBoundaryEdges(); + vector > inner_edges = model->getUniqueInnerEdges(); + for (size_t ki=0; kiface(); + shared_ptr face2_tmp = model->fetchAsSharedPtr(face); + shared_ptr face2 = face2_tmp; + + // Previous sampling point at start vertex + shared_ptr start_vx = bd_edges[ki]->getVertex(true); + int start_ix = getVertexIndex(vertices, start_vx); + ftSamplePoint *prev = (*points)[start_ix]; + + // Compute number of sampling points + double len = bd_edges[ki]->estimatedCurveLength(); + int num_sample = (int)(len/density) - 1; + num_sample = std::min(max_num, std::max(num_sample, min_num)); + + double tmin = bd_edges[ki]->tMin(); + double tmax = bd_edges[ki]->tMax(); + double tdel = (tmax - tmin)/(double)(num_sample + 1); + double tpar = tmin + tdel; + for (int ka=0; kapoint(tpar); + Vector3D pos3d(pos[0], pos[1], pos[2]); + shared_ptr sfsetpnt(new ftSurfaceSetPoint(pos3d, outer)); + + // Associated parameter value + Point par = bd_edges[ki]->faceParameter(tpar);; + sfsetpnt->addPair(face2, Vector2D(par[0], par[1])); + points->addEntry(sfsetpnt); + + sfsetpnt->addNeighbour(prev); + prev->addNeighbour(sfsetpnt.get()); + + if (ka == num_sample-1) + { + // Next sampling point at end vertex + shared_ptr end_vx = bd_edges[ki]->getVertex(false); + int end_ix = getVertexIndex(vertices, end_vx); + ftSamplePoint *next = (*points)[end_ix]; + sfsetpnt->addNeighbour(next); + next->addNeighbour(sfsetpnt.get()); + } + + prev = sfsetpnt.get(); + } + } + + for (size_t ki=0; kiface(); + shared_ptr face_tmp = model->fetchAsSharedPtr(face); + shared_ptr face_shr = face_tmp; + + // Previous sampling point at start vertex + shared_ptr start_vx = inner_edges[ki]->getVertex(true); + int start_ix = getVertexIndex(vertices, start_vx); + ftSamplePoint *prev = (*points)[start_ix]; + + // Compute number of sampling points + double len = inner_edges[ki]->estimatedCurveLength(); + int num_sample = (int)(len/density) - 1; + num_sample = std::min(max_num, std::max(num_sample, min_num)); + + // Adjacent face + ftEdgeBase *twin = inner_edges[ki]->twin(); + ftEdge *twin2 = twin->geomEdge(); + bool turned = (twin2->getVertex(false).get() == start_vx.get()); + ftFaceBase *face2 = twin->face(); + shared_ptr face2_tmp = model->fetchAsSharedPtr(face2); + shared_ptr face2_shr = face2_tmp; + + double tmin = inner_edges[ki]->tMin(); + double tmax = inner_edges[ki]->tMax(); + double tmin2 = twin->tMin(); + double tmax2 = twin->tMax(); + double tdel = (tmax - tmin)/(double)(num_sample + 1); + double tpar = tmin + tdel; + double tpar2; + for (int ka=0; kapoint(tpar); + Vector3D pos3d(pos[0], pos[1], pos[2]); + shared_ptr sfsetpnt(new ftSurfaceSetPoint(pos3d, inner)); + // Associated parameter value + Point par = inner_edges[ki]->faceParameter(tpar);; + sfsetpnt->addPair(face_shr, Vector2D(par[0], par[1])); + + // Parameter value for adjacent surface + double td2 = (tpar - tmin)*(tmax2 - tmin2)/(tmax - tmin); + tpar2 = (turned) ? tmax2 - td2 : tmin2 + td2; + Point par2 = twin2->faceParameter(tpar2); + sfsetpnt->addPair(face2_shr, Vector2D(par2[0], par2[1])); + + points->addEntry(sfsetpnt); + + sfsetpnt->addNeighbour(prev); + prev->addNeighbour(sfsetpnt.get()); + + if (ka == num_sample-1) + { + // Next sampling point at end vertex + shared_ptr end_vx = inner_edges[ki]->getVertex(false); + int end_ix = getVertexIndex(vertices, end_vx); + ftSamplePoint *next = (*points)[end_ix]; + sfsetpnt->addNeighbour(next); + next->addNeighbour(sfsetpnt.get()); + } + + prev = sfsetpnt.get(); + } + + } +} + +ftFaceBase* commonFace(vector& pol_pts) +{ + if (pol_pts.size() < 1) + return 0; + + int num_face = pol_pts[0]->nmbFaces(); + vector faces(num_face); + for (int ka=0; kaface(ka).get(); + + for (size_t ki=1; ki=0; --ka) + { + bool contained = pol_pts[ki]->containsFace(faces[ka]); + if (!contained) + faces.erase(faces.begin()+ka); + } + } + + if (faces.size() == 1) + return faces[0]; + else + return 0; +} + +//=========================================================================== +void SurfaceModelUtils::performModelTriangulate(shared_ptr& model, + shared_ptr& points, + double density) +//=========================================================================== +{ + double eps = 1.0e-9; + double angtol = 0.01; //model->getTolerances().kink; + double tol = model->getTolerances().gap; + double len_lim = 2.0*density; + int num_faces = model->nmbEntities(); + vector > faces_points(num_faces); + + // Fetch the sample points belonging to the individual faces + for (int ki = 0; ki < num_faces; ++ki) + { + // For each surface we rescale parameter domain of spline/underlying surface. + // Parameter values are then mapped accordingly. + shared_ptr face = model->getFace(ki); + shared_ptr surf = face->surface(); + RectDomain rect_dom = surf->containingDomain(); + RectDomain new_dom = cmUtils::geometricParamDomain(surf.get()); + for (size_t kj = 0; kj < points->size(); ++kj) + { + ftSurfaceSetPoint* sspnt = dynamic_cast((*points)[kj]); + if (sspnt == 0) + continue; + + for (int km = 0; km < sspnt->nmbFaces(); ++km) + { + if (face.get() == sspnt->face(km).get()) + { + double u = sspnt->parValue(km)[0]; + double v = sspnt->parValue(km)[1]; + double new_u = (new_dom.umax() - new_dom.umin())*(u - rect_dom.umin())/ + (rect_dom.umax() - rect_dom.umin()) + new_dom.umin(); + double new_v = (new_dom.vmax() - new_dom.vmin())*(v - rect_dom.vmin())/ + (rect_dom.vmax() - rect_dom.vmin()) + new_dom.vmin(); + faces_points[ki].push_back(new ttlPoint(sspnt, new_u, new_v)); + } + } + } + } + + vector triang(faces_points.size()); + bool missing_face_points = false; + for (int ki = 0; ki < num_faces; ++ki) + { + if (faces_points[ki].size() == 0) + { + missing_face_points = true; + break; + } + triang[ki].createDelaunay(faces_points[ki].begin(), + faces_points[ki].end()); +#ifdef DEBUG + bool ok = triang[ki].checkDelaunay(); + std::cout << "Delaunay: " << ok << std::endl; +#endif + } + + if (missing_face_points) + { + // Cleanup + for (int ki = 0; ki < num_faces; ++ki) + { + for (size_t kj = 0; kj < faces_points[ki].size(); ++kj) + { + if (faces_points[ki][kj]) + delete faces_points[ki][kj]; // No more need for object. + } + } + return; + } + +#ifdef DEBUG + std::ofstream of("tri_edgs.g2"); + std::ofstream ofn("not_transferred.g2"); +#endif + // We run through the vector, updating structure for each face. + vector l_bdtri; + for (size_t kr = 0; kr < triang.size(); ++kr) + { + const list& l_edges = triang[kr].getLeadingEdges(); + + // For triangulation of current face, we run through all triangles. + list::const_iterator leading_edge_it = l_edges.begin(); + int num = 0; + while (leading_edge_it != l_edges.end()) + { + // Pre check + hetriang::Edge* tmp_edge = *leading_edge_it; + int num_bd = 0; + shared_ptr tmp_node = tmp_edge->getSourceNode(); + bool tmp_bd = (tmp_node->pointIter()->isOnBoundary() || + tmp_node->pointIter()->isOnSubSurfaceBoundary()); + if (tmp_bd) + num_bd++; + for (int km = 1; km < 3; ++km) + { + tmp_edge = tmp_edge->getNextEdgeInFace(); + tmp_node = tmp_edge->getSourceNode(); + tmp_bd = (tmp_node->pointIter()->isOnBoundary() || + tmp_node->pointIter()->isOnSubSurfaceBoundary()); + if (tmp_bd) + num_bd++; + } + + bool inner = true; + if (num_bd == 3) + { +// #ifdef DEBUG +// std::cout << "Possible outside triangle" << std::endl; +// #endif + double fac = 1.0/3.0; + hetriang::Edge* tmp_edge = *leading_edge_it; + vector tri_pts(3); + shared_ptr tmp_node = tmp_edge->getSourceNode(); + tri_pts[0] = tmp_node->pointIter()->asSurfaceSetPoint(); + vector par(3); + par[0] = Point(tmp_node->x(), tmp_node->y()); + Vector2D mid_par = fac*Vector2D(par[0][0], par[0][1]); + for (int km = 1; km < 3; ++km) + { + tmp_edge = tmp_edge->getNextEdgeInFace(); + tmp_node = tmp_edge->getSourceNode(); + par[km] = Point(tmp_node->x(), tmp_node->y()); + tri_pts[km] = tmp_node->pointIter()->asSurfaceSetPoint(); + mid_par += fac*Vector2D(par[km][0], par[km][1]); + } + double min_ang = M_PI; + double max_len = 0.0; + for (int km=0; km<3; ++km) + { + int kn1 = (km + 1)%3; + int kn2 = (km + 2)%3; + Point vec1 = par[kn1] - par[km]; + Point vec2 = par[kn2] - par[km]; + double ang = vec1.angle(vec2); + ang = std::min(ang, M_PI-ang); + min_ang = std::min(min_ang, ang); + max_len = std::max(max_len, vec1.length()); + } + if (min_ang < angtol || max_len > len_lim) + inner = false; + // else + // { + // ftFaceBase *face = commonFace(tri_pts); + // if (face) + // { + // shared_ptr surf = face->surface(); + // const Domain& dom = surf->parameterDomain(); + // int inside = dom.isInDomain2(mid_par, eps); + // if (!inside) + // inner = false; + // } + // else + // inner = false; + // } + } + + if (num_bd == 3 && inner) + { + l_bdtri.push_back(*leading_edge_it); + ++leading_edge_it; // We iterate to next triangle. + continue; + } + + num++; + hetriang::Edge* curr_edge = *leading_edge_it; + shared_ptr source_node = curr_edge->getSourceNode(); + shared_ptr target_node = curr_edge->getTargetNode(); + bool bd1 = (source_node->pointIter()->isOnBoundary() || + source_node->pointIter()->isOnSubSurfaceBoundary()); + bool bd2 = (target_node->pointIter()->isOnBoundary() || + target_node->pointIter()->isOnSubSurfaceBoundary()); +#ifdef DEBUG + vector nodept; + Point ptn1(source_node->x(), source_node->y(), source_node->z()); + Point ptn2(target_node->x(), target_node->y(), target_node->z()); + nodept.push_back(ptn1); + nodept.push_back(ptn2); +#endif + for (int km = 0; km < 3; ++km) + { + std::vector neighbours = + source_node->pointIter()->getNeighbours(); + size_t kj; + for (kj = 0; kj < neighbours.size(); ++kj) + if (target_node->pointIter() == neighbours[kj]) + break; + +#ifdef DEBUG + if (kj == neighbours.size()) + { + std::ofstream ofc("curr_edg.g2"); + ofc << "410 0 0 4 155 100 0 255" << std::endl; + ofc << "1" << std::endl; + ofc << ptn1 << " " << ptn2 << std::endl; + + if (bd1 && bd2 && inner) + { + ofn << "410 0 0 4 155 0 100 255" << std::endl; + ofn << "1" << std::endl; + ofn << source_node->pointIter()->getPoint(); + ofn << " " << target_node->pointIter()->getPoint() << std::endl; + } + int stop_edg = 1; + } +#endif + // If break was executed, connection already exists. + // If both nodes are on boundary we dont't make the points they + // are referring to neighbours (as all boundary points already + // have got their maximum of two boundary neighbours). + // We could have allowed two points on a subsurfaceboundary + // to be neighbours, but it would reault in a conflict when + // topology is to be used in the context of a graph. + if (kj == neighbours.size() && (!(bd1 && bd2))) //inner) //(!(bd1 && bd2) || num_bd < 3)) + { + // Add neighbour + (source_node->pointIter())-> + addNeighbour(target_node->pointIter()); + target_node->pointIter()-> + addNeighbour(source_node->pointIter()); + } + + if (km == 2) + break; + curr_edge = curr_edge->getNextEdgeInFace(); + source_node = curr_edge->getSourceNode(); + target_node = curr_edge->getTargetNode(); +#ifdef DEBUG + Point ptn1(source_node->x(), source_node->y(), source_node->z()); + Point ptn2(target_node->x(), target_node->y(), target_node->z()); + nodept.push_back(ptn1); + nodept.push_back(ptn2); +#endif + } +#ifdef DEBUG + of << "410 1 0 4 0 0 0 255" << std::endl; + of << nodept.size()/2 << std::endl; + for (size_t kj=0; kj curr_node = curr_edge->getSourceNode(); + ftSamplePoint *curr_point = curr_node->pointIter(); + if (curr_point->isOnInnerBoundary()) + num_inner_edge++; + curr_edge = curr_edge->getNextEdgeInFace(); + } + + curr_edge = l_bdtri[kj]; + for (int ka=0; ka<3; ++ka) + { + // Compare number of associated triangles to the current + // vertex with the number of neigbours + shared_ptr curr_node = curr_edge->getSourceNode(); + ftSamplePoint *curr_point = curr_node->pointIter(); + int num_neighbour = curr_point->getNmbNeighbour(); + vector > curr_tri; + curr_point->getAttachedTriangles(curr_tri); + int out = (curr_point->isOnBoundary() && num_inner_edge<3); + if ((int)curr_tri.size() < num_neighbour-out) + num_miss++; + + curr_edge = curr_edge->getNextEdgeInFace(); + } + if (num_miss == 3) + { + // An associated face with more than three edges + hetriang::Edge* curr_edge = l_bdtri[kj]; + shared_ptr source_node = curr_edge->getSourceNode(); + shared_ptr target_node = curr_edge->getTargetNode(); + for (int km = 0; km < 3; ++km) + { + std::vector neighbours = + source_node->pointIter()->getNeighbours(); + size_t kj; + for (kj = 0; kj < neighbours.size(); ++kj) + if (target_node->pointIter() == neighbours[kj]) + break; + if (kj == neighbours.size()) + { + // Add neighbour + (source_node->pointIter())-> + addNeighbour(target_node->pointIter()); + target_node->pointIter()-> + addNeighbour(source_node->pointIter()); + } + + if (km == 2) + break; + curr_edge = curr_edge->getNextEdgeInFace(); + source_node = curr_edge->getSourceNode(); + target_node = curr_edge->getTargetNode(); + } + } + } + + // Cleanup + for (int ki = 0; ki < num_faces; ++ki) + { + for (size_t kj = 0; kj < faces_points[ki].size(); ++kj) + { + if (faces_points[ki][kj]) + delete faces_points[ki][kj]; // No more need for object. + } + } +} + + + //=========================================================================== void SurfaceModelUtils::reduceUnderlyingSurface(shared_ptr& bd_sf, diff --git a/compositemodel/src/ftPointSet.C b/compositemodel/src/ftPointSet.C index 0290521c0..f9b128185 100644 --- a/compositemodel/src/ftPointSet.C +++ b/compositemodel/src/ftPointSet.C @@ -200,25 +200,64 @@ bool ftSamplePoint::isNeighbour(ftSamplePoint* other) const pnt = next_[ki]->next_[kj]; size_t kh = 0; for (kh=0; khnext_.size(); ++kh) - if (pnt->next_[kh] == this) - { - // A triangle is found. Arrange points after increasing - // index - vector index(3); - index[0] = index_; - index[1] = next_[ki]->index_; - index[2] = pnt->index_; - std::sort(index.begin(), index.end()); - auto it = std::find(triangles.begin(), triangles.end(), index); - - if (it == triangles.end()) - triangles.push_back(index); - break; - } + { + size_t curr_size = triangles.size(); + if (pnt->next_[kh] == this) + { + // A triangle is found. + vector index(3); + index[0] = index_; + index[1] = next_[ki]->index_; + index[2] = pnt->index_; + + // Arrange points after increasing index + std::sort(index.begin(), index.end()); + auto it = std::find(triangles.begin(), triangles.end(), index); + + if (it == triangles.end()) + { + // Check if the triangle contains an inner node + bool tri_OK = true; + vector adj1(next_.begin(), next_.end()); + vector adj2(pnt->next_.begin(), pnt->next_.end()); + vector adj3(pnt->next_[kh]->next_.begin(), + pnt->next_[kh]->next_.end()); + std::sort(adj1.begin(), adj1.end()); + std::sort(adj2.begin(), adj2.end()); + vector common1_2; + std::set_intersection(adj1.begin(), adj1.end(), + adj2.begin(), adj2.end(), + std::back_inserter(common1_2)); + if (common1_2.size() > 0) + { + std::sort(adj3.begin(), adj3.end()); + vector common; + std::set_intersection(common1_2.begin(), common1_2.end(), + adj3.begin(), adj3.end(), + std::back_inserter(common)); + for (size_t kr=0; krgetNmbNeighbour(); + if (num_adj == 3) + { + tri_OK = false; + break; + } + } + } + + if (tri_OK) + { + triangles.push_back(index); + break; + } + } + + } + } // if (khnext_.size()) // break; } - } //=========================================================================== diff --git a/gotools-core/app/creators/approxPointSeq.C b/gotools-core/app/creators/approxPointSeq.C new file mode 100644 index 000000000..3caef42ff --- /dev/null +++ b/gotools-core/app/creators/approxPointSeq.C @@ -0,0 +1,93 @@ +/* + * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT, + * Applied Mathematics, Norway. + * + * Contact information: E-mail: tor.dokken@sintef.no + * SINTEF ICT, Department of Applied Mathematics, + * P.O. Box 124 Blindern, + * 0314 Oslo, Norway. + * + * This file is part of GoTools. + * + * GoTools is free software: you can redistribute it and/or modify + * it under the terms of the GNU Affero General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * GoTools 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 Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public + * License along with GoTools. If not, see + * . + * + * In accordance with Section 7(b) of the GNU Affero General Public + * License, a covered work must retain the producer line in every data + * file that is created or manipulated using GoTools. + * + * Other Usage + * You can be released from the requirements of the license by purchasing + * a commercial license. Buying such a license is mandatory as soon as you + * develop commercial activities involving the GoTools library without + * disclosing the source code of your own applications. + * + * This file may be used in accordance with the terms contained in a + * written agreement between you and SINTEF ICT. + */ + +#include "GoTools/creators/ApproxCurve.h" +#include "GoTools/geometry/SplineCurve.h" +#include "GoTools/geometry/ObjectHeader.h" +#include "GoTools/geometry/PointCloud.h" +#include "GoTools/geometry/Utils.h" + +#include + + +using namespace Go; +using std::vector; + + +int main(int argc, char* argv[]) +{ + if (argc != 6) { + MESSAGE("Usage: input point file, approximation tolerance, degree, initial no of coefs, output curve file."); + return 0; + } + + // Read input arguments + std::ifstream filein(argv[1]); + ALWAYS_ERROR_IF(filein.bad(), "Input file not found or file corrupt"); + double tol = atof(argv[2]); + int degree = atoi(argv[3]); + int num_coef = atoi(argv[4]); + std::ofstream fileout(argv[5]); + int dim = 3; + + // Input surface should be a a PointCloud. + ObjectHeader header; + header.read(filein); + PointCloud3D pointseq; + pointseq.read(filein); + double* raw_data = pointseq.rawData(); + int num_pts = pointseq.numPoints(); + vector points(raw_data, raw_data+3*num_pts); + + // Parameterize + vector param(num_pts); + param[0] = 0.0; + for (int ka=1; ka crv = approx.getApproxCurve(maxdist, avdist, max_iter); + std::cout << "Maxdist: " << maxdist << ", avdist: " << avdist << std::endl; + + crv->writeStandardHeader(fileout); + crv->write(fileout); +} + diff --git a/gotools-core/include/GoTools/creators/ApproxSurf.h b/gotools-core/include/GoTools/creators/ApproxSurf.h index ca6188f3f..370813545 100644 --- a/gotools-core/include/GoTools/creators/ApproxSurf.h +++ b/gotools-core/include/GoTools/creators/ApproxSurf.h @@ -146,6 +146,46 @@ class ApproxSurf bool repar=true); + /// Constructor where the user specifies a spline surface that should be + /// modified, the points to approximate and their parameter values, as well + /// as the geometric tolerance. The surface that is given as argument is not + /// copied internally, only pointed to, so it \em will be modified. + /// \param srf the surface that will be modified to approximate the points. + /// Assumed to contain k-regular knots. + /// \param points vector containing the coordinates of the points that this + /// surface should interpolate. They are stored in + /// "xyzxyz...-fashion". + /// \param parvals vector containing the parameter values of the points given in + /// the 'points' vector. They are stored in "uvuv...-fashion". + /// \param pointwgts vector containing individual approximation weights for each + /// point + /// \param dim spatial dimension of the points (usually 3). + /// \param aepsge geometric tolerance to use internally + /// \param constdir The points \em will be reparameterized internally according + /// to their spatial position with respect to the surface that + /// shall be generated. However, they might be reparameterized + /// in both their u and v parameters, only in their u parameters + /// or only in their v parameters. The user can specify this + /// with 'constdir'. If 'constdir' is set to 0, the points will + /// be reparameterized in both u and v. If 'constdir' is set to + /// 1, they will only be reparameterized in the v parameter. + /// If 'constdir' is set to 2, they will only be reparameterized in + /// the u parameter. + ApproxSurf(shared_ptr& srf, + const std::vector& points, + const std::vector& parvals, + const std::vector& pointwgts, + int dim, double aepsge, int constdir = 0, + bool approx_orig = false, + bool repar=true); + + + ApproxSurf(const std::vector& points, + const std::vector& parvals, + int order1, int order2, int num_coef1, int num_coef2, + int dim, double aepsge, bool repar=true); + + /// Destructor ~ApproxSurf(); @@ -170,6 +210,11 @@ class ApproxSurf edge_derivs_[0] = edge_derivs_[1] = edge_derivs_[2] = edge_derivs_[3] = fix; } + void setFixCorners(bool fix_corners) + { + corner_fix_ = fix_corners; + } + /// Decide whether specific edges of the surface's boundary should be kept fixed /// (i.e. unchanged by approximation process), as well as a certain number of cross- /// derivatives across these curves. @@ -289,6 +334,7 @@ class ApproxSurf int dim_; std::vector points_; std::vector parvals_; + std::vector pt_weight_; int pts_stabil_; std::vector norm_points_; std::vector norm_parvals_; @@ -298,6 +344,7 @@ class ApproxSurf double c1fac1_, c1fac2_; int acc_criter_; double acc_frac_; + bool corner_fix_; /// Generate an initial curve representing the spline space int makeInitSurf(std::vector > &crvs, diff --git a/gotools-core/include/GoTools/geometry/SplineUtils.h b/gotools-core/include/GoTools/geometry/SplineUtils.h index e16ec9e76..6633e8378 100644 --- a/gotools-core/include/GoTools/geometry/SplineUtils.h +++ b/gotools-core/include/GoTools/geometry/SplineUtils.h @@ -290,6 +290,11 @@ namespace SplineUtils { const std::vector new_knots_u, const std::vector new_knots_v); + void extractMissingKnots(std::vector& union_vec, + std::vector& vec, + double tol, int order, + std::vector& resvec); + } // End of namespace SplineUtils } // End of namespace Go diff --git a/gotools-core/src/creators/ApproxSurf.C b/gotools-core/src/creators/ApproxSurf.C index cf38194a5..5b96e0fbd 100644 --- a/gotools-core/src/creators/ApproxSurf.C +++ b/gotools-core/src/creators/ApproxSurf.C @@ -93,6 +93,7 @@ ApproxSurf::ApproxSurf() c1fac1_ = 0.0; c1fac2_ = 0.0; acc_criter_ = ACCURACY_MAXDIST; + corner_fix_ = false; } //*************************************************************************** @@ -132,6 +133,8 @@ ApproxSurf::ApproxSurf(std::vector >& crvs, repar_ = repar; refine_ = true; mba_ = false; + vector tmp_weight(parvals.size()/2, 1.0); + pt_weight_ = tmp_weight; c1fac1_ = 0.0; c1fac2_ = 0.0; acc_criter_ = ACCURACY_MAXDIST; @@ -142,6 +145,7 @@ ApproxSurf::ApproxSurf(std::vector >& crvs, smoothfac_ = 1.0/((curr_srf_->endparam_u() - curr_srf_->startparam_u()) + (curr_srf_->endparam_v() - curr_srf_->startparam_v())); + corner_fix_ = false; } //*************************************************************************** @@ -178,6 +182,7 @@ ApproxSurf::ApproxSurf(shared_ptr& srf, orig_ = approx_orig; repar_ = repar; refine_ = true; + pt_weight_ = std::vector(parvals.size()/2, 1.0); mba_ = false; c1fac1_ = 0.0; c1fac2_ = 0.0; @@ -191,10 +196,148 @@ ApproxSurf::ApproxSurf(shared_ptr& srf, curr_srf_ = srf; init_srf_ = shared_ptr(srf->clone()); + corner_fix_ = false; } //*************************************************************************** +ApproxSurf::ApproxSurf(shared_ptr& srf, + const std::vector& points, + const std::vector& parvals, + const std::vector& pointwgts, + int dim, double aepsge, int constdir, + bool approx_orig, + bool repar) + //-------------------------------------------------------------------------- + // Constructor for class ApproxSurf. + // + // Purpose : Initialize class variables + // + // Calls : + // + //-------------------------------------------------------------------------- +{ + prevdist_ = maxdist_ = -10000.0; + prevav_ = avdist_ = 0; + outsideeps_ = 0; + dim_ = dim; + aepsge_ = aepsge; + smoothweight_ = 1.0e-3; // 1.0e-9; + constdir_ = constdir; + use_normals_ = false; + close_belt_ = false; + edge_derivs_[0] = edge_derivs_[1] = edge_derivs_[2] = edge_derivs_[3] = 1; + pts_stabil_ = 0; + norm_stabil_ = 0; + orig_ = approx_orig; + repar_ = repar; + refine_ = true; + mba_ = false; + c1fac1_ = 0.0; + c1fac2_ = 0.0; + acc_criter_ = ACCURACY_MAXDIST; + + points_ = points; + parvals_ = parvals; + pt_weight_ = pointwgts; + + smoothfac_ = 1.0/((srf->endparam_u() - srf->startparam_u()) + + (srf->endparam_v() - srf->startparam_v())); + + curr_srf_ = srf; + init_srf_ = shared_ptr(srf->clone()); + corner_fix_ = false; +} + +//*************************************************************************** + + +//*************************************************************************** + +ApproxSurf::ApproxSurf(const std::vector& points, + const std::vector& parvals, + int order1, int order2, int num_coef1, int num_coef2, + int dim, double aepsge, bool repar) + //-------------------------------------------------------------------------- + // Constructor for class ApproxSurf. + // + // Purpose : Initialize class variables + // + // Calls : + // + + //-------------------------------------------------------------------------- +{ + prevdist_ = maxdist_ = -10000.0; + prevav_ = avdist_ = 0; + outsideeps_ = 0; + dim_ = dim; + aepsge_ = aepsge; + smoothweight_ = 1.0e-3; // 1.0e-9; + constdir_ = -1; + use_normals_ = false; + close_belt_ = false; + edge_derivs_[0] = edge_derivs_[1] = edge_derivs_[2] = edge_derivs_[3] = 0; + pts_stabil_ = 0; + norm_stabil_ = 0; + orig_ = false; + repar_ = repar; + refine_ = true; + mba_ = false; + vector tmp_weight(parvals.size()/2, 1.0); + pt_weight_ = tmp_weight; + c1fac1_ = 0.0; + c1fac2_ = 0.0; + acc_criter_ = ACCURACY_MAXDIST; + + points_ = points; + parvals_ = parvals; + + double umin, umax, vmin, vmax; + umin = umax = parvals_[0]; + vmin = vmax = parvals_[1]; + for (size_t kr=2; kr knots_u(order1+num_coef1); + vector knots_v(order2+num_coef2); + + int ki, kj; + for (kj=0; kj coefs(num_coef1*num_coef2*dim, 0.0); + curr_srf_ = shared_ptr(new SplineSurface(num_coef1, num_coef2, + order1, order2, &knots_u[0], + &knots_v[0], &coefs[0], + dim)); + init_srf_ = shared_ptr(curr_srf_->clone()); + + smoothfac_ = 1.0/((umax - umin) + (vmax - vmin)); + corner_fix_ = false; +} + +//*************************************************************************** ApproxSurf::~ApproxSurf() //-------------------------------------------------------------------------- @@ -422,7 +565,7 @@ int ApproxSurf::makeSmoothSurf() normweight /= weight_sum; } double approxweight = 1.0 - wgt1 - wgt2 - wgt3 - normweight; - std::vector pt_weight(parvals_.size()/2, 1.0); + //std::vector pt_weight(parvals_.size()/2, 1.0); double wgt_orig = 0.0; if (orig_) wgt_orig = 0.1*approxweight; @@ -434,11 +577,11 @@ int ApproxSurf::makeSmoothSurf() srfgen.setOptimize(wgt1, wgt2, wgt3); srfgen.setLeastSquares(points_, parvals_, - pt_weight, approxweight); + pt_weight_, approxweight); if (use_normals_) { stat = srfgen.setNormalCond(norm_points_, norm_parvals_, - pt_weight, normweight); + pt_weight_, normweight); if (stat < 0) return stat; } @@ -991,6 +1134,11 @@ void ApproxSurf::setCoefKnown() for (kj = 0; kj < kn2; ++kj) coef_known_[kj*kn1+ki] = 1; + if (corner_fix_) + { + coef_known_[0] = coef_known_[kn1-1] = coef_known_[(kn2-1)*kn1] = + coef_known_[kn1*kn2-1] = 1; + } } diff --git a/gotools-core/src/geometry/SplineUtils.C b/gotools-core/src/geometry/SplineUtils.C index 403916d7c..bc9c4b2dd 100644 --- a/gotools-core/src/geometry/SplineUtils.C +++ b/gotools-core/src/geometry/SplineUtils.C @@ -662,6 +662,33 @@ shared_ptr GO_API SplineUtils::insertKnots(const Go::SplineSurfac } +//============================================================================== +void SplineUtils::extractMissingKnots(vector& union_vec, + vector& vec, + double tol, int order, + vector& resvec) +//============================================================================== +{ + int ki, kj; + int size1 = (int)vec.size() - order; + int size2 = (int)union_vec.size() - order; + for (ki=order, kj=order; ki