1 #include <iostream>
2 #include <algorithm>
3 #include "volm_io_tools.h"
4 //:
5 // \file
6
7
8 #include "vul/vul_file.h"
9 #include "vul/vul_file_iterator.h"
10 #include "vpgl/vpgl_utm.h"
11 #include <bkml/bkml_write.h>
12 #include "vil/vil_load.h"
13 #include "vil/vil_image_view.h"
14 #include "vil/vil_crop.h"
15 #include "vpgl/vpgl_lvcs.h"
16 #include "vpgl/vpgl_lvcs_sptr.h"
17 #ifdef _MSC_VER
18 # include "vcl_msvc_warnings.h"
19 #endif
20 #include <vgl/algo/vgl_fit_lines_2d.h>
21 #include "vgl/vgl_intersection.h"
22 #include "vgl/vgl_line_segment_3d.h"
23
24
25 unsigned int volm_io_tools::northing = 0; // WARNING: north hard-coded
26
27 static double eps = 1.0e-5;
near_zero(double x)28 inline bool near_zero(double x) { return x < eps && x > -eps; }
near_equal(double x,double y)29 inline bool near_equal(double x, double y) { return near_zero(x-y); }
near_eq_pt(vgl_point_2d<double> a,vgl_point_2d<double> b)30 bool near_eq_pt(vgl_point_2d<double> a, vgl_point_2d<double> b)
31 {
32 return (near_equal(a.x(),b.x()) && near_equal(a.y(), b.y()));
33 }
34
35
save_box_kml(const std::string & out_name)36 void volm_img_info::save_box_kml(const std::string& out_name) {
37 std::ofstream ofs(out_name.c_str());
38 bkml_write::open_document(ofs);
39 bkml_write::write_box(ofs, name, "", bbox);
40 bkml_write::close_document(ofs);
41 }
42
read_box(const std::string & bbox_file,vgl_box_2d<double> & bbox)43 bool read_box(const std::string& bbox_file, vgl_box_2d<double>& bbox) {
44 char buffer[1000];
45 std::ifstream ifs(bbox_file.c_str());
46 if (!ifs.is_open()) {
47 std::cerr << " cannot open: " << bbox_file << "!\n";
48 return false;
49 }
50
51 std::string dummy; double top_lat, bottom_lat, left_lon, right_lon;
52 for (unsigned kk = 0; kk < 22; kk++)
53 ifs.getline(buffer, 1000);
54 // top
55 ifs.getline(buffer, 1000);
56 std::stringstream top_edge_line(buffer);
57 top_edge_line >> dummy; top_edge_line >> dummy; top_edge_line >> dummy;
58 top_edge_line >> top_lat;
59 // bottom
60 ifs.getline(buffer, 1000);
61 std::stringstream bot_edge_line(buffer);
62 bot_edge_line >> dummy; bot_edge_line >> dummy; bot_edge_line >> dummy;
63 bot_edge_line >> bottom_lat;
64 // left
65 ifs.getline(buffer, 1000);
66 std::stringstream left_edge_line(buffer);
67 left_edge_line >> dummy; left_edge_line >> dummy; left_edge_line >> dummy;
68 left_edge_line >> left_lon;
69 // right
70 ifs.getline(buffer, 1000);
71 std::stringstream right_edge_line(buffer);
72 right_edge_line >> dummy; right_edge_line >> dummy; right_edge_line >> dummy;
73 right_edge_line >> right_lon;
74
75 vgl_point_2d<double> lower_left(left_lon, bottom_lat);
76 vgl_point_2d<double> upper_right(right_lon, top_lat);
77 bbox = vgl_box_2d<double>(lower_left, upper_right);
78 //std::cout << "bbox: " << bbox << std::endl;
79 return true;
80 }
81
load_naip_img(std::string const & img_folder,std::string const & name,vpgl_lvcs_sptr & lvcs,volm_img_info & info,bool load_resource)82 bool volm_io_tools::load_naip_img(std::string const& img_folder, std::string const& name, vpgl_lvcs_sptr& lvcs, volm_img_info& info, bool load_resource)
83 {
84 std::string filename = img_folder + "\\" + name;
85 std::string img_name = filename + "\\" + name + ".tif";
86 std::string tfw_name = filename + "\\" + name + ".tfw";
87 if (!vul_file::exists(tfw_name) || !vul_file::exists(img_name))
88 return false;
89
90 info.name = name; info.img_name = img_name;
91 std::string bbox_file = filename + "\\output_parameters.txt";
92 if (!read_box(bbox_file, info.bbox)) {
93 std::cerr << " cannot find: " << bbox_file << std::endl;
94 return false;
95 }
96 //std::cout << "NAIP bbox: " << info.bbox << std::endl;
97 // figure out utm zone
98 vpgl_utm utm; int utm_zone, zone_max; double xx, yy;
99 utm.transform(info.bbox.min_point().y(), info.bbox.min_point().x(), xx, yy, utm_zone);
100 utm.transform(info.bbox.max_point().y(), info.bbox.max_point().x(), xx, yy, zone_max);
101 if (utm_zone != zone_max) {
102 std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!WARNING! img: " << img_name << " has min and max points in different UTM zones, using zone of min point!\n";
103
104 }
105
106 vpgl_geo_camera *cam = nullptr;
107 if (!vpgl_geo_camera::init_geo_camera(tfw_name, lvcs, utm_zone, northing, cam))
108 return false;
109 info.cam = cam;
110
111 if (load_resource) {
112 vil_image_resource_sptr img = vil_load_image_resource(img_name.c_str());
113 std::cout << "ni: " << img->ni() <<" nj: " << img->nj() <<std::endl;
114 info.ni = img->ni(); info.nj = img->nj();
115 }
116 return true;
117 }
118
load_naip_imgs(std::string const & img_folder,std::vector<volm_img_info> & imgs,bool load_resource)119 bool volm_io_tools::load_naip_imgs(std::string const& img_folder, std::vector<volm_img_info>& imgs, bool load_resource) {
120
121 vpgl_lvcs_sptr lvcs = new vpgl_lvcs; // just the default, no concept of local coordinate system here, so won't be used
122
123 std::string in_dir = img_folder + "*";
124 for (vul_file_iterator fn = in_dir.c_str(); fn; ++fn) {
125 std::string filename = fn();
126 //std::cout << "filename: " << filename << std::endl;
127 std::string file = vul_file::strip_directory(filename);
128
129 volm_img_info info;
130 if (load_naip_img(img_folder, file, lvcs, info, load_resource))
131 imgs.push_back(info);
132 }
133 return true;
134 }
135
load_lidar_img(const std::string & img_file,volm_img_info & info,bool load_image_resource,bool is_cam_global,bool load_cam_from_tfw,std::string const & cam_tfw_file)136 int volm_io_tools::load_lidar_img(const std::string& img_file, volm_img_info& info, bool load_image_resource,
137 bool is_cam_global,
138 bool load_cam_from_tfw, std::string const& cam_tfw_file)
139 {
140 vpgl_lvcs_sptr lvcs = new vpgl_lvcs; // just the default, no concept of local coordinate system here, so won't be used
141
142 if (load_image_resource) {
143 info.img_r = vil_load(img_file.c_str());
144 info.ni = info.img_r->ni(); info.nj = info.img_r->nj();
145 } else {
146 vil_image_view_base_sptr img_sptr = vil_load(img_file.c_str());
147 info.ni = img_sptr->ni(); info.nj = img_sptr->nj();
148 }
149
150 info.name = vul_file::strip_directory(vul_file::strip_extension(img_file));
151 info.img_name = img_file;
152
153 vpgl_geo_camera *cam = nullptr;
154 // try to load camera from image header first
155 vil_image_resource_sptr img_res = vil_load_image_resource(info.img_name.c_str());
156 vpgl_lvcs_sptr lvcs_dummy = new vpgl_lvcs;
157 vpgl_geo_camera::init_geo_camera(img_res, lvcs_dummy, cam);
158
159 if (!cam) {
160 if (!is_cam_global && !load_cam_from_tfw) {
161 vpgl_geo_camera::init_geo_camera(img_file, info.ni, info.nj, lvcs, cam);
162 info.cam = cam;
163 } else if (load_cam_from_tfw) {
164 vpgl_geo_camera::init_geo_camera(cam_tfw_file, lvcs, 0, 0, cam);
165 info.cam = cam;
166 } else if (is_cam_global) {
167 vpgl_geo_camera::init_geo_camera_from_filename(img_file, info.ni, info.nj, lvcs, cam);
168 info.cam = cam;
169 }
170 }
171 else
172 info.cam = cam;
173
174 // obtain the bounding box of current image
175 std::string name = vul_file::strip_directory(img_file);
176 name = name.substr(name.find_first_of('_')+1, name.size());
177
178 std::string n_coords = name.substr(0, name.find_first_of('_'));
179 std::string hemisphere, direction;
180 std::size_t n = n_coords.find('N');
181 if (n < n_coords.size())
182 hemisphere = "N";
183 else
184 hemisphere = "S";
185 n = n_coords.find('E');
186 if (n < n_coords.size())
187 direction = "E";
188 else
189 direction = "W";
190
191 double lat, lon;
192 cam->img_to_global(0.0, info.nj-1, lon, lat);
193 if (!is_cam_global) {
194 if (direction == "W")
195 lon = -lon;
196 if (hemisphere == "S")
197 lat = -lat;
198 }
199 vgl_point_2d<double> lower_left(lon, lat);
200
201 vpgl_utm utm; int utm_zone; double x,y;
202 utm.transform(lat, lon, x, y, utm_zone);
203 std::cout << " zone of lidar img: " << img_file << ": " << utm_zone << " from lower left corner: " << lon << ',' << lat << "!\n";
204
205 cam->img_to_global(info.ni-1, 0.0, lon, lat);
206 if (!is_cam_global) {
207 if (direction == "W")
208 lon = -lon;
209 if (hemisphere == "S")
210 lat = -lat;
211 }
212 vgl_point_2d<double> upper_right(lon, lat);
213 vgl_box_2d<double> bbox(lower_left, upper_right);
214 std::cout << "bbox: " << bbox << std::endl;
215 info.bbox = bbox;
216
217 return utm_zone;
218 }
219
load_lidar_imgs(std::string const & folder,std::vector<volm_img_info> & infos)220 void volm_io_tools::load_lidar_imgs(std::string const& folder, std::vector<volm_img_info>& infos) {
221 std::string in_dir = folder + "/lidar_*.tif";
222 for (vul_file_iterator fn = in_dir.c_str(); fn; ++fn) {
223 std::string filename = fn();
224 volm_img_info info;
225 volm_io_tools::load_lidar_img(filename, info, true, true, false);
226 infos.push_back(info);
227 }
228 }
229
load_nlcd_imgs(std::string const & folder,std::vector<volm_img_info> & infos)230 void volm_io_tools::load_nlcd_imgs(std::string const& folder, std::vector<volm_img_info>& infos)
231 {
232 std::string in_dir = folder + "*.tif*"; // sometimes .tif is written .tiff
233 for (vul_file_iterator fn = in_dir.c_str(); fn; ++fn) {
234 std::string filename = fn();
235 volm_img_info info;
236 load_geotiff_image(filename, info, true);
237 infos.push_back(info);
238 }
239 }
load_imgs(std::string const & folder,std::vector<volm_img_info> & infos,bool load_image_resource,bool is_cam_global,bool load_cam_from_tfw)240 void volm_io_tools::load_imgs(std::string const& folder, std::vector<volm_img_info>& infos, bool load_image_resource, bool is_cam_global, bool load_cam_from_tfw)
241 {
242 std::string in_dir = folder + "/*.tif"; // sometimes .tif is written .tiff
243 for (vul_file_iterator fn = in_dir.c_str(); fn; ++fn) {
244 std::string filename = fn();
245 std::string cam_tfw_file = vul_file::strip_extension(filename) + ".tfw";
246 volm_img_info info;
247 load_lidar_img(filename, info, load_image_resource, is_cam_global, load_cam_from_tfw, cam_tfw_file);
248 infos.push_back(info);
249 }
250 // also load .tiff
251 in_dir = folder + "/*.tiff";
252 for (vul_file_iterator fn = in_dir.c_str(); fn; ++fn) {
253 std::string filename = fn();
254 std::string cam_tfw_file = vul_file::strip_extension(filename) + ".tfw";
255 volm_img_info info;
256 load_lidar_img(filename, info, load_image_resource, is_cam_global, load_cam_from_tfw, cam_tfw_file);
257 infos.push_back(info);
258 }
259 }
260
get_location_nlcd(std::vector<volm_img_info> & NLCD_imgs,double lat,double lon,double elev,unsigned char & label)261 bool volm_io_tools::get_location_nlcd(std::vector<volm_img_info>& NLCD_imgs, double lat, double lon, double elev, unsigned char& label)
262 {
263 bool found_it = false;
264 for (auto & NLCD_img : NLCD_imgs) {
265 if (NLCD_img.bbox.contains(lon, lat)) {
266 vil_image_view<vxl_byte> img(NLCD_img.img_r);
267
268 // get the land type of the location
269 double u, v;
270 NLCD_img.cam->global_to_img(lon, lat, elev, u, v);
271 //NLCD_imgs[i].cam->global_to_img(-lon, lat, elev, u, v);
272 auto uu = (unsigned)std::floor(u + 0.5);
273 auto vv = (unsigned)std::floor(v + 0.5);
274 if (uu > 0 && vv > 0 && uu < NLCD_img.ni && vv < NLCD_img.nj) {
275 label = img(uu, vv);
276 found_it = true;
277 break;
278 }
279 }
280 }
281 return found_it;
282 }
283
expend_line(std::vector<vgl_point_2d<double>> line,double const & w,vgl_polygon<double> & poly)284 bool volm_io_tools::expend_line(std::vector<vgl_point_2d<double> > line, double const& w, vgl_polygon<double>& poly)
285 {
286 std::vector<vgl_point_2d<double> > sheet;
287 std::vector<vgl_point_2d<double> > pts_u;
288 std::vector<vgl_point_2d<double> > pts_d;
289 unsigned n_pts = line.size();
290 for (unsigned i = 0; i < n_pts-1; i++) {
291 vgl_point_2d<double> s, e;
292 s = line[i]; e = line[i+1];
293 if (near_eq_pt(s,e))
294 continue;
295 vgl_line_2d<double> seg(s, e);
296 vgl_vector_2d<double> n = seg.normal();
297 vgl_point_2d<double> su, sd, eu, ed;
298 su = s + 0.5*w*n; sd = s - 0.5*w*n;
299 pts_u.push_back(su); pts_d.push_back(sd);
300 if (i == n_pts-2) {
301 eu = e + 0.5*w*n; ed = e - 0.5*w*n;
302 pts_u.push_back(eu); pts_d.push_back(ed);
303 }
304 }
305 // rearrange the point list
306 sheet.reserve(pts_u.size());
307 for (auto i : pts_u)
308 sheet.push_back(i);
309 for (int i = pts_d.size()-1; i >=0; i--)
310 sheet.push_back(pts_d[i]);
311 poly.push_back(sheet);
312 return true;
313 }
314
315
316 #include <volm/volm_geo_index2.h>
317 #include <volm/volm_osm_objects.h>
318
319 // a method to read the binary osm object file and also contstruct the volm_geo_index2, the method returns the root of the tree
read_osm_data_and_tree(const std::string & geoindex_filename_pre,const std::string & osm_bin_filename,volm_osm_objects & osm_objs,double & min_size)320 volm_geo_index2_node_sptr volm_io_tools::read_osm_data_and_tree(const std::string& geoindex_filename_pre, const std::string& osm_bin_filename, volm_osm_objects& osm_objs, double& min_size)
321 {
322 std::string filename = geoindex_filename_pre + ".txt";
323 volm_geo_index2_node_sptr root = volm_geo_index2::read_and_construct<volm_osm_object_ids_sptr>(filename, min_size);
324 // obtain all leaves
325 std::vector<volm_geo_index2_node_sptr> leaves;
326 volm_geo_index2::get_leaves(root, leaves);
327 // load the content for valid leaves
328 for (auto & leave : leaves) {
329 std::string bin_file = leave->get_label_name(geoindex_filename_pre, "osm");
330 if (!vul_file::exists(bin_file))
331 continue;
332 auto* ptr = dynamic_cast<volm_geo_index2_node<volm_osm_object_ids_sptr>* >(leave.ptr());
333 ptr->contents_ = new volm_osm_object_ids(bin_file);
334 }
335
336 // load the osm bin file to get real location date with unit of lon and lat, associated with ids stored in geo_index2
337 if (!vul_file::exists(osm_bin_filename)) {
338 std::cout << "ERROR: can not find osm binary: " << osm_bin_filename << std::endl;
339 return nullptr;
340 }
341
342 vsl_b_ifstream is(osm_bin_filename.c_str());
343 if (!is) {
344 std::cerr << "In volm_osm_object::volm_osm_object() -- cannot open: " << osm_bin_filename << std::endl;
345 return nullptr;
346 }
347 osm_objs.b_read(is);
348 is.close();
349
350 #if 0
351 if (is_kml()) {
352 std::stringstream kml_pts, kml_roads, kml_regions;
353 kml_pts << out_pre() << "/p1b_wr" << world_id() << "_tile_" << tile_id() << "_osm_pts.kml";
354 kml_roads << out_pre() << "/p1b_wr" << world_id() << "_tile_" << tile_id() << "_osm_roads.kml";
355 kml_regions << out_pre() << "/p1b_wr" << world_id() << "_tile_" << tile_id() << "_osm_regions.kml";
356 osm_objs.write_pts_to_kml(kml_pts.str());
357 osm_objs.write_lines_to_kml(kml_roads.str());
358 osm_objs.write_polys_to_kml(kml_regions.str());
359 }
360 #endif
361
362 return root;
363 }
364
365 //: load a geotiff file with .tif file and read its ortho camera info from its header, puts a dummy lvcs to vpgl_geo_cam object so lvcs is not valid
366 // even though it reads the camera from filename with N/S and W/E distinction, it constructs a camera in global WGS84 coordinates, so the global coordinates should be used to fetch pixels, (i.e. no need to make them always positive)
load_geotiff_image(const std::string & filename,volm_img_info & info,bool load_cam_from_name)367 void volm_io_tools::load_geotiff_image(const std::string& filename, volm_img_info& info, bool load_cam_from_name)
368 {
369 info.img_name = filename;
370 info.name = vul_file::strip_directory(info.img_name);
371 info.name = vul_file::strip_extension(info.name);
372
373 info.img_r = vil_load(info.img_name.c_str());
374 info.ni = info.img_r->ni(); info.nj = info.img_r->nj();
375
376 vpgl_geo_camera *cam;
377 vpgl_lvcs_sptr lvcs_dummy = new vpgl_lvcs;
378 if (load_cam_from_name) {
379 vpgl_geo_camera::init_geo_camera_from_filename(filename, info.ni, info.nj, lvcs_dummy, cam); // constructs in global WGS84 (no distinction of N/S or W/E)
380 std::cout << cam->trans_matrix() << std::endl;
381 } else {
382 vil_image_resource_sptr img_res = vil_load_image_resource(info.img_name.c_str());
383 vpgl_geo_camera::init_geo_camera(img_res, lvcs_dummy, cam);
384 }
385
386 info.cam = cam;
387
388 double lat, lon;
389 cam->img_to_global(0.0, info.nj-1, lon, lat);
390 vgl_point_2d<double> lower_left(lon, lat);
391
392 vpgl_utm utm; int utm_zone; double x,y;
393 utm.transform(lat, lon, x, y, utm_zone);
394 cam->img_to_global(info.ni-1, 0.0, lon, lat);
395 vgl_point_2d<double> upper_right(lon, lat);
396 vgl_box_2d<double> bbox(lower_left, upper_right);
397 info.bbox = bbox;
398 }
399
load_aster_dem_imgs(std::string const & folder,std::vector<volm_img_info> & infos)400 void volm_io_tools::load_aster_dem_imgs(std::string const& folder, std::vector<volm_img_info>& infos)
401 {
402 std::string file_glob = folder + "//*.tif";
403 for (vul_file_iterator fn = file_glob.c_str(); fn; ++fn) {
404 volm_img_info info;
405 volm_io_tools::load_geotiff_image(fn(), info);
406 infos.push_back(info);
407 //volm_img_info info;
408 //std::string folder = fn();
409 //std::string file_glob2 = folder + "//" + "*_dem*.tif";
410 //for (vul_file_iterator fn2 = file_glob2.c_str(); fn2; ++fn2) {
411 // volm_io_tools::load_geotiff_image(fn2(), info);
412 // infos.push_back(info);
413 //}
414 }
415 }
416
load_satellite_height_map(std::string const & filename,volm_img_info & info,bool const & load_cam_from_file)417 bool volm_io_tools::load_satellite_height_map(std::string const& filename, volm_img_info& info, bool const& load_cam_from_file)
418 {
419 info.img_name = filename;
420 info.name = vul_file::strip_directory(info.img_name);
421 info.name = vul_file::strip_extension(info.name);
422
423 info.img_r = vil_load(info.img_name.c_str());
424 info.ni = info.img_r->ni(); info.nj = info.img_r->nj();
425 std::cout << "satellite height image ni: " << info.ni << " nj: " << info.nj << std::endl;
426 // load the camera (either from a tfw file or from image header)
427 vpgl_geo_camera *cam;
428 vpgl_lvcs_sptr lvcs_dummy = new vpgl_lvcs;
429 if (load_cam_from_file) {
430 std::string cam_file = vul_file::dirname(filename) + "/" + info.name + "_geo.tfw";
431 if (!vul_file::exists(cam_file)) {
432 std::cerr << " can not find camera file: " << cam_file << "!\n";
433 return false;
434 }
435 vpgl_geo_camera::init_geo_camera(cam_file, lvcs_dummy, 0, 0, cam);
436 }
437 else { // load from image header
438 vil_image_resource_sptr img_res = vil_load_image_resource(info.img_name.c_str());
439 vpgl_geo_camera::init_geo_camera(img_res, lvcs_dummy, cam);
440 }
441
442 // defined the bbox of the image, use lower_left corner to define utm and northing of camera
443 double lat, lon;
444 cam->img_to_global(0.0, info.nj-1, lon, lat);
445 vgl_point_2d<double> lower_left(lon, lat);
446 vpgl_utm utm;
447 int utm_zone; double x, y;
448 utm.transform(lat, lon, x, y, utm_zone);
449 std::cout << " zone of satellite height image " << info.name << ": " << utm_zone << " from lower left corner!\n";
450 cam->img_to_global(info.ni-1,0.0, lon, lat);
451 vgl_point_2d<double> upper_right(lon, lat);
452 #if 0
453 int utm_zone_ur;
454 utm.transform(lat, lon, x, y, utm_zone_ur);
455 if (utm_zone_ur != utm_zone)
456 std::cout << " warning: satellite height image " << info.name << " crosses two utm zone: " << utm_zone << ", " << utm_zone_ur << std::endl;
457 unsigned northing = 1;
458 if (upper_right.x() < 0 && lower_left.x() < 0)
459 northing = 0;
460 if (upper_right.x() * lower_left.x() < 0)
461 std::cout << " warning: satellite height image " << info.name << " crosses the Equator (set it to be northing)" << std::endl;
462 cam->set_utm(utm_zone, northing);
463 #endif
464 info.cam = cam;
465 vgl_box_2d<double> bbox(lower_left, upper_right);
466 std::cout << "bbox: " << bbox << std::endl;
467 info.bbox = bbox;
468
469 return true;
470 }
471
load_satellite_height_imgs(std::string const & folder,std::vector<volm_img_info> & infos,bool const & load_cam_from_file,std::string const & keywords)472 bool volm_io_tools::load_satellite_height_imgs(std::string const& folder, std::vector<volm_img_info>& infos, bool const& load_cam_from_file, std::string const& keywords)
473 {
474 std::string file_glob;
475 if (keywords.compare("") == 0)
476 file_glob = folder + "//scene_*.tif";
477 else
478 file_glob = folder + "//scene_*" + keywords + ".tif";
479 for (vul_file_iterator fn = file_glob.c_str(); fn; ++fn) {
480 volm_img_info info;
481 if (!volm_io_tools::load_satellite_height_map(fn(), info, load_cam_from_file))
482 return false;
483 infos.push_back(info);
484 }
485 return true;
486 }
487
load_geocover_imgs(std::string const & folder,std::vector<volm_img_info> & infos)488 void volm_io_tools::load_geocover_imgs(std::string const& folder, std::vector<volm_img_info>& infos)
489 {
490 std::string file_glob = folder + "//*.tif";
491 for (vul_file_iterator fn = file_glob.c_str(); fn; ++fn) {
492 volm_img_info info;
493 volm_io_tools::load_geotiff_image(fn(), info, false); // last argument true so load camera from the file name
494 infos.push_back(info);
495 }
496 }
497
load_urban_imgs(std::string const & folder,std::vector<volm_img_info> & infos)498 void volm_io_tools::load_urban_imgs(std::string const& folder, std::vector<volm_img_info>& infos)
499 {
500 std::string file_glob = folder + "//Urextent_*.tif";
501 for (vul_file_iterator fn = file_glob.c_str(); fn; ++fn) {
502 volm_img_info info;
503 volm_io_tools::load_geotiff_image(fn(), info, true);
504 infos.push_back(info);
505 }
506 }
507
crop_and_find_min_max(std::vector<volm_img_info> & infos,unsigned img_id,int i0,int j0,int crop_ni,int crop_nj,double & min,double & max)508 void crop_and_find_min_max(std::vector<volm_img_info>& infos, unsigned img_id, int i0, int j0, int crop_ni, int crop_nj, double& min, double& max)
509 {
510 #if 0
511 vil_image_view<vxl_int_16> img(infos[img_id].img_r);
512 vil_image_view<vxl_int_16> img_crop = vil_crop(img, i0, crop_ni, j0, crop_nj);
513 for (unsigned ii = 0; ii < img_crop.ni(); ii++)
514 for (unsigned jj = 0; jj < img_crop.nj(); jj++) {
515 if (min > img_crop(ii, jj)) min = img_crop(ii, jj);
516 if (max < img_crop(ii, jj)) max = img_crop(ii, jj);
517 }
518 #endif
519 if (auto* img = dynamic_cast<vil_image_view<vxl_int_16>*>(infos[img_id].img_r.ptr())) {
520 vil_image_view<vxl_int_16> img_crop = vil_crop(*img, i0, crop_ni, j0, crop_nj);
521 for (unsigned ii = 0; ii < img_crop.ni(); ii++)
522 for (unsigned jj = 0; jj < img_crop.nj(); jj++) {
523 if (min > img_crop(ii, jj)) min = img_crop(ii, jj);
524 if (max < img_crop(ii, jj)) max = img_crop(ii, jj);
525 }
526 }
527 else if (auto* img = dynamic_cast<vil_image_view<float>*>(infos[img_id].img_r.ptr()) ) {
528 vil_image_view<float> img_crop = vil_crop(*img, i0, crop_ni, j0, crop_nj);
529 for (unsigned ii = 0; ii < img_crop.ni(); ii++)
530 for (unsigned jj = 0; jj < img_crop.nj(); jj++) {
531 if (min > img_crop(ii, jj)) min = img_crop(ii, jj);
532 if (max < img_crop(ii, jj)) max = img_crop(ii, jj);
533 }
534 }
535
536
537 }
538
find_min_max_height(vgl_point_2d<double> & lower_left,vgl_point_2d<double> & upper_right,std::vector<volm_img_info> & infos,double & min,double & max)539 bool volm_io_tools::find_min_max_height(vgl_point_2d<double>& lower_left, vgl_point_2d<double>& upper_right, std::vector<volm_img_info>& infos, double& min, double& max)
540 {
541 // find the image of all four corners
542 std::vector<std::pair<unsigned, std::pair<int, int> > > corners;
543 std::vector<vgl_point_2d<double> > pts;
544 pts.emplace_back(lower_left.x(), upper_right.y());
545 pts.emplace_back(upper_right.x(), lower_left.y());
546 pts.push_back(lower_left);
547 pts.push_back(upper_right);
548
549 for (auto & pt : pts) {
550 // find the image
551 for (unsigned j = 0; j < (unsigned)infos.size(); j++) {
552 double u, v;
553 infos[j].cam->global_to_img(pt.x(), pt.y(), 0, u, v);
554 int uu = (int)std::floor(u+0.5);
555 int vv = (int)std::floor(v+0.5);
556 if (uu < 0 || vv < 0 || uu >= (int)infos[j].ni || vv >= (int)infos[j].nj)
557 continue;
558 std::pair<unsigned, std::pair<int, int> > pp(j, std::pair<int, int>(uu, vv));
559 corners.push_back(pp);
560 break;
561 }
562 }
563 if (corners.size() != 4) {
564 std::cerr << "Cannot locate all 4 corners among these DEM tiles!\n";
565 return false;
566 }
567 // case 1: all corners are in the same image
568 if (corners[0].first == corners[1].first) {
569 // crop the image
570 int i0 = corners[0].second.first;
571 int j0 = corners[0].second.second;
572 int crop_ni = corners[1].second.first-corners[0].second.first+1;
573 int crop_nj = corners[1].second.second-corners[0].second.second+1;
574 crop_and_find_min_max(infos, corners[0].first, i0, j0, crop_ni, crop_nj, min, max);
575 return true;
576 }
577 // case 2: two corners are in the same image
578 if (corners[0].first == corners[2].first && corners[1].first == corners[3].first) {
579 // crop the first image
580 int i0 = corners[0].second.first;
581 int j0 = corners[0].second.second;
582 int crop_ni = infos[corners[0].first].ni - corners[0].second.first;
583 int crop_nj = corners[2].second.second-corners[0].second.second+1;
584 crop_and_find_min_max(infos, corners[0].first, i0, j0, crop_ni, crop_nj, min, max);
585
586 // crop the second image
587 i0 = 0;
588 j0 = corners[3].second.second;
589 crop_ni = corners[3].second.first + 1;
590 crop_nj = corners[1].second.second-corners[3].second.second+1;
591 crop_and_find_min_max(infos, corners[1].first, i0, j0, crop_ni, crop_nj, min, max);
592 return true;
593 }
594 // case 3: two corners are in the same image
595 if (corners[0].first == corners[3].first && corners[1].first == corners[2].first) {
596 // crop the first image
597 int i0 = corners[0].second.first;
598 int j0 = corners[0].second.second;
599 int crop_ni = corners[3].second.first - corners[0].second.first + 1;
600 int crop_nj = infos[corners[0].first].nj - corners[0].second.second;
601 crop_and_find_min_max(infos, corners[0].first, i0, j0, crop_ni, crop_nj, min, max);
602
603 // crop the second image
604 i0 = corners[2].second.first;
605 j0 = 0;
606 crop_ni = corners[1].second.first - corners[2].second.first + 1;
607 crop_nj = corners[2].second.second + 1;
608 crop_and_find_min_max(infos, corners[1].first, i0, j0, crop_ni, crop_nj, min, max);
609 return true;
610 }
611 // case 4: all corners are in a different image
612 // crop the first image, image of corner 0
613 int i0 = corners[0].second.first;
614 int j0 = corners[0].second.second;
615 int crop_ni = infos[corners[0].first].ni - corners[0].second.first;
616 int crop_nj = infos[corners[0].first].nj - corners[0].second.second;
617 crop_and_find_min_max(infos, corners[0].first, i0, j0, crop_ni, crop_nj, min, max);
618
619 // crop the second image, image of corner 1
620 i0 = 0;
621 j0 = 0;
622 crop_ni = corners[1].second.first + 1;
623 crop_nj = corners[1].second.second + 1;
624 crop_and_find_min_max(infos, corners[1].first, i0, j0, crop_ni, crop_nj, min, max);
625
626 // crop the third image, image of corner 2
627 i0 = corners[2].second.first;
628 j0 = 0;
629 crop_ni = infos[corners[2].first].ni - corners[2].second.first;
630 crop_nj = corners[2].second.second + 1;
631 crop_and_find_min_max(infos, corners[2].first, i0, j0, crop_ni, crop_nj, min, max);
632
633 // crop the fourth image, image of corner 3
634 i0 = 0;
635 j0 = corners[3].second.second;
636 crop_ni = corners[3].second.first + 1;
637 crop_nj = infos[corners[3].first].nj - corners[3].second.second;
638 crop_and_find_min_max(infos, corners[3].first, i0, j0, crop_ni, crop_nj, min, max);
639
640 return true;
641 }
642
643 //: use the following method to get the multiplier for conversion of meters to degrees, uses vpgl_lvcs internally
meter_to_seconds(double lat,double lon)644 double volm_io_tools::meter_to_seconds(double lat, double lon)
645 {
646 vpgl_lvcs_sptr lvcs = new vpgl_lvcs(lat, lon, 0, vpgl_lvcs::wgs84, vpgl_lvcs::DEG, vpgl_lvcs::METERS);
647 double lon1, lat1, lon2, lat2, gz;
648 lvcs->local_to_global(0, 0, 0, vpgl_lvcs::wgs84, lon1, lat1, gz);
649 lvcs->local_to_global(1, 0, 0, vpgl_lvcs::wgs84, lon2, lat2, gz);
650 double dif_lon = lon2-lon1;
651 double dif_lat = lat2-lat1;
652 double dif = std::sqrt(dif_lon*dif_lon + dif_lat*dif_lat); // 1 meter is this many degrees in this area
653 return dif;
654 }
655
find_intersect(vgl_box_2d<double> const & bbox,vgl_point_2d<double> const & s,vgl_point_2d<double> e,vgl_point_2d<double> & pt)656 bool find_intersect(vgl_box_2d<double> const& bbox, vgl_point_2d<double> const& s, vgl_point_2d<double> e, vgl_point_2d<double>& pt)
657 {
658 vgl_line_2d<double> line(s, e);
659 vgl_point_2d<double> pi0, pi1;
660 if (!vgl_intersection(bbox, line, pi0, pi1))
661 return false;
662 double x1 = s.x(), y1 = s.y();
663 double x2 = e.x(), y2 = e.y();
664 double xp = pi0.x(), yp = pi0.y();
665 double d1p = (xp-x1)*(xp-x1) + (yp-y1)*(yp-y1);
666 double d2p = (xp-x2)*(xp-x2) + (yp-y2)*(yp-y2);
667 double d12 = (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1);
668 double diff = std::sqrt(d1p) + std::sqrt(d2p) - std::sqrt(d12);
669 if (diff < 1E-5) {
670 pt = pi0;
671 return true;
672 }
673 xp = pi1.x(); yp = pi1.y();
674 d1p = (xp-x1)*(xp-x1) + (yp-y1)*(yp-y1);
675 d2p = (xp-x2)*(xp-x2) + (yp-y2)*(yp-y2);
676 d12 = (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1);
677 diff = std::sqrt(d1p) + std::sqrt(d2p) - std::sqrt(d12);
678 if (diff < 1E-5) {
679 pt = pi1;
680 return true;
681 }
682 return false;
683 }
684
line_inside_the_box(vgl_box_2d<double> const & bbox,std::vector<vgl_point_2d<double>> const & line,std::vector<vgl_point_2d<double>> & road)685 bool volm_io_tools::line_inside_the_box(vgl_box_2d<double> const& bbox, std::vector<vgl_point_2d<double> > const& line, std::vector<vgl_point_2d<double> >& road)
686 {
687 // obtain points that lie inside the bounding box
688 std::vector<vgl_point_2d<double> > line_in = vgl_intersection(line, bbox);
689 if (line_in.empty())
690 return false;
691 for (auto i : line_in)
692 road.push_back(i);
693
694 // find the intersection points
695 for (auto curr_pt : line_in) {
696 auto vit = std::find(line.begin(), line.end(), curr_pt);
697 if (vit == line.begin() ) {
698 vgl_point_2d<double> next = *(vit+1);
699 if (bbox.contains(next))
700 continue;
701 else {
702 vgl_point_2d<double> intersect_pt;
703 if (!find_intersect(bbox, curr_pt, next, intersect_pt))
704 return false;
705 // insert the intersect after current point
706 auto it = std::find(road.begin(), road.end(), curr_pt);
707 road.insert(it+1, intersect_pt);
708 }
709 }
710 else if (vit == line.end()-1) {
711 vgl_point_2d<double> prev = *(vit-1);
712 if (bbox.contains(prev))
713 continue;
714 else {
715 vgl_point_2d<double> intersect_pt;
716 if (!find_intersect(bbox, prev, curr_pt,intersect_pt))
717 return false;
718 // insert the intersect point before current point
719 auto it = std::find(road.begin(), road.end(), curr_pt);
720 road.insert(it, intersect_pt);
721 }
722 }
723 else if (vit != line.end()) {
724 vgl_point_2d<double> prev = *(vit-1);
725 vgl_point_2d<double> next = *(vit+1);
726 if (bbox.contains(next) && bbox.contains(prev))
727 continue;
728 else if (bbox.contains(next)) {
729 vgl_point_2d<double> intersect_pt;
730 if (!find_intersect(bbox, prev, curr_pt, intersect_pt))
731 return false;
732 auto it = std::find(road.begin(), road.end(), curr_pt);
733 road.insert(it, intersect_pt);
734 }
735 else if (bbox.contains(prev)) {
736 vgl_point_2d<double> intersect_pt;
737 if (!find_intersect(bbox, curr_pt, next, intersect_pt))
738 return false;
739 // insert the intersect after current point
740 auto it = std::find(road.begin(), road.end(), curr_pt);
741 road.insert(it+1, intersect_pt);
742 }
743 else {
744 vgl_point_2d<double> intersect_pt;
745 if (!find_intersect(bbox, curr_pt, next, intersect_pt))
746 return false;
747 // find and insert the intersection after current point
748 auto it = std::find(road.begin(), road.end(), curr_pt);
749 road.insert(it+1, intersect_pt);
750 // find and insert the intersection before current point
751 if (!find_intersect(bbox, prev, curr_pt,intersect_pt))
752 return false;
753 it = std::find(road.begin(), road.end(), curr_pt);
754 road.insert(it, intersect_pt);
755 }
756 }
757 else
758 return false;
759 }
760 return true;
761 }
762
form_line_segment_from_pts(std::vector<vgl_point_2d<double>> const & road,std::vector<vgl_line_segment_2d<double>> & road_seg)763 void form_line_segment_from_pts(std::vector<vgl_point_2d<double> > const& road, std::vector<vgl_line_segment_2d<double> >& road_seg)
764 {
765 unsigned num_pts = road.size();
766 unsigned num_segs = num_pts - 1;
767 for (unsigned i = 0; i < num_segs; i++) {
768 vgl_point_2d<double> s = road[i]; vgl_point_2d<double> e = road[i+1];
769 if (near_eq_pt(s,e))
770 continue;
771 road_seg.emplace_back(s, e);
772 }
773 #if 0
774 // define a 2d line fit
775 unsigned min_pts = 3;
776 double tol = 5.0; // in pixel unit
777 vgl_fit_lines_2d<double> fitter(min_pts, tol);
778
779 unsigned num_pts = road.size();
780 unsigned num_segs = num_pts - 1;
781 if (num_pts <= min_pts) {
782 for (unsigned i = 0; i < num_segs; i++) {
783 vgl_point_2d<double> s = road[i]; vgl_point_2d<double> e = road[i+1];
784 road_seg.push_back(vgl_line_segment_2d<double>(s, e));
785 }
786 }
787 else {
788 fitter.add_curve(road);
789 fitter.fit();
790 road_seg = fitter.get_line_segs();
791 if (road_seg.empty()) {
792 // fitting failed form the segment directly
793 road_seg.clear();
794 for (unsigned i = 0; i < num_segs; i++) {
795 vgl_point_2d<double> s = road[i]; vgl_point_2d<double> e = road[i+1];
796 road_seg.push_back(vgl_line_segment_2d<double>(s, e));
797 }
798 }
799 }
800 fitter.clear();
801 #endif
802 }
803
find_junctions(vgl_line_segment_2d<double> const & seg,volm_land_layer const & seg_prop,std::vector<vgl_line_segment_2d<double>> const & lines,volm_land_layer const & line_prop,std::vector<vgl_point_2d<double>> & cross_pts,std::vector<volm_land_layer> & cross_prop)804 void find_junctions(vgl_line_segment_2d<double> const& seg,
805 volm_land_layer const& seg_prop,
806 std::vector<vgl_line_segment_2d<double> > const& lines,
807 volm_land_layer const& line_prop,
808 std::vector<vgl_point_2d<double> >& cross_pts,
809 std::vector<volm_land_layer>& cross_prop)
810 {
811 vgl_line_segment_3d<double> l1(vgl_point_3d<double>(seg.point1().x(), seg.point1().y(), 0.0), vgl_point_3d<double>(seg.point2().x(), seg.point2().y(), 0.0));
812
813 unsigned n_seg = lines.size();
814 for (unsigned i = 0; i < n_seg; i++) {
815 vgl_line_segment_3d<double> l2(vgl_point_3d<double>(lines[i].point1().x(), lines[i].point1().y(), 0.0),
816 vgl_point_3d<double>(lines[i].point2().x(), lines[i].point2().y(), 0.0));
817 vgl_point_3d<double> pt;
818 // check intersection
819 if (!vgl_intersection(l1, l2, pt)) {
820 // check whether connect at end point
821 if (near_eq_pt(seg.point1(), lines[i].point1()))
822 pt.set(seg.point1().x(), seg.point1().y(), 0.0);
823 else if (near_eq_pt(seg.point1(), lines[i].point2()))
824 pt.set(seg.point1().x(), seg.point1().y(), 0.0);
825 else if (near_eq_pt(seg.point2(), lines[i].point1()))
826 pt.set(seg.point2().x(), seg.point2().y(), 0.0);
827 else if (near_eq_pt(seg.point2(), lines[i].point2()))
828 pt.set(seg.point2().x(), seg.point2().y(), 0.0);
829 else
830 continue;
831 }
832 cross_pts.emplace_back(pt.x(), pt.y());
833 std::pair<int,int> key(seg_prop.id_, line_prop.id_);
834 cross_prop.push_back(volm_osm_category_io::road_junction_table[key]);
835 }
836 }
837
count_line_start_from_cross(vgl_point_2d<double> const & cross_pt,std::vector<vgl_point_2d<double>> const & rd,std::vector<std::vector<vgl_point_2d<double>>> const & net)838 unsigned count_line_start_from_cross(vgl_point_2d<double> const& cross_pt,
839 std::vector<vgl_point_2d<double> > const& rd,
840 std::vector<std::vector<vgl_point_2d<double> > > const& net)
841 {
842 unsigned cnt = 0;
843 // check current road first
844 vgl_point_2d<double> s = *(rd.begin()); vgl_point_2d<double> e = *(rd.end()-1);
845 if ( near_eq_pt(cross_pt, s) || near_eq_pt(cross_pt, e))
846 cnt++;
847 for (const auto & i : net) {
848 s = *(i.begin()); e = *(i.end()-1);
849 if ( near_eq_pt(cross_pt, s) || near_eq_pt(cross_pt, e))
850 cnt++;
851 }
852 return cnt;
853 }
854
search_junctions(std::vector<vgl_point_2d<double>> const & road,volm_land_layer const & road_prop,std::vector<std::vector<vgl_point_2d<double>>> net,std::vector<volm_land_layer> net_props,std::vector<vgl_point_2d<double>> & cross_pts,std::vector<volm_land_layer> & cross_props,std::vector<volm_land_layer> & cross_geo_props)855 bool volm_io_tools::search_junctions(std::vector<vgl_point_2d<double> > const& road, volm_land_layer const& road_prop,
856 std::vector<std::vector<vgl_point_2d<double> > > net, std::vector<volm_land_layer> net_props,
857 std::vector<vgl_point_2d<double> >& cross_pts, std::vector<volm_land_layer>& cross_props,
858 std::vector<volm_land_layer>& cross_geo_props)
859 {
860 unsigned n_rds = net.size();
861
862 // form the line segment for each road in the network
863 std::vector<vgl_line_segment_2d<double> > road_seg;
864 form_line_segment_from_pts(road, road_seg);
865
866 std::vector<std::vector<vgl_line_segment_2d<double> > > net_segs;
867 for (unsigned r_idx = 0; r_idx < n_rds; r_idx++) {
868 std::vector<vgl_line_segment_2d<double> > seg;
869 form_line_segment_from_pts(net[r_idx], seg);
870 net_segs.push_back(seg);
871 }
872
873 // find the cross for each segment
874 unsigned n_seg = road_seg.size();
875 for (unsigned s_idx = 0; s_idx < n_seg; s_idx++) {
876 vgl_line_segment_2d<double> curr_seg = road_seg[s_idx];
877 for (unsigned r_idx = 0; r_idx < n_rds; r_idx++) {
878 std::vector<vgl_line_segment_2d<double> > curr_net_seg = net_segs[r_idx];
879 std::vector<vgl_point_2d<double> > pt;
880 std::vector<volm_land_layer> prop;
881 find_junctions(curr_seg, road_prop, curr_net_seg, net_props[r_idx], pt, prop);
882 if (pt.empty())
883 continue;
884 if (prop.size() != pt.size())
885 return false;
886 for (unsigned p_idx = 0; p_idx < pt.size(); p_idx++)
887 if (std::find(cross_pts.begin(), cross_pts.end(), pt[p_idx]) == cross_pts.end()) {
888 cross_pts.push_back(pt[p_idx]); cross_props.push_back(prop[p_idx]); cross_geo_props.push_back(volm_osm_category_io::volm_land_table_name["4_way"]);
889 }
890 }
891 }
892
893 //check whether the cross pt are T_section
894 //principle, count the number of roads whose end points are on the cross_pt
895 //if there is 1 or 3 lines start from this cross pt, the cross is a T_section
896 for (unsigned c_idx = 0; c_idx < cross_pts.size(); c_idx++) {
897 unsigned num_lines = count_line_start_from_cross(cross_pts[c_idx], road, net);
898 if (num_lines == 1 || num_lines == 3) {
899 cross_props[c_idx] = volm_osm_category_io::volm_land_table_name["T_section"];
900 }
901 }
902 // define the geometric properties of intersections
903 for (unsigned c_idx = 0; c_idx < cross_pts.size(); c_idx++) {
904 unsigned num_lines = count_line_start_from_cross(cross_pts[c_idx], road, net);
905 if (num_lines == 1 || num_lines == 3) {
906 cross_geo_props[c_idx] = volm_osm_category_io::volm_land_table_name["T_section"];
907 }
908 if (num_lines == 2) {
909 cross_geo_props[c_idx] = volm_osm_category_io::volm_land_table_name["2_way"];
910 }
911 }
912
913 return true;
914 }
915