From 28c4843567a99117da0c38bb17ae3032ce4dda00 Mon Sep 17 00:00:00 2001 From: cignoni Date: Thu, 8 Jan 2009 11:26:49 +0000 Subject: [PATCH] Heavily restructured. Now more robust and allow also the creation of discretized surfaces --- vcg/complex/trimesh/create/resampler.h | 158 ++++++++++++++++++------- 1 file changed, 117 insertions(+), 41 deletions(-) diff --git a/vcg/complex/trimesh/create/resampler.h b/vcg/complex/trimesh/create/resampler.h index 622cddde..78f80a57 100644 --- a/vcg/complex/trimesh/create/resampler.h +++ b/vcg/complex/trimesh/create/resampler.h @@ -89,15 +89,9 @@ template *Vo,float in,const Box3i &bbox,vcg::Point3i &resolution) - {*/ - /* init=in; - Vol=Vo;*/ - - - + float max_dim; // the limit value of the search (that takes into account of the offset) + float offset; // an offset value that is always added to the returned value. Useful for extrarting isosurface at a different threshold + bool DiscretizeFlag; // if the extracted surface should be discretized or not. Walker(const Box3f &_bbox, Point3i _siz ) { this->bbox= _bbox; @@ -107,6 +101,7 @@ template siz.X()+1)*(this->siz.Z()+1); CurrentSlice = 0; offset=0; + DiscretizeFlag=false; _x_cs = new VertexIndex[ SliceSize ]; _y_cs = new VertexIndex[ SliceSize ]; @@ -125,10 +120,11 @@ template VV(int x,int y,int z) { assert ((y==CurrentSlice)||(y==(CurrentSlice+1))); @@ -138,14 +134,14 @@ template IPiToPf(Point3i(x,y,z),testPt); - vcg::Point3f closestNorm; + vcg::Point3f closestNormV,closestNormF; vcg::Point3f closestPt; - vcg::Point3f pip; + vcg::Point3f pip(-1,-1,-1); // Note that PointDistanceBaseFunctor does not require the edge and plane precomptued. // while the PointDistanceFunctor requires them. @@ -167,23 +163,38 @@ template V(0)->cN())*pip[0]+ (f->V(1)->cN())*pip[1] + (f->V(2)->cN())*pip[2] ; - assert(!f->IsD()); - Point3f dir=(testPt-closestPt); - /* dist=dir.Norm();*/ + bool retIP; + + // To compute the interpolated normal we use the more robust function that require to know what is the most orhogonal direction of the face. + if((*f).Flags() & Old_Mesh::FaceType::NORMX) retIP=InterpolationParameters(*f,0,closestPt, pip); + else if((*f).Flags() & Old_Mesh::FaceType::NORMY) retIP=InterpolationParameters(*f,1,closestPt, pip); + else if((*f).Flags() & Old_Mesh::FaceType::NORMZ) retIP=InterpolationParameters(*f,2,closestPt, pip); + else assert(0); + assert(retIP); // this should happen only if the starting mesh has degenerate faces. + + closestNormV = (f->V(0)->cN())*pip[0] + (f->V(1)->cN())*pip[1] + (f->V(2)->cN())*pip[2] ; + closestNormF = f->cN() ; - dir.Normalize(); - //direction of normal inside the mesh - if ((dir.dot(closestNorm))<0) - dist=-dist; - //the intersection exist - return true; + Point3f dir=(testPt-closestPt).Normalize(); + + // Compute test if the point see the surface normal from inside or outside + // Surface normal for improved robustness is computed both by face and interpolated from vertices. + float signV = dir.dot(closestNormV) ; + float signF = dir.dot(closestNormF) ; + + // Note that the two sings could be discordant. + // Always choose the best one according to the magnitude. + float signBest; + if(fabs(signV) > fabs(signF)) signBest = signV; + else signBest = signF; + + if(signBest<0) dist=-dist; + + return true; } - ///compute the values if an entire slice (per y) distances>dig of a cell are signed with double of + /// compute the values if an entire slice (per y) distances>dig of a cell are signed with double of /// the distance of the bb void ComputeSliceValues(int slice,field_value *slice_values) { @@ -200,11 +211,60 @@ template voxel[0],this->voxel[1]),this->voxel[2]); + int flippedCnt=0; + int flippedTot=0; + int flippedTimes=0; + do + { + flippedCnt=0; + for (int i=0; i<=this->siz.X(); i++) + { + for (int k=0; k<=this->siz.Z(); k++) + { + int goodCnt=0; + int badCnt=0; + int index=GetSliceIndex(i,k); + int index_l,index_r,index_u,index_d; + if(slice_values[index].first) + { + float curVal= slice_values[index].second; + if(i > 0 ) index_l=GetSliceIndex(i-1,k); else index_l = index; + if(i < this->siz.X() ) index_r=GetSliceIndex(i+1,k); else index_r = index; + if(k > 0 ) index_d=GetSliceIndex(i,k-1); else index_d = index; + if(k < this->siz.Z() ) index_u=GetSliceIndex(i,k+1); else index_u = index; + + if(slice_values[index_l].first) { goodCnt++; if(fabs(slice_values[index_l].second - curVal) > max_dist) badCnt++; } + if(slice_values[index_r].first) { goodCnt++; if(fabs(slice_values[index_r].second - curVal) > max_dist) badCnt++; } + if(slice_values[index_u].first) { goodCnt++; if(fabs(slice_values[index_u].second - curVal) > max_dist) badCnt++; } + if(slice_values[index_d].first) { goodCnt++; if(fabs(slice_values[index_d].second - curVal) > max_dist) badCnt++; } + + if(badCnt >= goodCnt) { + slice_values[index].second *=-1.0f; + //slice_values[index].first = false; + flippedCnt++; + } + } + } + } + flippedTot+=flippedCnt; + flippedTimes++; + } while(flippedCnt>0); + + if(flippedTot>0) + qDebug("Flipped %i values in %i times",flippedTot,flippedTimes); + } template void ProcessSlice(EXTRACTOR_TYPE &extractor) { @@ -212,9 +272,15 @@ template siz.Z(); k++) { + bool goodCell=true; Point3i p1(i,CurrentSlice,k); Point3i p2=p1+Point3i(1,1,1); - extractor.ProcessCell(p1, p2); + for(int ii=0;ii<2;++ii) + for(int jj=0;jj<2;++jj) + for(int kk=0;kk<2;++kk) + goodCell &= VV(p1[0]+ii,p1[1]+jj,p1[2]+kk).first; + + if(goodCell) extractor.ProcessCell(p1, p2); } } } @@ -227,24 +293,33 @@ template ::PerFaceNormalized(old_mesh); + tri::UpdateNormals::PerVertexNormalizedPerFaceNormalized(old_mesh); tri::UpdateFlags::FaceProjection(old_mesh); + int _size=(int)old_mesh.fn*100; - _g.Set(_oldM->face.begin(),_oldM->face.end()); + _g.Set(_oldM->face.begin(),_oldM->face.end(),_size); markerFunctor.SetMesh(&old_mesh); _newM->Clear(); Begin(); extractor.Initialize(); - + int computeTime =0; + int extractTime =0; + int t0,t1,t2; for (int j=0; j<=this->siz.Y(); j++) { cb((100*j)/this->siz.Y(),"Marching "); + t0 = clock(); ProcessSlice(extractor);//find cells where there is the isosurface and examine it + t1 = clock(); NextSlice(); + t2 = clock(); + computeTime += t1-t0; + extractTime += t2-t1; } extractor.Finalize(); + qDebug("Extract %i, Compute %i",t1-t0,t2-t1); typename New_Mesh::VertexIterator vi; for(vi=new_mesh.vert.begin();vi!=new_mesh.vert.end();++vi) @@ -499,7 +574,7 @@ typedef Walker /*< Old_Mesh,New_Mesh>*/ MyWalker; typedef vcg::tri::MarchingCubes MyMarchingCubes; ///resample the mesh using marching cube algorithm ,the accuracy is the dimension of one cell the parameter -static void Resample(Old_Mesh &old_mesh,New_Mesh &new_mesh,vcg::Point3 accuracy,float max_dist, float thr=0, vcg::CallBackPos *cb=0 ) +static void Resample(Old_Mesh &old_mesh,New_Mesh &new_mesh,vcg::Point3 accuracy,float max_dist, float thr=0, bool DiscretizeFlag=false, vcg::CallBackPos *cb=0 ) { ///be sure that the bounding box is updated vcg::tri::UpdateBounding::Box(old_mesh); @@ -510,6 +585,7 @@ static void Resample(Old_Mesh &old_mesh,New_Mesh &new_mesh,vcg::Point3 accu walker.max_dim=max_dist+fabs(thr); walker.offset = - thr; + walker.DiscretizeFlag = DiscretizeFlag; MyMarchingCubes mc(new_mesh, walker); walker.BuildMesh(old_mesh,new_mesh,mc,cb); }