00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifndef __VCG_MESH_RESAMPLER
00024 #define __VCG_MESH_RESAMPLER
00025
00026 #include <vcg/complex/trimesh/update/normal.h>
00027 #include <vcg/complex/trimesh/update/flag.h>
00028 #include <vcg/complex/trimesh/update/bounding.h>
00029 #include <vcg/complex/trimesh/update/edges.h>
00030 #include <vcg/complex/trimesh/create/marching_cubes.h>
00031 #include <vcg/space/index/grid_static_ptr.h>
00032 #include <vcg/complex/trimesh/closest.h>
00033 #include <vcg/space/box3.h>
00034
00035 namespace vcg {
00036 namespace tri {
00037
00038
00041
00048 template <class OLD_MESH_TYPE,class NEW_MESH_TYPE, class FLT, class DISTFUNCTOR = vcg::face::PointDistanceBaseFunctor<typename OLD_MESH_TYPE::ScalarType > >
00049 class Resampler : public BasicGrid<FLT>
00050 {
00051 typedef OLD_MESH_TYPE Old_Mesh;
00052 typedef NEW_MESH_TYPE New_Mesh;
00053
00054
00055 class Walker : BasicGrid<float>
00056 {
00057 private:
00058 typedef int VertexIndex;
00059 typedef OLD_MESH_TYPE Old_Mesh;
00060 typedef NEW_MESH_TYPE New_Mesh;
00061 typedef typename New_Mesh::CoordType NewCoordType;
00062 typedef typename New_Mesh::VertexType* VertexPointer;
00063 typedef typename Old_Mesh::FaceContainer FaceCont;
00064 typedef typename vcg::GridStaticPtr<typename Old_Mesh::FaceType> GridType;
00065
00066 protected:
00067
00068 int SliceSize;
00069 int CurrentSlice;
00070 typedef tri::FaceTmark<Old_Mesh> MarkerFace;
00071 MarkerFace markerFunctor;
00072
00073
00074 VertexIndex *_x_cs;
00075 VertexIndex *_y_cs;
00076 VertexIndex *_z_cs;
00077 VertexIndex *_x_ns;
00078 VertexIndex *_z_ns;
00079
00080
00081
00082
00083 typedef typename std::pair<bool,float> field_value;
00084 field_value* _v_cs;
00085 field_value* _v_ns;
00086
00087 New_Mesh *_newM;
00088 Old_Mesh *_oldM;
00089 GridType _g;
00090
00091 public:
00092 float max_dim;
00093 float offset;
00094 bool DiscretizeFlag;
00095 bool MultiSampleFlag;
00096 bool AbsDistFlag;
00097 Walker(const Box3f &_bbox, Point3i _siz )
00098 {
00099 this->bbox= _bbox;
00100 this->siz=_siz;
00101 ComputeDimAndVoxel();
00102
00103 SliceSize = (this->siz.X()+1)*(this->siz.Z()+1);
00104 CurrentSlice = 0;
00105 offset=0;
00106 DiscretizeFlag=false;
00107 MultiSampleFlag=false;
00108 AbsDistFlag=false;
00109
00110 _x_cs = new VertexIndex[ SliceSize ];
00111 _y_cs = new VertexIndex[ SliceSize ];
00112 _z_cs = new VertexIndex[ SliceSize ];
00113 _x_ns = new VertexIndex[ SliceSize ];
00114 _z_ns = new VertexIndex[ SliceSize ];
00115
00116 _v_cs= new field_value[(this->siz.X()+1)*(this->siz.Z()+1)];
00117 _v_ns= new field_value[(this->siz.X()+1)*(this->siz.Z()+1)];
00118
00119 };
00120
00121 ~Walker()
00122 {}
00123
00124
00125 float V(const Point3i &p)
00126 {
00127 return V(p.V(0),p.V(1),p.V(2));
00128 }
00129
00130
00131 std::pair<bool,float> VV(int x,int y,int z)
00132 {
00133 assert ((y==CurrentSlice)||(y==(CurrentSlice+1)));
00134
00135
00136
00137
00138
00139 int index=GetSliceIndex(x,z);
00140
00141 if (y==CurrentSlice) return _v_cs[index];
00142 else return _v_ns[index];
00143 }
00144
00145 float V(int x,int y,int z)
00146 {
00147 if(DiscretizeFlag) return VV(x,y,z).second+offset<0?-1:1;
00148 return VV(x,y,z).second+offset;
00149 }
00151 field_value DistanceFromMesh(Point3f &pp,Old_Mesh *)
00152 {
00153 float dist;
00154 typename Old_Mesh::FaceType *f=NULL;
00155 const float max_dist = max_dim;
00156 vcg::Point3f testPt;
00157 this->IPfToPf(pp,testPt);
00158
00159 vcg::Point3f closestNormV,closestNormF;
00160 vcg::Point3f closestPt;
00161 vcg::Point3f pip(-1,-1,-1);
00162
00163
00164
00165
00166 DISTFUNCTOR PDistFunct;
00167 f = _g.GetClosest(PDistFunct,markerFunctor,testPt,max_dist,dist,closestPt);
00168 if (f==NULL) return field_value(false,0);
00169 if(AbsDistFlag) return field_value(true,dist);
00170 assert(!f->IsD());
00171 bool retIP;
00172
00173
00174 if((*f).Flags() & Old_Mesh::FaceType::NORMX) retIP=InterpolationParameters(*f,0,closestPt, pip);
00175 else if((*f).Flags() & Old_Mesh::FaceType::NORMY) retIP=InterpolationParameters(*f,1,closestPt, pip);
00176 else if((*f).Flags() & Old_Mesh::FaceType::NORMZ) retIP=InterpolationParameters(*f,2,closestPt, pip);
00177 else assert(0);
00178 assert(retIP);
00179
00180 const float InterpolationEpsilon = 0.00001f;
00181 int zeroCnt=0;
00182 if(pip[0]<InterpolationEpsilon) ++zeroCnt;
00183 if(pip[1]<InterpolationEpsilon) ++zeroCnt;
00184 if(pip[2]<InterpolationEpsilon) ++zeroCnt;
00185 assert(zeroCnt<3);
00186
00187 Point3f dir=(testPt-closestPt).Normalize();
00188
00189
00190
00191 float signBest;
00192
00193
00194
00195 if(zeroCnt>0)
00196 {
00197 closestNormV = (f->V(0)->cN())*pip[0] + (f->V(1)->cN())*pip[1] + (f->V(2)->cN())*pip[2] ;
00198 signBest = dir.dot(closestNormV) ;
00199 }
00200 else
00201 {
00202 closestNormF = f->cN() ;
00203 signBest = dir.dot(closestNormF) ;
00204 }
00205
00206 if(signBest<0) dist=-dist;
00207
00208 return field_value(true,dist);
00209 }
00210
00211 field_value MultiDistanceFromMesh(Point3f &pp, Old_Mesh *)
00212 {
00213 float distSum=0;
00214 int positiveCnt=0;
00215 const int MultiSample=7;
00216 const Point3f delta[7]={Point3f(0,0,0),
00217 Point3f( 0.2, -0.01, -0.02),
00218 Point3f(-0.2, 0.01, 0.02),
00219 Point3f( 0.01, 0.2, 0.01),
00220 Point3f( 0.03, -0.2, -0.03),
00221 Point3f(-0.02, -0.03, 0.2 ),
00222 Point3f(-0.01, 0.01, -0.2 )};
00223
00224 for(int qq=0;qq<MultiSample;++qq)
00225 {
00226 Point3f pp2=pp+delta[qq];
00227 field_value ff= DistanceFromMesh(pp2,_oldM);
00228 if(ff.first==false) return field_value(false,0);
00229 distSum += fabs(ff.second);
00230 if(ff.second>0) positiveCnt ++;
00231 }
00232 if(positiveCnt<=MultiSample/2) distSum = -distSum;
00233 return field_value(true, distSum/MultiSample);
00234 }
00235
00238 void ComputeSliceValues(int slice,field_value *slice_values)
00239 {
00240 for (int i=0; i<=this->siz.X(); i++)
00241 {
00242 for (int k=0; k<=this->siz.Z(); k++)
00243 {
00244 int index=GetSliceIndex(i,k);
00245 Point3f pp(i,slice,k);
00246 if(this->MultiSampleFlag) slice_values[index] = MultiDistanceFromMesh(pp,_oldM);
00247 else slice_values[index] = DistanceFromMesh(pp,_oldM);
00248 }
00249 }
00250
00251 }
00252
00253
00254
00255
00256
00257 void ComputeConsensus(int slice, field_value *slice_values)
00258 {
00259 float max_dist = min(min(this->voxel[0],this->voxel[1]),this->voxel[2]);
00260 int flippedCnt=0;
00261 int flippedTot=0;
00262 int flippedTimes=0;
00263 do
00264 {
00265 flippedCnt=0;
00266 for (int i=0; i<=this->siz.X(); i++)
00267 {
00268 for (int k=0; k<=this->siz.Z(); k++)
00269 {
00270 int goodCnt=0;
00271 int badCnt=0;
00272 int index=GetSliceIndex(i,k);
00273 int index_l,index_r,index_u,index_d;
00274 if(slice_values[index].first)
00275 {
00276 float curVal= slice_values[index].second;
00277 if(i > 0 ) index_l=GetSliceIndex(i-1,k); else index_l = index;
00278 if(i < this->siz.X() ) index_r=GetSliceIndex(i+1,k); else index_r = index;
00279 if(k > 0 ) index_d=GetSliceIndex(i,k-1); else index_d = index;
00280 if(k < this->siz.Z() ) index_u=GetSliceIndex(i,k+1); else index_u = index;
00281
00282 if(slice_values[index_l].first) { goodCnt++; if(fabs(slice_values[index_l].second - curVal) > max_dist) badCnt++; }
00283 if(slice_values[index_r].first) { goodCnt++; if(fabs(slice_values[index_r].second - curVal) > max_dist) badCnt++; }
00284 if(slice_values[index_u].first) { goodCnt++; if(fabs(slice_values[index_u].second - curVal) > max_dist) badCnt++; }
00285 if(slice_values[index_d].first) { goodCnt++; if(fabs(slice_values[index_d].second - curVal) > max_dist) badCnt++; }
00286
00287 if(badCnt >= goodCnt) {
00288 slice_values[index].second *=-1.0f;
00289
00290 flippedCnt++;
00291 }
00292 }
00293 }
00294 }
00295 flippedTot+=flippedCnt;
00296 flippedTimes++;
00297 } while(flippedCnt>0);
00298
00299
00300 #ifndef NO_QT
00301 if(flippedTot>0)
00302 qDebug("Flipped %i values in %i times",flippedTot,flippedTimes);
00303 #endif
00304 }
00305 template<class EXTRACTOR_TYPE>
00306 void ProcessSlice(EXTRACTOR_TYPE &extractor)
00307 {
00308 for (int i=0; i<this->siz.X(); i++)
00309 {
00310 for (int k=0; k<this->siz.Z(); k++)
00311 {
00312 bool goodCell=true;
00313 Point3i p1(i,CurrentSlice,k);
00314 Point3i p2=p1+Point3i(1,1,1);
00315 for(int ii=0;ii<2;++ii)
00316 for(int jj=0;jj<2;++jj)
00317 for(int kk=0;kk<2;++kk)
00318 goodCell &= VV(p1[0]+ii,p1[1]+jj,p1[2]+kk).first;
00319
00320 if(goodCell) extractor.ProcessCell(p1, p2);
00321 }
00322 }
00323 }
00324
00325
00326 template<class EXTRACTOR_TYPE>
00327 void BuildMesh(Old_Mesh &old_mesh,New_Mesh &new_mesh,EXTRACTOR_TYPE &extractor,vcg::CallBackPos *cb)
00328 {
00329 _newM=&new_mesh;
00330 _oldM=&old_mesh;
00331
00332
00333 tri::UpdateNormals<Old_Mesh>::PerFaceNormalized(old_mesh);
00334 tri::UpdateNormals<Old_Mesh>::PerVertexAngleWeighted(old_mesh);
00335 tri::UpdateFlags<Old_Mesh>::FaceProjection(old_mesh);
00336 int _size=(int)old_mesh.fn*100;
00337
00338 _g.Set(_oldM->face.begin(),_oldM->face.end(),_size);
00339 markerFunctor.SetMesh(&old_mesh);
00340
00341 _newM->Clear();
00342
00343 Begin();
00344 extractor.Initialize();
00345 for (int j=0; j<=this->siz.Y(); j++)
00346 {
00347 cb((100*j)/this->siz.Y(),"Marching ");
00348 ProcessSlice<EXTRACTOR_TYPE>(extractor);
00349 NextSlice();
00350 }
00351 extractor.Finalize();
00352 typename New_Mesh::VertexIterator vi;
00353 for(vi=new_mesh.vert.begin();vi!=new_mesh.vert.end();++vi)
00354 if(!(*vi).IsD())
00355 {
00356 IPfToPf((*vi).cP(),(*vi).P());
00357 }
00358 }
00359
00360
00361 int GetSliceIndex(int x,int z)
00362 {
00363 VertexIndex index = x+z*(this->siz.X()+1);
00364 return (index);
00365 }
00366
00367
00368 void NextSlice()
00369 {
00370
00371 memset(_x_cs, -1, SliceSize*sizeof(VertexIndex));
00372 memset(_y_cs, -1, SliceSize*sizeof(VertexIndex));
00373 memset(_z_cs, -1, SliceSize*sizeof(VertexIndex));
00374
00375
00376 std::swap(_x_cs, _x_ns);
00377 std::swap(_z_cs, _z_ns);
00378
00379 std::swap(_v_cs, _v_ns);
00380
00381 CurrentSlice ++;
00382
00383 ComputeSliceValues(CurrentSlice + 1,_v_ns);
00384 }
00385
00386
00387 void Begin()
00388 {
00389
00390 CurrentSlice = 0;
00391
00392 memset(_x_cs, -1, SliceSize*sizeof(VertexIndex));
00393 memset(_y_cs, -1, SliceSize*sizeof(VertexIndex));
00394 memset(_z_cs, -1, SliceSize*sizeof(VertexIndex));
00395 memset(_x_ns, -1, SliceSize*sizeof(VertexIndex));
00396 memset(_z_ns, -1, SliceSize*sizeof(VertexIndex));
00397
00398 ComputeSliceValues(CurrentSlice,_v_cs);
00399 ComputeSliceValues(CurrentSlice+1,_v_ns);
00400 }
00401
00402
00403
00404
00405 bool Exist(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
00406 {
00407 int i = p1.X();
00408 int z = p1.Z();
00409 VertexIndex index = i+z*this->siz.X();
00410
00411
00412 int v_ind = 0;
00413 if (p1.X()!=p2.X())
00414 {
00415 if (p1.Y()==CurrentSlice)
00416 {
00417 if (_x_cs[index]!=-1)
00418 {
00419 v_ind = _x_cs[index];
00420 v = &_newM->vert[v_ind];
00421 assert(!v->IsD());
00422 return true;
00423 }
00424
00425 }
00426 else
00427 {
00428 if (_x_ns[index]!=-1)
00429 {
00430 v_ind = _x_ns[index];
00431 v = &_newM->vert[v_ind];
00432 assert(!v->IsD());
00433 return true;
00434 }
00435 }
00436 v = NULL;
00437 return false;
00438 }
00439 else if (p1.Y()!=p2.Y())
00440 {
00441 if (_y_cs[index]!=-1)
00442 {
00443 v_ind =_y_cs[index];
00444 v = &_newM->vert[v_ind];
00445 assert(!v->IsD());
00446 return true;
00447 }
00448 else
00449 {
00450 v = NULL;
00451 return false;
00452 }
00453
00454 }
00455 else if (p1.Z()!=p2.Z())
00456
00457 {
00458 if (p1.Y()==CurrentSlice)
00459 {
00460 if ( _z_cs[index]!=-1)
00461 {
00462 v_ind = _z_cs[index];
00463 v = &_newM->vert[v_ind];
00464 assert(!v->IsD());
00465 return true;
00466 }
00467
00468 }
00469 else
00470 {
00471 if (_z_ns[index]!=-1)
00472 {
00473 v_ind = _z_ns[index];
00474 v = &_newM->vert[v_ind];
00475 assert(!v->IsD());
00476 return true;
00477 }
00478 }
00479 v = NULL;
00480 return false;
00481 }
00482 assert (0);
00483 return false;
00484 }
00485
00487 NewCoordType Interpolate(const vcg::Point3i &p1, const vcg::Point3i &p2,int dir)
00488 {
00489 float f1 = (float)V(p1);
00490 float f2 = (float)V(p2);
00491 float u = (float) f1/(f1-f2);
00492 NewCoordType ret=vcg::Point3f((float)p1.V(0),(float)p1.V(1),(float)p1.V(2));
00493 ret.V(dir) = (float) p1.V(dir)*(1.f-u) + u*(float)p2.V(dir);
00494 return (ret);
00495 }
00496
00498 void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
00499 {
00500 assert(p1.X()+1 == p2.X());
00501 assert(p1.Y() == p2.Y());
00502 assert(p1.Z() == p2.Z());
00503
00504 int i = p1.X();
00505 int z = p1.Z();
00506 VertexIndex index = i+z*this->siz.X();
00507 VertexIndex pos=-1;
00508 if (p1.Y()==CurrentSlice)
00509 {
00510 if ((pos=_x_cs[index])==-1)
00511 {
00512 _x_cs[index] = (VertexIndex) _newM->vert.size();
00513 pos = _x_cs[index];
00514 Allocator<New_Mesh>::AddVertices( *_newM, 1 );
00515 v = &_newM->vert[pos];
00516 v->P()=Interpolate(p1,p2,0);
00517 return;
00518 }
00519 }
00520 if (p1.Y()==CurrentSlice+1)
00521 {
00522 if ((pos=_x_ns[index])==-1)
00523 {
00524 _x_ns[index] = (VertexIndex) _newM->vert.size();
00525 pos = _x_ns[index];
00526 Allocator<New_Mesh>::AddVertices( *_newM, 1 );
00527 v = &_newM->vert[pos];
00528 v->P()=Interpolate(p1,p2,0);
00529 return;
00530 }
00531 }
00532 assert(pos>=0);
00533 v = &_newM->vert[pos];
00534 }
00535
00537 void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
00538 {
00539 assert(p1.X() == p2.X());
00540 assert(p1.Y()+1 == p2.Y());
00541 assert(p1.Z() == p2.Z());
00542
00543 int i = p1.X();
00544 int z = p1.Z();
00545 VertexIndex index = i+z*this->siz.X();
00546 VertexIndex pos=-1;
00547 if ((pos=_y_cs[index])==-1)
00548 {
00549 _y_cs[index] = (VertexIndex) _newM->vert.size();
00550 pos = _y_cs[index];
00551 Allocator<New_Mesh>::AddVertices( *_newM, 1);
00552 v = &_newM->vert[ pos ];
00553 v->P()=Interpolate(p1,p2,1);
00554 }
00555 assert(pos>=0);
00556 v = &_newM->vert[pos];
00557 }
00558
00560 void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
00561 {
00562 assert(p1.X() == p2.X());
00563 assert(p1.Y() == p2.Y());
00564 assert(p1.Z()+1 == p2.Z());
00565
00566 int i = p1.X();
00567 int z = p1.Z();
00568 VertexIndex index = i+z*this->siz.X();
00569
00570 VertexIndex pos=-1;
00571 if (p1.Y()==CurrentSlice)
00572 {
00573 if ((pos=_z_cs[index])==-1)
00574 {
00575 _z_cs[index] = (VertexIndex) _newM->vert.size();
00576 pos = _z_cs[index];
00577 Allocator<New_Mesh>::AddVertices( *_newM, 1 );
00578 v = &_newM->vert[pos];
00579 v->P()=Interpolate(p1,p2,2);
00580 return;
00581 }
00582 }
00583 if (p1.Y()==CurrentSlice+1)
00584 {
00585 if ((pos=_z_ns[index])==-1)
00586 {
00587 _z_ns[index] = (VertexIndex) _newM->vert.size();
00588 pos = _z_ns[index];
00589 Allocator<New_Mesh>::AddVertices( *_newM, 1 );
00590 v = &_newM->vert[pos];
00591 v->P()=Interpolate(p1,p2,2);
00592 return;
00593 }
00594 }
00595 assert(pos>=0);
00596 v = &_newM->vert[pos];
00597 }
00598
00599 };
00600
00601 public:
00602
00603 typedef Walker MyWalker;
00604
00605 typedef vcg::tri::MarchingCubes<New_Mesh, MyWalker> MyMarchingCubes;
00606
00608 static void Resample(Old_Mesh &old_mesh,New_Mesh &new_mesh, Box3f volumeBox, vcg::Point3<int> accuracy,float max_dist, float thr=0, bool DiscretizeFlag=false, bool MultiSampleFlag=false, bool AbsDistFlag=false, vcg::CallBackPos *cb=0 )
00609 {
00611 vcg::tri::UpdateBounding<Old_Mesh>::Box(old_mesh);
00612
00613 MyWalker walker(volumeBox,accuracy);
00614
00615 walker.max_dim=max_dist+fabs(thr);
00616 walker.offset = - thr;
00617 walker.DiscretizeFlag = DiscretizeFlag;
00618 walker.MultiSampleFlag = MultiSampleFlag;
00619 walker.AbsDistFlag = AbsDistFlag;
00620 MyMarchingCubes mc(new_mesh, walker);
00621 walker.BuildMesh(old_mesh,new_mesh,mc,cb);
00622 }
00623
00624
00625 };
00626
00627 };
00628 };
00629 #endif
