diff --git a/Apps/Examples/viewer.cpp b/Apps/Examples/viewer.cpp index c6c8aae3b1a339686ed6ef2864ba6ef59e10d5b4..880beaca2c32efccce288fcca92f8296d1deaf58 100644 --- a/Apps/Examples/viewer.cpp +++ b/Apps/Examples/viewer.cpp @@ -186,7 +186,6 @@ void Viewer::importMesh(std::string& filename) m_render->initPrimitives(myMap, allDarts, Algo::Render::GL2::TRIANGLES) ; bb = Algo::Geometry::computeBoundingBox(myMap, position) ; - gPosObj = bb.center() ; normalBaseSize = bb.diagSize() / 100.0f ; vertexBaseSize = normalBaseSize * 2.0f ; @@ -198,7 +197,7 @@ void Viewer::importMesh(std::string& filename) m_positionVBO->updateData(position) ; m_normalVBO->updateData(normal) ; - setParamObject(bb.maxSize(), gPosObj.data()) ; + setParamObject(bb.maxSize(), bb.center().data()) ; updateGLMatrices() ; } diff --git a/Apps/Examples/viewer.h b/Apps/Examples/viewer.h index ec11e7bf3db6cbe07ce468222d2e751d1de5b7c0..f36f90b1e65eda5b8c1ac34c29ce4c0a3168da4b 100644 --- a/Apps/Examples/viewer.h +++ b/Apps/Examples/viewer.h @@ -78,7 +78,6 @@ public: float shininess ; Geom::BoundingBox bb ; - Geom::Vec3f gPosObj ; float normalBaseSize ; float normalScaleFactor ; float vertexBaseSize ; diff --git a/include/Algo/Decimation/edgeSelector.h b/include/Algo/Decimation/edgeSelector.h index 8ec430ba6fda1d16be4cb9216778bd31776422c4..4e01c5ff4ab0a1964a1a8dbf032056d12d87a14d 100644 --- a/include/Algo/Decimation/edgeSelector.h +++ b/include/Algo/Decimation/edgeSelector.h @@ -27,6 +27,11 @@ #include "Algo/Decimation/selector.h" +#include "Container/fakeAttribute.h" +#include "Utils/qem.h" +#include "Utils/quadricRGBfunctions.h" +#include "Algo/Geometry/curvature.h" + namespace CGoGN { @@ -246,12 +251,17 @@ private: } CurvatureEdgeInfo ; typedef NoMathIOAttribute EdgeInfo ; + Geom::BoundingBox bb ; + REAL radius ; + typename PFP::TVEC3 normal ; AttributeHandler edgeInfo ; - typename PFP::TREAL k1 ; - typename PFP::TREAL k2 ; - typename PFP::TVEC3 K1 ; - typename PFP::TVEC3 K2 ; + typename PFP::TREAL edgeangle ; + typename PFP::TREAL kmax ; + typename PFP::TREAL kmin ; + typename PFP::TVEC3 Kmax ; + typename PFP::TVEC3 Kmin ; + typename PFP::TVEC3 Knormal ; std::multimap edges ; typename std::multimap::iterator cur ; @@ -266,6 +276,9 @@ public: EdgeSelector_Curvature(MAP& m, typename PFP::TVEC3& pos, std::vector*>& approx, const FunctorSelect& select = SelectorTrue()) : EdgeSelector(m, pos, approx, select) { + bb = Algo::Geometry::computeBoundingBox(m, pos) ; + radius = bb.diagSize() * 0.003 ; + normal = m.template getAttribute(VERTEX, "normal") ; if(!normal.isValid()) { @@ -273,29 +286,40 @@ public: Algo::Geometry::computeNormalVertices(m, pos, normal) ; } - k1 = m.template getAttribute(VERTEX, "k1") ; - k2 = m.template getAttribute(VERTEX, "k2") ; - K1 = m.template getAttribute(VERTEX, "K1") ; - K2 = m.template getAttribute(VERTEX, "K2") ; + edgeangle = m.template getAttribute(VERTEX, "edgeangle") ; + if(!edgeangle.isValid()) + { + edgeangle = m.template addAttribute(EDGE, "edgeangle") ; + Algo::Geometry::computeAnglesBetweenNormalsOnEdges(m, pos, edgeangle) ; + } + + kmax = m.template getAttribute(VERTEX, "kmax") ; + kmin = m.template getAttribute(VERTEX, "kmin") ; + Kmax = m.template getAttribute(VERTEX, "Kmax") ; + Kmin = m.template getAttribute(VERTEX, "Kmin") ; + Knormal = m.template getAttribute(VERTEX, "Knormal") ; // as all these attributes are computed simultaneously by computeCurvatureVertices // one can assume that if one of them is not valid, the others must be created too - if(!k1.isValid()) + if(!kmax.isValid()) { - k1 = m.template addAttribute(VERTEX, "k1") ; - k2 = m.template addAttribute(VERTEX, "k2") ; - K1 = m.template addAttribute(VERTEX, "K1") ; - K2 = m.template addAttribute(VERTEX, "K2") ; - Algo::Geometry::computeCurvatureVertices(m, this->m_position, normal, k1, k2, K1, K2) ; + kmax = m.template addAttribute(VERTEX, "kmax") ; + kmin = m.template addAttribute(VERTEX, "kmin") ; + Kmax = m.template addAttribute(VERTEX, "Kmax") ; + Kmin = m.template addAttribute(VERTEX, "Kmin") ; + Knormal = m.template addAttribute(VERTEX, "Knormal") ; + Algo::Geometry::computeCurvatureVertices_NormalCycles(m, radius, pos, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ; } edgeInfo = m.template addAttribute(EDGE, "edgeInfo") ; } ~EdgeSelector_Curvature() { - this->m_map.removeAttribute(k1) ; - this->m_map.removeAttribute(k2) ; - this->m_map.removeAttribute(K1) ; - this->m_map.removeAttribute(K2) ; + this->m_map.removeAttribute(edgeangle) ; + this->m_map.removeAttribute(kmax) ; + this->m_map.removeAttribute(kmin) ; + this->m_map.removeAttribute(Kmax) ; + this->m_map.removeAttribute(Kmin) ; + this->m_map.removeAttribute(Knormal) ; this->m_map.removeAttribute(edgeInfo) ; } SelectorType getType() { return S_Curvature ; } diff --git a/include/Algo/Decimation/edgeSelector.hpp b/include/Algo/Decimation/edgeSelector.hpp index dce20173231ce5b485dac6125aae3207da43e6e9..83d6c31dc5b0f3bb2e4fca053dab07837a3eb573 100644 --- a/include/Algo/Decimation/edgeSelector.hpp +++ b/include/Algo/Decimation/edgeSelector.hpp @@ -795,14 +795,14 @@ void EdgeSelector_Curvature::updateAfterCollapse(Dart d2, Dart dd2) MAP& m = this->m_map ; normal[d2] = Algo::Geometry::vertexNormal(m, d2, this->m_position) ; - Algo::Geometry::computeCurvatureVertex(m, d2, this->m_position, normal, k1, k2, K1, K2) ; + Algo::Geometry::computeCurvatureVertex_NormalCycles(m, d2, radius, this->m_position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ; Dart vit = d2 ; do { Dart nVert = m.phi1(vit) ; normal[nVert] = Algo::Geometry::vertexNormal(m, nVert, this->m_position) ; - Algo::Geometry::computeCurvatureVertex(m, nVert, this->m_position, normal, k1, k2, K1, K2) ; + Algo::Geometry::computeCurvatureVertex_NormalCycles(m, nVert, radius, this->m_position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ; updateEdgeInfo(m.phi1(vit), false) ; // must recompute some edge infos in the if(vit == d2 || vit == dd2) // neighborhood of the collapsed edge @@ -891,11 +891,11 @@ void EdgeSelector_Curvature::computeEdgeInfo(Dart d, EdgeInfo& einfo) // compute things on the coarse version of the mesh normal[newV] = Algo::Geometry::vertexNormal(m, d2, this->m_position) ; - Algo::Geometry::computeCurvatureVertex(m, d2, this->m_position, normal, k1, k2, K1, K2) ; + Algo::Geometry::computeCurvatureVertex_NormalCycles(m, d2, radius, this->m_position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ; - VEC3 norm = normal[newV] ; - REAL mCurv = ( k1[newV] + k2[newV] ) / REAL(2) ; - VEC3 cDir1 = K1[newV] ; +// VEC3 norm = normal[newV] ; + REAL mCurv = (kmax[newV] + kmin[newV]) / REAL(2) ; +// VEC3 cDir1 = Kmax[newV] ; // vertex split to reset the initial connectivity and embeddings m.insertTrianglePair(d, d2, dd2) ; @@ -904,17 +904,17 @@ void EdgeSelector_Curvature::computeEdgeInfo(Dart d, EdgeInfo& einfo) REAL err = 0 ; - REAL norm_deviation_1 = REAL(1) / abs(norm * normal[v1]) ; - REAL norm_deviation_2 = REAL(1) / abs(norm * normal[v2]) ; - err += norm_deviation_1 + norm_deviation_2 ; +// REAL norm_deviation_1 = REAL(1) / abs(norm * normal[v1]) ; +// REAL norm_deviation_2 = REAL(1) / abs(norm * normal[v2]) ; +// err += norm_deviation_1 + norm_deviation_2 ; - REAL mCurv_deviation_1 = abs(mCurv - (k1[v1] + k2[v1] / REAL(2))) ; - REAL mCurv_deviation_2 = abs(mCurv - (k1[v2] + k2[v2] / REAL(2))) ; + REAL mCurv_deviation_1 = abs(mCurv - (kmax[v1] + kmin[v1] / REAL(2))) ; + REAL mCurv_deviation_2 = abs(mCurv - (kmax[v2] + kmin[v2] / REAL(2))) ; err += mCurv_deviation_1 + mCurv_deviation_2 ; - REAL cDir1_deviation_1 = REAL(1) / abs(cDir1 * K1[v1]) ; - REAL cDir1_deviation_2 = REAL(1) / abs(cDir1 * K1[v2]) ; - err += cDir1_deviation_1 + cDir1_deviation_2 ; +// REAL cDir1_deviation_1 = REAL(1) / abs(cDir1 * Kmax[v1]) ; +// REAL cDir1_deviation_2 = REAL(1) / abs(cDir1 * Kmax[v2]) ; +// err += cDir1_deviation_1 + cDir1_deviation_2 ; einfo.it = edges.insert(std::make_pair(err, d)) ; einfo.valid = true ; diff --git a/include/Algo/Decimation/selector.h b/include/Algo/Decimation/selector.h index ee0f8484a39440d31eb7301ef296c2805543957b..83ccb277bc0030f23e733d354fe717d9146b7c1e 100644 --- a/include/Algo/Decimation/selector.h +++ b/include/Algo/Decimation/selector.h @@ -25,11 +25,6 @@ #ifndef __SELECTOR_H__ #define __SELECTOR_H__ -#include "Utils/qem.h" -#include "Utils/quadricRGBfunctions.h" -#include "Container/fakeAttribute.h" -#include "Algo/Geometry/curvature.h" - namespace CGoGN { @@ -39,7 +34,6 @@ namespace Algo namespace Decimation { - enum SelectorType { S_MapOrder, @@ -53,11 +47,9 @@ enum SelectorType S_hLightfield } ; - template class ApproximatorGen ; template class Approximator ; - template class EdgeSelector { @@ -91,5 +83,4 @@ public: } // namespace CGoGN - #endif diff --git a/include/Algo/Geometry/centroid.hpp b/include/Algo/Geometry/centroid.hpp index 7fbe08aade4443d62fbe267e54da6ef9521d431b..ae7af359300f2908bc7d21a7c6eeaff335c7abd3 100644 --- a/include/Algo/Geometry/centroid.hpp +++ b/include/Algo/Geometry/centroid.hpp @@ -67,11 +67,10 @@ EMB volumeCentroidGen(typename PFP::MAP& map, Dart d, const EMBV& attributs) std::vector visitedFaces; // Faces that are traversed visitedFaces.reserve(100); visitedFaces.push_back(d); // Start with the face of d - std::vector::iterator face; mark.markOrbit(VERTEX, d) ; - for(face = visitedFaces.begin(); face != visitedFaces.end(); ++face) + for(std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) { Dart e = *face ; @@ -90,7 +89,6 @@ EMB volumeCentroidGen(typename PFP::MAP& map, Dart d, const EMBV& attributs) } while(e != *face) ; } - center /= double(count) ; return center ; } diff --git a/include/Algo/Geometry/curvature.h b/include/Algo/Geometry/curvature.h index ba077137493fff0902498859ee2aea3c6e9e7dc7..ea9e4e8e08b98650c09568a0e5a1ac68dc19b2b7 100644 --- a/include/Algo/Geometry/curvature.h +++ b/include/Algo/Geometry/curvature.h @@ -95,7 +95,7 @@ void computeCurvatureVertices_NormalCycles( typename PFP::REAL radius, const typename PFP::TVEC3& position, const typename PFP::TVEC3& normal, - const typename PFP::TREAL& angles, + const typename PFP::TREAL& edgeangle, typename PFP::TREAL& kmax, typename PFP::TREAL& kmin, typename PFP::TVEC3& Kmax, @@ -110,7 +110,7 @@ void computeCurvatureVertex_NormalCycles( typename PFP::REAL radius, const typename PFP::TVEC3& position, const typename PFP::TVEC3& normal, - const typename PFP::TREAL& angles, + const typename PFP::TREAL& edgeangle, typename PFP::TREAL& kmax, typename PFP::TREAL& kmin, typename PFP::TVEC3& Kmax, diff --git a/include/Algo/Geometry/curvature.hpp b/include/Algo/Geometry/curvature.hpp index bea074209de6bc3ae1a22f80b1fe9709caa647f3..815b1c4ccc08d8ad3620a2f3298ed9b846f05c7f 100644 --- a/include/Algo/Geometry/curvature.hpp +++ b/include/Algo/Geometry/curvature.hpp @@ -304,7 +304,7 @@ void computeCurvatureVertices_NormalCycles( typename PFP::REAL radius, const typename PFP::TVEC3& position, const typename PFP::TVEC3& normal, - const typename PFP::TREAL& angles, + const typename PFP::TREAL& edgeangle, typename PFP::TREAL& kmax, typename PFP::TREAL& kmin, typename PFP::TVEC3& Kmax, @@ -318,7 +318,7 @@ void computeCurvatureVertices_NormalCycles( if(select(d) && !marker.isMarked(d)) { marker.mark(d); - computeCurvatureVertex_NormalCycles(map, d, radius, position, normal, angles, kmax, kmin, Kmax, Kmin, Knormal) ; + computeCurvatureVertex_NormalCycles(map, d, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ; } } } @@ -330,7 +330,7 @@ void computeCurvatureVertex_NormalCycles( typename PFP::REAL radius, const typename PFP::TVEC3& position, const typename PFP::TVEC3& normal, - const typename PFP::TREAL& angles, + const typename PFP::TREAL& edgeangle, typename PFP::TREAL& kmax, typename PFP::TREAL& kmin, typename PFP::TVEC3& Kmax, @@ -354,16 +354,15 @@ void computeCurvatureVertex_NormalCycles( for (std::vector::const_iterator it = vd1.begin(); it != vd1.end(); ++it) { const VEC3 e = position[map.phi2(*it)] - position[*it] ; - tensor += Geom::transposed_vectors_mult(e,e) * angles[*it] * (1 / e.norm()) ; + tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) ; } - // border const std::vector& vd2 = neigh.getBorder() ; for (std::vector::const_iterator it = vd2.begin(); it != vd2.end(); ++it) { const VEC3 e = position[map.phi2(*it)] - position[*it] ; const REAL alpha = neigh.intersect_SphereEdge(*it, map.phi2(*it)) ; - tensor += Geom::transposed_vectors_mult(e,e) * angles[*it] * (1 / e.norm()) * alpha ; + tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) * alpha ; } tensor /= neigh.getArea() ; diff --git a/include/Algo/Geometry/intersection.h b/include/Algo/Geometry/intersection.h index 87178381a6c00073849c393444eba86fc910bb01..5c6dd7e2f087c1e6c002ee786d2261f1fb5b0586 100644 --- a/include/Algo/Geometry/intersection.h +++ b/include/Algo/Geometry/intersection.h @@ -27,6 +27,7 @@ #include "Geometry/basic.h" #include "Geometry/intersection.h" +#include "Geometry/inclusion.h" namespace CGoGN { @@ -75,11 +76,17 @@ bool intersectionSegmentConvexFace(typename PFP::MAP& map, Dart d, const typenam template bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, const typename PFP::TVEC3& positions) ; -} +/** + * + */ +template +bool intersectionSphereEdge(typename PFP::MAP& map, typename PFP::VEC3& center, typename PFP::REAL radius, Dart d, const typename PFP::TVEC3& positions, typename PFP::REAL& alpha) ; + +} // namespace Geometry -} +} // namespace Algo -} +} // namespace CGoGN #include "Algo/Geometry/intersection.hpp" diff --git a/include/Algo/Geometry/intersection.hpp b/include/Algo/Geometry/intersection.hpp index b66a9aa33645d652e3618ba211ab4a0ce32f634a..6127b4dfd6f3f84f9882cf8264785ba3ddb186ce 100644 --- a/include/Algo/Geometry/intersection.hpp +++ b/include/Algo/Geometry/intersection.hpp @@ -46,10 +46,10 @@ bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const typename P const float SMALL_NUM = std::numeric_limits::min() * 5.0f; VEC3 p1 = positions[d]; - VEC3 n = faceNormal(map,d,positions); + VEC3 n = faceNormal(map, d, positions); VEC3 w0 = P - p1; - float a = -(n*w0); - float b = n*Dir; + float a = -(n * w0); + float b = n * Dir; if (fabs(b) < SMALL_NUM) return false; @@ -225,8 +225,25 @@ bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, co return intersection; } +template +bool intersectionSphereEdge(typename PFP::MAP& map, typename PFP::VEC3& center, typename PFP::REAL radius, Dart d, const typename PFP::TVEC3& positions, typename PFP::REAL& alpha) +{ + typename PFP::VEC3& p1 = position[d]; + typename PFP::VEC3& p2 = position[map.phi1(d)]; + if(Geom::isPointInSphere(p1, center, radius) && !Geom::isPointInSphere(p2, center, radius)) + { + VEC3 p = p1 - center; + VEC3 qminusp = p2 - center - p; + REAL s = p * qminusp; + REAL n2 = qminusp.norm2(); + alpha = (- s + sqrt(s*s + n2 * (radius*radius - p.norm2()))) / n2; + return true ; + } + return false ; } -} +} // namespace Geometry -} +} // namespace Algo + +} // namespace CGoGN diff --git a/include/Algo/Selection/collector.h b/include/Algo/Selection/collector.h index e3afe986d03fbff070197b9ba9b2e0e1bd9bd352..5cbcd9df6e058962a47830f06adb5649e7b30959 100644 --- a/include/Algo/Selection/collector.h +++ b/include/Algo/Selection/collector.h @@ -42,8 +42,9 @@ namespace Selection { /********************************************************* - * Collector + * Generic Collector *********************************************************/ + template class Collector { @@ -69,6 +70,8 @@ public: virtual void collect() = 0; + bool apply(FunctorType& f); + void sort() { std::sort(insideVertices.begin(), insideVertices.end()); @@ -135,6 +138,7 @@ protected: typename PFP::REAL radius_2; typename PFP::VEC3 centerPosition; typename PFP::REAL area; + public: Collector_WithinSphere(typename PFP::MAP& mymap, const typename PFP::TVEC3& myposition) : Collector(mymap,myposition) diff --git a/include/Algo/Selection/collector.hpp b/include/Algo/Selection/collector.hpp index c6abf0ec673051884aa4f07cdee79496ac72881f..ad4f05e805e0d9c1890817ed9aedae38ad3f7718 100644 --- a/include/Algo/Selection/collector.hpp +++ b/include/Algo/Selection/collector.hpp @@ -31,15 +31,25 @@ namespace Algo namespace Selection { -/******************************************** - * GENERIC COLLECTOR - ********************************************/ +/********************************************************* + * Generic Collector + *********************************************************/ + template Collector::Collector(typename PFP::MAP& mymap, const typename PFP::TVEC3& myposition): map(mymap), position(myposition) {} +template +bool Collector::apply(FunctorType& f) +{ + for(std::vector::iterator iv = insideVertices.begin(); iv != insideVertices.end(); ++iv) + if(f(*iv)) + return true ; + return false ; +} + template std::ostream& operator<<(std::ostream &out, const Collector& c) { @@ -121,6 +131,7 @@ void Collector_WithinSphere::collect() vm.mark(this->centerDart); unsigned int i = 0; +// for(std::vector::iterator iv = this->insideVertices.begin(); iv != this->insideVertices.end(); ++iv) while (i < this->insideVertices.size()) { Dart end = this->insideVertices[i]; @@ -132,7 +143,7 @@ void Collector_WithinSphere::collect() const Dart f = this->map.phi1(e); const Dart g = this->map.phi1(f); - if (!this->isInside(f)) + if (! this->isInside(f)) { this->border.push_back(e); // add to border em.mark(e); @@ -163,17 +174,17 @@ void Collector_WithinSphere::collect() } } -template -typename PFP::REAL Collector_WithinSphere::intersect_SphereEdge(const Dart din, const Dart dout) -{ - typedef typename PFP::VEC3 VEC3; - typedef typename PFP::REAL REAL; - VEC3 p = this->position[din] - this->centerPosition; - VEC3 qminusp = this->position[dout] - this->centerPosition - p; - REAL s = p*qminusp; - REAL n2 = qminusp.norm2(); - return (- s + sqrt (s*s + n2*(this->radius_2 - p.norm2())))/(n2); -} +//template +//typename PFP::REAL Collector_WithinSphere::intersect_SphereEdge(const Dart din, const Dart dout) +//{ +// typedef typename PFP::VEC3 VEC3; +// typedef typename PFP::REAL REAL; +// VEC3 p = this->position[din] - this->centerPosition; +// VEC3 qminusp = this->position[dout] - this->centerPosition - p; +// REAL s = p * qminusp; +// REAL n2 = qminusp.norm2(); +// return (- s + sqrt(s * s + n2 * (this->radius_2 - p.norm2()))) / n2; +//} template void Collector_WithinSphere::computeArea() @@ -189,14 +200,14 @@ void Collector_WithinSphere::computeArea() const Dart g = this->map.phi1(f); if (this->isInside(g)) { // only f is outside - typename PFP::REAL alpha = this->intersect_SphereEdge (*it, f); - typename PFP::REAL beta = this->intersect_SphereEdge (g, f); + typename PFP::REAL alpha = this->intersect_SphereEdge(*it, f); + typename PFP::REAL beta = this->intersect_SphereEdge(g, f); this->area += (alpha+beta - alpha*beta) * Algo::Geometry::triangleArea(this->map, *it, this->position); } else { // f and g are outside - typename PFP::REAL alpha = this->intersect_SphereEdge (*it, f); - typename PFP::REAL beta = this->intersect_SphereEdge (*it, g); + typename PFP::REAL alpha = this->intersect_SphereEdge(*it, f); + typename PFP::REAL beta = this->intersect_SphereEdge(*it, g); this->area += alpha * beta * Algo::Geometry::triangleArea(this->map, *it, this->position); } } diff --git a/include/Algo/Selection/raySelectFunctor.hpp b/include/Algo/Selection/raySelectFunctor.hpp index 62c18f96f0fc74dcc588dc09004f5c3eaa0b0f13..8a735e0fcd0c310b8ace568ae7caadf6d7d50743 100644 --- a/include/Algo/Selection/raySelectFunctor.hpp +++ b/include/Algo/Selection/raySelectFunctor.hpp @@ -44,13 +44,13 @@ class FuncFaceInter: public FunctorMap { typedef typename PFP::MAP MAP; - protected: std::vector& m_faces; std::vector& m_Ipoints ; const typename PFP::VEC3& m_A; const typename PFP::VEC3& m_AB; const typename PFP::TVEC3& m_positions; + public: /** * @param map the map @@ -64,7 +64,7 @@ public: bool operator()(Dart d) { - const typename PFP::VEC3& Ta = m_positions[d];//this->m_map.getVertexEmb(d)->getPosition(); + const typename PFP::VEC3& Ta = m_positions[d]; Dart dd = this->m_map.phi1(d); Dart ddd = this->m_map.phi1(dd); @@ -72,13 +72,13 @@ public: do { // get back position of triangle Ta,Tb,Tc - const typename PFP::VEC3& Tb = m_positions[dd]; //this->m_map.getVertexEmb(dd)->getPosition(); - const typename PFP::VEC3& Tc = m_positions[ddd]; //this->m_map.getVertexEmb(ddd)->getPosition(); + const typename PFP::VEC3& Tb = m_positions[dd]; + const typename PFP::VEC3& Tc = m_positions[ddd]; typename PFP::VEC3 I; if (Geom::intersectionLineTriangle(m_A, m_AB, Ta, Tb, Tc, I)) { m_faces.push_back(d); - m_Ipoints.push_back(I) ; + m_Ipoints.push_back(I); notfound = false; } // next triangle if we are in polygon @@ -95,7 +95,6 @@ class FuncEdgeInter: public FunctorMap { typedef typename PFP::MAP MAP; - protected: std::vector& m_edges; const typename PFP::VEC3& m_A; @@ -103,6 +102,7 @@ protected: float m_AB2; float m_distMax; const typename PFP::TVEC3& m_positions; + public: /** * @param map the map @@ -119,9 +119,9 @@ public: bool operator()(Dart d) { // get back position of segment PQ - const typename PFP::VEC3& P = m_positions[d];//this->m_map.getVertexEmb(d)->getPosition(); + const typename PFP::VEC3& P = m_positions[d]; Dart dd = this->m_map.phi1(d); - const typename PFP::VEC3& Q = m_positions[dd];//this->m_map.getVertexEmb(dd)->getPosition(); + const typename PFP::VEC3& Q = m_positions[dd]; // the three distance to P, Q and (PQ) not used here float dist = Geom::squaredDistanceLine2Seg(m_A, m_AB, m_AB2, P, Q); @@ -161,7 +161,7 @@ public: bool operator()(Dart d) { - const typename PFP::VEC3& P = m_positions[d];//this->m_map.getVertexEmb(d)->getPosition(); + const typename PFP::VEC3& P = m_positions[d]; float dist = Geom::squaredDistanceLine2Point(m_A, m_AB, m_AB2, P); if (dist < m_distMax) { @@ -172,7 +172,7 @@ public: }; /** - * Fonctor which store the dart that correspond to the subpart of face + * Functor which store the dart that correspond to the subpart of face * that is intersected * Must be called in foreachface */ @@ -373,7 +373,6 @@ bool distnintOrdering(const std::pair& e1, con return (e1.first < e2.first); } - } //namespace Selection } //namespace Algo diff --git a/include/Algo/Selection/raySelector.h b/include/Algo/Selection/raySelector.h index 820319022e9f64c18bd6a58d635ebe406af37a7d..e3e7203c8e27966cf7714dcdb2019a446bf54bf3 100644 --- a/include/Algo/Selection/raySelector.h +++ b/include/Algo/Selection/raySelector.h @@ -37,43 +37,44 @@ namespace Algo namespace Selection { + /** * Function that does the selection of faces, returned darts are sorted from closest to farthest * @param map the map we want to test * @param good a dart selector - * @param rayA first point of ray (user side) - * @param rayAB vector of ray (directed ot the scene) - * @param vecFaces (out) vector to store dart of intersected faces + * @param rayA first point of ray (user side) + * @param rayAB direction of ray (directed to the scene) + * @param vecFaces (out) vector to store the darts of intersected faces */ template void facesRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& position, const FunctorSelect& good, const typename PFP::VEC3& rayA, const typename PFP::VEC3& rayAB, std::vector& vecFaces) { - std::vector iPoints ; + std::vector iPoints; // get back intersected faces vecFaces.clear(); Algo::Selection::FuncFaceInter ffi(map, position, vecFaces, iPoints, rayA, rayAB); map.foreach_orbit(FACE, ffi, good); - // compute all distances to observer for each dart of intersected face + // compute all distances to observer for each intersected face // and put them in a vector for sorting typedef std::pair DartDist; std::vector distndart; unsigned int nbi = vecFaces.size(); distndart.resize(nbi); - for (unsigned int i=0; i< nbi; ++i) + for (unsigned int i = 0; i < nbi; ++i) { distndart[i].second = vecFaces[i]; typename PFP::VEC3 V = iPoints[i] - rayA; - distndart[i].first = V*V; + distndart[i].first = V.norm2(); } // sort the vector of pair dist/dart std::sort(distndart.begin(), distndart.end(), distndartOrdering); // store sorted darts in returned vector - for (unsigned int i=0; i< nbi; ++i) + for (unsigned int i = 0; i < nbi; ++i) vecFaces[i] = distndart[i].second; } @@ -83,13 +84,13 @@ void facesRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& positi * @param rayA first point of ray (user side) * @param rayAB vector of ray (directed ot the scene) * @param vecEdges (out) vector to store dart of intersected edges - * @param dist radius of cylinder modelling the edge + * @param dist radius of the cylinder of selection */ template void edgesRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& position, const FunctorSelect& good, const typename PFP::VEC3& rayA, const typename PFP::VEC3& rayAB, std::vector& vecEdges, float dist) { - typename PFP::REAL dist2 = dist*dist; - typename PFP::REAL AB2 =rayAB*rayAB; + typename PFP::REAL dist2 = dist * dist; + typename PFP::REAL AB2 = rayAB * rayAB; // recuperation des aretes intersectees vecEdges.clear(); @@ -102,22 +103,22 @@ void edgesRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& positi unsigned int nbi = vecEdges.size(); distndart.resize(nbi); - // compute all distances to observer for each dart of middle of edge + // compute all distances to observer for each middle of intersected edge // and put them in a vector for sorting for (unsigned int i = 0; i < nbi; ++i) { Dart d = vecEdges[i]; distndart[i].second = d; - typename PFP::VEC3 V = (position[d] + position[map.phi1(d)])/2.0; + typename PFP::VEC3 V = (position[d] + position[map.phi1(d)]) / typename PFP::REAL(2); V -= rayA; - distndart[i].first = V*V; + distndart[i].first = V.norm2(); } // sort the vector of pair dist/dart std::sort(distndart.begin(), distndart.end(), distndartOrdering); // store sorted darts in returned vector - for (unsigned int i=0; i< nbi; ++i) + for (unsigned int i = 0; i < nbi; ++i) vecEdges[i] = distndart[i].second; } @@ -126,14 +127,15 @@ void edgesRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& positi * @param map the map we want to test * @param rayA first point of ray (user side) * @param rayAB vector of ray (directed ot the scene) - * @param vecEdges (out) vector to store dart of intersected vertices - * @param dist radius of sphere modelling the vertex + * @param vecVertices (out) vector to store dart of intersected vertices + * @param dist radius of the cylinder of selection */ template void verticesRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& position, const typename PFP::VEC3& rayA, const typename PFP::VEC3& rayAB, std::vector& vecVertices, float dist, const FunctorSelect& good= SelectorTrue()) { - typename PFP::REAL dist2 = dist*dist; - typename PFP::REAL AB2 = rayAB*rayAB; + typename PFP::REAL dist2 = dist * dist; + typename PFP::REAL AB2 = rayAB * rayAB; + // recuperation des sommets intersectes vecVertices.clear(); Algo::Selection::FuncVertexInter ffi(map, position, vecVertices, rayA, rayAB, AB2, dist2); @@ -145,21 +147,21 @@ void verticesRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& pos unsigned int nbi = vecVertices.size(); distndart.resize(nbi); - // compute all distances to observer for each dart of middle of edge + // compute all distances to observer for each intersected vertex // and put them in a vector for sorting for (unsigned int i = 0; i < nbi; ++i) { Dart d = vecVertices[i]; distndart[i].second = d; typename PFP::VEC3 V = position[d] - rayA; - distndart[i].first = V*V; + distndart[i].first = V.norm2(); } // sort the vector of pair dist/dart std::sort(distndart.begin(), distndart.end(), distndartOrdering); // store sorted darts in returned vector - for (unsigned int i=0; i< nbi; ++i) + for (unsigned int i = 0; i < nbi; ++i) vecVertices[i] = distndart[i].second; } @@ -173,11 +175,13 @@ void verticesRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& pos template void vertexRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& position, const typename PFP::VEC3& rayA, const typename PFP::VEC3& rayAB, Dart& vertex, const FunctorSelect& good = SelectorTrue()) { + std::vector vecFaces; + std::vector iPoints; + // recuperation des faces intersectes - std::vector vecFaces ; - std::vector iPoints ; - Algo::Selection::FuncFaceInter ffi(map, position, vecFaces, iPoints, rayA, rayAB) ; - map.foreach_orbit(FACE, ffi, good) ; + Algo::Selection::FuncFaceInter ffi(map, position, vecFaces, iPoints, rayA, rayAB); + map.foreach_orbit(FACE, ffi, good); + if(vecFaces.size() > 0) { typedef std::pair IndexedDist; @@ -189,38 +193,39 @@ void vertexRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& posit { distnint[i].second = i; typename PFP::VEC3 V = iPoints[i] - rayA; - distnint[i].first = V*V; + distnint[i].first = V.norm2(); } // sort the vector of pair dist/dart std::sort(distnint.begin(), distnint.end(), distnintOrdering); - // recuperation du poitns d'intersection sur la face la plus proche + // recuperation du point d'intersection sur la face la plus proche unsigned int first = distnint[0].second; - typename PFP::VEC3 ip = iPoints[first] ; + typename PFP::VEC3 ip = iPoints[first]; + // recuperation du sommet le plus proche du point d'intersection - Dart d = vecFaces[first] ; - Dart it = d ; - typename PFP::REAL minDist = 1e10 ; // ?? - do + Dart d = vecFaces[first]; + Dart it = d; + typename PFP::REAL minDist = (ip - position[it]).norm2(); + it = map.phi1(it); + while(it != d) { - typename PFP::VEC3 vec = ip - position[it]; - typename PFP::REAL dist = vec*vec ; + typename PFP::REAL dist = (ip - position[it]).norm2(); if(dist < minDist) { - minDist = dist ; - vertex = it ; + minDist = dist; + vertex = it; } - it = map.phi1(it) ; - } while(it != d) ; + it = map.phi1(it); + } } else - vertex = map.end() ; + vertex = map.end(); } /** - * Fonction that do the selection of darts, returned darts are sorted from closer to farther - * Dart is here considered as a small triangle from vertex to middle distance to center + * Fonction that do the selection of darts, returned darts are sorted from closest to farthest + * Dart is here considered as a triangle formed by the 2 end vertices of the edge and the face centroid * @param map the map we want to test * @param rayA first point of ray (user side) * @param rayAB vector of ray (directed ot the scene) @@ -246,9 +251,9 @@ void dartsRaySelection(typename PFP::MAP& map, const typename PFP::TVEC3& positi { Dart d = vecDarts[i]; distndart[i].second = d; - typename PFP::VEC3 V = (position[d] + position[map.phi1(d)])/2.0; + typename PFP::VEC3 V = (position[d] + position[map.phi1(d)]) / typename PFP::REAL(2); V -= rayA; - distndart[i].first = V*V; + distndart[i].first = V.norm2(); } // sort the vector of pair dist/dart diff --git a/include/Geometry/inclusion.h b/include/Geometry/inclusion.h index 4475fe0acbd102297cc1627232c98470a7a6de21..78094cff7b7a4ced0b2042b58e08549ae22a2afd 100644 --- a/include/Geometry/inclusion.h +++ b/include/Geometry/inclusion.h @@ -53,6 +53,9 @@ enum Inclusion template Inclusion isPointInTriangle(const VEC3& point, const VEC3& Ta, const VEC3& Tb, const VEC3& Tc) ; +template +Inclusion isPointInSphere(const VEC3& point, const VEC3& center, const typename VEC3::DATA_TYPE& radius) ; + /** * test if a segment is inside a triangle, the segment MUST be in the plane of the triangle * TODO to test @@ -93,9 +96,9 @@ bool isEdgeInOrIntersectingTetrahedron(VEC3 points[4], VEC3& point1, VEC3& point template bool arePointsEquals(const VEC3& point1,const VEC3& point2) ; -} +} // namespace Geom -} +} // namespace CGoGN #include "inclusion.hpp" diff --git a/include/Geometry/inclusion.hpp b/include/Geometry/inclusion.hpp index c87de80e23e99b3735ec790b1c480f46e4c77166..ed47a754a65e862ccc7eae31888fe65dea47ad76 100644 --- a/include/Geometry/inclusion.hpp +++ b/include/Geometry/inclusion.hpp @@ -81,6 +81,12 @@ Inclusion isPointInTriangle(const VEC3& point,const VEC3& Ta,const VEC3& Tb,cons return NO_INCLUSION ; } +template +bool isPointInSphere(const VEC3& point, const VEC3& center, const typename VEC3::DATA_TYPE& radius) +{ + return (point - center).norm() < radius; +} + template Inclusion isSegmentInTriangle2D(const VEC3& P1, const VEC3& P2, const VEC3& Ta, const VEC3& Tb, const VEC3& Tc, const VEC3& N) { diff --git a/include/Geometry/intersection.h b/include/Geometry/intersection.h index b50f7e4fbbd87b7221ec25b792d4f9e5359a17ae..8917bf802aaef03ca2ddb768d28fa31f95b249f7 100644 --- a/include/Geometry/intersection.h +++ b/include/Geometry/intersection.h @@ -43,6 +43,9 @@ enum Intersection FACE_INTERSECTION = 3 } ; +/** + * test the intersection between a ray and a triangle + */ template Intersection intersectionRayTriangle(const VEC3& P, const VEC3& Dir, const VEC3& Ta, const VEC3& Tb, const VEC3& Tc, VEC3& Inter) ; diff --git a/include/Geometry/intersection.hpp b/include/Geometry/intersection.hpp index 4dda3e7ad7749b594ab4b86aad94f0ccb7462fd8..b0b98b35f8c12285dc0443965c00347e8835bf7f 100644 --- a/include/Geometry/intersection.hpp +++ b/include/Geometry/intersection.hpp @@ -78,7 +78,7 @@ Intersection intersectionLineTriangle(const VEC3& P, const VEC3& Dir, const VEC3 #define PRECISION 1e-20 if(fabs(b) < PRECISION) //ray parallel to triangle - return NO_INTERSECTION ; + return NO_INTERSECTION ; #undef PRECISION T r = a / b ; @@ -282,35 +282,33 @@ Intersection intersectionSegmentTriangle(const VEC3& PA, const VEC3& PB, const V // } template -Intersection intersectionPlaneRay(const PLANE3D& pl,const VEC3& p1,const VEC3& dir, VEC3& Inter) +Intersection intersectionPlaneRay(const PLANE3D& pl, const VEC3& p1, const VEC3& dir, VEC3& Inter) { - typename VEC3::DATA_TYPE denom = pl.normal()*dir; + typename VEC3::DATA_TYPE denom = pl.normal() * dir; - if(denom==0) + if(denom == 0) { - if(pl.distance(p1)==0) + if(pl.distance(p1) == 0) { - Inter = p1; - return FACE_INTERSECTION; + Inter = p1; + return FACE_INTERSECTION; } else return NO_INTERSECTION; } - typename VEC3::DATA_TYPE isect = ( pl.normal() * (pl.normal()*-1.0f*pl.d()-p1) ) / denom; + typename VEC3::DATA_TYPE isect = ( pl.normal() * (pl.normal() * -1.0f * pl.d() - p1) ) / denom; Inter = p1 + dir * isect; if(0.0f <= isect) - { return FACE_INTERSECTION; - } return NO_INTERSECTION; } template -Intersection intersection2DSegmentSegment(const VEC3& PA, const VEC3& PB, const VEC3& QA, const VEC3& QB, VEC3& Inter) +Intersection intersection2DSegmentSegment(const VEC3& PA, const VEC3& PB, const VEC3& QA, const VEC3& QB, VEC3& Inter) { typedef typename VEC3::DATA_TYPE T ; @@ -320,21 +318,26 @@ Intersection intersection2DSegmentSegment(const VEC3& PA, const VEC3& PB, const vq1q2 -= QA; VEC3 vp1q1(QA); vp1q1 -= PA; - T delta = vp1p2[0]*vq1q2[1]- vp1p2[1]*vq1q2[0]; - T coeff = vp1q1[0]*vq1q2[1]- vp1q1[1]*vq1q2[0]; - Inter = VEC3((PA[0]*delta+vp1p2[0]*coeff)/delta,(PA[1]*delta+vp1p2[1]*coeff)/delta,(PA[2]*delta+vp1p2[2]*coeff)/delta); + T delta = vp1p2[0] * vq1q2[1] - vp1p2[1] * vq1q2[0]; + T coeff = vp1q1[0] * vq1q2[1] - vp1q1[1] * vq1q2[0]; + Inter = VEC3( + (PA[0] * delta + vp1p2[0] * coeff) / delta, + (PA[1] * delta + vp1p2[1] * coeff) / delta, + (PA[2] * delta + vp1p2[2] * coeff) / delta + ); //test if inter point is outside the edges - if( (Inter[0]PA[0] && Inter[0]>PB[0]) - || (Inter[0]QA[0] && Inter[0]>QB[0]) - || (Inter[1]PA[1] && Inter[1]>PB[1]) - || (Inter[1]QA[1] && Inter[1]>QB[1]) + if( + (Inter[0] < PA[0] && Inter[0] < PB[0]) || (Inter[0] > PA[0] && Inter[0] > PB[0]) || + (Inter[0] < QA[0] && Inter[0] < QB[0]) || (Inter[0] > QA[0] && Inter[0] > QB[0]) || + (Inter[1] < PA[1] && Inter[1] < PB[1]) || (Inter[1] > PA[1] && Inter[1] > PB[1]) || + (Inter[1] < QA[1] && Inter[1] < QB[1]) || (Inter[1] > QA[1] && Inter[1] > QB[1]) ) return NO_INTERSECTION; return EDGE_INTERSECTION; } -} +} // namespace Geom -} +} // namespace CGoGN diff --git a/include/Utils/qtSimple.h b/include/Utils/qtSimple.h index 7ad2c18068671f199517772ace7a799e54c16f21..de402af77eb445824b80fa34c160544342ce00a5 100644 --- a/include/Utils/qtSimple.h +++ b/include/Utils/qtSimple.h @@ -216,7 +216,6 @@ public: const glm::mat4& projectionMatrix() const { return m_projection_matrix; } glm::mat4& projectionMatrix() { return m_projection_matrix; } - float* curquat() { return m_curquat; } float* lastquat() { return m_lastquat; } diff --git a/src/Topology/generic/genericmap.cpp b/src/Topology/generic/genericmap.cpp index f6ef947ab4d38619a193be59ec3730c319a40b47..bc9eeb425d8bea1fb3163c3687b91e9820af1afa 100644 --- a/src/Topology/generic/genericmap.cpp +++ b/src/Topology/generic/genericmap.cpp @@ -550,8 +550,8 @@ bool GenericMap::foreach_dart_of_orbit(unsigned int orbit, Dart d, FunctorType& bool GenericMap::foreach_orbit(unsigned int orbit, FunctorType& fonct, const FunctorSelect& good, unsigned int thread) { - DartMarker marker(*this,thread); // Lock a marker - bool found = false; // Store the result + DartMarker marker(*this, thread); // Lock a marker + bool found = false; // Store the result // Scan all darts of the map for (Dart d = begin(); !found && d != end(); next(d))