manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

updated voronoi processing stuff to manage float/double and to use the new kdtree

This commit is contained in:
Paolo Cignoni 2014-08-23 01:32:50 +00:00
parent 7285fadd53
commit e032901b7f
2 changed files with 95 additions and 95 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
vcg/complex/algorithms/voronoi_processing.h
View File

@ -154,7 +154,7 @@ static void ComputePerVertexSources(MeshType &m, std::vector<VertexType *> &seed
PerVertexPointerHandle vertexSources = tri::Allocator<MeshType>:: template AddPerVertexAttribute<VertexPointer> (m,"sources"); PerVertexPointerHandle vertexSources = tri::Allocator<MeshType>:: template AddPerVertexAttribute<VertexPointer> (m,"sources");
tri::Allocator<MeshType>::DeletePerFaceAttribute(m,"sources"); // delete any conflicting handle regardless of the type... tri::Allocator<MeshType>::DeletePerFaceAttribute(m,"sources"); // delete any conflicting handle regardless of the type...
PerFacePointerHandle faceSources = tri::Allocator<MeshType>:: template AddPerFaceAttribute<VertexPointer> (m,"sources"); tri::Allocator<MeshType>::template AddPerFaceAttribute<VertexPointer> (m,"sources");
assert(tri::Allocator<MeshType>::IsValidHandle(m,vertexSources)); assert(tri::Allocator<MeshType>::IsValidHandle(m,vertexSources));

188
vcg/complex/algorithms/voronoi_volume_sampling.h
View File

@ -36,41 +36,45 @@ class VoronoiVolumeSampling
{ {
public: public:
typedef typename tri::VoronoiProcessing<MeshType>::QuadricSumDistance QuadricSumDistance; typedef typename tri::VoronoiProcessing<MeshType>::QuadricSumDistance QuadricSumDistance;
typedef SimpleVolume<SimpleVoxel> MyVolume; typedef typename MeshType::ScalarType ScalarType;
typedef typename vcg::tri::TrivialWalker<MeshType,MyVolume> MyWalker; typedef typename MeshType::BoxType BoxType;
typedef typename vcg::tri::MarchingCubes<MeshType, MyWalker> MyMarchingCubes;
typedef typename vcg::GridStaticPtr<typename MeshType::FaceType> GridType;
typedef typename MeshType::VertexIterator VertexIterator; typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::VertexPointer VertexPointer; typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::CoordType CoordType; typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::FacePointer FacePointer; typedef typename MeshType::FacePointer FacePointer;
typedef typename vcg::GridStaticPtr<typename MeshType::FaceType, ScalarType> GridType;
typedef SimpleVolume<SimpleVoxel<ScalarType> > MyVolume;
typedef typename vcg::tri::TrivialWalker<MeshType,MyVolume> MyWalker;
typedef typename vcg::tri::MarchingCubes<MeshType, MyWalker> MyMarchingCubes;
VoronoiVolumeSampling(MeshType &_baseMesh, MeshType &_seedMesh) VoronoiVolumeSampling(MeshType &_baseMesh, MeshType &_seedMesh)
:baseMesh(_baseMesh),seedMesh(_seedMesh),seedTree(0),surfTree(0) :seedTree(0),surfTree(0),baseMesh(_baseMesh),seedMesh(_seedMesh)
{ {
} }
KdTree<float> *seedTree; KdTree<ScalarType> *seedTree;
KdTree<float> *surfTree; KdTree<ScalarType> *surfTree;
typename KdTree<ScalarType>::PriorityQueue pq;
GridType surfGrid; GridType surfGrid;
typedef FaceTmark<MeshType> MarkerFace; typedef FaceTmark<MeshType> MarkerFace;
MarkerFace mf; MarkerFace mf;
vcg::face::PointDistanceBaseFunctor<float> PDistFunct; vcg::face::PointDistanceBaseFunctor<ScalarType> PDistFunct;
MeshType &baseMesh; MeshType &baseMesh;
MeshType &seedMesh; MeshType &seedMesh;
MeshType poissonSurfaceMesh; MeshType poissonSurfaceMesh;
float poissonRadiusSurface; ScalarType poissonRadiusSurface;
MeshType montecarloVolumeMesh; MeshType montecarloVolumeMesh;
void Init(float radius=0) void Init(ScalarType radius=0)
{ {
MeshType montecarloSurfaceMesh; MeshType montecarloSurfaceMesh;
if(radius==0) poissonRadiusSurface = baseMesh.bbox.Diag()/50.0f; if(radius==0) poissonRadiusSurface = baseMesh.bbox.Diag()/50.0f;
else poissonRadiusSurface = radius; else poissonRadiusSurface = radius;
float meshArea = Stat<MeshType>::ComputeMeshArea(baseMesh); ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(baseMesh);
int MontecarloSampleNum = 10 * meshArea / (radius*radius); int MontecarloSampleNum = 10 * meshArea / (radius*radius);
tri::MeshSampler<MeshType> sampler(montecarloSurfaceMesh); tri::MeshSampler<MeshType> sampler(montecarloSurfaceMesh);
tri::SurfaceSampling<MeshType,tri::MeshSampler<CMeshO> >::Montecarlo(baseMesh, sampler, MontecarloSampleNum); tri::SurfaceSampling<MeshType,tri::MeshSampler<CMeshO> >::Montecarlo(baseMesh, sampler, MontecarloSampleNum);
@ -85,61 +89,58 @@ public:
qDebug("Surface Sampling radius %f - montecarlo %ivn - Poisson %ivn",poissonRadiusSurface,montecarloSurfaceMesh.vn,poissonSurfaceMesh.vn); qDebug("Surface Sampling radius %f - montecarlo %ivn - Poisson %ivn",poissonRadiusSurface,montecarloSurfaceMesh.vn,poissonSurfaceMesh.vn);
VertexConstDataWrapper<MeshType> ww(poissonSurfaceMesh); VertexConstDataWrapper<MeshType> ww(poissonSurfaceMesh);
if(surfTree) delete surfTree; if(surfTree) delete surfTree;
surfTree = new KdTree<float>(ww); surfTree = new KdTree<ScalarType>(ww);
surfTree->setMaxNofNeighbors(1);
surfGrid.SetWithRadius(baseMesh.face.begin(),baseMesh.face.end(),poissonRadiusSurface); surfGrid.SetWithRadius(baseMesh.face.begin(),baseMesh.face.end(),poissonRadiusSurface);
mf.SetMesh(&baseMesh); mf.SetMesh(&baseMesh);
} }
// Compute the signed distance from the surface // Compute the signed distance from the surface
float DistanceFromSurface(Point3f &p) ScalarType DistanceFromSurface(CoordType &p)
{ {
surfTree->doQueryK(p); ScalarType squaredDist;
float dist = sqrtf(surfTree->getNeighborSquaredDistance(0)); unsigned int ind;
surfTree->doQueryClosest(p,ind,squaredDist);
ScalarType dist = sqrt(squaredDist);
if( dist > 3.0f*poissonRadiusSurface) if( dist > 3.0f*poissonRadiusSurface)
{ {
Point3f dir = surfTree->getNeighbor(0) - p; // CoordType dir = surfTree->getNeighbor(0) - p;
const Point3f &surfN = this->poissonSurfaceMesh.vert[surfTree->getNeighborId(0)].N(); CoordType dir = this->poissonSurfaceMesh.vert[ind].P() - p;
const CoordType &surfN = this->poissonSurfaceMesh.vert[ind].N();
if(dir* surfN > 0) dist= -dist; if(dir* surfN > 0) dist= -dist;
return dist; return dist;
} }
float _maxDist = this->poissonRadiusSurface*3.0f; ScalarType _maxDist = this->poissonRadiusSurface*3.0f;
dist=_maxDist; dist=_maxDist;
Point3f _closestPt; CoordType _closestPt;
FacePointer f=surfGrid.GetClosest(PDistFunct,mf,p,_maxDist,dist,_closestPt); FacePointer f=surfGrid.GetClosest(PDistFunct,mf,p,_maxDist,dist,_closestPt);
assert(f); assert(f);
assert (dist >=0); assert (dist >=0);
Point3f dir = _closestPt - p; CoordType dir = _closestPt - p;
if(dir*f->cN() > 0) dist = -dist; if(dir*f->cN() > 0) dist = -dist;
return dist; return dist;
} }
float DistanceFromVoronoiSeed(Point3f p_point) ScalarType DistanceFromVoronoiSeed(CoordType p_point)
{ {
// Calculating the closest point to p_point ScalarType squaredDist;
seedTree->doQueryK(p_point); unsigned int ind;
float minD = seedTree->getNeighborSquaredDistance(0); surfTree->doQueryClosest(p_point,ind,squaredDist);
for(int i=1;i<seedTree->getNofFoundNeighbors();++i) return math::Sqrt(squaredDist);
minD=std::min(minD,seedTree->getNeighborSquaredDistance(i));
return sqrtf(minD);
} }
float DistanceFromVoronoiFace(Point3f p_point) ScalarType DistanceFromVoronoiFace(CoordType p_point)
{ {
seedTree->doQueryK(p_point);
std::vector<std::pair<float, Point3f> > closeSeedVec; seedTree->doQueryK(p_point,2,pq);
for(int i=0;i<seedTree->getNofFoundNeighbors();++i)
closeSeedVec.push_back(std::make_pair(seedTree->getNeighborSquaredDistance(i),seedTree->getNeighbor(i)));
std::sort(closeSeedVec.begin(),closeSeedVec.end()); std::vector<std::pair<ScalarType, CoordType> > closeSeedVec;
Point3f p0=closeSeedVec[0].second; CoordType p0= this->seedMesh.vert[pq.getIndex(0)].P();
Point3f p1=closeSeedVec[1].second; CoordType p1= this->seedMesh.vert[pq.getIndex(1)].P();
Plane3f pl; pl.Init((p0+p1)/2.0f,p0-p1); Plane3<ScalarType> pl; pl.Init((p0+p1)/2.0f,p0-p1);
return fabs(SignedDistancePlanePoint(pl,p_point)); return fabs(SignedDistancePlanePoint(pl,p_point));
} }
@ -156,29 +157,24 @@ float DistanceFromVoronoiFace(Point3f p_point)
* returns: distance between the point P and the line R * returns: distance between the point P and the line R
*/ */
float DistanceFromVoronoiEdge(Point3f p_point) ScalarType DistanceFromVoronoiEdge(CoordType p_point)
{ {
seedTree->doQueryK(p_point); seedTree->doQueryK(p_point,3,pq);
std::vector<std::pair<float, Point3f> > closeSeedVec; std::vector<std::pair<ScalarType, CoordType> > closeSeedVec;
for(int i=0;i<seedTree->getNofFoundNeighbors();++i) CoordType p0= this->seedMesh.vert[pq.getIndex(0)].P();
closeSeedVec.push_back(std::make_pair(seedTree->getNeighborSquaredDistance(i),seedTree->getNeighbor(i))); CoordType p1= this->seedMesh.vert[pq.getIndex(1)].P();
CoordType p2= this->seedMesh.vert[pq.getIndex(2)].P();
std::sort(closeSeedVec.begin(),closeSeedVec.end()); Plane3<ScalarType> pl01; pl01.Init((p0+p1)/2.0f,p0-p1);
Point3f p0=closeSeedVec[0].second; Plane3<ScalarType> pl02; pl02.Init((p0+p2)/2.0f,p0-p2);
Point3f p1=closeSeedVec[1].second; Line3<ScalarType> voroLine;
Point3f p2=closeSeedVec[2].second;
// Calculating the line R that intersect the planes pl01 and pl02
Plane3f pl01; pl01.Init((p0+p1)/2.0f,p0-p1);
Plane3f pl02; pl02.Init((p0+p2)/2.0f,p0-p2);
Line3f voroLine;
// Calculating the line R that intersect the planes po1 and p02
vcg::IntersectionPlanePlane(pl01,pl02,voroLine); vcg::IntersectionPlanePlane(pl01,pl02,voroLine);
// Calculating the distance k between the point p_point and the line R. // Calculating the distance k between the point p_point and the line R.
Point3f closestPt; CoordType closestPt;
float closestDist; ScalarType closestDist;
vcg::LinePointDistance(voroLine,p_point,closestPt, closestDist); vcg::LinePointDistance(voroLine,p_point,closestPt, closestDist);
return closestDist; return closestDist;
@ -192,43 +188,47 @@ float DistanceFromVoronoiFace(Point3f p_point)
int i; int i;
for(i=0;i<relaxStep;++i) for(i=0;i<relaxStep;++i)
{ {
seedTree->setMaxNofNeighbors(1);
QuadricSumDistance dz; QuadricSumDistance dz;
std::vector<QuadricSumDistance> dVec(montecarloVolumeMesh.vert.size(),dz); std::vector<QuadricSumDistance> dVec(montecarloVolumeMesh.vert.size(),dz);
for(typename MeshType::VertexIterator vi=montecarloVolumeMesh.vert.begin();vi!=montecarloVolumeMesh.vert.end();++vi) for(typename MeshType::VertexIterator vi=montecarloVolumeMesh.vert.begin();vi!=montecarloVolumeMesh.vert.end();++vi)
{ {
seedTree->doQueryK(vi->P()); unsigned int seedInd;
int seedIndex = seedTree->getNeighborId(0); ScalarType sqdist;
dVec[seedIndex].AddPoint(vi->P()); seedTree->doQueryClosest(vi->P(),seedInd,sqdist);
dVec[seedInd].AddPoint(vi->P());
} }
// Search the local maxima for each region and use them as new seeds // Search the local maxima for each region and use them as new seeds
std::vector< std::pair<float,int> > seedMaximaVec(seedMesh.vert.size(),std::make_pair(std::numeric_limits<float>::max(),-1 )); std::vector< std::pair<ScalarType,int> > seedMaximaVec(seedMesh.vert.size(),std::make_pair(std::numeric_limits<ScalarType>::max(),-1 ));
for(typename MeshType::VertexIterator vi=montecarloVolumeMesh.vert.begin();vi!=montecarloVolumeMesh.vert.end();++vi) for(typename MeshType::VertexIterator vi=montecarloVolumeMesh.vert.begin();vi!=montecarloVolumeMesh.vert.end();++vi)
{ {
seedTree->doQueryK(vi->P()); unsigned int seedInd;
int seedIndex = seedTree->getNeighborId(0); ScalarType sqdist;
float val = dVec[seedIndex].Eval(vi->P()); seedTree->doQueryClosest(vi->P(),seedInd,sqdist);
if(val < seedMaximaVec[seedIndex].first)
ScalarType val = dVec[seedInd].Eval(vi->P());
if(val < seedMaximaVec[seedInd].first)
{ {
seedMaximaVec[seedIndex].first = val; seedMaximaVec[seedInd].first = val;
seedMaximaVec[seedIndex].second = tri::Index(montecarloVolumeMesh,*vi); seedMaximaVec[seedInd].second = tri::Index(montecarloVolumeMesh,*vi);
} }
} }
changed=false; changed=false;
for(int i=0;i<seedMesh.vert.size();++i) for(int i=0;i<seedMesh.vert.size();++i)
{ {
Point3f prevP = seedMesh.vert[i].P() ; CoordType prevP = seedMesh.vert[i].P() ;
if(seedMaximaVec[i].second == -1) tri::Allocator<MeshType>::DeleteVertex(seedMesh,seedMesh.vert[i]);
seedMesh.vert[i].P() = montecarloVolumeMesh.vert[seedMaximaVec[i].second].P(); seedMesh.vert[i].P() = montecarloVolumeMesh.vert[seedMaximaVec[i].second].P();
if(prevP != seedMesh.vert[i].P()) changed = true; if(prevP != seedMesh.vert[i].P()) changed = true;
} }
tri::Allocator<MeshType>::CompactVertexVector(seedMesh);
// Kdtree must be rebuilt at the end of each step; // Kdtree for the seeds must be rebuilt at the end of each step;
VertexConstDataWrapper<MeshType> vdw(seedMesh); VertexConstDataWrapper<MeshType> vdw(seedMesh);
delete seedTree; delete seedTree;
seedTree = new KdTree<float>(vdw); seedTree = new KdTree<ScalarType>(vdw);
if(!changed) if(!changed)
break; break;
} }
@ -247,38 +247,36 @@ float DistanceFromVoronoiFace(Point3f p_point)
* PruningPoisson: mesh of inside and surface points, it's the voronoi3d diagram * PruningPoisson: mesh of inside and surface points, it's the voronoi3d diagram
* n_voxel: number of voxels for the greater side * n_voxel: number of voxels for the greater side
*/ */
void BuildScaffoldingMesh(MeshType &scaffoldingMesh, int volumeSide, float isoThr,int elemEnum, bool surfFlag) void BuildScaffoldingMesh(MeshType &scaffoldingMesh, int volumeSide, ScalarType isoThr,int elemEnum, bool surfFlag)
{ {
printf("Scaffolding of the mesh \n"); printf("Scaffolding of the mesh \n");
MyVolume volume; MyVolume volume;
float max = math::Max(baseMesh.bbox.DimX(),baseMesh.bbox.DimY(),baseMesh.bbox.DimZ()); ScalarType max = math::Max(baseMesh.bbox.DimX(),baseMesh.bbox.DimY(),baseMesh.bbox.DimZ());
float voxel = max / volumeSide; ScalarType voxel = max / volumeSide;
int sizeX = (baseMesh.bbox.DimX() / voxel)+1; int sizeX = (baseMesh.bbox.DimX() / voxel)+1;
int sizeY = (baseMesh.bbox.DimY() / voxel)+1; int sizeY = (baseMesh.bbox.DimY() / voxel)+1;
int sizeZ = (baseMesh.bbox.DimZ() / voxel)+1; int sizeZ = (baseMesh.bbox.DimZ() / voxel)+1;
// Kdtree // Kdtree
seedTree->setMaxNofNeighbors(4); // seedTree->setMaxNofNeighbors(4);
volume.bbox=baseMesh.bbox; BoxType bb = BoxType::Construct(baseMesh.bbox);
volume.bbox.Offset(baseMesh.bbox.Diag()*0.04f); bb.Offset(baseMesh.bbox.Diag()*0.04f);
volume.siz = Point3i(sizeX,sizeY,sizeZ); volume.Init(Point3i(sizeX,sizeY,sizeZ),bb);
volume.ComputeDimAndVoxel();
volume.Init(Point3i(sizeX,sizeY,sizeZ));
qDebug("Init Volume of %i %i %i",sizeX,sizeY,sizeZ); qDebug("Init Volume of %i %i %i",sizeX,sizeY,sizeZ);
int cnt=0; int cnt=0;
float offset= volume.voxel.Norm()*isoThr; ScalarType offset= volume.voxel.Norm()*isoThr;
for(float i=0;i<sizeX;i++) for(ScalarType i=0;i<sizeX;i++)
for(float j=0;j<sizeY;j++) for(ScalarType j=0;j<sizeY;j++)
for(float k=0;k<sizeZ;k++) for(ScalarType k=0;k<sizeZ;k++)
{ {
// check if the point is inside the mesh // check if the point is inside the mesh
Point3f p; CoordType p;
volume.IPiToPf(Point3i(i,j,k),p); volume.IPiToPf(Point3i(i,j,k),p);
float surfDist = this->DistanceFromSurface(p); ScalarType surfDist = this->DistanceFromSurface(p);
float elemDist; ScalarType elemDist;
switch(elemEnum) switch(elemEnum)
{ {
case 0: elemDist = DistanceFromVoronoiSeed(p) - offset; break; case 0: elemDist = DistanceFromVoronoiSeed(p) - offset; break;
@ -287,7 +285,7 @@ float DistanceFromVoronoiFace(Point3f p_point)
default: assert(0); default: assert(0);
} }
float val; ScalarType val;
if(surfFlag) if(surfFlag)
val = std::max(-elemDist,surfDist); val = std::max(-elemDist,surfDist);
else else
@ -311,13 +309,15 @@ float DistanceFromVoronoiFace(Point3f p_point)
*/ */
void ThicknessEvaluator() void ThicknessEvaluator()
{ {
surfTree->setMaxNofNeighbors(1); // surfTree->setMaxNofNeighbors(1);
tri::UpdateQuality<MeshType>::VertexConstant(poissonSurfaceMesh,0); tri::UpdateQuality<MeshType>::VertexConstant(poissonSurfaceMesh,0);
for(VertexIterator vi=montecarloVolumeMesh.vert.begin(); vi!=montecarloVolumeMesh.vert.end(); ++vi) for(VertexIterator vi=montecarloVolumeMesh.vert.begin(); vi!=montecarloVolumeMesh.vert.end(); ++vi)
{ {
this->surfTree->doQueryK(vi->P()); unsigned int ind;
VertexPointer vp = &poissonSurfaceMesh.vert[surfTree->getNeighborId(0)]; ScalarType sqdist;
float dist = sqrt(surfTree->getNeighborSquaredDistance(0)); this->surfTree->doQueryClosest(vi->P(),ind,sqdist);
VertexPointer vp = &poissonSurfaceMesh.vert[ind];
ScalarType dist = math::Sqrt(sqdist);
if(vp->Q() < dist) vp->Q()=dist; if(vp->Q() < dist) vp->Q()=dist;
} }
tri::UpdateColor<MeshType>::PerVertexQualityRamp(poissonSurfaceMesh); tri::UpdateColor<MeshType>::PerVertexQualityRamp(poissonSurfaceMesh);
@ -331,16 +331,16 @@ float DistanceFromVoronoiFace(Point3f p_point)
* Build a Poisson-Disk Point cloud that cover all the space of the original mesh m * Build a Poisson-Disk Point cloud that cover all the space of the original mesh m
* *
*/ */
void BuildVolumeSampling(int montecarloSampleNum, int seedNum, float &poissonRadius, vcg::CallBackPos *cb=0) void BuildVolumeSampling(int montecarloSampleNum, int seedNum, ScalarType &poissonRadius, vcg::CallBackPos *cb=0)
{ {
montecarloVolumeMesh.Clear(); montecarloVolumeMesh.Clear();
math::SubtractiveRingRNG rng; math::SubtractiveRingRNG rng;
surfTree->setMaxNofNeighbors(1); // surfTree->setMaxNofNeighbors(1);
while(montecarloVolumeMesh.vn < montecarloSampleNum) while(montecarloVolumeMesh.vn < montecarloSampleNum)
{ {
Point3f point = math::GeneratePointInBox3Uniform(rng,baseMesh.bbox); CoordType point = math::GeneratePointInBox3Uniform(rng,baseMesh.bbox);
float d = this->DistanceFromSurface(point); ScalarType d = this->DistanceFromSurface(point);
if(d<0){ if(d<0){
vcg::tri::Allocator<MeshType>::AddVertex(montecarloVolumeMesh,point); vcg::tri::Allocator<MeshType>::AddVertex(montecarloVolumeMesh,point);
montecarloVolumeMesh.vert.back().Q() = fabs(d); montecarloVolumeMesh.vert.back().Q() = fabs(d);
@ -363,7 +363,7 @@ float DistanceFromVoronoiFace(Point3f p_point)
// Kdtree must be rebuilt at the end of each step; // Kdtree must be rebuilt at the end of each step;
VertexConstDataWrapper<MeshType> vdw(seedMesh); VertexConstDataWrapper<MeshType> vdw(seedMesh);
if(seedTree) delete seedTree; if(seedTree) delete seedTree;
seedTree = new KdTree<float>(vdw); seedTree = new KdTree<ScalarType>(vdw);
} }
}; // end class }; // end class