diff --git a/src/main.cpp b/src/main.cpp index b99c734..eec315d 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -40,9 +40,9 @@ Scalar pointRadius = 0.005; /// < display radius of the point cloud bool useKnnGraph = false; /// < use k-neighbor graph instead of kdtree // Slicer -int slice; +float slice; int axis; - +bool isHDSlicer=false; /// Convenience function measuring and printing the processing time of F template @@ -230,6 +230,39 @@ void mlsDryRun() { }); } +Scalar eval_scalar_field_ASO(const VectorType& input_pos) +{ + using FitASO = Ponca::BasketDiff< + Ponca::Basket< + PPAdapter, + SmoothWeightFunc, + Ponca::OrientedSphereFit>, + Ponca::DiffType::FitSpaceDer, + Ponca::OrientedSphereDer, + Ponca::MlsSphereFitDer, + Ponca::CurvatureEstimatorBase, + Ponca::NormalDerivativesCurvatureEstimator>; + VectorType current_pos = input_pos; + Scalar current_value = std::numeric_limits::max(); + for(int mm = 0; mm < mlsIter; ++mm) + { + FitASO fit; + fit.setWeightFunc(SmoothWeightFunc(NSize)); + fit.init(current_pos); // weighting function using current pos (not input pos) + for(int j : tree.range_neighbors(current_pos, NSize)) { + fit.addNeighbor(tree.point_data()[j]); + } + if(fit.finalize() == Ponca::STABLE) { + current_pos = fit.project(input_pos); // always project input pos + current_value = fit.potential(input_pos); + // current_gradient = fit.primitiveGradient(input_pos); + } else { + // not enough neighbors (if far from the point cloud) + return .0;//std::numeric_limits::max(); + } + } + return current_value; +} /// Define Polyscope callbacks void callback() { @@ -261,52 +294,16 @@ void callback() { if (ImGui::Button("ASO")) estimateDifferentialQuantitiesWithASO(); ImGui::Separator(); - ImGui::SliderInt("Slice / 1024", &slice, 0, 1024); + ImGui::SliderFloat("Slice", &slice, 0, 1.0); ImGui::SameLine(); + ImGui::Checkbox("HD", &isHDSlicer); ImGui::RadioButton("X axis", &axis, 0); ImGui::SameLine(); ImGui::RadioButton("Y axis", &axis, 1); ImGui::SameLine(); ImGui::RadioButton("Z axis", &axis, 2); if (ImGui::Button("Update")) { - //tmp - auto eval_scalar_field_ASO = [&](const VectorType& input_pos) -> Scalar - { - using FitASO = Ponca::BasketDiff< - Ponca::Basket< - PPAdapter, - SmoothWeightFunc, - Ponca::OrientedSphereFit>, - Ponca::DiffType::FitSpaceDer, - Ponca::OrientedSphereDer, - Ponca::MlsSphereFitDer, - Ponca::CurvatureEstimatorBase, - Ponca::NormalDerivativesCurvatureEstimator>; - VectorType current_pos = input_pos; - Scalar current_value = std::numeric_limits::max(); - for(int mm = 0; mm < mlsIter; ++mm) - { - FitASO fit; - fit.setWeightFunc(SmoothWeightFunc(NSize)); - fit.init(current_pos); // weighting function using current pos (not input pos) - for(int j : tree.range_neighbors(current_pos, NSize)) { - fit.addNeighbor(tree.point_data()[j]); - } - //std::cout<<"Final : "<::max(); - } - } - return current_value; - }; - VectorType lower(-2,-2,-2),upper(2,2,2); - auto distX= [](const VectorType &v){ return v[0];}; - auto mySlicer = makeSlicer("slicer",eval_scalar_field_ASO,lower, upper, 1024,axis); - mySlicer.regiterSlicer(slice); + auto mySlicer = registerRegularSlicer("slicer", eval_scalar_field_ASO,lower, upper, + isHDSlicer?1024:256, axis, slice); } ImGui::SameLine(); ImGui::PopItemWidth(); diff --git a/src/polyscopeSlicer.hpp b/src/polyscopeSlicer.hpp index b3e7d97..cf13bb5 100644 --- a/src/polyscopeSlicer.hpp +++ b/src/polyscopeSlicer.hpp @@ -1,127 +1,83 @@ #pragma once #include -#include +#include #include #include + +/** + * Create and register a polyscope surface mesh that slices a given implicit function. + * This function uses a regular grid and evaluates the implicit function at the grid vertex positions. + * + * The implicit function could be a C/C++ function or any lambda taking as input a Point and returning a double. + * + * Note: if openmp is available, the implicit function is evaluated in parallel. + * + * @param name the name of the slicer + * @param implicit the implicit function that maps points (type Point) to scalars. + * @param lowerBound the bbox lower bound of the implicit function domain. + * @param upperBound the bbox upper bound of the implicit function domain. + * @param nbSteps step size for the regular grid construction (eg 256) + * @param axis the axis to slice {0,1,2}. + * @param slice the relative position of the slice in [0,1] + * + * @return a pointer to the polyscope surface mesh object + */ template -struct PolyscopeSlicer +polyscope::SurfaceMesh* registerRegularSlicer(std::string name, + Functor &implicit, + const Point &lowerBound, + const Point &upperBound, + size_t nbSteps, + size_t axis, + float slice=0.0) { + size_t sliceid = static_cast(std::floor(slice*nbSteps)); - Functor _implicit; - Point _lowerBound; - Point _upperBound; - size_t _nbSteps; - size_t _axis; - size_t _slice; - std::string _name; + auto dim1 = (axis+1)%3; + auto dim2 = (axis+2)%3; - PolyscopeSlicer(std::string name, - Functor &implicit, - const Point &lowerBound, - const Point &upperBound, - size_t nbSteps, - size_t axis): _name(name),_implicit(implicit),_lowerBound(lowerBound), - _upperBound(upperBound),_nbSteps(nbSteps),_axis(axis) - { - } - - polyscope::SurfaceMesh* regiterSlicer(size_t slice=0) - { - _slice = slice; - auto dim1 = (_axis+1)%3; - auto dim2 = (_axis+2)%3; - - double du = (_upperBound[dim1]-_lowerBound[dim1])/(double)_nbSteps; - double dv = (_upperBound[dim2]-_lowerBound[dim2])/(double)_nbSteps; - double dw = (_upperBound[_axis]-_lowerBound[_axis])/(double)_nbSteps; - - double u = _lowerBound[dim1]; - double v = _lowerBound[dim2]; - double w = _lowerBound[_axis] + _slice*dw; - - Point p; - Point vu,vv; - switch (_axis) { - case 0: p=Point(w,u,v); vu=Point(0,du,0); vv=Point(0,0,dv);break; - case 1: p=Point(u,w,v); vu=Point(du,0,0); vv=Point(0,0,dv);break; - case 2: p=Point(u,v,w); vu=Point(du,0,0); vv=Point(0,dv,0);break; - } - - std::vector vertices(_nbSteps*_nbSteps); - std::vector values(_nbSteps*_nbSteps); - std::vector> faces; - std::array face; - - for(size_t id=0; id < _nbSteps*_nbSteps; ++id) - { - auto i = id % _nbSteps; - auto j = id / _nbSteps; - p = _lowerBound + i*vu + j*vv; - p[_axis] += _slice*dw; - - vertices[id] = p; - values[id] = _implicit(p); - face = { id, id+1, id+1+_nbSteps, id+_nbSteps }; - if (((i+1) < _nbSteps) && ((j+1)<_nbSteps)) - faces.push_back(face); - } - std::cout<<"vertices: "<addVertexScalarQuantity("values",values)->setEnabled(true); - return psm; + double du = (upperBound[dim1]-lowerBound[dim1])/(double)nbSteps; + double dv = (upperBound[dim2]-lowerBound[dim2])/(double)nbSteps; + double dw = (upperBound[axis]-lowerBound[axis])/(double)nbSteps; + + double u = lowerBound[dim1]; + double v = lowerBound[dim2]; + double w = lowerBound[axis] + sliceid*dw; + + Point p; + Point vu,vv; + switch (axis) { + case 0: p=Point(w,u,v); vu=Point(0,du,0); vv=Point(0,0,dv);break; + case 1: p=Point(u,w,v); vu=Point(du,0,0); vv=Point(0,0,dv);break; + case 2: p=Point(u,v,w); vu=Point(du,0,0); vv=Point(0,dv,0);break; } + std::vector vertices(nbSteps*nbSteps); + std::vector values(nbSteps*nbSteps); + std::vector> faces; + faces.reserve(nbSteps*nbSteps); + std::array face; - void updateSlice(size_t slice) + //Regular grid construction + for(size_t id=0; id < nbSteps*nbSteps; ++id) { - auto dim1 = (_axis+1)%3; - auto dim2 = (_axis+2)%3; - double du = (_upperBound[dim1]-_lowerBound[dim1])/(double)_nbSteps; - double dv = (_upperBound[dim2]-_lowerBound[dim2])/(double)_nbSteps; - double dw = (_upperBound[_axis]-_lowerBound[_axis])/(double)_nbSteps; - double u = _lowerBound[dim1]; - double v = _lowerBound[dim2]; - double w = _lowerBound[_axis] + slice*dw; - - Point p; - Point vu,vv; - switch (_axis) { - case 0: p=Point(w,u,v); vu=Point(0,du,0); vv=Point(0,0,dv);break; - case 1: p=Point(u,w,v); vu=Point(du,0,0); vv=Point(0,0,dv);break; - case 2: p=Point(u,v,w); vu=Point(du,0,0); vv=Point(0,dv,0);break; - } - - std::vector vertices(_nbSteps*_nbSteps); - std::vector values(_nbSteps*_nbSteps); - for(size_t id=0; id < _nbSteps*_nbSteps; ++id) - { - auto i = id % _nbSteps; - auto j = id / _nbSteps; - p = _lowerBound + i*vu + j*vv; - vertices[id] = p; - values[id] = _implicit(p); - } - //auto psm = polyscope::registerSurfaceMesh("Slice "+std::to_string(axis), vertices,faces); - polyscope::getSurfaceMesh(_name)->updateVertexPositions(vertices); - polyscope::getSurfaceMesh(_name)->addVertexDistanceQuantity("values",values); + auto i = id % nbSteps; + auto j = id / nbSteps; + p = lowerBound + i*vu + j*vv; + p[axis] += sliceid*dw; + vertices[id] = p; + face = { id, id+1, id+1+nbSteps, id+nbSteps }; + if (((i+1) < nbSteps) && ((j+1) -PolyscopeSlicer makeSlicer(std::string name, - Functor & implicit, - const Point &lowerBound, - const Point &upperBound, - size_t nbSteps, - size_t axis) -{ - return PolyscopeSlicer(name, - implicit, - lowerBound, - upperBound, - nbSteps, - axis); + auto psm = polyscope::registerSurfaceMesh(name, vertices,faces); + psm->addVertexScalarQuantity("values",values)->setEnabled(true); + return psm; }