#include <HiReconstruction.hpp>
Public Member Functions | |
HiReconstruction (Core *impl, ParallelComm *comm=0, EntityHandle meshIn=0, int minpnts=5, bool recwhole=true) | |
~HiReconstruction () | |
ErrorCode | initialize (bool recwhole) |
ErrorCode | reconstruct3D_surf_geom (int degree, bool interp, bool safeguard, bool reset=false) |
Reconstruct a high order surface on given surface mesh. More... | |
ErrorCode | reconstruct3D_surf_geom (size_t npts, int *degrees, bool *interps, bool safeguard, bool reset=false) |
Reconstruct a high order surface on given surface mesh. More... | |
ErrorCode | reconstruct3D_curve_geom (int degree, bool interp, bool safeguard, bool reset=false) |
Reconstruct a high order curve on given curve mesh. More... | |
ErrorCode | reconstruct3D_curve_geom (size_t npts, int *degrees, bool *interps, bool safeguard, bool reset=false) |
Reconstruct a high order curve on given curve mesh. More... | |
ErrorCode | polyfit3d_walf_surf_vertex (const EntityHandle vid, const bool interp, int degree, int minpnts, const bool safeguard, const int ncoords, double *coords, int *degree_out, const int ncoeffs, double *coeffs) |
Construct vertex based polynomial fitting on a surface mesh. More... | |
ErrorCode | polyfit3d_walf_curve_vertex (const EntityHandle vid, const bool interp, int degree, int minpnts, const bool safeguard, const int ncoords, double *coords, int *degree_out, const int ncoeffs, double *coeffs) |
Construct vertex based polynomial fitting on a curve mesh. More... | |
ErrorCode | hiproj_walf_in_element (EntityHandle elem, const int nvpe, const int npts2fit, const double *naturalcoords2fit, double *newcoords) |
Perform high order projection of points in an element, using estimated geometry by reconstruction class. More... | |
ErrorCode | hiproj_walf_around_vertex (EntityHandle vid, const int npts2fit, const double *coords2fit, double *hiproj_new) |
Perform high order projection of points around a vertex, using estimated geometry by reconstruction class. More... | |
void | walf3d_surf_vertex_eval (const double *local_origin, const double *local_coords, const int local_deg, const double *local_coeffs, const bool interp, const int npts2fit, const double *coords2fit, double *hiproj_new) |
Perform high order projection of points around a center vertex, assume geometry is surface. More... | |
void | walf3d_curve_vertex_eval (const double *local_origin, const double *local_coords, const int local_deg, const double *local_coeffs, const bool interp, const int npts2fit, const double *coords2fit, double *hiproj_new) |
Perform high order projection of points around a center vertex, assume geometry is curve. More... | |
bool | get_fittings_data (EntityHandle vid, GEOMTYPE &geomtype, std::vector< double > &coords, int °ree_out, std::vector< double > &coeffs, bool &interp) |
Get interally stored fitting results. More... | |
Static Public Member Functions | |
static int | estimate_num_ghost_layers (int degree, bool interp=false) |
Protected Member Functions | |
HiReconstruction (const HiReconstruction &source) | |
HiReconstruction & | operator= (const HiReconstruction &right) |
int | estimate_num_rings (int degree, bool interp) |
ErrorCode | vertex_get_incident_elements (const EntityHandle &vid, const int elemdim, std::vector< EntityHandle > &adjents) |
Given a vertex, return the incident elements with dimension elemdim. More... | |
ErrorCode | obtain_nring_ngbvs (const EntityHandle vid, int ring, const int minpnts, Range &ngbvs) |
Get n-ring neighbor vertices, assuming curve/surface mesh, not volume mesh. More... | |
void | initialize_surf_geom (const int degree) |
void | initialize_surf_geom (const size_t npts, const int *degrees) |
void | initialize_3Dcurve_geom (const int degree) |
void | initialize_3Dcurve_geom (const size_t npts, const int *degrees) |
ErrorCode | average_vertex_normal (const EntityHandle vid, double *nrm) |
ErrorCode | compute_average_vertex_normals_surf () |
ErrorCode | get_normals_surf (const Range &vertsh, double *nrms) |
ErrorCode | average_vertex_tangent (const EntityHandle vid, double *tang) |
ErrorCode | compute_average_vertex_tangents_curve () |
ErrorCode | get_tangents_curve (const Range &vertsh, double *tangs) |
void | polyfit3d_surf_get_coeff (const int nverts, const double *ngbcors, const double *ngbnrms, int degree, const bool interp, const bool safeguard, const int ncoords, double *coords, const int ncoeffs, double *coeffs, int *degree_out, int *degree_pnt, int *degree_qr) |
Compute local coordinates system of a vertex, and perform vertex based polynomial fittings of local height function. More... | |
void | eval_vander_bivar_cmf (const int npts2fit, const double *us, const int ndim, double *bs, int degree, const double *ws, const bool interp, const bool safeguard, int *degree_out, int *degree_pnt, int *degree_qr) |
Form and solve Vandermonde system of bi-variables. More... | |
void | polyfit3d_curve_get_coeff (const int nverts, const double *ngbcors, const double *ngbtangs, int degree, const bool interp, const bool safeguard, const int ncoords, double *coords, const int ncoeffs, double *coeffs, int *degree_out) |
Compute local single variable coordinate system of a vertex, and perform vertex based polynomial fittings of three global coordinates axis. More... | |
void | eval_vander_univar_cmf (const int npts2fit, const double *us, const int ndim, double *bs, int degree, const double *ws, const bool interp, const bool safeguard, int *degree_out) |
Form and solve Vandermonde system of single-variables. More... | |
int | compute_weights (const int nrows, const int ncols, const double *us, const int nngbs, const double *ngbnrms, const int degree, const double toler, double *ws) |
Compute weights for points selected in weighted least square fittigns. More... | |
bool | check_barycentric_coords (const int nws, const double *naturalcoords) |
Check the correctness of barycentric coordination, wi>=0 and sum(wi)=1. More... | |
Protected Attributes | |
Core * | mbImpl |
ParallelComm * | pcomm |
HalfFacetRep * | ahf |
EntityHandle | _mesh2rec |
Range | _verts2rec |
Range | _inverts |
Range | _inedges |
Range | _infaces |
Range | _incells |
size_t | _nv2rec |
int | _MAXPNTS |
int | _MINPNTS |
double | _MINEPS |
bool | _hasderiv |
GEOMTYPE | _geom |
int | _dim |
bool | _hasfittings |
bool | _initfittings |
std::vector< double > | _local_coords |
std::vector< double > | _local_fit_coeffs |
std::vector< size_t > | _vertID2coeffID |
std::vector< int > | _degrees_out |
std::vector< bool > | _interps |
Definition at line 39 of file HiReconstruction.hpp.
moab::HiReconstruction::HiReconstruction | ( | Core * | impl, |
ParallelComm * | comm = 0 , |
||
EntityHandle | meshIn = 0 , |
||
int | minpnts = 5 , |
||
bool | recwhole = true |
||
) |
Definition at line 18 of file HiReconstruction.cpp.
19 : mbImpl( impl ), pcomm( comm ), _mesh2rec( meshIn ), _MINPNTS( minpnts )
20 {
21 assert( NULL != impl );
22 ErrorCode error;
23 _MINEPS = 1e-12;
24 _dim = 0;
25 _hasfittings = false;
26 _hasderiv = false;
27 #ifdef MOAB_HAVE_MPI
28 if( !pcomm )
29 {
30 pcomm = moab::ParallelComm::get_pcomm( mbImpl, 0 );
31 }
32 #endif
33
34 error = initialize( recwhole );
35 if( MB_SUCCESS != error )
36 {
37 std::cout << "Error initializing HiReconstruction\n" << std::endl;
38 exit( 1 );
39 }
40 }
References _dim, _hasderiv, _hasfittings, _MINEPS, moab::error(), ErrorCode, moab::ParallelComm::get_pcomm(), initialize(), MB_SUCCESS, mbImpl, and pcomm.
moab::HiReconstruction::~HiReconstruction | ( | ) |
Definition at line 42 of file HiReconstruction.cpp.
43 {
44 #ifdef MOAB_HAVE_AHF
45 ahf = NULL;
46 #else
47 delete ahf;
48 #endif
49 }
References ahf.
|
protected |
|
protected |
Save fitting results of a vertex into internal storage
vid | EntityHandle, a vertex in _mesh2rec, in _verts2rec |
coords | Pointer to double array, global coordinates of local uvw coordinate system axis directions |
degree_out | Integer, order of polynomial fittings for vid |
coeffs | Pointer to double array, coefficients of local polynomial fittings in monomial basis |
interp | Boolean, true = interpolation Compute area weighted average vertex normals for given vertex, assuming surface mesh For arbitrary polygon mesh, use incident two edges of each incident polygon of this vertex to form a triangle, then use these incident "triangles" to compute area weighted average vertex normals |
vid | EntityHandle, vertex in _mesh2rec, might be ghost vertex |
nrm | Pointer to 3-doubles array, preallocated by user |
Definition at line 841 of file HiReconstruction.cpp.
842 {
843 ErrorCode error;
844 std::vector< EntityHandle > adjfaces;
845 error = vertex_get_incident_elements( vid, 2, adjfaces );MB_CHK_ERR( error );
846 int npolys = adjfaces.size();
847 if( !npolys )
848 {
849 MB_SET_ERR( MB_FAILURE, "Vertex has no incident 2D entities" );
850 }
851 else
852 {
853 double v1[3], v2[3], v3[3], a[3], b[3], c[3];
854 nrm[0] = nrm[1] = nrm[2] = 0;
855 for( int i = 0; i < npolys; ++i )
856 {
857 // get incident "triangles"
858 std::vector< EntityHandle > elemconn;
859 error = mbImpl->get_connectivity( &adjfaces[i], 1, elemconn );MB_CHK_ERR( error );
860 EntityHandle pre, nxt;
861 int nvpe = elemconn.size();
862 for( int j = 0; j < nvpe; ++j )
863 {
864 if( vid == elemconn[j] )
865 {
866 pre = j == 0 ? elemconn[nvpe - 1] : elemconn[j - 1];
867 nxt = elemconn[( j + 1 ) % nvpe];
868 break;
869 }
870 }
871 // compute area weighted normals
872 error = mbImpl->get_coords( &pre, 1, a );MB_CHK_ERR( error );
873 error = mbImpl->get_coords( &vid, 1, b );MB_CHK_ERR( error );
874 error = mbImpl->get_coords( &nxt, 1, c );MB_CHK_ERR( error );
875 DGMSolver::vec_linear_operation( 3, 1, c, -1, b, v1 );
876 DGMSolver::vec_linear_operation( 3, 1, a, -1, b, v2 );
877 DGMSolver::vec_crossprod( v1, v2, v3 );
878 DGMSolver::vec_linear_operation( 3, 1, nrm, 1, v3, nrm );
879 }
880 #ifndef NDEBUG
881 assert( DGMSolver::vec_normalize( 3, nrm, nrm ) );
882 #endif
883 }
884 return error;
885 }
References moab::error(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), MB_CHK_ERR, MB_SET_ERR, mbImpl, moab::DGMSolver::vec_crossprod(), moab::DGMSolver::vec_linear_operation(), moab::DGMSolver::vec_normalize(), and vertex_get_incident_elements().
Referenced by compute_average_vertex_normals_surf(), and get_normals_surf().
|
protected |
Compute area weighted average vertex tangent vector for given vertex, assuming curve mesh Use incident two edges of vertex as estimatation of tangent vectors, weighted by length
vid | EntityHandle, vertex in _mesh2rec, might be ghost vertex |
tang | Pointer to 3-doubles array, preallocated by user |
Definition at line 943 of file HiReconstruction.cpp.
944 {
945 ErrorCode error;
946 std::vector< EntityHandle > adjedges;
947 error = vertex_get_incident_elements( vid, 1, adjedges );MB_CHK_ERR( error );
948 int nedges = adjedges.size();
949 if( !nedges )
950 {
951 MB_SET_ERR( MB_FAILURE, "Vertex has no incident edges" );
952 }
953 else
954 {
955 assert( nedges <= 2 );
956 tang[0] = tang[1] = tang[2] = 0;
957 for( int i = 0; i < nedges; ++i )
958 {
959 std::vector< EntityHandle > edgeconn;
960 error = mbImpl->get_connectivity( &adjedges[i], 1, edgeconn );
961 double istr[3], iend[3], t[3];
962 error = mbImpl->get_coords( &( edgeconn[0] ), 1, istr );
963 error = mbImpl->get_coords( &( edgeconn[1] ), 1, iend );
964 DGMSolver::vec_linear_operation( 3, 1, iend, -1, istr, t );
965 DGMSolver::vec_linear_operation( 3, 1, tang, 1, t, tang );
966 }
967 #ifndef NDEBUG
968 assert( DGMSolver::vec_normalize( 3, tang, tang ) );
969 #endif
970 }
971 return error;
972 }
References moab::error(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), MB_CHK_ERR, MB_SET_ERR, mbImpl, moab::DGMSolver::vec_linear_operation(), moab::DGMSolver::vec_normalize(), and vertex_get_incident_elements().
Referenced by compute_average_vertex_tangents_curve(), and get_tangents_curve().
|
protected |
Check the correctness of barycentric coordination, wi>=0 and sum(wi)=1.
Definition at line 1593 of file HiReconstruction.cpp.
1594 {
1595 double sum = 0;
1596 for( int i = 0; i < nws; ++i )
1597 {
1598 if( naturalcoords[i] < -_MINEPS )
1599 {
1600 return false;
1601 }
1602 sum += naturalcoords[i];
1603 }
1604 if( fabs( 1 - sum ) > _MINEPS )
1605 {
1606 return false;
1607 }
1608 else
1609 {
1610 return true;
1611 }
1612 }
References _MINEPS, and moab::sum().
Referenced by hiproj_walf_in_element().
|
protected |
Compute weighted average vertex normals for all vertices in _verts2rec, not including ghost vertices, results are stored interally in _local_coords
Definition at line 887 of file HiReconstruction.cpp.
888 {
889 if( _hasderiv )
890 {
891 return MB_SUCCESS;
892 }
893 ErrorCode error;
894 _local_coords.assign( 9 * _nv2rec, 0 );
895 size_t index = 0;
896 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
897 {
898 error = average_vertex_normal( *ivert, &( _local_coords[9 * index + 6] ) );MB_CHK_ERR( error );
899 }
900 return error;
901 }
References _hasderiv, _local_coords, _nv2rec, _verts2rec, average_vertex_normal(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, MB_CHK_ERR, and MB_SUCCESS.
Referenced by initialize(), and initialize_surf_geom().
|
protected |
Compute weighted average vertex tangent vectors for all vertices in _verts2rec, not including ghost vertices, results are stored interally in _local_coords
Definition at line 974 of file HiReconstruction.cpp.
975 {
976 if( _hasderiv )
977 {
978 return MB_SUCCESS;
979 }
980 ErrorCode error;
981 _local_coords.assign( 3 * _nv2rec, 0 );
982 size_t index = 0;
983 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
984 {
985 error = average_vertex_tangent( *ivert, &( _local_coords[3 * index] ) );MB_CHK_ERR( error );
986 }
987 return error;
988 }
References _hasderiv, _local_coords, _nv2rec, _verts2rec, average_vertex_tangent(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, MB_CHK_ERR, and MB_SUCCESS.
Referenced by initialize(), and initialize_3Dcurve_geom().
|
protected |
Compute weights for points selected in weighted least square fittigns.
Definition at line 1544 of file HiReconstruction.cpp.
1552 {
1553 assert( nrows <= _MAXPNTS && ws );
1554 bool interp = false;
1555 if( nngbs != nrows )
1556 {
1557 assert( nngbs == nrows + 1 );
1558 interp = true;
1559 }
1560 double epsilon = 1e-2;
1561
1562 // First, compute squared distance from each input piont to the center
1563 for( int i = 0; i < nrows; ++i )
1564 {
1565 ws[i] = DGMSolver::vec_innerprod( ncols, us + i * ncols, us + i * ncols );
1566 }
1567
1568 // Second, compute a small correction termt o guard against zero
1569 double h = 0;
1570 for( int i = 0; i < nrows; ++i )
1571 {
1572 h += ws[i];
1573 }
1574 h /= (double)nrows;
1575
1576 // Finally, compute the weights for each vertex
1577 int nzeros = 0;
1578 for( int i = 0; i < nrows; ++i )
1579 {
1580 double costheta = DGMSolver::vec_innerprod( 3, ngbnrms, ngbnrms + 3 * ( i + interp ) );
1581 if( costheta > toler )
1582 {
1583 ws[i] = costheta * pow( ws[i] / h + epsilon, -1 * (double)degree / 2.0 );
1584 }
1585 else
1586 {
1587 ws[i] = 0;
1588 ++nzeros;
1589 }
1590 }
1591 return nzeros;
1592 }
References _MAXPNTS, and moab::DGMSolver::vec_innerprod().
Referenced by polyfit3d_curve_get_coeff(), and polyfit3d_surf_get_coeff().
|
inlinestatic |
Definition at line 277 of file HiReconstruction.hpp.
278 {
279 return 1 + ( interp ? ( ( degree + 1 ) >> 1 ) + ( ( degree + 1 ) & 1 )
280 : ( ( degree + 2 ) >> 1 ) + ( ( degree + 2 ) & 1 ) );
281 };
|
protected |
Definition at line 653 of file HiReconstruction.cpp.
654 {
655 return interp ? ( ( degree + 1 ) >> 1 ) + ( ( degree + 1 ) & 1 ) : ( ( degree + 2 ) >> 1 ) + ( ( degree + 2 ) & 1 );
656 }
Referenced by polyfit3d_walf_curve_vertex(), and polyfit3d_walf_surf_vertex().
|
protected |
Form and solve Vandermonde system of bi-variables.
Definition at line 1176 of file HiReconstruction.cpp.
1187 {
1188 // step 1. adjust the degree according to number of points to fit
1189 int ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1190 while( 1 < degree && npts2fit < ncols )
1191 {
1192 --degree;
1193 ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1194 }
1195 *degree_pnt = degree;
1196
1197 // std::cout << "degree_pnt: " << *degree_pnt << std::endl;
1198
1199 // step 2. construct Vandermonde matrix, stored in columnwise
1200 std::vector< double > V; // V(npts2fit*(ncols+interp)); //double *V_init = new double[npts2fit*(ncols+interp)];
1201 DGMSolver::gen_vander_multivar( npts2fit, 2, us, degree, V );
1202 // remove the first column of 1s if interpolation
1203 if( interp )
1204 {
1205 V.erase( V.begin(), V.begin() + npts2fit );
1206 }
1207 /*double* V;
1208 if(interp){
1209 V = new double[npts2fit*ncols];
1210 std::memcpy(V,V_init+npts2fit,ncols*npts2fit*sizeof(double));
1211 delete [] V_init; V_init = 0;
1212 }else{
1213 V = V_init;
1214 }*/
1215
1216 // step 3. Scale rows to assign different weights to different points
1217 for( int i = 0; i < npts2fit; ++i )
1218 {
1219 for( int j = 0; j < ncols; ++j )
1220 {
1221 V[j * npts2fit + i] *= ws[i];
1222 }
1223 for( int k = 0; k < ndim; ++k )
1224 {
1225 bs[k * npts2fit + i] *= ws[i];
1226 }
1227 }
1228
1229 // step 4. scale columns to reduce condition number
1230 std::vector< double > ts( ncols ); // double *ts = new double[ncols];
1231 DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
1232
1233 // step 5. Perform Householder QR factorization
1234 std::vector< double > D( ncols ); // double *D = new double[ncols];
1235 int rank;
1236 DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
1237
1238 // step 6. adjust degree of fitting according to rank of Vandermonde matrix
1239 int ncols_sub = ncols;
1240 while( rank < ncols_sub )
1241 {
1242 --degree;
1243 if( degree == 0 )
1244 {
1245 // surface is flat, return 0
1246 *degree_out = *degree_qr = degree;
1247 for( int i = 0; i < npts2fit; ++i )
1248 {
1249 for( int k = 0; k < ndim; ++k )
1250 {
1251 bs[k * npts2fit + i] = 0;
1252 }
1253 }
1254 return;
1255 }
1256 else
1257 {
1258 ncols_sub = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1259 }
1260 }
1261 *degree_qr = degree;
1262
1263 // std::cout << "degree_qr: " << *degree_qr << std::endl;
1264
1265 /* DBG
1266 * std::cout<<"before Qtb"<<std::endl;
1267 std::cout<<std::endl;
1268 std::cout<<"bs = "<<std::endl;
1269 std::cout<<std::endl;
1270 for (int k=0; k< ndim; k++){
1271 for (int j=0; j<npts2fit; ++j){
1272 std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
1273 }
1274 }
1275 std::cout<<std::endl;*/
1276
1277 // step 7. compute Q'b
1278 DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
1279
1280 /* DBG
1281 * std::cout<<"after Qtb"<<std::endl;
1282 std::cout<<"bs = "<<std::endl;
1283 std::cout<<std::endl;
1284 for (int k=0; k< ndim; k++){
1285 for (int j=0; j<npts2fit; ++j){
1286 std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
1287 }
1288 }
1289 std::cout<<std::endl;*/
1290
1291 // step 8. perform backward substitution and scale the solution
1292 // assign diagonals of V
1293 for( int i = 0; i < ncols_sub; ++i )
1294 {
1295 V[i * npts2fit + i] = D[i];
1296 }
1297
1298 // backsolve
1299 if( safeguard )
1300 {
1301 // for debug
1302 // std::cout << "ts size " << ts.size() << std::endl;
1303 DGMSolver::backsolve_polyfit_safeguarded( 2, degree, interp, npts2fit, ncols_sub, &( V[0] ), ndim, bs,
1304 &( ts[0] ), degree_out );
1305 }
1306 else
1307 {
1308 DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), 1, bs, &( ts[0] ) );
1309 *degree_out = degree;
1310 }
1311 /*if(V_init){
1312 delete [] V_init;
1313 }else{
1314 delete [] V;
1315 }*/
1316 }
References moab::DGMSolver::backsolve(), moab::DGMSolver::backsolve_polyfit_safeguarded(), moab::DGMSolver::compute_qtransposeB(), moab::DGMSolver::gen_vander_multivar(), moab::DGMSolver::qr_polyfit_safeguarded(), and moab::DGMSolver::rescale_matrix().
Referenced by polyfit3d_surf_get_coeff().
|
protected |
Form and solve Vandermonde system of single-variables.
Definition at line 1434 of file HiReconstruction.cpp.
1443 {
1444 // step 1. determine degree of polynomials to fit according to number of points
1445 int ncols = degree + 1 - interp;
1446 while( npts2fit < ncols && degree >= 1 )
1447 {
1448 --degree;
1449 ncols = degree + 1 - interp;
1450 }
1451 if( !degree )
1452 {
1453 if( interp )
1454 {
1455 for( int icol = 0; icol < ndim; ++icol )
1456 {
1457 bs[icol * npts2fit] = 0;
1458 }
1459 }
1460 for( int irow = 1; irow < npts2fit; ++irow )
1461 {
1462 for( int icol = 0; icol < ndim; ++icol )
1463 {
1464 bs[icol * npts2fit + irow] = 0;
1465 }
1466 }
1467 *degree_out = 0;
1468 return;
1469 }
1470 // step 2. construct Vandermonde matrix
1471 std::vector< double > V; // V(npts2fit*(ncols+interp));
1472 DGMSolver::gen_vander_multivar( npts2fit, 1, us, degree, V );
1473
1474 if( interp )
1475 {
1476 V.erase( V.begin(), V.begin() + npts2fit );
1477 }
1478
1479 // step 3. scale rows with respect to weights
1480 for( int i = 0; i < npts2fit; ++i )
1481 {
1482 for( int j = 0; j < ncols; ++j )
1483 {
1484 V[j * npts2fit + i] *= ws[i];
1485 }
1486 for( int k = 0; k < ndim; ++k )
1487 {
1488 bs[k * npts2fit + i] *= ws[i];
1489 }
1490 }
1491
1492 // step 4. scale columns to reduce condition number
1493 std::vector< double > ts( ncols );
1494 DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
1495
1496 // step 5. perform Householder QR factorization
1497 std::vector< double > D( ncols );
1498 int rank;
1499 DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
1500
1501 // step 6. adjust degree of fitting
1502 int ncols_sub = ncols;
1503 while( rank < ncols_sub )
1504 {
1505 --degree;
1506 if( !degree )
1507 {
1508 // matrix is singular, consider curve on flat plane passing center and orthogonal to
1509 // tangent line
1510 *degree_out = 0;
1511 for( int i = 0; i < npts2fit; ++i )
1512 {
1513 for( int k = 0; k < ndim; ++k )
1514 {
1515 bs[k * npts2fit + i] = 0;
1516 }
1517 }
1518 }
1519 ncols_sub = degree + 1 - interp;
1520 }
1521
1522 // step 7. compute Q'*bs
1523 DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
1524
1525 // step 8. perform backward substitution and scale solutions
1526 // assign diagonals of V
1527 for( int i = 0; i < ncols_sub; ++i )
1528 {
1529 V[i * npts2fit + i] = D[i];
1530 }
1531 // backsolve
1532 if( safeguard )
1533 {
1534 DGMSolver::backsolve_polyfit_safeguarded( 1, degree, interp, npts2fit, ncols, &( V[0] ), ndim, bs, ws,
1535 degree_out );
1536 }
1537 else
1538 {
1539 DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), ndim, bs, &( ts[0] ) );
1540 *degree_out = degree;
1541 }
1542 }
References moab::DGMSolver::backsolve(), moab::DGMSolver::backsolve_polyfit_safeguarded(), moab::DGMSolver::compute_qtransposeB(), moab::DGMSolver::gen_vander_multivar(), moab::DGMSolver::qr_polyfit_safeguarded(), and moab::DGMSolver::rescale_matrix().
Referenced by polyfit3d_curve_get_coeff().
bool moab::HiReconstruction::get_fittings_data | ( | EntityHandle | vid, |
GEOMTYPE & | geomtype, | ||
std::vector< double > & | coords, | ||
int & | degree_out, | ||
std::vector< double > & | coeffs, | ||
bool & | interp | ||
) |
Get interally stored fitting results.
Get fittings results of a vertex, stored internally, results will be writtend to user provided memory
vid | EntityHandle, a vertex in _verts2rec |
geomtype | GEOMTYPE, one of HISURFACE,HI3DCURVE,HI2DCURVE |
coords | vector, global coordinates of local uvw coordinate system axis directions will be appended to the end of coords |
degree_out | Reference to Integer, order of polynomial fittings for vid |
coeffs | vector, coefficients of local polynomial fittings in monomial basis will be appended to the end of coeffs |
interp | Reference to Boolean, true = interpolation |
Definition at line 607 of file HiReconstruction.cpp.
613 {
614 if( !_hasfittings )
615 {
616 return false;
617 }
618 else
619 {
620 int index = _verts2rec.index( vid );
621 if( -1 == index )
622 {
623 std::cout << "Input vertex is not locally hosted vertex in this mesh set" << std::endl;
624 return false;
625 }
626 geomtype = _geom;
627 if( HISURFACE == _geom )
628 {
629 coords.insert( coords.end(), _local_coords.begin() + 9 * index, _local_coords.begin() + 9 * index + 9 );
630 degree_out = _degrees_out[index];
631 interp = _interps[index];
632 int ncoeffs = ( degree_out + 2 ) * ( degree_out + 1 ) >> 1;
633 size_t istr = _vertID2coeffID[index];
634 coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
635 }
636 else if( HI3DCURVE == _geom )
637 {
638 coords.insert( coords.end(), _local_coords.begin() + 3 * index, _local_coords.begin() + 3 * index + 3 );
639 degree_out = _degrees_out[index];
640 interp = _interps[index];
641 int ncoeffs = 3 * ( degree_out + 1 );
642 size_t istr = _vertID2coeffID[index];
643 coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
644 }
645 return true;
646 }
647 }
References _degrees_out, _geom, _hasfittings, _interps, _local_coords, _local_fit_coeffs, _vertID2coeffID, _verts2rec, moab::HI3DCURVE, moab::HISURFACE, and moab::Range::index().
|
protected |
Return the normals of given vertices in a Range, writing to preallocated memory If normals have been computed and stored, just access them If not, compute on the fly
vertsh | Range, EntityHandles of vertices |
nrms | Pointer of array of doubles, size = 3*vertsh.size() |
Definition at line 903 of file HiReconstruction.cpp.
904 {
905 ErrorCode error = MB_SUCCESS;
906 if( _hasderiv )
907 {
908 size_t id = 0;
909 for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
910 {
911 int index = _verts2rec.index( *ivert );
912 #ifdef MOAB_HAVE_MPI
913 if( -1 == index )
914 {
915 // ghost vertex
916 error = average_vertex_normal( *ivert, nrms + 3 * id );MB_CHK_ERR( error );
917 }
918 else
919 {
920 nrms[3 * id] = _local_coords[9 * index + 6];
921 nrms[3 * id + 1] = _local_coords[9 * index + 7];
922 nrms[3 * id + 2] = _local_coords[9 * index + 8];
923 }
924 #else
925 assert( -1 != index );
926 nrms[3 * id] = _local_coords[9 * index + 6];
927 nrms[3 * id + 1] = _local_coords[9 * index + 7];
928 nrms[3 * id + 2] = _local_coords[9 * index + 8];
929 #endif
930 }
931 }
932 else
933 {
934 size_t id = 0;
935 for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
936 {
937 error = average_vertex_normal( *ivert, nrms + 3 * id );MB_CHK_ERR( error );
938 }
939 }
940 return error;
941 }
References _hasderiv, _local_coords, _verts2rec, average_vertex_normal(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::Range::index(), MB_CHK_ERR, and MB_SUCCESS.
Referenced by polyfit3d_walf_surf_vertex().
|
protected |
Return the tangent vectors of given vertices in a Range, writing to preallocated memory If tangent vectors have been computed and stored, just access them If not, compute on the fly
vertsh | Range, EntityHandles of vertices |
tangs | Pointer of array of doubles, size = 3*vertsh.size() |
Definition at line 990 of file HiReconstruction.cpp.
991 {
992 ErrorCode error = MB_SUCCESS;
993 if( _hasderiv )
994 {
995 size_t id = 0;
996 for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
997 {
998 int index = _verts2rec.index( *ivert );
999 #ifdef MOAB_HAVE_MPI
1000 if( -1 != index )
1001 {
1002 tangs[3 * id] = _local_coords[3 * index];
1003 tangs[3 * id + 1] = _local_coords[3 * index + 1];
1004 tangs[3 * id + 2] = _local_coords[3 * index + 2];
1005 }
1006 else
1007 {
1008 error = average_vertex_tangent( *ivert, tangs + 3 * id );MB_CHK_ERR( error );
1009 }
1010 #else
1011 assert( -1 != index );
1012 tangs[3 * id] = _local_coords[3 * index];
1013 tangs[3 * id + 1] = _local_coords[3 * index + 1];
1014 tangs[3 * id + 2] = _local_coords[3 * index + 2];
1015 #endif
1016 }
1017 }
1018 else
1019 {
1020 size_t id = 0;
1021 for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
1022 {
1023 error = average_vertex_tangent( *ivert, tangs + 3 * id );MB_CHK_ERR( error );
1024 }
1025 }
1026 return error;
1027 }
References _hasderiv, _local_coords, _verts2rec, average_vertex_tangent(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::Range::index(), MB_CHK_ERR, and MB_SUCCESS.
Referenced by polyfit3d_walf_curve_vertex().
ErrorCode moab::HiReconstruction::hiproj_walf_around_vertex | ( | EntityHandle | vid, |
const int | npts2fit, | ||
const double * | coords2fit, | ||
double * | hiproj_new | ||
) |
Perform high order projection of points around a vertex, using estimated geometry by reconstruction class.
Given an vertex on the input mesh, and new points around this vertex, estimate their position in surface. This is done by first projecting input points onto the local uv-plane around this vertex and use the precomputed local fitting to estimate the ideal position of input points. The result will be returned to the user preallocated memory
vid | EntityHandle, the vertex around which to perform high order projection. |
npts2fit | Integer, number of points lying around vid to be fitted. |
coords2fit | Pointer to array of doubles, size=3*npts2fit, current coordinates of points to be projected. |
newcoords | Pointer to array of doubles, preallocated by user, size=3*npts2fit, estimated positions of input points. |
Definition at line 468 of file HiReconstruction.cpp.
472 {
473 if( !_hasfittings )
474 {
475 MB_SET_ERR( MB_FAILURE, "There is no existing fitting results" );
476 }
477 else if( -1 == _verts2rec.index( vid ) )
478 {
479 std::ostringstream convert;
480 convert << vid;
481 std::string VID = convert.str();
482 MB_SET_ERR( MB_FAILURE, "There is no existing fitting results for vertex " + VID );
483 }
484 ErrorCode error;
485 // get center of local coordinates system
486 double local_origin[3];
487 error = mbImpl->get_coords( &vid, 1, local_origin );MB_CHK_ERR( error );
488 // get local fitting parameters
489 int index = _verts2rec.index( vid );
490 bool interp = _interps[index];
491 int local_deg = _degrees_out[index];
492 double *uvw_coords, *local_coeffs;
493 if( _geom == HISURFACE )
494 {
495 uvw_coords = &( _local_coords[9 * index] );
496 // int ncoeffs = (local_deg+2)*(local_deg+1)>>1;
497 size_t istr = _vertID2coeffID[index];
498 local_coeffs = &( _local_fit_coeffs[istr] );
499 walf3d_surf_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
500 hiproj_new );
501 }
502 else if( _geom == HI3DCURVE )
503 {
504 uvw_coords = &( _local_coords[3 * index] );
505 size_t istr = _vertID2coeffID[index];
506 local_coeffs = &( _local_fit_coeffs[istr] );
507 walf3d_curve_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
508 hiproj_new );
509 }
510 return error;
511 }
References _degrees_out, _geom, _hasfittings, _interps, _local_coords, _local_fit_coeffs, _vertID2coeffID, _verts2rec, moab::error(), ErrorCode, moab::Core::get_coords(), moab::HI3DCURVE, moab::HISURFACE, moab::Range::index(), MB_CHK_ERR, MB_SET_ERR, mbImpl, walf3d_curve_vertex_eval(), and walf3d_surf_vertex_eval().
Referenced by hiproj_walf_in_element().
ErrorCode moab::HiReconstruction::hiproj_walf_in_element | ( | EntityHandle | elem, |
const int | nvpe, | ||
const int | npts2fit, | ||
const double * | naturalcoords2fit, | ||
double * | newcoords | ||
) |
Perform high order projection of points in an element, using estimated geometry by reconstruction class.
Given an element on the input mesh, and new points in this element, represented as natural coordinates in element, estimate their position in surface. This is done by weighted averaging of local fittings: for each vertex of this elment, a fitting has been computed and the new points could be projected by this fitting. The final result of projection is the weighted average of these projections, weights are chosen as the barycentric coordinates of the point in this element. The result will be returned to the user preallocated memory
elem | EntityHandle, the element on which to perform high order projection. |
nvpe | Integer, number of nodes of this element, triangle is 3, quad is four. |
npts2fit | Integer, number of points lying in elem to be projected. |
naturalcoords2fit | Pointer to array of doubles, size=nvpe*npts2fit, natural coordinates in elem of points to be projected. |
newcoords | Pointer to array of doubles, preallocated by user, size=3*npts2fit, estimated positions of input points. |
Definition at line 389 of file HiReconstruction.cpp.
394 {
395 assert( newcoords );
396 ErrorCode error;
397 // get connectivity table
398 std::vector< EntityHandle > elemconn;
399 error = mbImpl->get_connectivity( &elem, 1, elemconn );MB_CHK_ERR( error );
400 if( nvpe != (int)elemconn.size() )
401 {
402 MB_SET_ERR( MB_FAILURE, "element connectivity table size doesn't match input size" );
403 }
404
405 if( !_hasfittings )
406 {
407 MB_SET_ERR( MB_FAILURE, "There is no existing fitting results" );
408 }
409 else
410 {
411 std::ostringstream convert;
412 convert << elem;
413 std::string ID = convert.str();
414 for( int i = 0; i < nvpe; ++i )
415 {
416 if( -1 == _verts2rec.index( elemconn[i] ) )
417 {
418 MB_SET_ERR( MB_FAILURE, "There is no existing fitting results for element " + ID );
419 }
420 }
421 }
422 // check correctness of input
423 for( int i = 0; i < npts2fit; ++i )
424 {
425 if( !check_barycentric_coords( nvpe, naturalcoords2fit + i * nvpe ) )
426 {
427 MB_SET_ERR( MB_FAILURE, "Wrong barycentric coordinates" );
428 }
429 }
430
431 double* elemcoords = new double[nvpe * 3];
432 error = mbImpl->get_coords( &( elemconn[0] ), nvpe, elemcoords );MB_CHK_ERR( error );
433
434 double* coords2fit = new double[3 * npts2fit]();
435 for( int i = 0; i < npts2fit; ++i )
436 {
437 for( int j = 0; j < nvpe; ++j )
438 {
439 coords2fit[3 * i] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j];
440 coords2fit[3 * i + 1] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j + 1];
441 coords2fit[3 * i + 2] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j + 2];
442 }
443 }
444
445 double* hiproj_new = new double[3 * npts2fit];
446 // initialize output
447 for( int i = 0; i < npts2fit; ++i )
448 {
449 newcoords[3 * i] = newcoords[3 * i + 1] = newcoords[3 * i + 2] = 0;
450 }
451 // for each input vertex, call nvpe fittings and take average
452 for( int j = 0; j < nvpe; ++j )
453 {
454 error = hiproj_walf_around_vertex( elemconn[j], npts2fit, coords2fit, hiproj_new );MB_CHK_ERR( error );
455 for( int i = 0; i < npts2fit; ++i )
456 {
457 newcoords[3 * i] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i];
458 newcoords[3 * i + 1] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i + 1];
459 newcoords[3 * i + 2] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i + 2];
460 }
461 }
462 delete[] elemcoords;
463 delete[] coords2fit;
464 delete[] hiproj_new;
465 return error;
466 }
References _hasfittings, _verts2rec, check_barycentric_coords(), moab::error(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), hiproj_walf_around_vertex(), moab::Range::index(), MB_CHK_ERR, MB_SET_ERR, and mbImpl.
ErrorCode moab::HiReconstruction::initialize | ( | bool | recwhole | ) |
Definition at line 51 of file HiReconstruction.cpp.
52 {
53 ErrorCode error;
54
55 #ifdef HIREC_USE_AHF
56 std::cout << "HIREC_USE_AHF: Initializing" << std::endl;
57 ahf = new HalfFacetRep( mbImpl, pcomm, _mesh2rec, false );
58 if( !ahf )
59 {
60 return MB_MEMORY_ALLOCATION_FAILED;
61 }
62 error = ahf->initialize();MB_CHK_ERR( error );
63 #else
64 ahf = NULL;
65 #endif
66
67 // error = ahf->get_entity_ranges(_inverts,_inedges,_infaces,_incells); MB_CHK_ERR(error);
68 error = mbImpl->get_entities_by_dimension( _mesh2rec, 0, _inverts );MB_CHK_ERR( error );
69 error = mbImpl->get_entities_by_dimension( _mesh2rec, 1, _inedges );MB_CHK_ERR( error );
70 error = mbImpl->get_entities_by_dimension( _mesh2rec, 2, _infaces );MB_CHK_ERR( error );
71 error = mbImpl->get_entities_by_dimension( _mesh2rec, 3, _incells );MB_CHK_ERR( error );
72 if( _inedges.size() && _infaces.empty() && _incells.empty() )
73 {
74 _dim = 1;
75 _MAXPNTS = 13;
76 }
77 else if( _infaces.size() && _incells.empty() )
78 {
79 _dim = 2;
80 _MAXPNTS = 128;
81 }
82 else
83 {
84 MB_SET_ERR( MB_FAILURE, "Encountered a non-manifold mesh or a mesh with volume elements" );
85 }
86
87 // get locally hosted vertices by filtering pstatus
88 #ifdef MOAB_HAVE_MPI
89 if( pcomm )
90 {
91 error = pcomm->filter_pstatus( _inverts, PSTATUS_GHOST, PSTATUS_NOT, -1, &_verts2rec );MB_CHK_ERR( error );
92 }
93 else
94 {
95 _verts2rec = _inverts;
96 }
97 #else
98 _verts2rec = _inverts;
99 #endif
100 _nv2rec = _verts2rec.size();
101
102 if( recwhole )
103 {
104 // compute normals(surface) or tangent vector(curve) for all locally hosted vertices
105 if( 2 == _dim )
106 {
107 compute_average_vertex_normals_surf();
108 }
109 else if( 1 == _dim )
110 {
111 compute_average_vertex_tangents_curve();
112 }
113 else
114 {
115 MB_SET_ERR( MB_FAILURE, "Unknow space dimension" );
116 }
117 _hasderiv = true;
118 }
119 return error;
120 }
References _dim, _hasderiv, _incells, _inedges, _infaces, _inverts, _MAXPNTS, _mesh2rec, _nv2rec, _verts2rec, ahf, compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), moab::Range::empty(), moab::error(), ErrorCode, moab::ParallelComm::filter_pstatus(), moab::Core::get_entities_by_dimension(), moab::HalfFacetRep::initialize(), MB_CHK_ERR, MB_MEMORY_ALLOCATION_FAILED, MB_SET_ERR, mbImpl, pcomm, PSTATUS_GHOST, PSTATUS_NOT, and moab::Range::size().
Referenced by HiReconstruction().
|
protected |
Definition at line 780 of file HiReconstruction.cpp.
781 {
782 if( !_hasderiv )
783 {
784 compute_average_vertex_tangents_curve();
785 _hasderiv = true;
786 }
787 if( !_initfittings )
788 {
789 int ncoeffspvpd = degree + 1;
790 _degrees_out.assign( _nv2rec, 0 );
791 _interps.assign( _nv2rec, false );
792 _vertID2coeffID.resize( _nv2rec );
793 _local_fit_coeffs.assign( _nv2rec * ncoeffspvpd * 3, 0 );
794 for( size_t i = 0; i < _nv2rec; ++i )
795 {
796 _vertID2coeffID[i] = i * ncoeffspvpd * 3;
797 }
798 _initfittings = true;
799 }
800 }
References _degrees_out, _hasderiv, _initfittings, _interps, _local_fit_coeffs, _nv2rec, _vertID2coeffID, and compute_average_vertex_tangents_curve().
Referenced by reconstruct3D_curve_geom().
|
protected |
Definition at line 802 of file HiReconstruction.cpp.
803 {
804 if( !_hasderiv )
805 {
806 compute_average_vertex_tangents_curve();
807 _hasderiv = true;
808 }
809 if( !_hasfittings )
810 {
811 assert( _nv2rec == npts );
812 _degrees_out.assign( _nv2rec, 0 );
813 _interps.assign( _nv2rec, false );
814 _vertID2coeffID.reserve( _nv2rec );
815 size_t index = 0;
816 for( size_t i = 0; i < _nv2rec; ++i )
817 {
818 _vertID2coeffID[i] = index;
819 index += 3 * ( degrees[i] + 1 );
820 }
821 _local_fit_coeffs.assign( index, 0 );
822 _initfittings = true;
823 }
824 }
References _degrees_out, _hasderiv, _hasfittings, _initfittings, _interps, _local_fit_coeffs, _nv2rec, _vertID2coeffID, and compute_average_vertex_tangents_curve().
|
protected |
Initialize the storage for fitting results over _mesh2rec, curve/surface mesh Two options are provided: a) use uniform degree for all vertices b) use customized degrees for different vertices After calling of initializing functions, _initfitting is set to be true, the fitting result could be stored internally
Definition at line 734 of file HiReconstruction.cpp.
735 {
736 if( !_hasderiv )
737 {
738 compute_average_vertex_normals_surf();
739 _hasderiv = true;
740 }
741 if( !_initfittings )
742 {
743 int ncoeffspv = ( degree + 2 ) * ( degree + 1 ) / 2;
744 _degrees_out.assign( _nv2rec, 0 );
745 _interps.assign( _nv2rec, false );
746 _vertID2coeffID.resize( _nv2rec );
747 _local_fit_coeffs.assign( _nv2rec * ncoeffspv, 0 );
748 for( size_t i = 0; i < _nv2rec; ++i )
749 {
750 _vertID2coeffID[i] = i * ncoeffspv;
751 }
752 _initfittings = true;
753 }
754 }
References _degrees_out, _hasderiv, _initfittings, _interps, _local_fit_coeffs, _nv2rec, _vertID2coeffID, and compute_average_vertex_normals_surf().
Referenced by reconstruct3D_surf_geom().
|
protected |
Definition at line 756 of file HiReconstruction.cpp.
757 {
758 if( !_hasderiv )
759 {
760 compute_average_vertex_normals_surf();
761 _hasderiv = true;
762 }
763 if( !_initfittings )
764 {
765 assert( _nv2rec == npts );
766 _degrees_out.assign( _nv2rec, 0 );
767 _interps.assign( _nv2rec, false );
768 _vertID2coeffID.resize( _nv2rec );
769 size_t index = 0;
770 for( size_t i = 0; i < _nv2rec; ++i )
771 {
772 _vertID2coeffID[i] = index;
773 index += ( degrees[i] + 2 ) * ( degrees[i] + 1 ) / 2;
774 }
775 _local_fit_coeffs.assign( index, 0 );
776 _initfittings = true;
777 }
778 }
References _degrees_out, _hasderiv, _initfittings, _interps, _local_fit_coeffs, _nv2rec, _vertID2coeffID, and compute_average_vertex_normals_surf().
|
protected |
Get n-ring neighbor vertices, assuming curve/surface mesh, not volume mesh.
Given a vertex, find its n-ring neighbor vertices including itself in _mesrh2rec. 1-ring neighbor vertices of a vertex are the vertices connected with this vertex with an edge n-ring vertices are obtained first get the 1-ring vertices and then get the 1-ring of these vertices, and so on
vid | EntityHandle, vertex around which to get n-ring vertices |
ring | Integer, number of rings |
minpnts | Integer, number of minimum vertices to obtain, if the input ring could not provide enough vertices, i.e. more than minpnts, then expand the number of rings |
ngbvs | Range, the n-ring vertices of vid, including vid. If too many points found, i.e. more than _MAXPNTS, then terminate early. |
Definition at line 672 of file HiReconstruction.cpp.
673 {
674 ErrorCode error;
675 std::deque< EntityHandle > todo;
676 todo.push_back( vid );
677 ngbvs.insert( vid );
678 EntityHandle pre, nxt;
679 for( int i = 1; i <= ring; ++i )
680 {
681 int count = todo.size();
682 while( count )
683 {
684 EntityHandle center = todo.front();
685 todo.pop_front();
686 --count;
687 std::vector< EntityHandle > adjents;
688 error = vertex_get_incident_elements( center, _dim, adjents );MB_CHK_ERR( error );
689 for( size_t j = 0; j < adjents.size(); ++j )
690 {
691 std::vector< EntityHandle > elemconn;
692 error = mbImpl->get_connectivity( &adjents[j], 1, elemconn );MB_CHK_ERR( error );
693 int nvpe = elemconn.size();
694 for( int k = 0; k < nvpe; ++k )
695 {
696 if( elemconn[k] == center )
697 {
698 pre = k == 0 ? elemconn[nvpe - 1] : elemconn[k - 1];
699 nxt = elemconn[( k + 1 ) % nvpe];
700 if( ngbvs.find( pre ) == ngbvs.end() )
701 {
702 ngbvs.insert( pre );
703 todo.push_back( pre );
704 }
705 if( ngbvs.find( nxt ) == ngbvs.end() )
706 {
707 ngbvs.insert( nxt );
708 todo.push_back( nxt );
709 }
710 break;
711 }
712 }
713 }
714 }
715 if( _MAXPNTS <= (int)ngbvs.size() )
716 {
717 // obtain enough points
718 return error;
719 }
720 if( !todo.size() )
721 {
722 // current ring cannot introduce any points, return incase deadlock
723 return error;
724 }
725 if( ( i == ring ) && ( minpnts > (int)ngbvs.size() ) )
726 {
727 // reach maximum ring but not enough points
728 ++ring;
729 }
730 }
731 return error;
732 }
References _dim, _MAXPNTS, center(), moab::Range::end(), moab::error(), ErrorCode, moab::Range::find(), moab::Core::get_connectivity(), moab::Range::insert(), MB_CHK_ERR, mbImpl, moab::Range::size(), and vertex_get_incident_elements().
Referenced by polyfit3d_walf_curve_vertex(), and polyfit3d_walf_surf_vertex().
|
protected |
|
protected |
Compute local single variable coordinate system of a vertex, and perform vertex based polynomial fittings of three global coordinates axis.
This function take the first vertex of input as center, and build local u-line by estimating tangent vector Then other vertices form vectors originating from center and vectors then are projectd onto this u-plane to form three local height functions, one for each coordinates axis. Local fitting of these local height functions are performed in WLS sense, according if interpolation required or not.
nverts | Integer, number of vertices of input |
ngbcors | Pointer to array of doubles, size = 3*nverts, coordinates of input vertices, first will be center |
ngbtangs | Pointer to array of doubles, size = 3*nverts, vertex tangent vectors of input vertices |
degree | Integer, user specified fitting degree |
interp | Boolean, user input, interpolation or not |
safeguard | Boolean, true = use safeguarded numerical method in computing |
ncoords | Integer, size of *coords, should be 3, used for check |
coords | Pointer to array of doubles, preallocated memory for storing the glocal coordinates of local u axis direction |
ncoeffs | Integer, size of coeffs, should be 3(degree+1), used for check |
coeffs | Pointer to array of doubles, preallocated memory for storing coefficients of local fittings in monomial basis |
degree_out | Pointer to integer, order of resulted polynomial of fitting, could be downgraded due to numerical issues |
Definition at line 1318 of file HiReconstruction.cpp.
1329 {
1330 if( !nverts )
1331 {
1332 return;
1333 }
1334 // step 1. compute local coordinates system
1335 double tang[3] = { ngbtangs[0], ngbtangs[1], ngbtangs[2] };
1336 if( coords && ncoords > 2 )
1337 {
1338 coords[0] = tang[0];
1339 coords[1] = tang[1];
1340 coords[2] = tang[2];
1341 }
1342 if( !coeffs || !ncoeffs )
1343 {
1344 return;
1345 }
1346 else
1347 {
1348 assert( ncoeffs >= 3 * ( degree + 1 ) );
1349 }
1350 // step 2. project onto local coordinate system
1351 int npts2fit = nverts - interp;
1352 if( !npts2fit )
1353 {
1354 *degree_out = 0;
1355 for( int i = 0; i < ncoeffs; ++i )
1356 {
1357 coeffs[0] = 0;
1358 }
1359 return;
1360 }
1361 std::vector< double > us( npts2fit ); // double *us = new double[npts2fit];
1362 std::vector< double > bs( npts2fit * 3 ); // double *bs = new double[npts2fit*3];
1363 double uu[3];
1364 for( int i = interp; i < nverts; ++i )
1365 {
1366 int k = i - interp;
1367 DGMSolver::vec_linear_operation( 3, 1, ngbcors + 3 * i, -1, ngbcors, uu );
1368 us[k] = DGMSolver::vec_innerprod( 3, uu, tang );
1369 bs[k] = uu[0];
1370 bs[npts2fit + k] = uu[1];
1371 bs[2 * npts2fit + k] = uu[2];
1372 }
1373
1374 // step 3. copmute weights
1375 std::vector< double > ws( npts2fit );
1376 int nzeros = compute_weights( npts2fit, 1, &( us[0] ), nverts, ngbtangs, degree, _MINEPS, &( ws[0] ) );
1377 assert( nzeros <= npts2fit );
1378
1379 // step 4. adjust according to number of zero-weights
1380 if( nzeros )
1381 {
1382 if( nzeros == npts2fit )
1383 {
1384 // singular case
1385 *degree_out = 0;
1386 for( int i = 0; i < ncoeffs; ++i )
1387 {
1388 coeffs[i] = 0;
1389 }
1390 return;
1391 }
1392 int npts_new = npts2fit - nzeros;
1393 std::vector< double > bs_temp( 3 * npts_new );
1394 int index = 0;
1395 for( int i = 0; i < npts2fit; ++i )
1396 {
1397 if( ws[i] )
1398 {
1399 if( i > index )
1400 {
1401 us[index] = us[i];
1402 ws[index] = ws[i];
1403 }
1404 bs_temp[index] = bs[i];
1405 bs_temp[index + npts_new] = bs[i + npts2fit];
1406 bs_temp[index + 2 * npts_new] = bs[i + 2 * npts2fit];
1407 ++index;
1408 }
1409 }
1410 assert( index == npts_new );
1411 npts2fit = npts_new;
1412 us.resize( index );
1413 ws.resize( index );
1414 bs = bs_temp;
1415 // destroy bs_temp;
1416 std::vector< double >().swap( bs_temp );
1417 }
1418
1419 // step 5. fitting
1420 eval_vander_univar_cmf( npts2fit, &( us[0] ), 3, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out );
1421 // step 6. write results to output pointers
1422 int ncoeffs_out_pvpd = *degree_out + 1;
1423 assert( ncoeffs >= 3 * ncoeffs_out_pvpd );
1424 // write to coeffs consecutively, bs's size is not changed by eval_vander_univar_cmf
1425 coeffs[0] = coeffs[ncoeffs_out_pvpd] = coeffs[2 * ncoeffs_out_pvpd] = 0;
1426 for( int i = 0; i < ncoeffs_out_pvpd - interp; ++i )
1427 {
1428 coeffs[i + interp] = bs[i];
1429 coeffs[i + interp + ncoeffs_out_pvpd] = bs[i + npts2fit];
1430 coeffs[i + interp + 2 * ncoeffs_out_pvpd] = bs[i + 2 * npts2fit];
1431 }
1432 }
References _MINEPS, compute_weights(), eval_vander_univar_cmf(), moab::DGMSolver::vec_innerprod(), and moab::DGMSolver::vec_linear_operation().
Referenced by polyfit3d_walf_curve_vertex().
|
protected |
Compute local coordinates system of a vertex, and perform vertex based polynomial fittings of local height function.
This function take the first vertex of input as center, and build local uv-plane by estimating vertex normals and tangent planes Then other vertices forms vectors starting from center and then are projectd onto this uv-plane to form a local height function. Local fitting of this local height function is performed in WLS sense, according if interpolation required or not.
nverts | Integer, number of vertices of input |
ngbcors | Pointer to array of doubles, size = 3*nverts, coordinates of input vertices, first will be center |
ngbnrms | Pointer to array of doubles, size = 3*nverts, vertex normals of input vertices |
degree | Integer, user specified fitting degree |
interp | Boolean, user input, interpolation or not |
safeguard | Boolean, true = use safeguarded numerical method in computing |
ncoords | Integer, size of *coords, should be 9, used for check |
coords | Pointer to array of doubles, preallocated memory for storing the glocal coordinates of local uvw axis directions |
ncoeffs | Integer, size of *coeffs, should be (degree+2)(degree+1)/2, used for check |
coeffs | Pointer to array of doubles, preallocated memory for storing coefficients of local fittings in monomial basis |
degree_out | Pointer to integer, order of resulted polynomial of fitting, could be downgraded due to numerical issues |
degree_pnt | Pointer to integer, polynomial fitting order determined by stencil size/number of points |
degree_qr | Pointer to integer, polynomial fitting order determined by Vandermonde system condition number |
Definition at line 1033 of file HiReconstruction.cpp.
1046 {
1047 if( nverts <= 0 )
1048 {
1049 return;
1050 }
1051
1052 // std::cout << "npnts in initial stencil = " << nverts << std::endl;
1053 // std::cout << "centered at (" << ngbcoords[0] << "," << ngbcoords[1] << "," << ngbcoords[2] <<
1054 // ")" << std::endl;
1055
1056 // step 1. copmute local coordinate system
1057 double nrm[3] = { ngbnrms[0], ngbnrms[1], ngbnrms[2] }, tang1[3] = { 0, 0, 0 }, tang2[3] = { 0, 0, 0 };
1058 if( fabs( nrm[0] ) > fabs( nrm[1] ) && fabs( nrm[0] ) > fabs( nrm[2] ) )
1059 {
1060 tang1[1] = 1.0;
1061 }
1062 else
1063 {
1064 tang1[0] = 1.0;
1065 }
1066
1067 DGMSolver::vec_projoff( 3, tang1, nrm, tang1 );
1068 #ifndef NDEBUG
1069 assert( DGMSolver::vec_normalize( 3, tang1, tang1 ) );
1070 #endif
1071 DGMSolver::vec_crossprod( nrm, tang1, tang2 );
1072 if( 9 <= ncoords && coords )
1073 {
1074 coords[0] = tang1[0];
1075 coords[1] = tang1[1];
1076 coords[2] = tang1[2];
1077 coords[3] = tang2[0];
1078 coords[4] = tang2[1];
1079 coords[5] = tang2[2];
1080 coords[6] = nrm[0];
1081 coords[7] = nrm[1];
1082 coords[8] = nrm[2];
1083 }
1084 if( !ncoeffs || !coeffs )
1085 {
1086 return;
1087 }
1088 else
1089 {
1090 assert( ncoeffs >= ( degree + 2 ) * ( degree + 1 ) / 2 );
1091 }
1092
1093 // step 2. project onto local coordinates system
1094 int npts2fit = nverts - interp;
1095 if( 0 == npts2fit )
1096 {
1097 *degree_out = *degree_pnt = *degree_qr = 0;
1098 for( int i = 0; i < ncoeffs; ++i )
1099 {
1100 coeffs[i] = 0;
1101 }
1102 // coeffs[0] = 0;
1103 return;
1104 }
1105 std::vector< double > us( npts2fit * 2 ); // double *us = new double[npts2fit*2];
1106 std::vector< double > bs( npts2fit ); // double *bs = new double[npts2fit];
1107 for( int i = interp; i < nverts; ++i )
1108 {
1109 int k = i - interp;
1110 double uu[3];
1111 DGMSolver::vec_linear_operation( 3, 1, ngbcoords + 3 * i, -1, ngbcoords, uu );
1112 us[k * 2] = DGMSolver::vec_innerprod( 3, tang1, uu );
1113 us[k * 2 + 1] = DGMSolver::vec_innerprod( 3, tang2, uu );
1114 bs[k] = DGMSolver::vec_innerprod( 3, nrm, uu );
1115 }
1116
1117 // step 3. compute weights
1118 std::vector< double > ws( npts2fit ); // double *ws = new double[npts2fit];
1119 int nzeros = compute_weights( npts2fit, 2, &( us[0] ), nverts, ngbnrms, degree, _MINEPS, &( ws[0] ) );
1120
1121 // step 4. adjust according to zero-weights
1122 if( nzeros )
1123 {
1124 if( nzeros == npts2fit )
1125 {
1126 *degree_out = *degree_pnt = *degree_qr = 0;
1127 for( int i = 0; i < ncoeffs; ++i )
1128 {
1129 coeffs[i] = 0;
1130 }
1131 // coeffs[0] = 0;
1132 return;
1133 }
1134 int index = 0;
1135 for( int i = 0; i < npts2fit; ++i )
1136 {
1137 if( ws[i] )
1138 {
1139 if( i > index )
1140 {
1141 us[index * 2] = us[i * 2];
1142 us[index * 2 + 1] = us[i * 2 + 1];
1143 bs[index] = bs[i];
1144 ws[index] = ws[i];
1145 }
1146 ++index;
1147 }
1148 }
1149 npts2fit -= nzeros;
1150 assert( index == npts2fit );
1151 us.resize( npts2fit * 2 );
1152 bs.resize( npts2fit );
1153 ws.resize( npts2fit );
1154 /*us = realloc(us,npts2fit*2*sizeof(double));
1155 bs = realloc(bs,npts2fit*sizeof(double));
1156 ws = realloc(ws,npts2fit*sizeof(double));*/
1157 }
1158
1159 // std::cout<<"npnts after weighting = "<<npts2fit<<std::endl;
1160
1161 // step 5. fitting
1162 eval_vander_bivar_cmf( npts2fit, &( us[0] ), 1, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out,
1163 degree_pnt, degree_qr );
1164
1165 // step 6. organize output
1166 int ncoeffs_out = ( *degree_out + 2 ) * ( *degree_out + 1 ) / 2;
1167 assert( ncoeffs_out <= ncoeffs );
1168 coeffs[0] = 0;
1169 for( int j = 0; j < ncoeffs_out - interp; ++j )
1170 {
1171 coeffs[j + interp] = bs[j];
1172 }
1173 // delete [] us; delete [] bs; delete [] ws;
1174 }
References _MINEPS, compute_weights(), eval_vander_bivar_cmf(), moab::DGMSolver::vec_crossprod(), moab::DGMSolver::vec_innerprod(), moab::DGMSolver::vec_linear_operation(), moab::DGMSolver::vec_normalize(), and moab::DGMSolver::vec_projoff().
Referenced by polyfit3d_walf_surf_vertex().
ErrorCode moab::HiReconstruction::polyfit3d_walf_curve_vertex | ( | const EntityHandle | vid, |
const bool | interp, | ||
int | degree, | ||
int | minpnts, | ||
const bool | safeguard, | ||
const int | ncoords, | ||
double * | coords, | ||
int * | degree_out, | ||
const int | ncoeffs, | ||
double * | coeffs | ||
) |
Construct vertex based polynomial fitting on a curve mesh.
Given a vertex on a curve mesh, construct three one-parameter local fittings for each coordinates axis around this vertex. Stencils around this vertex will be selected according to input degree and if data is noise. Local u-line, or the single parameter will be the estimated tangent line at this vertex. On each axis of xyz, a polynomial fitting will be performed according to user input. minpnts will be used to specify the minimum number allowed in the local stencil. The result will be returned to user by preallocated memory coords, degree_out, coeffs.
vid | EntityHandle, the fittings will be performed around this vertex. |
interp | Boolean, true=Interpolation, false=least square fitting. |
degree | Integer, order of polynomials used for local fittings. |
minpnts | Integer, the allowed minimum number of vertices in local stencil. If too small, the resulted fitting might be low order accurate. If too large, it may introduce overfitting. |
safeguard | Boolean, true=using safe guarded method in numerical computing. |
coords | Pointer to double, preallocated memory by user, should have at least 3 doubles; stores the global coordinates of local coordinate system u direction. |
degree_out | Pointer to integer, used to store the degree of resulted fitting |
coeffs,Pointer | to double, preallocated memory for coefficients of local fittings, should have at least 3*(degree+1) doubles. |
Definition at line 344 of file HiReconstruction.cpp.
354 {
355 ErrorCode error;
356 int ring = estimate_num_rings( degree, interp );
357 // get n-ring neighbors
358 Range ngbvs;
359 error = obtain_nring_ngbvs( vid, ring, minpnts, ngbvs );MB_CHK_ERR( error );
360 // get coordinates
361 size_t nverts = ngbvs.size();
362 assert( nverts );
363 double* ngbcoords = new double[nverts * 3];
364 error = mbImpl->get_coords( ngbvs, ngbcoords );MB_CHK_ERR( error );
365 // get tangent vectors
366 double* ngbtangs = new double[nverts * 3];
367 error = get_tangents_curve( ngbvs, ngbtangs );MB_CHK_ERR( error );
368 // switch vid to first one
369 int index = ngbvs.index( vid );
370 assert( index != -1 );
371 std::swap( ngbcoords[0], ngbcoords[3 * index] );
372 std::swap( ngbcoords[1], ngbcoords[3 * index + 1] );
373 std::swap( ngbcoords[2], ngbcoords[3 * index + 2] );
374 std::swap( ngbtangs[0], ngbtangs[3 * index] );
375 std::swap( ngbtangs[1], ngbtangs[3 * index + 1] );
376 std::swap( ngbtangs[2], ngbtangs[3 * index + 2] );
377 // local WLS fittings
378 polyfit3d_curve_get_coeff( nverts, ngbcoords, ngbtangs, degree, interp, safeguard, ncoords, coords, ncoeffs, coeffs,
379 degree_out );
380 delete[] ngbcoords;
381 delete[] ngbtangs;
382 return error;
383 }
References moab::error(), ErrorCode, estimate_num_rings(), moab::Core::get_coords(), get_tangents_curve(), moab::Range::index(), MB_CHK_ERR, mbImpl, obtain_nring_ngbvs(), polyfit3d_curve_get_coeff(), and moab::Range::size().
Referenced by reconstruct3D_curve_geom().
ErrorCode moab::HiReconstruction::polyfit3d_walf_surf_vertex | ( | const EntityHandle | vid, |
const bool | interp, | ||
int | degree, | ||
int | minpnts, | ||
const bool | safeguard, | ||
const int | ncoords, | ||
double * | coords, | ||
int * | degree_out, | ||
const int | ncoeffs, | ||
double * | coeffs | ||
) |
Construct vertex based polynomial fitting on a surface mesh.
Given a vertex on a surface mesh, construct a local fitting around this vertex. Stencils around this vertex will be selected according to input degree and if data is noise. Local uv-plane will be the estimated tangent plane at this vertex. minpnts will be used to specify the minimum number allowed in the local stencil. The result will be returned to user by preallocated memory coords, degree_out, coeffs.
vid | EntityHandle, the fitting will be performed around this vertex for the local height function over the uv-plane. |
interp | Boolean, true=Interpolation, false=least square fitting. |
degree | Integer, order of polynomials used for local fittings. |
minpnts | Integer, the allowed minimum number of vertices in local stencil. If too small, the resulted fitting might be low order accurate. If too large, it may introduce overfitting. |
safeguard | Boolean, true=using safe guarded method in numerical computing. |
coords | Pointer to double, preallocated memory by user, should have at least 9 doubles; stores the global coordinates of local coordinates system uvw directions. |
degree_out | Pointer to integer, used to store the degree of resulted fitting |
coeffs,Pointer | to double, preallocated memory for coefficients of local fittings, should have at least (degree+2)(degree+1)/2 doubles. |
Definition at line 294 of file HiReconstruction.cpp.
304 {
305 assert( _dim == 2 );
306 ErrorCode error;
307 int ring = estimate_num_rings( degree, interp );
308 // std::cout<<"ring = "<<ring<<std::endl;
309 // get n-ring neighbors
310 Range ngbvs;
311 error = obtain_nring_ngbvs( vid, ring, minpnts, ngbvs );MB_CHK_ERR( error );
312 // for debug
313 /*if(_verts2rec.index(vid)==70){
314 for(Range::iterator ingb=ngbvs.begin();ingb!=ngbvs.end();++ingb) std::cerr <<
315 _verts2rec.index(*ingb) << " "; std::cout << std::endl;
316 }*/
317
318 // get coordinates;
319 size_t nverts = ngbvs.size();
320 assert( nverts );
321 double* ngbcoords = new double[nverts * 3];
322 error = mbImpl->get_coords( ngbvs, ngbcoords );MB_CHK_ERR( error );
323 // get normals
324 double* ngbnrms = new double[nverts * 3];
325 error = get_normals_surf( ngbvs, ngbnrms );MB_CHK_ERR( error );
326 // switch vid to first one
327 int index = ngbvs.index( vid );
328 assert( index != -1 );
329 std::swap( ngbcoords[0], ngbcoords[3 * index] );
330 std::swap( ngbcoords[1], ngbcoords[3 * index + 1] );
331 std::swap( ngbcoords[2], ngbcoords[3 * index + 2] );
332 std::swap( ngbnrms[0], ngbnrms[3 * index] );
333 std::swap( ngbnrms[1], ngbnrms[3 * index + 1] );
334 std::swap( ngbnrms[2], ngbnrms[3 * index + 2] );
335 // local WLS fitting
336 int degree_pnt, degree_qr;
337 polyfit3d_surf_get_coeff( nverts, ngbcoords, ngbnrms, degree, interp, safeguard, ncoords, coords, ncoeffs, coeffs,
338 degree_out, °ree_pnt, °ree_qr );
339 delete[] ngbcoords;
340 delete[] ngbnrms;
341 return error;
342 }
References _dim, moab::error(), ErrorCode, estimate_num_rings(), moab::Core::get_coords(), get_normals_surf(), moab::Range::index(), MB_CHK_ERR, mbImpl, obtain_nring_ngbvs(), polyfit3d_surf_get_coeff(), and moab::Range::size().
Referenced by reconstruct3D_surf_geom().
ErrorCode moab::HiReconstruction::reconstruct3D_curve_geom | ( | int | degree, |
bool | interp, | ||
bool | safeguard, | ||
bool | reset = false |
||
) |
Reconstruct a high order curve on given curve mesh.
Given a curve mesh, compute vertex based polynomail fittings for all vertices hosted by current processor. The vertex based fitting is done by perfoming three one-parameter fittings along each axis, i.e. x,y,z. The result will be stored interally for later usage of evalution.
degree | Integer, order of polynomials used for local fittings. |
interp | Boolean, true=Interpolation, false=least square fitting. |
safeguard | Boolean, true=using safe guarded method in numerical computing. |
reset | Boolean, reset=true will recompute the fittings based on user input and replace the existing one. |
Definition at line 221 of file HiReconstruction.cpp.
222 {
223 assert( _dim == 1 );
224 if( _hasfittings && !reset )
225 {
226 return MB_SUCCESS;
227 }
228 else
229 {
230 _initfittings = _hasfittings = false;
231 }
232 initialize_3Dcurve_geom( degree );
233 ErrorCode error;
234 double *coords = 0, *coeffs;
235 int* degree_out;
236 int ncoeffs = 3 * ( degree + 1 );
237 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert )
238 {
239 int index = _verts2rec.index( *ivert );
240 assert( index != -1 );
241 size_t istr = _vertID2coeffID[index];
242 coeffs = &( _local_fit_coeffs[istr] );
243 degree_out = &( _degrees_out[index] );
244 _interps[index] = interp;
245 error = polyfit3d_walf_curve_vertex( *ivert, interp, degree, _MINPNTS, safeguard, 0, coords, degree_out,
246 ncoeffs, coeffs );MB_CHK_ERR( error );
247 }
248 _geom = HI3DCURVE;
249 _hasfittings = true;
250 return error;
251 }
References _degrees_out, _dim, _geom, _hasfittings, _initfittings, _interps, _local_fit_coeffs, _MINPNTS, _vertID2coeffID, _verts2rec, moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::HI3DCURVE, moab::Range::index(), initialize_3Dcurve_geom(), MB_CHK_ERR, MB_SUCCESS, and polyfit3d_walf_curve_vertex().
ErrorCode moab::HiReconstruction::reconstruct3D_curve_geom | ( | size_t | npts, |
int * | degrees, | ||
bool * | interps, | ||
bool | safeguard, | ||
bool | reset = false |
||
) |
Reconstruct a high order curve on given curve mesh.
Given a curve mesh, compute vertex based polynomail fittings for all vertices hosted by current processor. The vertex based fitting is done by perfoming three one-parameter fittings along each axis, i.e. x,y,z. User could specify various degrees for different vertices. It assumes that the input degrees for vertices stored in the same order as that this class stores vertices: 1) reconstruction will be only performed at vertices hosted by current processor, thus input npts should match the number of hosted vertices. 2) all hosted vertices will be stored in a MOAB::Range object, degrees for all these vertices should be stored in degrees as the same order in the MOAB::Range object The result will be stored interally for later usage of evalution.
npts | Integer size of array pointed by degrees, used for check |
degrees | Integer arrray, order of polynomials for local fitting at all hosted vertices. |
interp | Boolean, true=Interpolation, false=least square fitting. |
safeguard | Boolean, true=using safe guarded method in numerical computing. |
reset | Boolean, reset=true will recompute the fittings based on user input and replace the existing one. |
Definition at line 253 of file HiReconstruction.cpp.
258 {
259 assert( _dim == 1 );
260 ErrorCode error;
261 if( npts != _nv2rec )
262 {
263 MB_SET_ERR( MB_FAILURE, "Input number of degrees doesn't match the number of vertices" );
264 }
265 if( _hasfittings && !reset )
266 {
267 return MB_SUCCESS;
268 }
269 else
270 {
271 _initfittings = _hasfittings = false;
272 }
273 // initialize
274 initialize_3Dcurve_geom( npts, degrees );
275 double *coords = 0, *coeffs;
276 int* degree_out;
277 size_t i = 0;
278 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++i )
279 {
280 int index = _verts2rec.index( *ivert );
281 size_t istr = _vertID2coeffID[index];
282 coeffs = &( _local_fit_coeffs[istr] );
283 degree_out = &( _degrees_out[index] );
284 _interps[index] = interps[i];
285 int ncoeffs = 3 * ( degrees[i] + 1 );
286 error = polyfit3d_walf_curve_vertex( *ivert, interps[i], degrees[i], _MINPNTS, safeguard, 0, coords, degree_out,
287 ncoeffs, coeffs );MB_CHK_ERR( error );
288 }
289 _geom = HI3DCURVE;
290 _hasfittings = true;
291 return error;
292 }
References _degrees_out, _dim, _geom, _hasfittings, _initfittings, _interps, _local_fit_coeffs, _MINPNTS, _nv2rec, _vertID2coeffID, _verts2rec, moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::HI3DCURVE, moab::Range::index(), initialize_3Dcurve_geom(), MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, and polyfit3d_walf_curve_vertex().
ErrorCode moab::HiReconstruction::reconstruct3D_surf_geom | ( | int | degree, |
bool | interp, | ||
bool | safeguard, | ||
bool | reset = false |
||
) |
Reconstruct a high order surface on given surface mesh.
Given a mesh, compute vertex based polynomial fittings for all vertices hosted by current processor. The result will be stored interally for later usage of evalution. The inputs are: a) degree, which is the order of polynomial used for vertex based fitting. b) interp, if it's true, then interpolation will be applied for local fitting, otherwise it's least square fitting. c) safeguard, specifies whether to use safeguarded numeric method. d) reset, if fittings have been computed and stored in current object, then reset=true will recompute the fittings based on user input and replace the existing one.
degree | Integer, order of polynomials used for local fittings. |
interp | Boolean, true=Interpolation, false=least square fitting. |
safeguard | Boolean, true=using safe guarded method in numerical computing. |
reset | Boolean, reset=true will recompute the fittings based on user input and replace the existing one. |
Definition at line 126 of file HiReconstruction.cpp.
127 {
128 assert( 2 == _dim );
129 if( _hasfittings && !reset )
130 {
131 // This object has precomputed fitting results and user don't want to reset
132 return MB_SUCCESS;
133 }
134 else
135 {
136 _initfittings = _hasfittings = false;
137 }
138 // initialize for geometric information
139 initialize_surf_geom( degree );
140 ErrorCode error;
141 double *coeffs, *coords;
142 int* degree_out;
143 int ncoeffs = ( degree + 2 ) * ( degree + 1 ) / 2;
144
145 // DBG
146 // int dcount = 0;
147
148 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert )
149 {
150 int index = _verts2rec.index( *ivert );
151 // for debug
152 /*if(index==70){
153 EntityHandle vid = *ivert;
154 double vertcoords[3];
155 error = mbImpl->get_coords(&vid,1,vertcoords);
156 }*/
157
158 size_t istr = _vertID2coeffID[index];
159 coords = &( _local_coords[9 * index] );
160 coeffs = &( _local_fit_coeffs[istr] );
161 degree_out = &( _degrees_out[index] );
162 _interps[index] = interp;
163 error = polyfit3d_walf_surf_vertex( *ivert, interp, degree, _MINPNTS, safeguard, 9, coords, degree_out, ncoeffs,
164 coeffs );MB_CHK_ERR( error );
165
166 // DBG
167 // if( degree_out[0] < degree ) dcount += 1;
168 }
169
170 // DBG
171 // std::cout<<"Total #points ="<<_verts2rec.size()<<", #degraded points = "<<dcount<<std::endl;
172
173 _geom = HISURFACE;
174 _hasfittings = true;
175 return error;
176 }
References _degrees_out, _dim, _geom, _hasfittings, _initfittings, _interps, _local_coords, _local_fit_coeffs, _MINPNTS, _vertID2coeffID, _verts2rec, moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::HISURFACE, moab::Range::index(), initialize_surf_geom(), MB_CHK_ERR, MB_SUCCESS, and polyfit3d_walf_surf_vertex().
ErrorCode moab::HiReconstruction::reconstruct3D_surf_geom | ( | size_t | npts, |
int * | degrees, | ||
bool * | interps, | ||
bool | safeguard, | ||
bool | reset = false |
||
) |
Reconstruct a high order surface on given surface mesh.
Given a mesh, compute vertex based polynomial fittings for all vertices hosted by current processor. User could specify various degrees for different vertices. It assumes that the input degrees for vertices stored in the same order as that this class stores vertices: 1) reconstruction will be only performed at vertices hosted by current processor, thus input npts should match the number of hosted vertices. 2) all hosted vertices will be stored in a MOAB::Range object, degrees for all these vertices should be stored in degrees as the same order in the MOAB::Range object The result will be stored interally for later usage of evalution.
npts | Integer size of array pointed by degrees, used for check |
degrees | Integer arrray, order of polynomials for local fitting at all hosted vertices |
interp | Boolean, true=Interpolation, false=least square fitting. |
safeguard | Boolean, true=using safe guarded method in numerical computing. |
reset | Boolean, reset=true will recompute the fittings based on user input and replace the existing one. |
Definition at line 178 of file HiReconstruction.cpp.
183 {
184 assert( _dim == 2 );
185 if( npts != _nv2rec )
186 {
187 MB_SET_ERR( MB_FAILURE, "Input number of degrees doesn't match number of vertices" );
188 }
189 if( _hasfittings && !reset )
190 {
191 return MB_SUCCESS;
192 }
193 else
194 {
195 _initfittings = _hasfittings = false;
196 }
197 ErrorCode error;
198 // initialization for fitting
199 initialize_surf_geom( npts, degrees );
200 double *coeffs, *coords;
201 int* degree_out;
202 size_t i = 0;
203 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++i )
204 {
205 int index = _verts2rec.index( *ivert );
206 assert( -1 != index );
207 size_t istr = _vertID2coeffID[index];
208 coords = &( _local_coords[9 * index] );
209 coeffs = &( _local_fit_coeffs[istr] );
210 degree_out = &( _degrees_out[index] );
211 _interps[index] = interps[i];
212 int ncoeffs = ( degrees[i] + 2 ) * ( degrees[i] + 1 ) / 2;
213 error = polyfit3d_walf_surf_vertex( *ivert, interps[i], degrees[i], _MINPNTS, safeguard, 9, coords, degree_out,
214 ncoeffs, coeffs );MB_CHK_ERR( error );
215 }
216 _geom = HISURFACE;
217 _hasfittings = true;
218 return error;
219 }
References _degrees_out, _dim, _geom, _hasfittings, _initfittings, _interps, _local_coords, _local_fit_coeffs, _MINPNTS, _nv2rec, _vertID2coeffID, _verts2rec, moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::HISURFACE, moab::Range::index(), initialize_surf_geom(), MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, and polyfit3d_walf_surf_vertex().
|
protected |
Given a vertex, return the incident elements with dimension elemdim.
Wrapper of MOAB Core->get_adjacencies and HalfRep->get_up_adjacencies, depends on if USE_AHF is defined
vid | EntityHandle of vertex |
elemdim | Integer, dimension of elements incidented in vid |
adjents | vector<EntityHandle>, container which push incident elements in |
Definition at line 658 of file HiReconstruction.cpp.
661 {
662 ErrorCode error;
663 assert( elemdim == _dim );
664 #ifdef HIREC_USE_AHF
665 error = ahf->get_up_adjacencies( vid, elemdim, adjents );MB_CHK_ERR( error );
666 #else
667 error = mbImpl->get_adjacencies( &vid, 1, elemdim, false, adjents );MB_CHK_ERR( error );
668 #endif
669 return error;
670 }
References _dim, ahf, moab::error(), ErrorCode, moab::Core::get_adjacencies(), moab::HalfFacetRep::get_up_adjacencies(), MB_CHK_ERR, and mbImpl.
Referenced by average_vertex_normal(), average_vertex_tangent(), and obtain_nring_ngbvs().
void moab::HiReconstruction::walf3d_curve_vertex_eval | ( | const double * | local_origin, |
const double * | local_coords, | ||
const int | local_deg, | ||
const double * | local_coeffs, | ||
const bool | interp, | ||
const int | npts2fit, | ||
const double * | coords2fit, | ||
double * | hiproj_new | ||
) |
Perform high order projection of points around a center vertex, assume geometry is curve.
Given a vertex position and the local one-parameter fittings parameter around this vertex, estimate the ideal position of input position according to the local fittings. This is done by first projecting input points onto the local u-direction at this vertex and then use the value u as parameter for the three fittings, one for each coordinates axis of xyz. The result will be returned to user preallocated memory
local_origin | Pointer to 3 doubles, coordinates of the center vertex |
local_coords | Pointer to 3 doubles, global coordinates of direction of local u coordinate axis at center vertex |
local_deg | Integer, order of local polynomial fitting |
local_coeffs | Pointer to array of doubles, size=3*(local_deg+1), coefficients of three local polynomial fittings, in monomial basis. For each fitting, local_deg+1 parameters. |
interp | Boolean, true=local fitting is interpolation, false=local fitting is least square fitting |
npts2fit | Integer, number of points to be estimated, around the center vertices |
coords2fit | Pointer to array of doubles, size=3*npts2fit, current coordinates of points to be estimated |
hiproj_new | Pointer to array of doubles, size=3*npts2fit, memory preallocated by user to store the fitting/estimated positions of input points. |
Definition at line 569 of file HiReconstruction.cpp.
577 {
578 assert( local_origin && local_coords && local_coeffs );
579 int ncoeffspvpd = local_deg + 1;
580 for( int i = 0; i < npts2fit; ++i )
581 {
582 // get the vector from center to current point, project to tangent line
583 double vec[3], ans[3] = { 0, 0, 0 };
584 DGMSolver::vec_linear_operation( 3, 1, coords2fit + 3 * i, -1, local_origin, vec );
585 double u = DGMSolver::vec_innerprod( 3, local_coords, vec );
586 // evaluate polynomials
587 if( !interp )
588 {
589 ans[0] = local_coeffs[0];
590 ans[1] = local_coeffs[ncoeffspvpd];
591 ans[2] = local_coeffs[2 * ncoeffspvpd];
592 }
593 double uk = 1; // degree_out and degree different, stored in columnwise contiguously
594 for( int j = 1; j < ncoeffspvpd; ++j )
595 {
596 uk *= u;
597 ans[0] += uk * local_coeffs[j];
598 ans[1] += uk * local_coeffs[j + ncoeffspvpd];
599 ans[2] += uk * local_coeffs[j + 2 * ncoeffspvpd];
600 }
601 hiproj_new[3 * i] = ans[0] + local_origin[0];
602 hiproj_new[3 * i + 1] = ans[1] + local_origin[1];
603 hiproj_new[3 * i + 2] = ans[2] + local_origin[2];
604 }
605 }
References moab::DGMSolver::vec_innerprod(), and moab::DGMSolver::vec_linear_operation().
Referenced by hiproj_walf_around_vertex().
void moab::HiReconstruction::walf3d_surf_vertex_eval | ( | const double * | local_origin, |
const double * | local_coords, | ||
const int | local_deg, | ||
const double * | local_coeffs, | ||
const bool | interp, | ||
const int | npts2fit, | ||
const double * | coords2fit, | ||
double * | hiproj_new | ||
) |
Perform high order projection of points around a center vertex, assume geometry is surface.
Given a vertex position and the local fitting parameter around this vertex, estimate the ideal position of input position according to the local fitting. This is done by first projecting input points onto the local uv-plane around this vertex and use the given fitting to estimate the ideal position of input points. The result will be returned to user preallocated memory
local_origin | Pointer to 3 doubles, coordinates of the center vertex |
local_coords | Pointer to 9 doubles, global coordinates of directions of local uvw coordinates axis at center vertex |
local_deg | Integer, order of local polynomial fitting |
local_coeffs | Pointer to array of doubles, size=(local_deg+2)(local_deg+1)/2, coefficients of local polynomial fittings, in monomial basis |
interp | Boolean, true=local fitting is interpolation, false=local fitting is least square fitting |
npts2fit | Integer, number of points to be estimated, around the center vertices |
coords2fit | Pointer to array of doubles, size=3*npts2fit, current coordinates of points to be estimated |
hiproj_new | Pointer to array of doubles, size=3*npts2fit, memory preallocated by user to store the fitting/estimated positions of input points. |
Definition at line 513 of file HiReconstruction.cpp.
521 {
522 double xaxis[3], yaxis[3], zaxis[3];
523 for( int i = 0; i < 3; ++i )
524 {
525 xaxis[i] = local_coords[i];
526 yaxis[i] = local_coords[3 + i];
527 zaxis[i] = local_coords[6 + i];
528 }
529 // double *basis = new double[(local_deg+2)*(local_deg+1)/2-1];
530 std::vector< double > basis( ( local_deg + 2 ) * ( local_deg + 1 ) / 2 - 1 );
531 for( int i = 0; i < npts2fit; ++i )
532 {
533 double local_pos[3];
534 for( int j = 0; j < 3; ++j )
535 {
536 local_pos[j] = coords2fit[3 * i + j] - local_origin[j];
537 }
538 double u, v, height = 0;
539 u = DGMSolver::vec_innerprod( 3, local_pos, xaxis );
540 v = DGMSolver::vec_innerprod( 3, local_pos, yaxis );
541 basis[0] = u;
542 basis[1] = v;
543 int l = 1;
544 for( int k = 2; k <= local_deg; ++k )
545 {
546 ++l;
547 basis[l] = u * basis[l - k];
548 for( int id = 0; id < k; ++id )
549 {
550 ++l;
551 basis[l] = basis[l - k - 1] * v;
552 }
553 }
554 if( !interp )
555 {
556 height = local_coeffs[0];
557 }
558 for( int p = 0; p <= l; ++p )
559 {
560 height += local_coeffs[p + 1] * basis[p];
561 }
562 hiproj_new[3 * i] = local_origin[0] + u * xaxis[0] + v * yaxis[0] + height * zaxis[0];
563 hiproj_new[3 * i + 1] = local_origin[1] + u * xaxis[1] + v * yaxis[1] + height * zaxis[1];
564 hiproj_new[3 * i + 2] = local_origin[2] + u * xaxis[2] + v * yaxis[2] + height * zaxis[2];
565 }
566 // delete [] basis;
567 }
References moab::DGMSolver::vec_innerprod().
Referenced by hiproj_walf_around_vertex().
|
protected |
Definition at line 312 of file HiReconstruction.hpp.
Referenced by get_fittings_data(), hiproj_walf_around_vertex(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 306 of file HiReconstruction.hpp.
Referenced by HiReconstruction(), initialize(), obtain_nring_ngbvs(), polyfit3d_walf_surf_vertex(), reconstruct3D_curve_geom(), reconstruct3D_surf_geom(), and vertex_get_incident_elements().
|
protected |
Definition at line 305 of file HiReconstruction.hpp.
Referenced by get_fittings_data(), hiproj_walf_around_vertex(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 303 of file HiReconstruction.hpp.
Referenced by compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), get_normals_surf(), get_tangents_curve(), HiReconstruction(), initialize(), initialize_3Dcurve_geom(), and initialize_surf_geom().
|
protected |
Definition at line 307 of file HiReconstruction.hpp.
Referenced by get_fittings_data(), hiproj_walf_around_vertex(), hiproj_walf_in_element(), HiReconstruction(), initialize_3Dcurve_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 295 of file HiReconstruction.hpp.
Referenced by initialize().
|
protected |
Definition at line 295 of file HiReconstruction.hpp.
Referenced by initialize().
|
protected |
Definition at line 295 of file HiReconstruction.hpp.
Referenced by initialize().
|
protected |
Definition at line 308 of file HiReconstruction.hpp.
Referenced by initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 313 of file HiReconstruction.hpp.
Referenced by get_fittings_data(), hiproj_walf_around_vertex(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 295 of file HiReconstruction.hpp.
Referenced by initialize().
|
protected |
Definition at line 309 of file HiReconstruction.hpp.
Referenced by compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), get_fittings_data(), get_normals_surf(), get_tangents_curve(), hiproj_walf_around_vertex(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 310 of file HiReconstruction.hpp.
Referenced by get_fittings_data(), hiproj_walf_around_vertex(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 298 of file HiReconstruction.hpp.
Referenced by compute_weights(), initialize(), and obtain_nring_ngbvs().
|
protected |
Definition at line 292 of file HiReconstruction.hpp.
Referenced by initialize().
|
protected |
Definition at line 299 of file HiReconstruction.hpp.
Referenced by check_barycentric_coords(), HiReconstruction(), polyfit3d_curve_get_coeff(), and polyfit3d_surf_get_coeff().
|
protected |
Definition at line 298 of file HiReconstruction.hpp.
Referenced by reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 296 of file HiReconstruction.hpp.
Referenced by compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), initialize(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 311 of file HiReconstruction.hpp.
Referenced by get_fittings_data(), hiproj_walf_around_vertex(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 295 of file HiReconstruction.hpp.
Referenced by compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), get_fittings_data(), get_normals_surf(), get_tangents_curve(), hiproj_walf_around_vertex(), hiproj_walf_in_element(), initialize(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().
|
protected |
Definition at line 286 of file HiReconstruction.hpp.
Referenced by initialize(), vertex_get_incident_elements(), and ~HiReconstruction().
|
protected |
Definition at line 284 of file HiReconstruction.hpp.
Referenced by average_vertex_normal(), average_vertex_tangent(), hiproj_walf_around_vertex(), hiproj_walf_in_element(), HiReconstruction(), initialize(), obtain_nring_ngbvs(), polyfit3d_walf_curve_vertex(), polyfit3d_walf_surf_vertex(), and vertex_get_incident_elements().
|
protected |
Definition at line 285 of file HiReconstruction.hpp.
Referenced by HiReconstruction(), and initialize().