Skip to content

Commit

Permalink
Merge pull request #14 from dcoeurjo/main
Browse files Browse the repository at this point in the history
Polyscope slicer and visualization of the ASO scalar field as example
  • Loading branch information
nmellado authored Dec 15, 2023
2 parents e9fec31 + 78687d0 commit 65f49ad
Show file tree
Hide file tree
Showing 3 changed files with 203 additions and 32 deletions.
1 change: 1 addition & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ endif()
# Create an executable
add_executable( poncascope
src/poncaAdapters.hpp
src/polyscopeSlicer.hpp
src/main.cpp
)

Expand Down
125 changes: 93 additions & 32 deletions src/main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,12 @@
#include <igl/readOBJ.h>
#include <igl/per_vertex_normals.h>

#include "polyscope/messages.h"
#include "polyscope/point_cloud.h"

#include <Ponca/Fitting>
#include <Ponca/SpatialPartitioning>
#include "poncaAdapters.hpp"
#include "polyscopeSlicer.hpp"

#include <iostream>
#include <utility>
Expand Down Expand Up @@ -38,6 +38,14 @@ Scalar pointRadius = 0.005; /// < display radius of the point cloud
bool useKnnGraph = false; /// < use k-neighbor graph instead of kdtree



// Slicer
float slice = 0.f;
int axis = 0;
bool isHDSlicer=false;
VectorType lower, upper;


