From 0ebd1f6f91f9574fbd25780186a41464a328b2f2 Mon Sep 17 00:00:00 2001 From: nicopietroni Date: Tue, 8 Feb 2005 17:50:41 +0000 Subject: [PATCH] optimized ( distance map are calculated 1 time each point) --- vcg/complex/trimesh/create/resampler.h | 507 ++++++++++++++++++------- 1 file changed, 378 insertions(+), 129 deletions(-) diff --git a/vcg/complex/trimesh/create/resampler.h b/vcg/complex/trimesh/create/resampler.h index 4b423548..f44b0b62 100644 --- a/vcg/complex/trimesh/create/resampler.h +++ b/vcg/complex/trimesh/create/resampler.h @@ -10,6 +10,8 @@ #include #include +//#include //debugghe + namespace vcg { namespace trimesh { @@ -48,6 +50,7 @@ class Resampler:RES typedef typename Old_Mesh::FaceContainer FaceCont; typedef typename GridStaticPtr Grid; typedef typename vcg::Box3 BoundingBox; + //typedef typename std::pair PointPair; typedef vcg::tri::Allocator< New_Mesh > Allocator; protected: @@ -56,6 +59,7 @@ class Resampler:RES vcg::Point3i _cell_size; float max_dim; + float dim_diag; int _slice_dimension; int _current_slice; @@ -67,23 +71,29 @@ class Resampler:RES VertexIndex *_x_ns; // indici dell'intersezioni della superficie lungo gli Xedge della prossima fetta VertexIndex *_z_ns; // indici dell'intersezioni della superficie lungo gli Zedge della prossima fetta + //float *_v_cs;///values of distance fields for each direction in current slice + //float *_v_ns;///values of distance fields for each direction in next slice + + typedef typename std::pair field_value; + field_value* _v_cs; + field_value* _v_ns; + New_Mesh *_newM; Old_Mesh *_oldM; Grid _g; - + public: - Walker(const BoundingBox &bbox,vcg::Point3i &resolution) - { - assert (resolution.V(0)<=bbox.DimX()); - assert (resolution.V(1)<=bbox.DimY()); - assert (resolution.V(2)<=bbox.DimZ()); + /*Walker(Volume_Dataset *Vo,float in,const BoundingBox &bbox,vcg::Point3i &resolution) + {*/ + /* init=in; + Vol=Vo;*/ - _bbox= bbox; - _resolution = resolution; - _cell_size.X() = _bbox.DimX()/_resolution.X(); - _cell_size.Y() = _bbox.DimY()/_resolution.Y(); - _cell_size.Z() = _bbox.DimZ()/_resolution.Z(); + void SetBBParameters() + { + _cell_size.X() =_bbox.DimX()/_resolution.X(); + _cell_size.Y() =_bbox.DimY()/_resolution.Y(); + _cell_size.Z() =_bbox.DimZ()/_resolution.Z(); ///extend bb until the box - resolution and cell matches while ((_bbox.DimX()%_cell_size.X())!=0) @@ -117,71 +127,204 @@ class Resampler:RES Point3f diag=Point3f((float)_cell_size.V(0),(float)_cell_size.V(1),(float)_cell_size.V(2)); max_dim=diag.Norm();///diagonal of a cell + + _current_slice = _bbox.min.Y(); + + Point3f minD=Point3f((float)_bbox.min.V(0),(float)_bbox.min.V(1),(float)_bbox.min.V(2)); + Point3f maxD=Point3f((float)_bbox.max.V(0),(float)_bbox.max.V(1),(float)_bbox.max.V(2)); + Point3f d=(maxD-minD); + dim_diag=d.Norm(); + } + + Walker(const BoundingBox &bbox,vcg::Point3i &resolution) + { + assert (resolution.V(0)<=bbox.DimX()); + assert (resolution.V(1)<=bbox.DimY()); + assert (resolution.V(2)<=bbox.DimZ()); + + _bbox= bbox; + + _resolution = resolution; + + SetBBParameters(); _x_cs = new VertexIndex[ _slice_dimension ]; _y_cs = new VertexIndex[ _slice_dimension ]; _z_cs = new VertexIndex[ _slice_dimension ]; _x_ns = new VertexIndex[ _slice_dimension ]; _z_ns = new VertexIndex[ _slice_dimension ]; + + _v_cs= new field_value[(_resolution.X()+1)*(_resolution.Z()+1)]; + _v_ns= new field_value[(_resolution.X()+1)*(_resolution.Z()+1)]; }; ~Walker() {} + + float V(Point3i p) + { + return (V(p.V(0),p.V(1),p.V(2))); + } + + float V(int x,int y,int z) + { + assert ((y==_current_slice)||(y==(_current_slice+_cell_size.Y()))); + + //test if it is outside the bb of the mesh + //vcg::Point3f test=vcg::Point3f((float)x,(float)y,(float)z); + /*if (!_oldM->bbox.IsIn(test)) + return (1.f);*/ + int index=GetSliceIndex(x,z); + + if (y==_current_slice) + { + //assert(_v_cs[index]bbox.IsIn(test)) + { + dist=1.f; + return true; + }*/ + + vcg::Point3f Norm; + vcg::Point3f Target; + vcg::Point3f pip; + + vcg::trimesh::Closest((*mesh),test,_g,dist,Norm,Target,f,pip); + + if (f==NULL) + return false; + else + { + Point3f dir=(test-Target); + //dist=dir.Norm(); + + dir.Normalize(); + //direction of normal inside the mesh + if ((dir*Norm)<0) + dist=-dist; + //the intersection exist + return true; + } + } + + ///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 CumputeSliceValues(int slice,field_value *slice_values) + { + float dist; + /*for (int i=_bbox.min.X(); i<=_bbox.max.X()-_cell_size.X(); i+=_cell_size.X()) + { + for (int k=_bbox.min.Z(); k<=_bbox.max.Z()-_cell_size.X(); k+=_cell_size.Z()) + { + */ + for (int i=_bbox.min.X(); i<=_bbox.max.X(); i+=_cell_size.X()) + { + for (int k=_bbox.min.Z(); k<=_bbox.max.Z(); k+=_cell_size.Z()) + { + + int index=GetSliceIndex(i,k); + if (DistanceFromMesh(i,slice,k,_oldM,dist))///compute the distance,inside volume of the mesh is negative + { + //put computed values in the slice values matrix + slice_values[index]=field_value(true,dist); + //end putting values + } + else + slice_values[index]=field_value(false,dist); + } + } + } + + template + void ProcessSlice(std::vector cells,EXTRACTOR_TYPE &extractor) + { + std::vector::iterator it; + + for (it=cells.begin();it BBf=vcg::Box3(min-cell,max+cell); + + _g.SetBBox(BBf); + + _g.Set(_oldM->face); + } + template void BuildMesh(Old_Mesh &old_mesh,New_Mesh &new_mesh,EXTRACTOR_TYPE &extractor) { _newM=&new_mesh; _oldM=&old_mesh; - - Point3f min=Point3f((float)_bbox.min.V(0),(float)_bbox.min.V(1),(float)_bbox.min.V(2)); - Point3f max=Point3f((float)_bbox.max.V(0),(float)_bbox.max.V(1),(float)_bbox.max.V(2)); - - vcg::Box3 BBf=vcg::Box3(min,max); - - _g.SetBBox(BBf); - - _g.Set(_oldM->face); + + SetUGrid(); _newM->Clear(); - vcg::Point3i p1, p2; + vcg::Point3i p1, p2; Begin(); extractor.Initialize(); - for (int j=_bbox.min.Y(); j<_bbox.max.Y()-_cell_size.Y(); j+=_cell_size.Y()) + for (int j=_bbox.min.Y(); j<=_bbox.max.Y()-_cell_size.Y(); j+=_cell_size.Y()) { - for (int i=_bbox.min.X(); i<_bbox.max.X()-_cell_size.X(); i+=_cell_size.X()) - { - for (int k=_bbox.min.Z(); k<_bbox.max.Z()-_cell_size.Z(); k+=_cell_size.Z()) - { - p1.X()=i; - p1.Y()=j; - p1.Z()=k; - p2.X()=i+_cell_size.X(); - p2.Y()=j+_cell_size.Y(); - p2.Z()=k+_cell_size.Z(); - if (ExistIntersection(p1,p2)) - extractor.ProcessCell(p1, p2); - } - } + ProcessSlice(FindCells(),extractor);//find cells where there is the isosurface and examine it NextSlice(); } extractor.Finalize(); - /*_newM= NULL;*/ - }; - - float V(Point3i p) - { - return (V(p.V(0),p.V(1),p.V(2))); } - - ///control if exist any intersection with the mesh using all points of the cell - bool ExistIntersection(Point3i min,Point3i max) + + //return the index of a vertex in slide as it was stored + int GetSliceIndex(int x,int z) { - vcg::Point3i _corners[8]; + int ii = (x - _bbox.min.X())/_cell_size.X(); + int zz = (z - _bbox.min.Z())/_cell_size.Z(); + VertexIndex index = ii+zz*(_resolution.X()+1); + return (index); + } + ///return true if exist in the cell one value <0 and another one >0 + bool FindMinMax(vcg::Point3i min,vcg::Point3i max) + { + assert((min.X()((*_oldM),test,_g,distm,Norm,Target,f,pip); - if (f==NULL) + //if one value is > that bbox.diag this value is not valid + //that is the mark + + if (_corners[i].Y()==_current_slice) + value=_v_cs[GetSliceIndex(_corners[i].X(),_corners[i].Z())]; + else + value=_v_ns[GetSliceIndex(_corners[i].X(),_corners[i].Z())]; + + if (value.first==false) return false; - } - return true; + + //assign new values of min and max + if (value.secondmax_value) + max_value=value.second; + } + + /////do not test with zero.. + if ((min_value<=0.f)&&(max_value>=0.f)) + return true; + + return false; + //return true; } - - ///si potrebbe ottimizzare sfruttando valori che in realta' ho gia' calcolato - float V(int pi, int pj, int pk) + + ///filter the cells from to_hexamine vector to the ones that + /// min and max of the cell are <0 and >0 + std::vector FindCells() { - vcg::Point3f Norm; - Old_Mesh::FaceType *f=NULL; - float dist=max_dim;; - vcg::Point3f test=vcg::Point3f((float)pi,(float)pj,(float)pk); - - Point3f pip; - Point3f Target; - vcg::trimesh::Closest((*_oldM),test,_g,dist,Norm,Target,f,pip); - - assert(f!=NULL); - - Point3f dir=(test-Target); - - //dist=dir.Norm(); - dir=dir.Normalize(); - - //direction of normal inside or outside the mesh - if ((f->N()*dir)>0) - return (dist); - else - return (-dist); + std::vector res; + for (int i=_bbox.min.X(); i<=_bbox.max.X()-_cell_size.X(); i+=_cell_size.X()) + { + for (int k=_bbox.min.Z(); k<=_bbox.max.Z()-_cell_size.Z(); k+=_cell_size.Z()) + { + int x0=i; + int y0=_current_slice; + int z0=k; + int x1=x0+_cell_size.X(); + int y1=y0+_cell_size.Y(); + int z1=z0+_cell_size.Z(); + if (FindMinMax(Point3i(x0,y0,z0),Point3i(x1,y1,z1))) + res.push_back(Point3i(x0,y0,z0)); + } + } + return res; } + //swap slices , the initial value of distance fields ids set as double of bbox of space + void NextSlice() + { + + memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); + + + std::swap(_x_cs, _x_ns); + std::swap(_z_cs, _z_ns); + + std::swap(_v_cs, _v_ns); + + _current_slice += _cell_size.Y(); + + //memset(_v_ns, dim_diag*2.f, _slice_dimension*sizeof(float)); + //memset(_v_ns, field_value(false,0.f), _slice_dimension*sizeof(field_value)); + + CumputeSliceValues(_current_slice+ _cell_size.Y(),_v_ns); + } + + //initialize data strucures , the initial value of distance fields ids set as double of bbox of space + void Begin() + { + + _current_slice = _bbox.min.Y(); + + memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_x_ns, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_z_ns, -1, _slice_dimension*sizeof(VertexIndex)); + + /*memset(_v_cs, dim_diag*2.f, _slice_dimension*sizeof(float)); + memset(_v_ns, dim_diag*2.f, _slice_dimension*sizeof(float));*/ + + /*memset(_v_cs, field_value(false,0.f), _slice_dimension*sizeof(field_value)); + memset(_v_ns, field_value(false,0.f), _slice_dimension*sizeof(field_value));*/ + + CumputeSliceValues(_current_slice,_v_cs); + CumputeSliceValues(_current_slice+_cell_size.Y(),_v_ns); + + } + + + + bool Exist(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) { - return false; - //int i_idx = p1.X()-_bbox.min.X(); - //int k_idx = p2.Z()-_bbox.min.Z(); - //int index = i_idx+k_idx*_resolution.X(); - //if (p1.X()!=p2.X()) //intersezione della superficie con un Xedge - // return (p1.Y()==_current_slice)? _x_cs[index]!=-1 : _x_ns[index]!=-1; - //else if (p1.Y()!=p2.Y()) //intersezione della superficie con un Yedge - // return _y_cs[index]!=-1; - //else if (p1.Z()!=p2.Z()) //intersezione della superficie con un Zedge - // return (p1.Y()==_current_slice)? _z_cs[index]!=-1 : _z_ns[index]!=-1; + int i = (p1.X() - _bbox.min.X())/_cell_size.X(); + int z = (p1.Z() - _bbox.min.Z())/_cell_size.Z(); + VertexIndex index = i+z*_resolution.X(); + + //VertexIndex index =GetSliceIndex(// + int v_ind = 0; + if (p1.X()!=p2.X()) //intersezione della superficie con un Xedge + { + if (p1.Y()==_current_slice) + { + if (_x_cs[index]!=-1) + { + v_ind = _x_cs[index]; + v = &_newM->vert[v_ind]; + assert(!v->IsD()); + return true; + } + + } + else + { + if (_x_ns[index]!=-1) + { + v_ind = _x_ns[index]; + v = &_newM->vert[v_ind]; + assert(!v->IsD()); + return true; + } + } + v = NULL; + return false; + } + else if (p1.Y()!=p2.Y()) //intersezione della superficie con un Yedge + { + if (_y_cs[index]!=-1) + { + v_ind =_y_cs[index]; + v = &_newM->vert[v_ind]; + assert(!v->IsD()); + return true; + } + else + { + v = NULL; + return false; + } + + } + else if (p1.Z()!=p2.Z()) + //intersezione della superficie con un Zedge + { + if (p1.Y()==_current_slice) + { + if ( _z_cs[index]!=-1) + { + v_ind = _z_cs[index]; + v = &_newM->vert[v_ind]; + assert(!v->IsD()); + return true; + } + + } + else + { + if (_z_ns[index]!=-1) + { + v_ind = _z_ns[index]; + v = &_newM->vert[v_ind]; + assert(!v->IsD()); + return true; + } + } + v = NULL; + return false; + } + assert (0); } - + ///interpolate NewCoordType Interpolate(const vcg::Point3i &p1, const vcg::Point3i &p2,int dir) { - float f1 = V(p1); - float f2 = V(p2); + float f1 = (float)V(p1); + float f2 = (float)V(p2); float u = (float) f1/(f1-f2); NewCoordType ret=vcg::Point3f((float)p1.V(0),(float)p1.V(1),(float)p1.V(2)); - ret.V(dir) = (float) p1.V(dir)*(1-u) + u*p2.V(dir); + ret.V(dir) = (float) p1.V(dir)*(1.f-u) + u*(float)p2.V(dir); return (ret); } + ///if there is a vertex in z axis of a cell return the vertex or create it void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) { + assert ((p1.Y()==_current_slice)||(p1.Y()==(_current_slice+_cell_size.Y()))); + int i = (p1.X() - _bbox.min.X())/_cell_size.X(); int z = (p1.Z() - _bbox.min.Z())/_cell_size.Z(); VertexIndex index = i+z*_resolution.X(); @@ -293,9 +561,12 @@ class Resampler:RES } v = &_newM->vert[pos]; } - + + ///if there is a vertex in y axis of a cell return the vertex or create it void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) { + assert ((p1.Y()==_current_slice)||(p1.Y()==(_current_slice+_cell_size.Y()))); + int i = (p1.X() - _bbox.min.X())/_cell_size.X(); int z = (p1.Z() - _bbox.min.Z())/_cell_size.Z(); VertexIndex index = i+z*_resolution.X(); @@ -310,9 +581,12 @@ class Resampler:RES } v = &_newM->vert[pos]; } - + + ///if there is a vertex in z axis of a cell return the vertex or create it void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) { + assert ((p1.Y()==_current_slice)||(p1.Y()==(_current_slice+_cell_size.Y()))); + int i = (p1.X() - _bbox.min.X())/_cell_size.X(); int z = (p1.Z() - _bbox.min.Z())/_cell_size.Z(); VertexIndex index = i+z*_resolution.X(); @@ -345,31 +619,6 @@ class Resampler:RES v = &_newM->vert[pos]; } - - - void NextSlice() - { - memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); - - std::swap(_x_cs, _x_ns); - std::swap(_z_cs, _z_ns); - - _current_slice += _cell_size.Y(); - } - - void Begin() - { - _current_slice = _bbox.min.Y(); - - memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_x_ns, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_z_ns, -1, _slice_dimension*sizeof(VertexIndex)); - - } };//end class walker public: @@ -402,11 +651,10 @@ static void Resample(Old_Mesh &old_mesh,New_Mesh &new_mesh,vcg::Point3 accu Point3i min=Point3i((int)ceil(old_mesh.bbox.min.V(0)),(int)ceil(old_mesh.bbox.min.V(1)),(int)ceil(old_mesh.bbox.min.V(2))); Point3i max=Point3i((int)ceil(old_mesh.bbox.max.V(0)),(int)ceil(old_mesh.bbox.max.V(1)),(int)ceil(old_mesh.bbox.max.V(2))); - ///exetend out to BB for resample the limits of the mesh - /*min-=Point3i(1,1,1); - max+=Point3i(1,1,1);*/ - vcg::Box3 boxInt=Box3(min,max); + vcg::Box3 boxInt=Box3(min,max); + + float rx=((float)boxInt.DimX())/(float)accuracy.X(); float ry=((float)boxInt.DimY())/(float)accuracy.Y(); float rz=((float)boxInt.DimZ())/(float)accuracy.Z(); @@ -418,6 +666,7 @@ static void Resample(Old_Mesh &old_mesh,New_Mesh &new_mesh,vcg::Point3 accu Point3i res=Point3i(rxi,ryi,rzi); MyWalker walker(boxInt,res); + if (mm==MMarchingCubes) { MarchingCubes mc(new_mesh, walker);