00001
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
00032
00033
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
00062
00063
00064
00065
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
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
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());
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
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
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;
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
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
00195 const imesh_vertex_array<2>& verts = mesh_one.vertices<2>();
00196 imesh_vertex_array<3>* verts3 = new imesh_vertex_array<3>();
00197
00198
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
00218
00219
00220
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
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
00252 void sdet_image_mesh::set_anchor_points(imesh_mesh& mesh, vil_image_view<float> dt_img)
00253 {
00254
00255 int ni = resc_->ni();
00256 int nj = resc_->nj();
00257 vil_image_view<float> tri_depth(ni, nj);
00258
00259
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
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
00290 float minz, maxz;
00291 vil_math_value_range(tri_depth, minz, maxz);
00292
00293
00294 vil_image_view<float> z_img = brip_vil_float_ops::convert_to_float(resc_);
00295
00296
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
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
00315 vil_image_view_base_sptr stretched = vil_convert_stretch_range( float(0.0f), resc_->get_view());
00316
00317
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
00324 vil_image_view_base_sptr converted = vil_convert_cast<vxl_byte>(0, dest_sptr);
00325
00326
00327 resc_ = vil_new_image_resource_of_view(*converted.ptr());
00328 }
00329 }
00330