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_TRIVIAL_WALKER
00024 #define __VCG_TRIVIAL_WALKER
00025 #include <wrap/callback.h>
00026
00027 namespace vcg {
00028
00029
00030
00031
00032 template <class VOX_TYPE>
00033 class SimpleVolume
00034 {
00035 public:
00036 typedef VOX_TYPE VoxelType;
00037 std::vector<VoxelType> Vol;
00038
00039 Point3i sz;
00040
00041 const Point3i &ISize() {return sz;};
00042
00043 void Init(Point3i _sz)
00044 {
00045 sz=_sz;
00046 Vol.resize(sz[0]*sz[1]*sz[2]);
00047 }
00048
00049 float Val(const int &x,const int &y,const int &z) const {
00050 return cV(x,y,z).V();
00051
00052 }
00053
00054 float &Val(const int &x,const int &y,const int &z) {
00055 return V(x,y,z).V();
00056
00057 }
00058
00059 VOX_TYPE &V(const int &x,const int &y,const int &z) {
00060 return Vol[x+y*sz[0]+z*sz[0]*sz[1]];
00061 }
00062
00063 const VOX_TYPE &cV(const int &x,const int &y,const int &z) const {
00064 return Vol[x+y*sz[0]+z*sz[0]*sz[1]];
00065 }
00066
00067
00068 typedef enum { XAxis=0,YAxis=1,ZAxis=2} VolumeAxis;
00069
00070 template < class VertexPointerType, VolumeAxis AxisVal >
00071 void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr)
00072 {
00073 float f1 = Val(p1.X(), p1.Y(), p1.Z())-thr;
00074 float f2 = Val(p2.X(), p2.Y(), p2.Z())-thr;
00075 float u = (float) f1/(f1-f2);
00076 if(AxisVal==XAxis) v->P().X() = (float) p1.X()*(1-u) + u*p2.X();
00077 else v->P().X() = (float) p1.X();
00078 if(AxisVal==YAxis) v->P().Y() = (float) p1.Y()*(1-u) + u*p2.Y();
00079 else v->P().Y() = (float) p1.Y();
00080 if(AxisVal==ZAxis) v->P().Z() = (float) p1.Z()*(1-u) + u*p2.Z();
00081 else v->P().Z() = (float) p1.Z();
00082 }
00083
00084 template < class VertexPointerType >
00085 void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr)
00086 { GetIntercept<VertexPointerType,XAxis>(p1,p2,v,thr); }
00087
00088 template < class VertexPointerType >
00089 void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr)
00090 { GetIntercept<VertexPointerType,YAxis>(p1,p2,v,thr); }
00091
00092 template < class VertexPointerType >
00093 void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr)
00094 { GetIntercept<VertexPointerType,ZAxis>(p1,p2,v,thr); }
00095 };
00096 template <class VolumeType>
00097 class RawVolumeImporter
00098 {
00099 public:
00100 enum DataType
00101 {
00102
00103 UNDEF=0,
00104 BYTE=1,
00105 SHORT=2,
00106 FLOAT=3
00107 };
00108
00109 static bool Open(const char *filename, VolumeType &V, Point3i sz, DataType d)
00110 {
00111 return true;
00112 }
00113 };
00114
00115 class SimpleVoxel
00116 {
00117 private:
00118 float _v;
00119 public:
00120 float &V() {return _v;};
00121 float V() const {return _v;};
00122 };
00123
00124
00125 namespace tri {
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136 template <class MeshType, class VolumeType>
00137 class TrivialWalker
00138 {
00139 private:
00140 typedef int VertexIndex;
00141 typedef typename MeshType::ScalarType ScalarType;
00142 typedef typename MeshType::VertexPointer VertexPointer;
00143 public:
00144
00145
00146
00147
00148
00149
00150 void Init(VolumeType &volume)
00151 {
00152 _bbox = Box3i(Point3i(0,0,0),volume.ISize());
00153 _slice_dimension = _bbox.DimX()*_bbox.DimZ();
00154
00155 _x_cs = new VertexIndex[ _slice_dimension ];
00156 _y_cs = new VertexIndex[ _slice_dimension ];
00157 _z_cs = new VertexIndex[ _slice_dimension ];
00158 _x_ns = new VertexIndex[ _slice_dimension ];
00159 _z_ns = new VertexIndex[ _slice_dimension ];
00160
00161 };
00162
00163 ~TrivialWalker()
00164 {_thr=0;}
00165
00166 template<class EXTRACTOR_TYPE>
00167 void BuildMesh(MeshType &mesh, VolumeType &volume, EXTRACTOR_TYPE &extractor, const float threshold, vcg::CallBackPos * cb=0)
00168 {
00169 Init(volume);
00170 _volume = &volume;
00171 _mesh = &mesh;
00172 _mesh->Clear();
00173 _thr=threshold;
00174 vcg::Point3i p1, p2;
00175
00176 Begin();
00177 extractor.Initialize();
00178 for (int j=_bbox.min.Y(); j<(_bbox.max.Y()-1)-1; j+=1)
00179 {
00180
00181 if(cb && ((j%10)==0) ) cb(j*_bbox.DimY()/100.0,"Marching volume");
00182
00183 for (int i=_bbox.min.X(); i<(_bbox.max.X()-1)-1; i+=1)
00184 {
00185 for (int k=_bbox.min.Z(); k<(_bbox.max.Z()-1)-1; k+=1)
00186 {
00187 p1.X()=i; p1.Y()=j; p1.Z()=k;
00188 p2.X()=i+1; p2.Y()=j+1; p2.Z()=k+1;
00189 extractor.ProcessCell(p1, p2);
00190 }
00191 }
00192 NextSlice();
00193 }
00194 extractor.Finalize();
00195 _volume = NULL;
00196 _mesh = NULL;
00197 };
00198
00199 float V(int pi, int pj, int pk)
00200 {
00201 return _volume->Val(pi, pj, pk)-_thr;
00202 }
00203
00204 bool Exist(const vcg::Point3i &p0, const vcg::Point3i &p1, VertexPointer &v)
00205 {
00206 int pos = p0.X()+p0.Z()*_bbox.max.X();
00207 int vidx;
00208
00209 if (p0.X()!=p1.X())
00210 vidx = (p0.Y()==_current_slice) ? _x_cs[pos] : _x_ns[pos];
00211 else if (p0.Y()!=p1.Y())
00212 vidx = _y_cs[pos];
00213 else if (p0.Z()!=p1.Z())
00214 vidx = (p0.Y()==_current_slice)? _z_cs[pos] : _z_ns[pos];
00215 else
00216 assert(false);
00217
00218 v = (vidx!=-1)? &_mesh->vert[vidx] : NULL;
00219 return v!=NULL;
00220 }
00221
00222 void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
00223 {
00224 int i = p1.X() - _bbox.min.X();
00225 int z = p1.Z() - _bbox.min.Z();
00226 VertexIndex index = i+z*_bbox.max.X();
00227 VertexIndex pos;
00228 if (p1.Y()==_current_slice)
00229 {
00230 if ((pos=_x_cs[index])==-1)
00231 {
00232 _x_cs[index] = (VertexIndex) _mesh->vert.size();
00233 pos = _x_cs[index];
00234 Allocator<MeshType>::AddVertices( *_mesh, 1 );
00235 v = &_mesh->vert[pos];
00236 _volume->GetXIntercept(p1, p2, v, _thr);
00237 return;
00238 }
00239 }
00240 if (p1.Y()==_current_slice+1)
00241 {
00242 if ((pos=_x_ns[index])==-1)
00243 {
00244 _x_ns[index] = (VertexIndex) _mesh->vert.size();
00245 pos = _x_ns[index];
00246 Allocator<MeshType>::AddVertices( *_mesh, 1 );
00247 v = &_mesh->vert[pos];
00248 _volume->GetXIntercept(p1, p2, v,_thr);
00249 return;
00250 }
00251 }
00252 assert(pos >=0 && size_t(pos)< _mesh->vert.size());
00253 v = &_mesh->vert[pos];
00254 }
00255 void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
00256 {
00257 int i = p1.X() - _bbox.min.X();
00258 int z = p1.Z() - _bbox.min.Z();
00259 VertexIndex index = i+z*_bbox.max.X();
00260 VertexIndex pos;
00261 if ((pos=_y_cs[index])==-1)
00262 {
00263 _y_cs[index] = (VertexIndex) _mesh->vert.size();
00264 pos = _y_cs[index];
00265 Allocator<MeshType>::AddVertices( *_mesh, 1);
00266 v = &_mesh->vert[ pos ];
00267 _volume->GetYIntercept(p1, p2, v,_thr);
00268 }
00269 v = &_mesh->vert[pos];
00270 }
00271 void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
00272 {
00273 int i = p1.X() - _bbox.min.X();
00274 int z = p1.Z() - _bbox.min.Z();
00275 VertexIndex index = i+z*_bbox.max.X();
00276 VertexIndex pos;
00277 if (p1.Y()==_current_slice)
00278 {
00279 if ((pos=_z_cs[index])==-1)
00280 {
00281 _z_cs[index] = (VertexIndex) _mesh->vert.size();
00282 pos = _z_cs[index];
00283 Allocator<MeshType>::AddVertices( *_mesh, 1 );
00284 v = &_mesh->vert[pos];
00285 _volume->GetZIntercept(p1, p2, v,_thr);
00286 return;
00287 }
00288 }
00289 if (p1.Y()==_current_slice+1)
00290 {
00291 if ((pos=_z_ns[index])==-1)
00292 {
00293 _z_ns[index] = (VertexIndex) _mesh->vert.size();
00294 pos = _z_ns[index];
00295 Allocator<MeshType>::AddVertices( *_mesh, 1 );
00296 v = &_mesh->vert[pos];
00297 _volume->GetZIntercept(p1, p2, v,_thr);
00298 return;
00299 }
00300 }
00301 v = &_mesh->vert[pos];
00302 }
00303
00304 protected:
00305 Box3i _bbox;
00306
00307 int _slice_dimension;
00308 int _current_slice;
00309
00310 VertexIndex *_x_cs;
00311 VertexIndex *_y_cs;
00312 VertexIndex *_z_cs;
00313 VertexIndex *_x_ns;
00314 VertexIndex *_z_ns;
00315
00316 MeshType *_mesh;
00317 VolumeType *_volume;
00318
00319 float _thr;
00320 void NextSlice()
00321 {
00322 memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex));
00323 memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex));
00324 memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex));
00325
00326 std::swap(_x_cs, _x_ns);
00327 std::swap(_z_cs, _z_ns);
00328
00329 _current_slice += 1;
00330 }
00331
00332 void Begin()
00333 {
00334 _current_slice = _bbox.min.Y();
00335
00336 memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex));
00337 memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex));
00338 memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex));
00339 memset(_x_ns, -1, _slice_dimension*sizeof(VertexIndex));
00340 memset(_z_ns, -1, _slice_dimension*sizeof(VertexIndex));
00341
00342 }
00343 };
00344 }
00345 }
00346 #endif // __VCGTEST_WALKER