/// Convenience function measuring and printing the processing time of F
template <typename Functor>
void measureTime( const std::string &actionName, Functor F ){
Expand Down Expand Up @@ -110,11 +118,8 @@ void colorizeKnn() {
/// \note Functor is called only if fit is stable
template<typename FitT, typename Functor>
void processPointCloud(const typename FitT::WeightFunction& w, Functor f){

int nvert = tree.index_data().size();

#pragma omp parallel for
for (int i = 0; i < nvert; ++i) {
for (int i = 0; i < tree.index_data().size(); ++i) {
VectorType pos = tree.point_data()[i].pos();

for( int mm = 0; mm < mlsIter; ++mm) {
Expand Down Expand Up @@ -181,49 +186,80 @@ void estimateDifferentialQuantities_impl(const std::string& name) {
});
}

using FitDry = Ponca::Basket<PPAdapter, SmoothWeightFunc, Ponca::DryFit>;

using FitPlane = Ponca::Basket<PPAdapter, SmoothWeightFunc, Ponca::CovariancePlaneFit>;
using FitPlaneDiff = Ponca::BasketDiff<
FitPlane,
Ponca::DiffType::FitSpaceDer,
Ponca::CovariancePlaneDer,
Ponca::CurvatureEstimatorBase, Ponca::NormalDerivativesCurvatureEstimator>;

using FitAPSS = Ponca::Basket<PPAdapter, SmoothWeightFunc, Ponca::OrientedSphereFit>;
using FitAPSSDiff = Ponca::BasketDiff<
FitAPSS,
Ponca::DiffType::FitSpaceDer,
Ponca::OrientedSphereDer,
Ponca::CurvatureEstimatorBase, Ponca::NormalDerivativesCurvatureEstimator>;

using FitASO = FitAPSS;
using FitASODiff = Ponca::BasketDiff<
FitASO,
Ponca::DiffType::FitSpaceDer,
Ponca::OrientedSphereDer, Ponca::MlsSphereFitDer,
Ponca::CurvatureEstimatorBase, Ponca::NormalDerivativesCurvatureEstimator>;

/// Compute curvature using Covariance Plane fitting
/// \see estimateDifferentialQuantities_impl
void estimateDifferentialQuantitiesWithPlane() {
estimateDifferentialQuantities_impl<
Ponca::BasketDiff<
Ponca::Basket<PPAdapter, SmoothWeightFunc, Ponca::CovariancePlaneFit>,
Ponca::DiffType::FitSpaceDer,
Ponca::CovariancePlaneDer,
Ponca::CurvatureEstimatorBase, Ponca::NormalDerivativesCurvatureEstimator>>("PSS");
estimateDifferentialQuantities_impl<FitPlaneDiff>("PSS");
}

/// Compute curvature using APSS
/// \see estimateDifferentialQuantities_impl
void estimateDifferentialQuantitiesWithAPSS() {
estimateDifferentialQuantities_impl<
Ponca::BasketDiff<
Ponca::Basket<PPAdapter, SmoothWeightFunc, Ponca::OrientedSphereFit>,
Ponca::DiffType::FitSpaceDer,
Ponca::OrientedSphereDer,
Ponca::CurvatureEstimatorBase, Ponca::NormalDerivativesCurvatureEstimator>>("APSS");
inline void estimateDifferentialQuantitiesWithAPSS() {
estimateDifferentialQuantities_impl<FitAPSSDiff>("APSS");
}

/// Compute curvature using Algebraic Shape Operator
/// \see estimateDifferentialQuantities_impl
void estimateDifferentialQuantitiesWithASO() {
estimateDifferentialQuantities_impl<
Ponca::BasketDiff<
Ponca::Basket<PPAdapter, SmoothWeightFunc, Ponca::OrientedSphereFit>,
Ponca::DiffType::FitSpaceDer,
Ponca::OrientedSphereDer, Ponca::MlsSphereFitDer,
Ponca::CurvatureEstimatorBase, Ponca::NormalDerivativesCurvatureEstimator>>("ASO");
inline void estimateDifferentialQuantitiesWithASO() {
estimateDifferentialQuantities_impl<FitASODiff>("ASO");
}

/// Dry run: loop over all vertices + run MLS loops without computation
/// This function is useful to monitor the KdTree performances
void mlsDryRun() {
using DryFit = Ponca::Basket<PPAdapter, SmoothWeightFunc, Ponca::DryFit>;
inline void mlsDryRun() {
measureTime( "[Ponca] Dry run MLS ", []() {
processPointCloud<DryFit>(
SmoothWeightFunc(NSize), [](int, const DryFit&, const VectorType& ){ });
processPointCloud<FitDry>(
SmoothWeightFunc(NSize), [](int, const FitDry&, const VectorType& ){ });
});
}

///Evaluate scalar field for generic FitType.
///// \tparam FitT Defines the type of estimator used for computation
template<typename FitT, bool isSigned = true>
Scalar evalScalarField_impl(const VectorType& input_pos)
{
VectorType current_pos = input_pos;
Scalar current_value = std::numeric_limits<Scalar>::max();
for(int mm = 0; mm < mlsIter; ++mm)
{
FitT fit;
fit.setWeightFunc(SmoothWeightFunc(NSize));
fit.init(current_pos); // weighting function using current pos (not input pos)
auto res = fit.computeWithIds(tree.range_neighbors(current_pos, NSize), tree.point_data());
if(res == Ponca::STABLE) {
current_pos = fit.project(input_pos); // always project input pos
current_value = isSigned ? fit.potential(input_pos) : std::abs(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<Scalar>::max();
}
}
return current_value;
}

/// Define Polyscope callbacks
void callback() {
Expand Down Expand Up @@ -253,13 +289,34 @@ void callback() {
if (ImGui::Button("APSS")) estimateDifferentialQuantitiesWithAPSS();
ImGui::SameLine();
if (ImGui::Button("ASO")) estimateDifferentialQuantitiesWithASO();


ImGui::Separator();

ImGui::Text("Implicit function slicer");
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);
const char* items[] = { "ASO", "APSS", "PSS"};
static int item_current = 0;
ImGui::Combo("Fit function", &item_current, items, IM_ARRAYSIZE(items));
if (ImGui::Button("Update"))
{
switch(item_current)
{
case 0: registerRegularSlicer("slicer", evalScalarField_impl<FitASO, true>,lower, upper, isHDSlicer?1024:256, axis, slice); break;
case 1: registerRegularSlicer("slicer", evalScalarField_impl<FitAPSS, true>,lower, upper, isHDSlicer?1024:256, axis, slice); break;
case 2: registerRegularSlicer("slicer", evalScalarField_impl<FitPlane, false>,lower, upper, isHDSlicer?1024:256, axis, slice); break;
}
}
ImGui::SameLine();
ImGui::PopItemWidth();
}

int main(int argc, char **argv) {
// Options
polyscope::options::autocenterStructures = true;
polyscope::options::autocenterStructures = false;
polyscope::options::programName = "poncascope";
polyscope::view::windowWidth = 1024;
polyscope::view::windowHeight = 1024;
Expand All @@ -284,6 +341,10 @@ int main(int argc, char **argv) {
return EXIT_FAILURE;
}

//Bounding Box (used in the slicer)
lower = cloudV.colwise().minCoeff();
upper = cloudV.colwise().maxCoeff();

// Build Ponca KdTree
measureTime( "[Ponca] Build KdTree", []() {
buildKdTree(cloudV, cloudN, tree);
Expand All @@ -302,4 +363,4 @@ int main(int argc, char **argv) {

delete knnGraph;
return EXIT_SUCCESS;
}
}
109 changes: 109 additions & 0 deletions src/polyscopeSlicer.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
//MIT License
//
//Copyright (c) 2023 Ponca Development Group
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
//of this software and associated documentation files (the "Software"), to deal
//in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
//furnished to do so, subject to the following conditions:
//
//The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
//IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
//OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
//SOFTWARE.

/// \file This file contains a slicer of the implicit function as a Polyscope surface mesh
/// \author David Coeurjolly <david.coeurjolly@cnrs.fr >

#pragma once
#include <vector>
#include <array>
#include <polyscope/surface_mesh.h>
#include <polyscope/polyscope.h>


/**
* 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<typename Point,typename Functor>
polyscope::SurfaceMesh* registerRegularSlicer(const std::string &name,
const Functor &implicit,
const Point &lowerBound,
const Point &upperBound,
size_t nbSteps,
size_t axis,
float slice=0.0)
{
size_t sliceid = static_cast<size_t>(std::floor(slice*nbSteps));

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] + 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<Point> vertices(nbSteps*nbSteps);
std::vector<double> values(nbSteps*nbSteps);
std::vector<std::array<size_t,4>> faces;
faces.reserve(nbSteps*nbSteps);
std::array<size_t,4> face;

//Regular grid construction
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] += sliceid*dw;
vertices[id] = p;
face = { id, id+1, id+1+nbSteps, id+nbSteps };
if (((i+1) < nbSteps) && ((j+1)<nbSteps))
faces.push_back(face);
}

//Evaluating the implicit function (in parallel using openmp)
#pragma omp parallel for default(none) shared(nbSteps,values,vertices,implicit)
for(int id=0; id < nbSteps*nbSteps; ++id)
values[id] = implicit(vertices[id]);

//Polyscope registration
auto psm = polyscope::registerSurfaceMesh(name, vertices,faces);
psm->addVertexScalarQuantity("values",values)->setEnabled(true);
return psm;
}

0 comments on commit 65f49ad

Please sign in to comment.