contrib/brl/bseg/sdet/sdet_image_mesh.cxx
Go to the documentation of this file.
00001 // This is brl/bseg/sdet/sdet_image_mesh.cxx
00002 #include "sdet_image_mesh.h"
00003 #include <vgl/vgl_point_2d.h>
00004 #include <vgl/vgl_line_segment_2d.h>
00005 #include <vil/vil_image_view.h>
00006 #include <vil/vil_math.h>
00007 #include <vil/vil_convert.h>
00008 #include <vil/vil_new.h>
00009 #include <vil/vil_pixel_format.h>
00010 #include <vil/vil_save.h>
00011 #include <vtol/vtol_edge_2d_sptr.h>
00012 #include <vsol/vsol_line_2d.h>
00013 #include <vtol/vtol_edge_2d.h>
00014 #include <vdgl/vdgl_digital_curve.h>
00015 #include <vdgl/vdgl_interpolator.h>
00016 #include <vdgl/vdgl_interpolator_sptr.h>
00017 #include <vdgl/vdgl_edgel_chain.h>
00018 #include <vdgl/vdgl_edgel_chain_sptr.h>
00019 
00020 #include <brip/brip_vil_float_ops.h>
00021 #include <brip/brip_line_generator.h>
00022 
00023 #include <sdet/sdet_detector_params.h>
00024 #include <sdet/sdet_detector.h>
00025 #include <sdet/sdet_fit_lines_params.h>
00026 #include <sdet/sdet_fit_lines.h>
00027 #include <bbas/bil/algo/bil_cedt.h>
00028 #include <bvgl/bvgl_triangle_interpolation_iterator.h>
00029 #include <imesh/algo/imesh_generate_mesh.h>
00030 
00031 //note: this method is somewhat of a hack and should be replaced by
00032 // a computed step function transition width, e.g. by a 2nd derivative
00033 // operator
00034 bool sdet_image_mesh:: step_boundary(vgl_line_segment_2d<double> const& parent,
00035                                      vgl_line_segment_2d<double>& child0,
00036                                      vgl_line_segment_2d<double>& child1)
00037 {
00038   if (! resc_ ) return false;
00039   unsigned ni = resc_->ni(), nj = resc_->nj();
00040   vgl_vector_2d<double> n = parent.normal();
00041   vgl_point_2d<double> p1 = parent.point1();
00042   vgl_point_2d<double> p2 = parent.point2();
00043   vgl_point_2d<double> p10 = p1 - step_half_width_*n;
00044   vgl_point_2d<double> p11 = p1 + step_half_width_*n;
00045   vgl_point_2d<double> p20 = p2 - step_half_width_*n;
00046   vgl_point_2d<double> p21 = p2 + step_half_width_*n;
00047   if (p10.x()<0 || p10.y()<0||p20.x()<0||p20.y()<0||
00048       p10.x()>=ni ||p10.y()>=nj||p20.x()>=ni||p20.y()>=nj)
00049     return false;
00050   else
00051     child0 = vgl_line_segment_2d<double>(p10, p20);
00052   if (p11.x()<0 || p11.y()<0||p21.x()<0||p21.y()<0||
00053       p11.x()>=ni ||p11.y()>=nj||p21.x()>=ni||p21.y()>=nj)
00054     return false;
00055   else
00056     child1 = vgl_line_segment_2d<double>(p11, p21);
00057   return true;
00058 }
00059 
00060 //---------------------------------------------------------------
00061 // Constructors
00062 //
00063 //----------------------------------------------------------------
00064 
00065 // constructor from a parameter block (the only way)
00066 sdet_image_mesh::sdet_image_mesh(sdet_image_mesh_params& imp)
00067   : sdet_image_mesh_params(imp), mesh_valid_(false), resc_(0)
00068 {
00069 }
00070 
00071 // Destructor
00072 sdet_image_mesh::~sdet_image_mesh()
00073 {
00074 }
00075 
00076 bool sdet_image_mesh::compute_line_segments(vil_image_resource_sptr const& resc,
00077                                             vcl_vector<vgl_line_segment_2d<double> > & segs)
00078 {
00079   if (!resc) return false;
00080   mesh_valid_ = false;
00081   sdet_detector_params dp;
00082   dp.smooth= smooth_;
00083   dp.noise_multiplier = thresh_;
00084   dp.aggressive_junction_closure=0;
00085   dp.junctionp=false;
00086   sdet_detector det(dp);
00087   det.SetImage(resc);
00088   det.DoContour();
00089   vcl_vector<vtol_edge_2d_sptr>* edges = det.GetEdges();
00090   if (!edges) {
00091     vcl_cout<<"sdet_image_mesh:: could not detect edges in buffer"<<vcl_endl;
00092     return false;
00093   }
00094 
00095   //fit lines on edges
00096   sdet_fit_lines_params flp;
00097   flp.min_fit_length_ = min_fit_length_;
00098   flp.rms_distance_ = rms_distance_;
00099   sdet_fit_lines fl(flp);
00100   fl.set_edges(*edges);
00101   bool fit_worked = fl.fit_lines();
00102   if (!fit_worked) {
00103     vcl_cout<<"sdet_image_mesh:: could not fit lines on edges"<<vcl_endl;
00104     return false;
00105   }
00106 
00107   fl.get_line_segs(segs);
00108   return true;
00109 }
00110 
00111 bool sdet_image_mesh::compute_mesh()
00112 {
00113   vil_image_view<unsigned char> line_img(resc_->ni(),resc_->nj());//some resource
00114   line_img.fill(255);
00115 
00116   vcl_vector<vgl_line_segment_2d<double> > segs, segs_pair;
00117   sdet_detector_params dp;
00118   dp.smooth= smooth_;
00119   dp.noise_multiplier = thresh_;
00120   dp.aggressive_junction_closure=0;
00121   dp.junctionp=false;
00122   sdet_detector det(dp);
00123   det.SetImage(resc_);
00124   det.DoContour();
00125   vcl_vector<vtol_edge_2d_sptr>* edges = det.GetEdges();
00126   for (vcl_vector<vtol_edge_2d_sptr>::iterator eit = edges->begin();
00127        eit != edges->end(); eit++)
00128   {
00129     vsol_curve_2d_sptr c = (*eit)->curve();
00130     vdgl_digital_curve_sptr dc = c->cast_to_vdgl_digital_curve();
00131     if (!dc)
00132       continue;
00133     vdgl_interpolator_sptr intp = dc->get_interpolator();
00134     vdgl_edgel_chain_sptr ec = intp->get_edgel_chain();
00135     int nedgl = ec->size();
00136     for (int i=0; i<nedgl; i++)
00137     {
00138       // anchor_points_.push_back(vgl_point_2d<double>((*ec)[i].x(),(*ec)[i].y()));
00139       line_img((int)(*ec)[i].x(),(int) (*ec)[i].y())=0;
00140     }
00141   }
00142 
00143  bil_cedt dt(line_img);
00144  if (!dt.compute_cedt())
00145   vcl_cout<<"Error in computing DT"<<vcl_endl;
00146 
00147  vil_image_view<float> dtimg=dt.cedtimg();
00148  vil_image_view<unsigned char> dtimg_threshed(dtimg.ni(),dtimg.nj());dtimg_threshed.fill(0);
00149  for (unsigned i=0;i<dtimg.ni();i++)
00150    for (unsigned j=0;j<dtimg.nj();j++)
00151      if (dtimg(i,j)<2.0)
00152        dtimg_threshed(i,j)=255;
00153  vil_save(dtimg_threshed,"F:/visdt/dt_thresh.png");
00154 
00155  //segs_pair.clear();
00156  vcl_vector<vgl_line_segment_2d<double> > lines;
00157  this->compute_line_segments(vil_new_image_resource_of_view(dtimg_threshed),lines);
00158  segs_pair.insert(segs_pair.end(),lines.begin(),lines.end());
00159  line_img.fill(255);
00160 
00161  for (unsigned i = 0;i<segs_pair.size();i++)
00162  {
00163    bool init = true;
00164    float xs= float(segs_pair[i].point1().x());
00165    float ys= float(segs_pair[i].point1().y());
00166    float xe= float(segs_pair[i].point2().x());
00167    float ye= float(segs_pair[i].point2().y());
00168    float x=0.0f;
00169    float y=0.0f;
00170    while (brip_line_generator::generate(init, xs, ys, xe, ye, x, y))
00171    {
00172      int xi = (int)x;
00173      int yi = (int)y; //convert the pixel location to integer
00174      line_img(xi, yi)=0;
00175    }
00176  }
00177  vil_save(line_img,"F:/visdt/line_img_more.png");
00178 
00179  bil_cedt dt1(line_img);
00180  if (!dt1.compute_cedt())
00181    vcl_cout<<"Error in computing DT"<<vcl_endl;
00182 
00183   //generate a 2d mesh based on the edges
00184   vgl_point_2d<double> ul(0.0, 0.0), ur(resc_->ni()-1,0.0);
00185   vgl_point_2d<double> lr(resc_->ni()-1, resc_->nj()-1), ll(0.0, resc_->nj()-1);
00186 
00187   vcl_vector<vgl_point_2d<double> > convex_hull;
00188   convex_hull.push_back(ul);   convex_hull.push_back(ur);
00189   convex_hull.push_back(lr);   convex_hull.push_back(ll);
00190   vcl_vector<vgl_point_2d<double> > cvexh = convex_hull;
00191   imesh_mesh mesh_one;
00192   imesh_generate_mesh_2d_2(convex_hull, segs_pair,anchor_points_, mesh_one);
00193 
00194   //lift vertices to 3-d
00195   const imesh_vertex_array<2>& verts = mesh_one.vertices<2>();
00196   imesh_vertex_array<3>* verts3 = new imesh_vertex_array<3>();
00197 
00198   // convert image to float
00199   vil_image_view<float> view = brip_vil_float_ops::convert_to_float(resc_);
00200   unsigned ni = view.ni(), nj = view.nj();
00201   float minv=0, maxv=0;
00202   vil_math_value_range(view, minv, maxv);
00203   unsigned nverts = mesh_one.num_verts();
00204   for (unsigned iv = 0; iv<nverts; ++iv)
00205   {
00206     unsigned i = static_cast<unsigned>(verts[iv][0]);
00207     unsigned j = static_cast<unsigned>(verts[iv][1]);
00208     double height =maxv;
00209     if (i<ni && j<nj)
00210       height = static_cast<double>(view(i,j));
00211     height = maxv-height;
00212     imesh_vertex<3> v3(verts[iv][0], verts[iv][1], height);
00213     verts3->push_back(v3);
00214   }
00215   vcl_auto_ptr<imesh_vertex_array_base> v3(verts3);
00216   mesh_one.set_vertices(v3);
00217   //mesh_valid_ = true;
00218 
00219   ///////////////////////////////////////////////////////
00220   //compute anchor points, and rerun generate_mesh_2d_2
00221   ///////////////////////////////////////////////////////
00222   this->set_anchor_points(mesh_one,dt1.cedtimg());
00223 
00224   vcl_cout<<"Number of anchor points: "<<anchor_points_.size()<<vcl_endl;
00225   imesh_generate_mesh_2d_2(cvexh, segs_pair, anchor_points_, mesh_);
00226   const imesh_vertex_array<2>& verts2 = mesh_.vertices<2>();
00227   imesh_vertex_array<3>* newVerts = new imesh_vertex_array<3>();
00228 
00229   // convert image to float
00230   ni = view.ni(), nj = view.nj();
00231   vil_math_value_range(view, minv, maxv);
00232   nverts = mesh_.num_verts();
00233   for (unsigned iv = 0; iv<nverts; ++iv)
00234   {
00235     unsigned i = static_cast<unsigned>(verts2[iv][0]);
00236     unsigned j = static_cast<unsigned>(verts2[iv][1]);
00237     double height =maxv;
00238     if (i<ni && j<nj)
00239       height = static_cast<double>(view(i,j));
00240     height = maxv-height;
00241     imesh_vertex<3> v3(verts2[iv][0], verts2[iv][1], height);
00242     newVerts->push_back(v3);
00243   }
00244   vcl_auto_ptr<imesh_vertex_array_base> v3_ptr(newVerts);
00245   mesh_.set_vertices(v3_ptr);
00246   mesh_valid_ = true;
00247 
00248   return true;
00249 }
00250 
00251 //takes the 3d mesh, calculates depth disparity map (from top)
00252 void sdet_image_mesh::set_anchor_points(imesh_mesh& mesh, vil_image_view<float> dt_img)
00253 {
00254   //create a tri_depth image
00255   int ni = resc_->ni();
00256   int nj = resc_->nj();
00257   vil_image_view<float> tri_depth(ni, nj);
00258 
00259   //find the range of triangles in the Z direction
00260   imesh_regular_face_array<3>& faces = (imesh_regular_face_array<3>&) mesh.faces();
00261   imesh_vertex_array<3>& verts = mesh.vertices<3>();
00262   unsigned nfaces = mesh.num_faces();
00263   for (unsigned iface = 0; iface<nfaces; ++iface)
00264   {
00265     unsigned v1 = faces[iface][0];
00266     unsigned v2 = faces[iface][1];
00267     unsigned v3 = faces[iface][2];
00268     double verts_x[] = { verts[v1][0], verts[v2][0], verts[v3][0] };
00269     double verts_y[] = { verts[v1][1], verts[v2][1], verts[v3][1] };
00270     double verts_z[] = { verts[v1][2], verts[v2][2], verts[v3][2] };
00271     bvgl_triangle_interpolation_iterator<double> tsi(verts_x, verts_y, verts_z);
00272 
00273     //scan the triangle, storing the depth
00274     for (tsi.reset(); tsi.next(); ) {
00275       int y = tsi.scany();
00276       if (y<0 || y>=nj) continue;
00277       int min_x = tsi.startx();
00278       int max_x = tsi.endx();
00279       if (min_x >= ni || max_x < 0)
00280         continue;
00281       if (min_x < 0) min_x = 0;
00282       if (max_x >= ni) max_x = ni-1;
00283       for (int x = min_x; x <= max_x; ++x) {
00284         tri_depth(x,y) = (float) tsi.value_at(x);
00285       }
00286     }
00287   }
00288 
00289   //get range of triangle depths
00290   float minz, maxz;
00291   vil_math_value_range(tri_depth, minz, maxz);
00292 
00293   //maximum z diff allowed for a triangle is 1/512 of the total z range
00294   vil_image_view<float> z_img = brip_vil_float_ops::convert_to_float(resc_);
00295 
00296   //for ( unsigned niter=1;niter<16;niter*=2)
00297   {
00298     double max_z_diff = (maxz-minz)/64.0;
00299     for (int i=0; i<ni; i+=4)
00300       for (int j=0; j<nj; j+=4)
00301         if ( vcl_fabs( tri_depth(i,j)- z_img(i,j) ) > max_z_diff && dt_img(i,j) >= 3.5 )
00302           anchor_points_.push_back(vgl_point_2d<double>(i,j));
00303   }
00304 }
00305 
00306 //ensure the image is a byte image (between 0 and 255)
00307 void sdet_image_mesh::set_image(vil_image_resource_sptr const& resource)
00308 {
00309   resc_ = resource;
00310   if (resc_->pixel_format() != VIL_PIXEL_FORMAT_BYTE )
00311   {
00312     vcl_cout<<"Converting image from "<<resc_->pixel_format()<<" to vxl_byte image"<<vcl_endl;
00313 
00314     //make the float image on the range of [0,1];
00315     vil_image_view_base_sptr stretched = vil_convert_stretch_range( float(0.0f), resc_->get_view());
00316 
00317     //makvil_image_view_base_sptr byte_img = vil_convert_cast<vxl_byte>(0, stretched*255.0f);
00318     vil_image_view_base_sptr dest_sptr = new vil_image_view<float>(stretched->ni(), stretched->nj());
00319     vil_image_view<float>* dest = (vil_image_view<float>*) dest_sptr.ptr();
00320     vil_image_view<float>* stf  = (vil_image_view<float>*) stretched.ptr();
00321     vil_convert_stretch_range<float>( *stf, *dest, 0.0f, 255.0f);
00322 
00323     //now turn em into bytes
00324     vil_image_view_base_sptr converted = vil_convert_cast<vxl_byte>(0, dest_sptr);
00325 
00326     //create new resource sptr
00327     resc_ = vil_new_image_resource_of_view(*converted.ptr());
00328   }
00329 }
00330