Commit edc040c7 authored by untereiner's avatar untereiner

Merge cgogn.u-strasbg.fr:~jund/CGoGN

parents cd6c1705 e8dbdda7
......@@ -48,7 +48,7 @@ template <typename PFP>
void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position, SubdivideType sType = S_QUAD);
template <typename PFP>
void subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position);
Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position);
template <typename PFP>
void subdivideLoop(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position);
......
......@@ -168,11 +168,13 @@ void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position
}
template <typename PFP>
void subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
{
assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ;
assert(!map.volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided volume") ;
std::cout << "SUUUUUUB marine" << std::endl;
unsigned int vLevel = map.volumeLevel(d) ;
Dart old = map.volumeOldestDart(d) ;
......@@ -362,6 +364,8 @@ void subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi
}
map.setCurrentLevel(cur) ;
return subdividedfaces.begin()->first;
}
......
#ifndef PARTCELL_H
#define PARTCELL_H
#include "particle_base.h"
#include "Algo/Geometry/inclusion.h"
#include "Geometry/intersection.h"
#include "Geometry/orientation.h"
#include "Geometry/plane_3d.h"
#include <iostream>
/* A particle cell is a particle base within a map, within a precise cell, the displacement function should indicate
after each displacement wherein lies the new position of the particle */
namespace CGoGN
{
namespace Algo
{
namespace MovingObjects
{
template <typename PFP>
class ParticleCell3D : public ParticleBase
{
public :
typedef typename PFP::MAP Map;
typedef typename PFP::VEC3 VEC3;
typedef typename PFP::TVEC3 TAB_POS;
Map& m;
const TAB_POS& position;
Dart d;
Dart lastCrossed;
VEC3 m_positionFace;
unsigned int state;
unsigned int crossCell ;
ParticleCell3D(Map& map) : m(map)
{}
ParticleCell3D(Map& map, Dart belonging_cell, VEC3 pos, const TAB_POS& tabPos) :
ParticleBase(pos), m(map), position(tabPos), d(belonging_cell), state(3)
{
m_positionFace = pointInFace(d);
}
void display();
Dart getCell() { return d; }
VEC3 pointInFace(Dart d);
Geom::Orientation3D isLeftENextVertex(VEC3 c, Dart d, VEC3 base);
bool isRightVertex(VEC3 c, Dart d, VEC3 base);
Geom::Orientation3D whichSideOfFace(VEC3 c, Dart d);
Geom::Orientation3D isLeftL1DVol(VEC3 c, Dart d, VEC3 base, VEC3 top);
Geom::Orientation3D isRightDVol(VEC3 c, Dart d, VEC3 base, VEC3 top);
Geom::Orientation3D isAbove(VEC3 c, Dart d, VEC3 top);
int isLeftL1DFace(VEC3 c, Dart d, VEC3 base, VEC3 normal);
bool isRightDFace(VEC3 c, Dart d, VEC3 base, VEC3 normal);
Dart nextDartOfVertexNotMarked(Dart d, CellMarker& mark);
Dart nextNonPlanar(Dart d);
Dart nextFaceNotMarked(Dart d,CellMarker& mark);
Geom::Orientation3D whichSideOfEdge(VEC3 c, Dart d);
bool isOnHalfEdge(VEC3 c, Dart d);
void vertexState(const VEC3& current);
void edgeState(const VEC3& current);
void faceState(const VEC3& current, Geom::Orientation3D sideOfFace=Geom::ON);
void volumeState(const VEC3& current);
void volumeSpecialCase(const VEC3& current);
void move(const VEC3& newCurrent)
{
if(!Geom::arePointsEquals(newCurrent, m_position))
{
switch(state) {
case VERTEX_ORBIT : vertexState(newCurrent); break;
case EDGE_ORBIT : edgeState(newCurrent); break;
case FACE_ORBIT : faceState(newCurrent); break;
case VOLUME_ORBIT : volumeState(newCurrent); break;
}
display();
}
}
};
#include "particle_cell_3D.hpp"
}
}
}
#endif
// #define DEBUG
#define DELTA 0.00001
//static const float DELTA=0.00001;
template <typename PFP>
void ParticleCell3D<PFP>::display()
{
// std::cout << "position : " << m_position << std::endl;
}
template <typename PFP>
typename PFP::VEC3 ParticleCell3D<PFP>::pointInFace(Dart d)
{
return Algo::Geometry::faceCentroid<PFP>(m,d,position);
// const VEC3& p1(m_positions[d]);
// Dart dd=m.phi1(d);
// const VEC3& p2(m_positions[dd]);
// dd=m.phi1(dd);
// VEC3& p3(m_positions[dd]);
//
// while(Geom::testOrientation2D(p3,p1,p2)==Geom::ALIGNED) {
// dd = m.phi1(dd);
// p3 = m_positions[dd];
// }
//
// CGoGNout << "pointInFace " << (p1+p3)*0.5f << CGoGNendl;
//
// return (p1+p3)*0.5f;
}
template <typename PFP>
Geom::Orientation3D ParticleCell3D<PFP>::isLeftENextVertex(VEC3 c, Dart d, VEC3 base)
{
VEC3 p = position[d];
VEC3 p1 = position[m.phi1(m.phi2(d))];
Geom::Plane3D<typename PFP::REAL> pl(p,base,p1);
return pl.orient(c);
}
template <typename PFP>
bool ParticleCell3D<PFP>::isRightVertex(VEC3 c, Dart d, VEC3 base)
{
VEC3 p = position[d];
VEC3 p1 = position[m.phi_1(d)];
Geom::Plane3D<typename PFP::REAL> pl(p,base,p1);
return pl.orient(c)==Geom::UNDER;
}
template <typename PFP>
Geom::Orientation3D ParticleCell3D<PFP>::whichSideOfFace(VEC3 c, Dart d)
{
Geom::Plane3D<typename PFP::REAL> pl = Algo::Geometry::facePlane<PFP>(m,d,position);
return pl.orient(c);
}
template <typename PFP>
Geom::Orientation3D ParticleCell3D<PFP>::isLeftL1DVol(VEC3 c, Dart d, VEC3 base, VEC3 top)
{
VEC3 p2 = position[d];
Geom::Plane3D<typename PFP::REAL> pl(top,p2,base);
return pl.orient(c);
}
template <typename PFP>
Geom::Orientation3D ParticleCell3D<PFP>::isRightDVol(VEC3 c, Dart d, VEC3 base, VEC3 top)
{
VEC3 p1= position[m.phi1(d)];
Geom::Plane3D<typename PFP::REAL> pl(top,p1,base);
return pl.orient(c);
}
template <typename PFP>
Geom::Orientation3D ParticleCell3D<PFP>::isAbove(VEC3 c, Dart d, VEC3 top)
{
VEC3 p1 = position[d];
VEC3 p2 = position[m.phi1(d)];
Geom::Plane3D<typename PFP::REAL> pl(top,p2,p1);
return pl.orient(c);
}
template <typename PFP>
int ParticleCell3D<PFP>::isLeftL1DFace(VEC3 c, Dart d, VEC3 base, VEC3 normal)
{
VEC3 p2 = position[d];
VEC3 v2(p2-base);
VEC3 np = normal ^ v2;
Geom::Plane3D<typename PFP::REAL> pl(np,base*np);
return pl.orient(c);
}
template <typename PFP>
bool ParticleCell3D<PFP>::isRightDFace(VEC3 c, Dart d, VEC3 base, VEC3 normal)
{
VEC3 p1 = position[m.phi1(d)];
VEC3 np = normal ^ VEC3(p1-base);
Geom::Plane3D<typename PFP::REAL> pl(np,base*np);
return pl.orient(c)==Geom::UNDER;
}
template <typename PFP>
Dart ParticleCell3D<PFP>::nextDartOfVertexNotMarked(Dart d, CellMarker& mark)
{
// lock a marker
Dart d1;
DartMarkerNoUnmark markCC(m);
// init algo with parameter dart
std::list<Dart> darts_list;
darts_list.push_back(d);
markCC.mark(d);
// use iterator for begin of not yet treated darts
std::list<Dart>::iterator beg = darts_list.begin();
// until all darts treated
while (beg != darts_list.end())
{
d1 = *beg;
// add phi1, phi2 and phi3 successor if they are not yet marked
Dart d2 = m.phi1(m.phi2(d1));
Dart d3 = m.phi1(m.phi3(d1));
if (!markCC.isMarked(d2)) {
darts_list.push_back(d2);
markCC.mark(d2);
}
if (!markCC.isMarked(d3)) {
darts_list.push_back(d3);
markCC.mark(d3);
}
beg++;
// apply functor
if (!mark.isMarked(d1)) {
for (std::list<Dart>::iterator it=darts_list.begin(); it!=darts_list.end(); ++it)
markCC.unmark(*it);
return d1;
}
}
// clear markers
for (std::list<Dart>::iterator it=darts_list.begin(); it!=darts_list.end(); ++it)
markCC.unmark(*it);
return d;
}
template <typename PFP>
Dart ParticleCell3D<PFP>::nextNonPlanar(Dart d)
{
// lock a marker
Dart d1;
DartMarkerNoUnmark markCC(m);
// init algo with parameter dart
std::list<Dart> darts_list;
darts_list.push_back(d);
markCC.mark(d);
// use iterator for begin of not yet treated darts
std::list<Dart>::iterator beg = darts_list.begin();
// until all darts treated
while (beg != darts_list.end()) {
d1 = *beg;
// add phi1, phi2 and phi3 successor if they are not yet marked
Dart d2 = m.phi1(d1);
Dart d3 = m.phi2(d1);
Dart d4 = m.phi_1(d1);
if (!markCC.isMarked(d2)) {
darts_list.push_back(d2);
markCC.mark(d2);
}
if (!markCC.isMarked(d3)) {
darts_list.push_back(d3);
markCC.mark(d3);
}
if (!markCC.isMarked(d4)) {
darts_list.push_back(d4);
markCC.mark(d4);
}
beg++;
// apply functor
Geom::Plane3D<typename PFP::REAL> pl1 = Algo::Geometry::facePlane<PFP>(m,d,position);
Geom::Plane3D<typename PFP::REAL> pl2 = Algo::Geometry::facePlane<PFP>(m,d1,position);
if ((pl1.normal()-pl2.normal()).norm2()>0.000001)
beg = darts_list.end();
// remove the head of the list
}
// clear markers
for (std::list<Dart>::iterator it=darts_list.begin(); it!=darts_list.end(); ++it)
markCC.unmark(*it);
return d1;
}
template <typename PFP>
Dart ParticleCell3D<PFP>::nextFaceNotMarked(Dart d,CellMarker& mark)
{
// lock a marker
Dart d1;
DartMarkerNoUnmark markCC(m);
// init algo with parameter dart
std::list<Dart> darts_list;
darts_list.push_back(d);
markCC.mark(d);
// use iterator for begin of not yet treated darts
std::list<Dart>::iterator beg = darts_list.begin();
// until all darts treated
while (beg != darts_list.end())
{
d1 = *beg;
if(!mark.isMarked(d1))
{
break;
}
// add phi1, phi2 and phi3 successor if they are not yet marked
Dart d2 = m.phi1(d1);
Dart d3 = m.phi2(d1);
Dart d4 = m.phi_1(d1);
if (!markCC.isMarked(d2))
{
darts_list.push_back(d2);
markCC.mark(d2);
}
if (!markCC.isMarked(d3))
{
darts_list.push_back(d3);
markCC.mark(d3);
}
if (!markCC.isMarked(d4))
{
darts_list.push_back(d4);
markCC.mark(d4);
}
beg++;
}
// clear markers
for (std::list<Dart>::iterator it=darts_list.begin(); it!=darts_list.end(); ++it)
{
markCC.unmark(*it);
}
if(beg==darts_list.end())
return d;
return d1;
}
template <typename PFP>
Geom::Orientation3D ParticleCell3D<PFP>::whichSideOfEdge(VEC3 c, Dart d)
{
VEC3 p1 = position[m.phi1(d)];
VEC3 p2 = position[d];
Geom::Plane3D<typename PFP::REAL> pl = Algo::Geometry::facePlane<PFP>(m,d,position);
VEC3 norm = pl.normal();
VEC3 n2 = norm ^ VEC3(p2-p1);
Geom::Plane3D<typename PFP::REAL> pl2(n2,p1*n2);
return pl2.orient(c);
}
template <typename PFP>
bool ParticleCell3D<PFP>::isOnHalfEdge(VEC3 c, Dart d)
{
VEC3 p1 = position[d];
VEC3 p2 = position[m.phi1(d)];
VEC3 norm(p2-p1);
norm.normalize();
Geom::Plane3D<typename PFP::REAL> pl(norm,p1*norm);
return pl.orient(c)==Geom::OVER && !Geom::arePointsEquals(c,p1);
}
/**when the ParticleCell3D trajectory go through a vertex
* searching the good volume "umbrella" where the ParticleCell3D is
* if the ParticleCell3D is on the vertex, do nothing */
template <typename PFP>
void ParticleCell3D<PFP>::vertexState(const VEC3& current)
{
#ifdef DEBUG
std::cout << "vertexState" << d << std::endl;
#endif
VEC3 som = position[d];
if(Geom::arePointsEquals(current, m_position)) {
m_position = m_positionFace = som;
state = VERTEX_ORBIT;
return;
}
Dart dd=d;
Geom::Orientation3D wsof;
CellMarkerStore mark(m, FACE_CELL);
do {
VEC3 dualsp = (som+ Algo::Geometry::vertexNormal<PFP>(m,d,position));
Dart ddd=d;
mark.mark(d);
//searching the good orientation in a volume
if(isLeftENextVertex(current,d,dualsp)!=Geom::UNDER) {
d=m.phi1(m.phi2(d));
while(isLeftENextVertex(current,d,dualsp)!=Geom::UNDER && ddd!=d)
d=m.phi1(m.phi2(d));
if(ddd==d) {
if(isnan(current[0]) || isnan(current[1]) || isnan(current[2])) {
std::cout << __FILE__ << " " << __LINE__ << " NaN !" << std::endl;
}
bool verif = true;
do {
if(whichSideOfFace(current,d)!=Geom::OVER)
verif = false;
else
d=m.alpha1(d);
} while(verif && d!=ddd);
if(verif) {
volumeState(current);
return;
}
}
}
else {
while(isRightVertex(current,d,dualsp) && dd!=m.alpha_1(d))
d=m.phi2(m.phi_1(d));
}
wsof=whichSideOfFace(current,d);
//if c is before the vertex on the edge, we have to change of umbrella, we are symetric to the good umbrella
if(wsof!=Geom::OVER) {
VEC3 p1=position[d];
VEC3 p2=position[m.phi1(d)];
VEC3 norm(p2-p1);
norm.normalize();
Geom::Plane3D<typename PFP::REAL> plane(norm,p1*norm);
wsof = plane.orient(current);
}
//if c is on the other side of the face, we have to change umbrella
//if the umbrella has already been tested, we just take another arbitrary one
if(wsof==1) {
mark.mark(d);
if(!mark.isMarked(m.alpha1(d)))
d=m.alpha1(d);
else {
Dart dtmp=d;
d = nextDartOfVertexNotMarked(d,mark);
if(dtmp==d) {
std::cout << "numerical rounding ?" << std::endl;
d = dd;
m_position = pointInFace(d);
m_positionFace = m_position;
volumeState(current);
return;
}
}
}
}while(wsof==1);
if(wsof!=0)
{
m_position = pointInFace(d);
d = nextNonPlanar(d);
m_positionFace = pointInFace(d);
volumeState(current);
}
else
{
Dart ddd=d;
edgeState(current);
}
}
template <typename PFP>
void ParticleCell3D<PFP>::edgeState(const VEC3& current)
{
#ifdef DEBUG
std::cout << "edgeState" << d << " " << m_position << std::endl;
#endif
bool onEdge=false;
Dart dd=d;
Geom::Orientation3D wsof = whichSideOfFace(current,m.alpha2(d));
if(wsof!=Geom::UNDER) {
do {
d = m.alpha2(d);
wsof = whichSideOfFace(current,m.alpha2(d));
}while(wsof!=1 && dd!=d);
if(d==dd)
onEdge = true;
if(wsof==0) {
switch(whichSideOfEdge(current,d)) {
case Geom::ON :
onEdge=true;
break;
case Geom::UNDER :
d = m.phi2(d);
break;
default :
break;
}
}
wsof = whichSideOfFace(current,d);
}
else {
wsof = whichSideOfFace(current,d);
while(wsof==1 && dd != m.alpha_2(d)) {
d = m.alpha_2(d);
wsof = whichSideOfFace(current,d);
}
switch(whichSideOfEdge(current,d)) {