1 #include "betr_edgel_factory.h"
2 #include <sdet/sdet_detector.h>
3 #include <vdgl/vdgl_digital_curve.h>
4 #include <vdgl/vdgl_interpolator.h>
5 #include <vdgl/vdgl_edgel_chain.h>
6 #include <vdgl/vdgl_edgel.h>
7 #include <brip/brip_vil_float_ops.h>
8 #include "vil/vil_image_resource.h"
9 #include "vil/vil_image_view_base.h"
10 #include "vil/vil_image_view.h"
11 #include "vil/vil_save.h"
12 #include "vil/vil_new.h"
13 #include "vil/vil_resample_bicub.h"
14 #include <bsol/bsol_algs.h>
15 #include <vsol/vsol_polygon_2d.h>
16 #include <vsol/vsol_point_2d.h>
17 #include "vsl/vsl_binary_io.h"
18 #include "vsl/vsl_vector_io.h"
19
add_image(std::string const & iname,vil_image_resource_sptr const & imgr)20 bool betr_edgel_factory::add_image(std::string const& iname, vil_image_resource_sptr const& imgr){
21 if (!imgr||!imgr->ni()||!imgr->nj())
22 {
23 std::cout << "empty image resource when adding image\n";
24 return false;
25 }
26 images_[iname]=imgr;
27 brip_roi_sptr broi = new brip_roi(imgr->ni(), imgr->nj());
28 rois_[iname] = broi;
29 return true;
30 }
add_region(std::string const & iname,std::string const & region_name,vsol_box_2d_sptr const & box)31 bool betr_edgel_factory::add_region(std::string const& iname, std::string const& region_name,
32 vsol_box_2d_sptr const& box){
33 brip_roi_sptr broi = rois_[iname];
34 if(!broi){
35 std::cout << "null roi for image " << iname << '\n';
36 return false;
37 }
38 unsigned n = broi->n_regions();
39 broi->add_region(box);
40 regions_[iname][region_name] = n;
41 return true;
42 }
add_region(std::string const & iname,std::string const & region_name,vsol_polygon_2d_sptr const & poly)43 bool betr_edgel_factory::add_region(std::string const& iname, std::string const& region_name, vsol_polygon_2d_sptr const& poly){
44 brip_roi_sptr broi = rois_[iname];
45 if(!broi){
46 std::cout << "null roi for image " << iname << '\n';
47 return false;
48 }
49 unsigned n = broi->n_regions();
50 vsol_box_2d_sptr box = poly->get_bounding_box();
51 bool good = this->add_region(iname, region_name, box);
52 if(good)
53 polys_[iname][n] = poly;
54 return good;
55 }
add_region_from_origin_and_size(std::string const & iname,std::string const & region_name,unsigned i0,unsigned j0,unsigned ni,unsigned nj)56 bool betr_edgel_factory::add_region_from_origin_and_size(std::string const& iname, std::string const& region_name,
57 unsigned i0, unsigned j0, unsigned ni, unsigned nj){
58 brip_roi_sptr broi = rois_[iname];
59 if(!broi){
60 std::cout << "null roi for image " << iname << '\n';
61 return false;
62 }
63 unsigned n = broi->n_regions();
64 broi->add_region(i0, j0, ni, nj);
65 regions_[iname][region_name] = n;
66 return true;
67 }
process(const std::string & iname,const std::string & region_name)68 bool betr_edgel_factory::process(const std::string& iname, const std::string& region_name){
69 vil_image_resource_sptr imgr = images_[iname];
70 if(!imgr){
71 std::cout << "image " << iname << " not found in map \n";
72 return false;
73 }
74 brip_roi_sptr roi = rois_[iname];
75 if(!roi){
76 std::cout << "roi for " << iname << " not found in map \n";
77 return false;
78 }
79 unsigned region_id = regions_[iname][region_name];
80 unsigned ni = roi->csize(region_id);
81 unsigned nj = roi->rsize(region_id);
82 if(ni < params_.min_region_edge_length_ || nj <params_.min_region_edge_length_ ){
83 std::cout << "roi " << region_name << " for " << iname << " is empty \n";
84 return false;
85 }
86 int imin = roi->cmin(region_id);
87 int jmin = roi->rmin(region_id);
88 //debug
89 int imax = imin + ni;
90 int jmax = jmin + nj;
91 std::cout << "region(" << iname << ':' << region_name << ")(" << imin << ',' << jmin << ")->(" << imax << ',' << jmax <<")\n";
92 //
93 brip_roi_sptr temp_roi = new brip_roi();
94 temp_roi->add_region(imin, jmin, ni, nj);
95 vil_image_resource_sptr clip_resc;
96 if(!brip_vil_float_ops::chip(imgr, temp_roi, clip_resc)){
97 std::cout << " clip function failed for (" << imin << ' ' << jmin << ' ' << ni << ' ' << nj << ") on " << iname << '\n';
98 return false;
99 }
100 if(clip_resc->pixel_format()==VIL_PIXEL_FORMAT_FLOAT){
101 vil_image_view<vxl_byte> view = brip_vil_float_ops::convert_to_byte(clip_resc);
102 clip_resc = vil_new_image_resource_of_view(view);
103 }
104 // check if chip needs to be upsampled
105 if(params_.upsample_factor_ != 1.0){
106 auto dni = static_cast<double>(clip_resc->ni()), dnj = static_cast<double>(clip_resc->nj());
107 dni *= params_.upsample_factor_; dnj *= params_.upsample_factor_;
108 auto ni = static_cast<unsigned>(dni), nj = static_cast<unsigned>(dnj);
109 if(clip_resc->pixel_format()==VIL_PIXEL_FORMAT_UINT_16){
110 vil_image_view<unsigned short> temp = clip_resc->get_view(), uptemp;
111 vil_resample_bicub(temp, uptemp, ni, nj);
112 clip_resc = vil_new_image_resource_of_view(uptemp);
113 }else if(clip_resc->pixel_format()==VIL_PIXEL_FORMAT_BYTE){
114 vil_image_view<unsigned char> temp = clip_resc->get_view(), uptemp;
115 vil_resample_bicub(temp, uptemp, ni, nj);
116 clip_resc = vil_new_image_resource_of_view(uptemp);
117 }
118 }
119 #if 0
120 //debug
121 std::string dir = "D:/tests/kandahar_test/";
122 std::string fname = dir + iname + "_" + region_name + ".tif";
123 vil_save_image_resource(clip_resc, fname.c_str());
124 #endif
125 sdet_detector det(params_.det_params_);
126
127 det.SetImage(clip_resc);
128
129 if(!det.DoContour()){
130 std::cout << "Edgel detection failed\n";
131 return false;
132 }
133 std::vector<vdgl_digital_curve_sptr> vd_edges;
134 if (!det.get_vdgl_edges(vd_edges)){
135 std::cout << "Detection worked but returned no edgels\n";
136 return false;
137 }
138
139 edgels_[iname][region_name] = vd_edges;
140
141 unsigned region_index = regions_[iname][region_name];
142 vsol_polygon_2d_sptr poly = polys_[iname][region_index];
143 std::vector<double> gmags;
144 if(!poly){
145 std::cout << "warning - can't find region poly - using roi instead\n";
146 this->grad_mags(iname, region_name, gmags);
147 }else
148 this->grad_mags(iname, region_name, poly, gmags);
149 bsta_histogram<double> h(params_.gradient_range_, params_.nbins_);
150 for(double & gmag : gmags)
151 h.upcount(gmag, (1.0 + gmag));//increase weight to favor high gradient values (small objects)
152
153 if(h.area()<3.0*params_.nbins_){
154 std::cout << "insufficient edges in region " << region_name << " - fatal" << std::endl;
155 return false;
156 }
157 grad_hists_[iname][region_name] = h;
158 return true;
159 }
grad_mags(const std::string & iname,const std::string & region_name,std::vector<double> & mags)160 bool betr_edgel_factory::grad_mags(const std::string& iname, const std::string& region_name, std::vector<double>& mags){
161 std::vector< vdgl_digital_curve_sptr > vd_edges = edgels_[iname][region_name];
162 if(!vd_edges.size()){
163 std::cout << "no edgels for the specified region " << iname << ':' << region_name << '\n';
164 return false;
165 }
166 for (auto & vd_edge : vd_edges)
167 {
168 //get the edgel chain
169 vdgl_interpolator_sptr itrp = vd_edge->get_interpolator();
170 vdgl_edgel_chain_sptr ech = itrp->get_edgel_chain();
171 unsigned int n = ech->size();
172 for (unsigned int i=0; i<n;i++)
173 {
174 vdgl_edgel ed = (*ech)[i];
175 double grad_mag = ed.get_grad();
176 mags.push_back(grad_mag);
177 }
178 }
179 return true;
180 }
grad_mags(const std::string & iname,const std::string & region_name,vsol_polygon_2d_sptr const & poly,std::vector<double> & mags)181 bool betr_edgel_factory::grad_mags(const std::string& iname, const std::string& region_name, vsol_polygon_2d_sptr const& poly,
182 std::vector<double>& mags){
183 vgl_polygon<double> vpoly = bsol_algs::vgl_from_poly(poly); // THIS IS WHERE YOU CAN CHECK IF IT IS INSIDE POLY
184 brip_roi_sptr broi = rois_[iname];
185 if(!broi){
186 std::cout << "no roi for " << iname << '\n';
187 return false;
188 }
189 unsigned region_index = regions_[iname][region_name];
190 std::vector< vdgl_digital_curve_sptr > vd_edges = edgels_[iname][region_name];
191 if(!vd_edges.size()){
192 std::cout << "no edgels for the specified region " << iname << ':' << region_name << '\n';
193 return false;
194 }
195 auto x0 = static_cast<double>(broi->cmin(region_index));
196 auto y0 = static_cast<double>(broi->rmin(region_index));
197 for (auto & vd_edge : vd_edges)
198 {
199 //get the edgel chain
200 vdgl_interpolator_sptr itrp = vd_edge->get_interpolator();
201 vdgl_edgel_chain_sptr ech = itrp->get_edgel_chain();
202 unsigned int n = ech->size();
203 for (unsigned int i=0; i<n;i++)
204 {
205 vdgl_edgel& ed = (*ech)[i];
206 double x = ed.get_x(), y = ed.get_y();
207 x/= params_.upsample_factor_; y/=params_.upsample_factor_;
208 x += x0, y += y0;
209 // ed.set_x(x); ed.set_y(y);
210 if(!vpoly.contains(x, y))
211 continue;
212 double grad_mag = ed.get_grad();
213 mags.push_back(grad_mag);
214 }
215 }
216 return true;
217 }
save_edgels(std::string const & dir) const218 bool betr_edgel_factory::save_edgels(std::string const& dir) const {
219 auto iit = edgels_.begin();
220 for(; iit != edgels_.end(); ++iit){
221 const std::pair<std::string, std::map<std::string, std::vector< vdgl_digital_curve_sptr > > >& emap = *iit;
222 for(const auto & eit : emap.second){
223 std::string region_name = eit.first;
224 std::vector< vdgl_digital_curve_sptr > edges = eit.second;
225 std::vector<vsol_digital_curve_2d_sptr> vsol_edges = sdet_detector::convert_vdgl_to_vsol(edges);
226 // convert to spatial object 2d
227 std::vector<vsol_spatial_object_2d_sptr> sos;
228 for(std::vector<vsol_digital_curve_2d_sptr>::const_iterator cit = vsol_edges.begin();
229 cit!= vsol_edges.end(); ++cit){
230 vsol_spatial_object_2d_sptr ec = dynamic_cast<vsol_spatial_object_2d*>(cit->ptr());
231 sos.push_back(ec);
232 }
233 std::string path = dir + region_name + "_edges.vsl";
234 vsl_b_ofstream ostr(path);
235 if(!ostr){
236 std::cout << "couldn't open binary stream for " << path << '\n';
237 return false;
238 }
239 vsl_b_write(ostr, sos);
240 ostr.close();
241 }
242 }
243 return true;
244 }
save_edgels_in_poly(std::string const & identifier,std::string const & dir)245 bool betr_edgel_factory::save_edgels_in_poly(std::string const& identifier, std::string const& dir){
246 std::map<std::string, std::map<std::string, std::vector< vdgl_digital_curve_sptr > > >::const_iterator iit = edgels_.begin();
247 for(; iit != edgels_.end(); ++iit){
248 const std::pair<std::string, std::map<std::string, std::vector< vdgl_digital_curve_sptr > > >& emap = *iit;
249 std::string iname = emap.first;
250 for(const auto & eit : emap.second){
251 std::string region_name = eit.first;
252 unsigned index = regions_[iname][region_name];
253 vsol_polygon_2d_sptr poly = polys_[iname][index];
254 if(!poly){
255 std::cout << "Null polygon in save_edgel_in_poly\n";
256 return false;
257 }
258 vgl_polygon<double> vpoly = bsol_algs::vgl_from_poly(poly);
259 std::vector< vdgl_digital_curve_sptr > edges = eit.second;
260 std::vector<vsol_digital_curve_2d_sptr> vsol_edges = sdet_detector::convert_vdgl_to_vsol(edges);
261 // convert to spatial object 2d
262 std::vector<vsol_spatial_object_2d_sptr> sos;
263 for(std::vector<vsol_digital_curve_2d_sptr>::const_iterator cit = vsol_edges.begin();
264 cit!= vsol_edges.end(); ++cit){
265 unsigned n = (*cit)->size();
266 for(unsigned i = 0; i<n; ++i){
267 vsol_point_2d_sptr p = (*cit)->point(i);
268 if(!vpoly.contains(p->x(), p->y()))
269 continue;
270 vsol_spatial_object_2d_sptr sp = dynamic_cast<vsol_point_2d*>(p.ptr());
271 sos.push_back(sp);
272 }
273 }
274 std::string path = dir + identifier + "_" + region_name + "_edges.vsl";
275 vsl_b_ofstream ostr(path);
276 if(!ostr){
277 std::cout << "couldn't open binary stream for " << path << '\n';
278 return false;
279 }
280 vsl_b_write(ostr, sos);
281 ostr.close();
282 }
283 }
284 return true;
285
286 }
287 vil_image_resource_sptr betr_edgel_factory::
edgel_image(const std::string & iname,const std::string & region_name,unsigned & i_offset,unsigned & j_offset)288 edgel_image(const std::string& iname, const std::string& region_name, unsigned& i_offset, unsigned& j_offset){
289 brip_roi_sptr roi = rois_[iname];
290 if(!roi){
291 std::cout << "roi for " << iname << " not found in map \n";
292 return nullptr;
293 }
294 unsigned region_id = regions_[iname][region_name];
295 unsigned ni = roi->csize(region_id);
296 unsigned nj = roi->rsize(region_id);
297 i_offset = static_cast<unsigned>(roi->cmin(region_id));
298 j_offset = static_cast<unsigned>(roi->rmin(region_id));
299 vil_image_view<vxl_byte> view(ni, nj);
300 view.fill(vxl_byte(0));
301 vsol_polygon_2d_sptr poly = polys_[iname][region_id];
302 if(!poly){
303 std::cout << "Null polygon in edgel_image for " << iname << ':' << region_name << "\n";
304 return nullptr;
305 }
306 vgl_polygon<double> vpoly = bsol_algs::vgl_from_poly(poly);
307 std::vector< vdgl_digital_curve_sptr > edges = edgels_[iname][region_name];
308 if(edges.size() == 0){
309 std::cout << "No edgels for " << iname << ':' << region_name << "\n";
310 return nullptr;
311 }
312 double grad_scale = 255.0/params_.gradient_range_;
313 for(auto & edge : edges){
314 vdgl_edgel_chain_sptr echain = edge->get_interpolator()->get_edgel_chain();
315 if(!echain){
316 std::cout << "Null edgel chain " << iname << ':' << region_name << "\n";
317 return nullptr;
318 }
319 unsigned n = echain->size();
320 unsigned bdr = 3;
321 for(unsigned k = 0; k<n; ++k){
322 vdgl_edgel& e = echain->edgel(k);
323 double x = e.get_x(), y = e.get_y();
324 x/= params_.upsample_factor_; y/=params_.upsample_factor_;
325 if(!vpoly.contains(x+i_offset, y+j_offset))
326 continue;
327 auto i = static_cast<unsigned>(x), j = static_cast<unsigned>(y);
328 // clean away border hash due to convolving and upsampling
329 if(i < bdr || j < bdr || i > (ni-bdr-1) || j > (nj-bdr-1))
330 continue;
331 double g = e.get_grad();
332 double v = g*grad_scale;
333 if(v>255.0)
334 v= 255;
335 view(i,j) = static_cast<vxl_byte>(v);
336 }
337 }
338 vil_image_resource_sptr ret = vil_new_image_resource_of_view(view);
339 return ret;
340 }
341