1 //:
2 // \file
3 // \brief executable to classify all the pixels in a single lidar image,
4 // finds all the aerial images which overlap this image to create filter responses
5 // run this exe in parallel to classify all lidar images
6 // \author Ozge C. Ozcanli
7 // \date Mar 01, 2013
8
9 #include <volm/volm_io.h>
10 #include <volm/volm_io_tools.h>
11 #include <volm/volm_tile.h>
12 #include <volm/volm_loc_hyp.h>
13 #include "vul/vul_arg.h"
14 #include "vul/vul_file.h"
15 #include "vul/vul_file_iterator.h"
16 #include "vil/vil_image_view.h"
17 #include "vil/vil_load.h"
18 #include "vil/vil_save.h"
19 #include "vil/vil_new.h"
20 #include "vil/vil_convert.h"
21 #include "vpgl/vpgl_utm.h"
22 #include "vpgl/vpgl_lvcs.h"
23 #include <bkml/bkml_parser.h>
24 #include "vul/vul_timer.h"
25 #include "vgl/vgl_intersection.h"
26 #include <bkml/bkml_write.h>
27 #include <volm/volm_geo_index.h>
28
29 #include <sdet/sdet_texture_classifier.h>
30 #include <sdet/sdet_texture_classifier_params.h>
31
32
33 // no need to check cause get data checks it anyway
get_block(int i,int j,int bb,std::vector<std::pair<int,int>> & pixels)34 void get_block(int i, int j, int bb, std::vector<std::pair<int, int> >& pixels)
35 {
36 for (int ii = i-bb; ii < i+bb; ii++)
37 for (int jj = j-bb; jj < j+bb; jj++) {
38 pixels.emplace_back(ii, jj);
39 }
40 }
41
42 // read the tiles of the region, create a geo index and write the hyps
main(int argc,char ** argv)43 int main(int argc, char** argv)
44 {
45 vul_timer t;
46 t.mark();
47
48 vul_arg<std::string> img_folder("-imgs", "input folder to read image files", "");
49 vul_arg<std::string> lidar_file("-lidar", "input file to read", "");
50 vul_arg<float> lidar_max("-max", "max lidar value to stretch lidar images", 80.0f);
51 vul_arg<std::string> filter_folder("-filter_dir", "folder to read/write filter responses", "");
52 vul_arg<std::string> out_folder("-out", "folder to write output to", "");
53 vul_arg<int> neighborhood_size("-n", "size of the neighborhood around a pixel to collect training examples, in pixels, e.g. 5", 5);
54 vul_arg<std::string> test("-test", "load the specified classifier instance and classify pixels in images in imgs folder specified in test_imgs.txt file", "");
55 vul_arg<int> block("-block", "size of the block around a pixel to compute its set of filter vectors", 5);
56 vul_arg_parse(argc, argv);
57 std::cout << "argc: " << argc << std::endl;
58
59 if (filter_folder().compare("") == 0) {
60 std::cout << " filter folder is not specified! exiting..\n"; return 0;
61 }
62
63 std::vector<volm_img_info> imgs;
64 volm_io_tools::load_naip_imgs(img_folder(), imgs);
65
66 volm_img_info lidar_info;
67 int utm_zone = volm_io_tools::load_lidar_img(lidar_file(), lidar_info);
68 lidar_info.save_box_kml(out_folder() + lidar_info.name + "_box.kml");
69
70 vpgl_lvcs_sptr lvcs;
71 if (utm_zone == 17)
72 lvcs = new vpgl_lvcs(30.371923, -81.539398, -8.0, vpgl_lvcs::utm, vpgl_lvcs::DEG, vpgl_lvcs::METERS);
73 else if (utm_zone == 18)
74 lvcs = new vpgl_lvcs(35.289477, -75.570596, -8.0, vpgl_lvcs::utm, vpgl_lvcs::DEG, vpgl_lvcs::METERS);
75
76 // find the images that it intersects
77 std::vector<volm_img_info> intersection_imgs;
78 std::cout << "checking intersection of bbox: " << lidar_info.bbox << std::endl;
79 for (auto & img : imgs) {
80 //std::cout << " \t with bbox: " << imgs[ii].bbox << std::endl;
81 if (lidar_info.intersects(img.bbox)) {
82 volm_img_info info;
83 volm_io_tools::load_naip_img(img_folder(), img.name, lvcs, info, true);
84 intersection_imgs.push_back(info);
85 std::cout << "intersects: " << img.name << std::endl;
86 }
87 }
88 std::cout << " !!!!!! lidar intersects: " << intersection_imgs.size() << " imgs!\n";
89
90 std::string name = vul_file::strip_extension(test()) + "_dict.vsl";
91
92 sdet_texture_classifier_params dummy;
93 sdet_texture_classifier tc(dummy);
94 tc.load_dictionary(name); // loads the params as well
95 unsigned ntextons = tc.get_number_of_textons();
96 std::cout << " testing using the dictionary: " << name << " with the number of textons: " << ntextons << "\n categories:\n";
97
98 std::vector<std::string> cats = tc.get_training_categories();
99 std::map<std::string, vil_rgb<vxl_byte> > cat_color_map;
100 for (const auto & cat : cats) {
101 unsigned char id = volm_label_table::get_id_closest_name(cat);
102 cat_color_map[cat] = volm_label_table::land_id[id].color_;
103 std::cout << "\t\t" << cat << " closest land type: " << volm_label_table::land_id[id].name_ << " color: " << volm_label_table::land_id[id].color_ << '\n';
104 std::cout.flush();
105 }
106
107 // process lidar img
108 sdet_texture_classifier lidar_c((sdet_texture_classifier_params)tc);
109 if (!lidar_c.compute_filter_bank_float_img(filter_folder(), lidar_file(), lidar_max())) {
110 std::cerr << " problems in computing filters of: " << lidar_file() << std::endl;
111 return 0;
112 }
113
114 // compute filters of intersection imgs
115 std::vector<sdet_texture_classifier*> intersection_imgs_naip_c;
116 for (auto & intersection_img : intersection_imgs) {
117 // process naip img
118 sdet_texture_classifier* naip_c = new sdet_texture_classifier((sdet_texture_classifier_params)tc);
119 if (!naip_c->compute_filter_bank_color_img(filter_folder(), intersection_img.img_name)) {
120 std::cerr << " problems in computing filters of: " << intersection_img.img_name << std::endl;
121 return 0;
122 }
123 intersection_imgs_naip_c.push_back(naip_c);
124 }
125
126 vil_image_view<float> lidar_orig = vil_load(lidar_file().c_str());
127
128 vil_image_view<float> out(lidar_info.ni, lidar_info.nj);
129 vil_image_view<vil_rgb<vxl_byte> > out_rgb(out.ni(), out.nj());
130 out.fill(0); out_rgb.fill(vil_rgb<vxl_byte>(0,0,0));
131
132 int bb = block();
133 unsigned larger_neigh = bb+neighborhood_size();
134 double lon, lat, u, v;
135
136 std::cout << "lidar ni: " << lidar_info.ni << " nj: " << lidar_info.nj << "\n";
137
138 // find the local coords of each pixel to retrieve corresponding img pixels
139 for (int i = bb+1; i < (int)lidar_info.ni-bb; i++) {
140 for (int j = bb+1; j < (int)lidar_info.nj-bb; j++) {
141
142 if (lidar_orig(i,j) <= 0)
143 continue;
144
145 lidar_info.cam->img_to_global(i, j, lon, lat);
146
147 // generate a block of pixels around lidar
148 std::vector<std::pair<int, int> > l_pixels;
149 get_block(i, j, bb, l_pixels);
150 std::vector<vnl_vector<double> > data;
151 lidar_c.compute_data(l_pixels, data);
152
153 for (unsigned kk = 0; kk < intersection_imgs.size(); kk++) {
154 int img_ni = intersection_imgs[kk].ni;
155 int img_nj = intersection_imgs[kk].nj;
156
157 //lvcs->global_to_local(-lon, lat, 0, vpgl_lvcs::wgs84, lx, ly, lz);
158 //intersection_imgs[kk].cam->project(lx, ly, 0.0, u, v);
159 intersection_imgs[kk].cam->global_to_img(-lon, lat, 0, u, v);
160 int uu = (int)std::floor(u + 0.5);
161 int vv = (int)std::floor(v + 0.5);
162 if (uu >= 0 && vv >= 0 && uu < img_ni && vv < img_nj) { // there is a correspondence
163 //std::cout << i << ",j "; std::cout.flush();
164 // now around img
165 std::vector<std::pair<int, int> > img_pixels;
166 get_block(uu, vv, larger_neigh, img_pixels); // enlarge the neighborhood to account for misregistration
167
168 // collect data from both images
169 std::vector<vnl_vector<double> > data2;
170 intersection_imgs_naip_c[kk]->compute_data(img_pixels, data2);
171
172 // create texton histogram and classify, use the same weight for all samples put into the histogram
173 float weight = 1.0f/ (data.size() + data2.size());
174 std::vector<float> hist(ntextons, 0.0f);
175 tc.update_hist(data, weight, hist);
176 tc.update_hist(data2, weight, hist);
177 std::pair<std::string, float> hc = tc.highest_prob_class(hist);
178 out_rgb(i,j) = cat_color_map[hc.first];
179 out(i,j) = hc.second;
180 break;
181 }
182 }
183
184 }
185 //std::cout << i << " ";
186 }
187 std::cout << '\n';
188
189 std::string out_name = out_folder() + lidar_info.name + "_outmap";
190 vil_save(out, (out_name + ".tif").c_str());
191 vil_save(out_rgb, (out_name + ".png").c_str());
192
193 std::cout << "total time: " << t.all()/1000 << " seconds = " << t.all()/(1000*60) << " mins.\n";
194 return volm_io::SUCCESS;
195 }
196