Main Page | Modules | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members

# Vol_GrabberT.hh

    00001 #ifndef VOL_GRABBERT_HH
    00002 #define VOL_GRABBERT_HH
    00003
    00004 #include <iostream>
    00005 #include <stdio.h>
    00006
    00007 #include <OpenMesh/Core/IO/MeshIO.hh>
    00008 #include <OpenMesh/Core/Mesh/TriMeshT.hh>
    00009
    00010 // ACG
    00011 #include <ACG/Geometry/Algorithms.hh>
    00012
    00013 #include "HelperVolume.hh"
    00014
    00015 #define Epsilon 0.001
    00016
    00017 template <class MyMesh,  class Volume>
    00018 class Vol_GrabberT
    00019 {
    00020 public:
    00021
    00022   // default constructor
    00023   Vol_GrabberT(MyMesh& _mesh) : epsilon_(1), mesh_(_mesh) { }
    00024
    00025   ~Vol_GrabberT() {}
    00026
    00027   void set_mesh(MyMesh& _mesh) { mesh_ = _mesh;}
    00028
    00029   inline void set_epsilon(unsigned int _e) { epsilon_ = _e;  }
    00030
    00031   void render_into_volume(Volume& _volume ,HelperVolume& _helper);
    00032
    00033 private:
    00034
    00035   typedef typename MyMesh::FaceVertexIter FaceVertexIter;
    00036   typedef typename MyMesh::VertexIter VertexIter;
    00037   typedef typename MyMesh::FaceIter FaceIter;
    00038   typedef typename MyMesh::VertexFaceIter VertexFaceIter;
    00039   typedef typename MyMesh::FaceHandle FaceHandle;
    00040   typedef typename MyMesh::VertexHandle VertexHandle;
    00041
    00042   OpenMesh::VPropHandleT<typename MyMesh::Normal> angle_normal;
    00043
    00044   unsigned int epsilon_;
    00045   MyMesh& mesh_;
    00046
    00047   unsigned int xres, yres, zres;
    00048
    00049   typedef typename Volume::Scalar Scalar;
    00050   typedef OpenMesh::VectorT<Scalar,3> Point;
    00051
    00052   //-----------------------------------------------------------------------
    00053   void update_angle_weighted_vertex_normals();
    00054
    00055   inline void get_triangle_bb(Point p0, Point p1, Point p2,
    00056                               OpenMesh::Vec3i &_bbmin,
    00057                               OpenMesh::Vec3i &_bbmax)
    00058   {
    00059
    00060     Point tmpbbmin(p0), tmpbbmax(p0);
    00061     tmpbbmin.minimize(p1); tmpbbmin.minimize(p2);
    00062     tmpbbmax.maximize(p1); tmpbbmax.maximize(p2);
    00063
    00064     _bbmin = OpenMesh::Vec3i((unsigned int)(floor(tmpbbmin[0])),
    00065                              (unsigned int)(floor(tmpbbmin[1])),
    00066                              (unsigned int)(floor(tmpbbmin[2])));
    00067     _bbmax = OpenMesh::Vec3i((unsigned int)(ceil(tmpbbmax[0])),
    00068                              (unsigned int)(ceil(tmpbbmax[1])),
    00069                              (unsigned int)(ceil(tmpbbmax[2])));
    00070   }
    00071
    00072
    00073   void bounding_box(Point &_BBMin, Point &_BBMax)
    00074   {
    00075     VertexIter v_it(mesh_.vertices_begin()), v_end(mesh_.vertices_end());
    00076
    00077     _BBMin = mesh_.point(v_it.handle());
    00078     _BBMax = mesh_.point(v_it.handle());
    00079
    00080     for (; v_it != v_end; ++v_it)
    00081     {
    00082       Point p = mesh_.point(v_it.handle());
    00083       _BBMin.minimize(p);
    00084       _BBMax.maximize(p);
    00085     }
    00086   }
    00087
    00088   inline Scalar maximum(Scalar a, Scalar b, Scalar c)
    00089   {
    00090     if (a > b && a > c) return a;
    00091     if (b > a && b > c) return b;
    00092     return c;
    00093   }
    00094
    00095   inline Point mat_mult(Point z1, Point z2, Point z3, Point x)
    00096   {
    00097     return Point(z1 | x,  z2 | x,  z3 | x);
    00098   };
    00099
    00100
    00101 };
    00102
    00103
    00104 //-- Implementation of special normals and main procedure --------------------
    00105
    00106
    00107 template <class MyMesh, class Volume>
    00108 void Vol_GrabberT< MyMesh, Volume >::
    00109 update_angle_weighted_vertex_normals()
    00110 {
    00111
    00112   VertexIter v_it, end = mesh_.vertices_end();
    00113   for (v_it = mesh_.vertices_begin(); v_it != end; ++v_it)
    00114   {
    00115     Point point = mesh_.point(v_it.handle());
    00116     Point final_normal(0,0,0);
    00117     Scalar angle_sum = 0;
    00118
    00119     VertexFaceIter vf_it;
    00120     for (vf_it = mesh_.vf_iter(v_it.handle()); vf_it; ++vf_it)
    00121     {
    00122       FaceVertexIter fv_it = mesh_.fv_iter(vf_it.handle());
    00123       Point p0 = mesh_.point(fv_it.handle()); ++fv_it;
    00124       Point p1 = mesh_.point(fv_it.handle()); ++fv_it;
    00125       Point p2 = mesh_.point(fv_it.handle());
    00126
    00127       Point normal = mesh_.normal(vf_it);
    00128
    00129       Scalar angle; Point e1, e2;
    00130
    00131       if (p0 == point) { e1 = p1-p0; e2 = p2-p0; }
    00132       else if (p1 == point) { e1 = p0-p1; e2 = p2-p1; }
    00133       else { e1 = p0-p2; e2 = p1-p2; }
    00134
    00135       angle = acos(std::min((Scalar)1.0,
    00136                             std::max((Scalar)-1.0,
    00137                                      (e1|e2)/(e1.norm()*e2.norm()))));
    00138       final_normal += angle*normal;
    00139       angle_sum += angle;
    00140     }
    00141     final_normal.normalize();
    00142     mesh_.property(angle_normal, v_it.handle()) = final_normal;
    00143   }
    00144 }
    00145
    00146 //-----------------------------------------------------------------------------
    00147
    00148 template < class MyMesh, class Volume>
    00149 void Vol_GrabberT< MyMesh ,Volume >::
    00150 render_into_volume(Volume& _volume, HelperVolume& _helper)
    00151 {
    00152
    00153   // register vertex property storing angle weighted vertex normals
    00154   mesh_.add_property( angle_normal );
    00155
    00156   Point z1 ,z2 ,z3, t; // Matrix-rows, Translation
    00157
    00158   // bbmin + x_vec * i + y_vec * j + z_vec * k ( vol_points   )
    00159   // -> invers:  compute basis transform
    00160
    00161   Point x_unit_dir = _volume.x_axis()/Scalar( _volume.x_resolution() -1);
    00162   Point y_unit_dir = _volume.y_axis()/Scalar( _volume.y_resolution() -1);
    00163   Point z_unit_dir = _volume.z_axis()/Scalar( _volume.z_resolution() -1);
    00164
    00165   // ?? -1
    00166
    00167   z1 =  x_unit_dir / (x_unit_dir | x_unit_dir );
    00168   z2 =  y_unit_dir / (y_unit_dir | y_unit_dir );
    00169   z3 =  z_unit_dir / (z_unit_dir | z_unit_dir );
    00170
    00171   t = mat_mult(z1, z2, z3 ,-_volume.origin());
    00172
    00173   Scalar voxel_diag_squared = (x_unit_dir | x_unit_dir );
    00174
    00175   //init volume
    00176   for (unsigned int i = 0; i < _volume.x_resolution(); ++i)
    00177     for (unsigned int j = 0; j < _volume.y_resolution(); ++j)
    00178       for (unsigned int k = 0; k < _volume.z_resolution(); ++k)
    00179         _volume(i,j,k) = FLT_MAX ;
    00180
    00181   // computing special normals for sign determination
    00182   update_angle_weighted_vertex_normals();
    00183
    00184   // compute distances
    00185
    00186   OpenMesh::Vec3i bb_tri_min, bb_tri_max;
    00187
    00188   FaceIter f_it, f_end = mesh_.faces_end();
    00189
    00190   for (f_it = mesh_.faces_begin(); f_it != f_end; ++f_it)
    00191   {
    00192     FaceHandle fh = f_it.handle();
    00193     Point tri_normal;
    00194     FaceVertexIter fv_it = mesh_.fv_iter(fh);
    00195     VertexHandle vh0, vh1, vh2;
    00196
    00197     vh0 = fv_it.handle();
    00198     Point p0 = mat_mult(z1, z2, z3, mesh_.point(vh0)) + t; ++fv_it;
    00199     Point p0_v = mesh_.point(vh0);
    00200     vh1 = fv_it.handle();
    00201     Point p1 = mat_mult(z1, z2, z3, mesh_.point(vh1)) + t; ++fv_it;
    00202     Point p1_v = mesh_.point(vh1);
    00203     vh2 = fv_it.handle();
    00204     Point p2 = mat_mult(z1, z2, z3, mesh_.point(vh2)) + t;
    00205     Point p2_v = mesh_.point(vh2);
    00206
    00207
    00208     // calculate triangle BB
    00209     get_triangle_bb(p0, p1, p2, bb_tri_min, bb_tri_max);
    00210
    00211     // compute max in unit-dimensions
    00212     int minx = bb_tri_min[0] - epsilon_;
    00213     int miny = bb_tri_min[1] - epsilon_;
    00214     int minz = bb_tri_min[2] - epsilon_;
    00215
    00216     int maxx = bb_tri_max[0] + epsilon_;
    00217     int maxy = bb_tri_max[1] + epsilon_;
    00218     int maxz = bb_tri_max[2] + epsilon_;
    00219
    00220     // check that epsilon BB does not extend beyond volume limits
    00221     if (minx < 0) minx = 0; if (miny < 0) miny = 0; if (minz < 0) minz = 0;
    00222     if (maxx >= (int)_volume.x_resolution()) maxx = _volume.x_resolution() - 1;
    00223     if (maxy >= (int)_volume.y_resolution()) maxy = _volume.y_resolution() - 1;
    00224     if (maxz >= (int)_volume.z_resolution()) maxz = _volume.z_resolution() - 1;
    00225
    00226
    00227     // for every voxel in the BB, calculate distance value to the triangle
    00228
    00229     Point nearest;  // nearest on triangle
    00230
    00231     // for all voxels in extended bounding box
    00232     for (int x = minx; x <= maxx; x++)
    00233       for (int y = miny; y <= maxy; y++)
    00234         for (int z = minz; z <= maxz; z++)
    00235         {
    00236           //Point voxel_center(x+0.5, y+0.5, z+0.5);
    00237           // Scalar Grid errechnet distances nicht im center sondern auf grid
    00238
    00239           Point voxel_center(x, y, z);
    00240           voxel_center = _volume.origin() +
    00241             Scalar(x) * _volume.x_axis() /Scalar(_volume.x_resolution()-1) +
    00242             Scalar(y) * _volume.y_axis() /Scalar(_volume.y_resolution()-1) +
    00243             Scalar(z) * _volume.z_axis() /Scalar(_volume.z_resolution()-1);
    00244
    00245
    00246           Scalar d = ACG::Geometry::distPointTriangleSquared (voxel_center,
    00247                                                               p0_v,p1_v,p2_v,
    00248                                                               nearest);
    00249
    00250           if (d < voxel_diag_squared) // > max dist in voxel
    00251           {
    00252             d = sqrt (d);
    00253
    00254             // test if inner or outer by using angle normals
    00255             if (fabs(d) < fabs (_volume(x,y,z)) )
    00256             {
    00257               // sign decision:
    00258               Point barys;
    00259               ACG::Geometry::baryCoord( nearest,
    00260                                         p0_v, p1_v, p2_v,
    00261                                         barys );
    00262
    00263               tri_normal = barys[0] * mesh_.property(angle_normal, vh0) +
    00264                            barys[1] * mesh_.property(angle_normal, vh1) +
    00265                            barys[2] * mesh_.property(angle_normal, vh2);
    00266
    00267               tri_normal.normalize();
    00268
    00269               Point diff_vec = nearest - voxel_center;
    00270
    00271               Scalar sign_ = diff_vec | tri_normal;
    00272
    00273               if (sign_ >= 0.0)
    00274                 _volume(x,y,z) = -d;
    00275               else _volume(x,y,z) = d;
    00276               _helper(x,y,z).set_conquered(true);
    00277             };
    00278           };
    00279         };
    00280   };
    00281   // initial values computed !!
    00282
    00283   mesh_.remove_property( angle_normal );
    00284 };
    00285
    00286
    00287 #endif // VOL_GRABBERT_HH