1 #include <iostream>
2 #include <iomanip>
3 #include <sstream>
4 #include "volm_tile.h"
5 #include "volm_io.h"
6 //:
7 // \file
8 #ifdef _MSC_VER
9 # include "vcl_msvc_warnings.h"
10 #endif
11 #include <cassert>
12
13 #include <brip/brip_vil_float_ops.h>
14 #include <bkml/bkml_write.h>
15 #include "vpgl/vpgl_utm.h"
16 #include "vul/vul_file.h"
17
18 //: specify lat lon in positive numbers but specify hemisphere ('N' or 'S') and direction ('W' or 'E')
volm_tile(float lat,float lon,char hemisphere,char direction,float scale_i,float scale_j,unsigned ni,unsigned nj)19 volm_tile::volm_tile(float lat, float lon, char hemisphere, char direction, float scale_i, float scale_j, unsigned ni, unsigned nj) :
20 lat_(lat), lon_(lon), hemisphere_(hemisphere), direction_(direction), scale_i_(scale_i), scale_j_(scale_j), ni_(ni), nj_(nj)
21 {
22 vnl_matrix<double> trans_matrix(4,4,0.0);
23 //divide by ni-1 to account for 1 pixel overlap with the next tile
24 if (direction_ == 'E') {
25 trans_matrix[0][3] = lon_ - 0.5/(ni-1.0);
26 trans_matrix[0][0] = scale_i_/(ni-1.0);
27 }
28 else {
29 trans_matrix[0][3] = lon_ + 0.5/(ni-1.0);
30 trans_matrix[0][0] = -scale_i_/(ni-1.0);
31 }
32 if (hemisphere_ == 'N') {
33 trans_matrix[1][1] = -scale_j_/(nj-1.0);
34 trans_matrix[1][3] = lat_+scale_j_+0.5/(nj-1.0);
35 }
36 else {
37 trans_matrix[1][1] = scale_j_/(nj-1.0);
38 trans_matrix[1][3] = lat_-scale_j_-0.5/(nj-1.0);
39 }
40 // just pass an empty lvcs, this geo cam will only be used to compute image pixels to global coords mappings
41 vpgl_lvcs_sptr lv = new vpgl_lvcs;
42 vpgl_geo_camera cam(trans_matrix, lv); cam.set_scale_format(true);
43 cam_ = cam;
44 }
45
46 //: parse the name string and construct tile
volm_tile(const std::string & file_name,unsigned ni,unsigned nj)47 volm_tile::volm_tile(const std::string& file_name, unsigned ni, unsigned nj) : ni_(ni), nj_(nj)
48 {
49 std::string name = vul_file::strip_directory(file_name);
50 name = name.substr(name.find_first_of('_')+1, name.size());
51
52 std::string n_coords = name.substr(0, name.find_first_of('_'));
53 std::string n_scale = name.substr(name.find_first_of('_')+1, name.find_last_of('_')-name.find_first_of('_')-1);
54
55 // determine the lat, lon, hemisphere (North or South) and direction (East or West)
56 std::size_t n = n_coords.find('N');
57 if (n < n_coords.size()) hemisphere_ = 'N';
58 else hemisphere_ = 'S';
59 n = n_coords.find('E');
60 if (n < n_coords.size()) direction_ = 'E';
61 else direction_ = 'W';
62
63 std::string n_str = n_coords.substr(n_coords.find_first_of(hemisphere_)+1,
64 n_coords.find_first_of(direction_)-n_coords.find_first_of(hemisphere_)-1);
65 std::stringstream str(n_str); str >> lat_;
66
67 n_str = n_coords.substr(n_coords.find_first_of(direction_)+1, n_coords.size());
68 std::stringstream str2(n_str); str2 >> lon_;
69
70 n_str = n_scale.substr(n_scale.find_first_of('S')+1, n_scale.find_first_of('x')-n_scale.find_first_of('S')-1);
71 std::stringstream str3(n_str); str3 >> scale_i_; scale_j_ = scale_i_;
72
73 #if 1
74 if (hemisphere_ == 'N') std::cout << " upper left corner in the image is: " << hemisphere_ << lat_+scale_i_ << direction_ << lon_ << std::endl;
75 else std::cout << " upper left corner in the image is: " << hemisphere_ << lat_-scale_i_ << direction_ << lon_ << std::endl;
76 if (direction_ == 'W') std::cout << " lower right corner in the image is: " << hemisphere_ << lat_ << direction_ << lon_-scale_j_ << std::endl;
77 else std::cout << " lower right corner in the image is: " << hemisphere_ << lat_ << direction_ << lon_+scale_j_ << std::endl;
78 #endif
79 vnl_matrix<double> trans_matrix(4,4,0.0);
80 //divide by ni-1 to account for 1 pixel overlap with the next tile
81 if (direction_ == 'E') {
82 trans_matrix[0][3] = lon_ - 0.5/(ni-1.0);
83 trans_matrix[0][0] = scale_i_/(ni-1.0);
84 }
85 else {
86 trans_matrix[0][3] = lon_ + 0.5/(ni-1.0);
87 trans_matrix[0][0] = -scale_i_/(ni-1.0);
88 }
89 if (hemisphere_ == 'N') {
90 trans_matrix[1][1] = -scale_j_/(nj-1.0);
91 trans_matrix[1][3] = lat_+scale_j_+0.5/(nj-1.0);
92 }
93 else {
94 trans_matrix[1][1] = scale_j_/(nj-1.0);
95 trans_matrix[1][3] = lat_-scale_j_-0.5/(nj-1.0);
96 }
97 vpgl_lvcs_sptr dummy_lvcs = new vpgl_lvcs;
98 cam_ = vpgl_geo_camera(trans_matrix, dummy_lvcs);
99 cam_.set_scale_format(true);
100
101 #if 0
102 #ifdef DEBUG
103 std::cout << "will determine transformation matrix from the file name: " << name << std::endl;
104 #endif
105 std::string n = name.substr(name.find_first_of('N')+1, name.find_first_of('W'));
106 assert(n.size() != 0); // for now there is no support for 'S' and/or 'E'
107 hemisphere_ = 'N';
108 direction_ = 'W';
109
110 std::stringstream str(n); str >> lat_;
111 n = name.substr(name.find_first_of('W')+1, name.find_first_of('_'));
112 std::stringstream str2(n); str2 >> lon_;
113
114 name = name.substr(name.find_first_of('_'), name.size());
115 n = name.substr(name.find_first_of('S')+1, name.find_first_of('x'));
116 std::stringstream str3(n); str3 >> scale_i_;
117 n = name.substr(name.find_first_of('x')+1, name.find_last_of('.'));
118 std::stringstream str4(n); str4 >> scale_j_;
119 #ifdef DEBUG
120 std::cout << " lat: " << lat_ << " lon: " << lon_ << " scale_i:" << scale_i_ << " scale_j: " << scale_j_ << std::endl;
121
122 // determine the upper left corner to use a vpgl_geo_cam, subtract from lat
123 std::cout << "upper left corner in the image is: " << lat_+scale_j_ << " N " << lon_ << " W\n"
124 << "lower right corner in the image is: " << lat_ << " N " << lon_-scale_i_ << " W" << std::endl;
125 #endif
126 vnl_matrix<double> trans_matrix(4,4,0.0);
127 //divide by ni-1 to account for 1 pixel overlap with the next tile
128 trans_matrix[0][3] = lon_;
129 if (direction_ == 'E')
130 trans_matrix[0][0] = scale_i_/(ni-1.0);
131 else
132 trans_matrix[0][0] = -scale_i_/(ni-1.0);
133 if (hemisphere_ == 'N') {
134 trans_matrix[1][1] = -scale_j_/(nj-1.0);
135 trans_matrix[1][3] = lat_+scale_j_+1/3600.0;
136 }
137 else {
138 trans_matrix[1][1] = scale_j_/(nj-1.0);
139 trans_matrix[1][3] = lat_-scale_j_-1/3600.0;
140 }
141 vpgl_lvcs_sptr dummy_lvcs = new vpgl_lvcs;
142 cam_ = vpgl_geo_camera(trans_matrix, dummy_lvcs);
143 cam_.set_scale_format(true);
144 #endif
145 }
146
147
148 //: specify lat lon as regular, e.g. negative lon for 'W'
volm_tile(float lat,float lon,float scale_i,float scale_j,unsigned ni,unsigned nj)149 volm_tile::volm_tile(float lat, float lon, float scale_i, float scale_j, unsigned ni, unsigned nj) :
150 scale_i_(scale_i), scale_j_(scale_j), ni_(ni), nj_(nj)
151 {
152 if (lat < 0) { lat_ = -lat; hemisphere_ = 'S'; } else { lat_ = lat; hemisphere_ = 'N'; }
153 if (lon < 0) { lon_ = -lon; direction_ = 'W'; } else { lon_ = lon; direction_ = 'E'; }
154
155 vnl_matrix<double> trans_matrix(4,4,0.0);
156 //divide by ni-1 to account for 1 pixel overlap with the next tile
157 if (direction_ == 'E') {
158 trans_matrix[0][3] = lon_ - 0.5/(ni-1.0);
159 trans_matrix[0][0] = scale_i_/(ni-1.0);
160 }
161 else {
162 trans_matrix[0][3] = lon_ + 0.5/(ni-1.0);
163 trans_matrix[0][0] = -scale_i_/(ni-1.0);
164 }
165 if (hemisphere_ == 'N') {
166 trans_matrix[1][1] = -scale_j_/(nj-1.0);
167 trans_matrix[1][3] = lat_+scale_j_+0.5/(nj-1.0);
168 }
169 else {
170 trans_matrix[1][1] = scale_j_/(nj-1.0);
171 trans_matrix[1][3] = lat_-scale_j_-0.5/(nj-1.0);
172 }
173 // just pass an empty lvcs, this geo cam will only be used to compute image pixels to global coords mappings
174 vpgl_lvcs_sptr lv = new vpgl_lvcs;
175 vpgl_geo_camera cam(trans_matrix, lv); cam.set_scale_format(true);
176 cam_ = cam;
177 }
178
generate_p1_tiles()179 std::vector<volm_tile> volm_tile::generate_p1_tiles()
180 {
181 std::vector<volm_tile> p1_tiles;
182 unsigned ni = 3601;
183 unsigned nj = 3601;
184 float scale_i = 1;
185 float scale_j = 1;
186 p1_tiles.emplace_back(37, 118, 'N', 'W', scale_i, scale_j, ni, nj);
187 p1_tiles.emplace_back(37, 119, 'N', 'W', scale_i, scale_j, ni, nj);
188 p1_tiles.emplace_back(38, 118, 'N', 'W', scale_i, scale_j, ni, nj);
189 p1_tiles.emplace_back(38, 119, 'N', 'W', scale_i, scale_j, ni, nj);
190 p1_tiles.emplace_back(30, 82, 'N', 'W', scale_i, scale_j, ni, nj);
191 p1_tiles.emplace_back(31, 81, 'N', 'W', scale_i, scale_j, ni, nj);
192 p1_tiles.emplace_back(31, 82, 'N', 'W', scale_i, scale_j, ni, nj);
193 p1_tiles.emplace_back(32, 80, 'N', 'W', scale_i, scale_j, ni, nj);
194 p1_tiles.emplace_back(32, 81, 'N', 'W', scale_i, scale_j, ni, nj);
195 p1_tiles.emplace_back(33, 78, 'N', 'W', scale_i, scale_j, ni, nj);
196 p1_tiles.emplace_back(33, 79, 'N', 'W', scale_i, scale_j, ni, nj);
197 p1_tiles.emplace_back(33, 80, 'N', 'W', scale_i, scale_j, ni, nj);
198 p1_tiles.emplace_back(34, 77, 'N', 'W', scale_i, scale_j, ni, nj);
199 p1_tiles.emplace_back(34, 78, 'N', 'W', scale_i, scale_j, ni, nj);
200 p1_tiles.emplace_back(34, 79, 'N', 'W', scale_i, scale_j, ni, nj);
201 p1_tiles.emplace_back(35, 76, 'N', 'W', scale_i, scale_j, ni, nj);
202 p1_tiles.emplace_back(35, 77, 'N', 'W', scale_i, scale_j, ni, nj);
203 p1_tiles.emplace_back(36, 76, 'N', 'W', scale_i, scale_j, ni, nj);
204 return p1_tiles;
205 }
206
generate_p1_wr1_tiles()207 std::vector<volm_tile> volm_tile::generate_p1_wr1_tiles()
208 {
209 std::vector<volm_tile> p1_tiles;
210 unsigned ni = 3601;
211 unsigned nj = 3601;
212 float scale_i = 1;
213 float scale_j = 1;
214 p1_tiles.emplace_back(37, 118, 'N', 'W', scale_i, scale_j, ni, nj);
215 p1_tiles.emplace_back(37, 119, 'N', 'W', scale_i, scale_j, ni, nj);
216 p1_tiles.emplace_back(38, 118, 'N', 'W', scale_i, scale_j, ni, nj);
217 p1_tiles.emplace_back(38, 119, 'N', 'W', scale_i, scale_j, ni, nj);
218 return p1_tiles;
219 }
220
generate_p1_wr2_tiles()221 std::vector<volm_tile> volm_tile::generate_p1_wr2_tiles()
222 {
223 std::vector<volm_tile> p1_tiles;
224 unsigned ni = 3601;
225 unsigned nj = 3601;
226 float scale_i = 1;
227 float scale_j = 1;
228 p1_tiles.emplace_back(30, 82, 'N', 'W', scale_i, scale_j, ni, nj);
229 p1_tiles.emplace_back(31, 81, 'N', 'W', scale_i, scale_j, ni, nj);
230 p1_tiles.emplace_back(31, 82, 'N', 'W', scale_i, scale_j, ni, nj);
231 p1_tiles.emplace_back(32, 80, 'N', 'W', scale_i, scale_j, ni, nj);
232 p1_tiles.emplace_back(32, 81, 'N', 'W', scale_i, scale_j, ni, nj);
233 p1_tiles.emplace_back(33, 78, 'N', 'W', scale_i, scale_j, ni, nj);
234 p1_tiles.emplace_back(33, 79, 'N', 'W', scale_i, scale_j, ni, nj);
235 p1_tiles.emplace_back(33, 80, 'N', 'W', scale_i, scale_j, ni, nj);
236 p1_tiles.emplace_back(34, 77, 'N', 'W', scale_i, scale_j, ni, nj);
237 p1_tiles.emplace_back(34, 78, 'N', 'W', scale_i, scale_j, ni, nj);
238 p1_tiles.emplace_back(34, 79, 'N', 'W', scale_i, scale_j, ni, nj);
239 p1_tiles.emplace_back(35, 76, 'N', 'W', scale_i, scale_j, ni, nj);
240 p1_tiles.emplace_back(35, 77, 'N', 'W', scale_i, scale_j, ni, nj);
241 p1_tiles.emplace_back(36, 76, 'N', 'W', scale_i, scale_j, ni, nj);
242 return p1_tiles;
243 }
244
generate_p1b_wr1_tiles()245 std::vector<volm_tile> volm_tile::generate_p1b_wr1_tiles()
246 {
247 std::vector<volm_tile> p1b_tiles;
248 unsigned ni = 3601;
249 unsigned nj = 3601;
250 float scale_i = 1;
251 float scale_j = 1;
252 //p1b_tiles.push_back(volm_tile(33, 71, 'S', 'W', scale_i, scale_j, ni, nj));
253 //p1b_tiles.push_back(volm_tile(33, 72, 'S', 'W', scale_i, scale_j, ni, nj));
254 p1b_tiles.emplace_back(34, 71, 'S', 'W', scale_i, scale_j, ni, nj);
255 p1b_tiles.emplace_back(34, 72, 'S', 'W', scale_i, scale_j, ni, nj);
256 p1b_tiles.emplace_back(35, 71, 'S', 'W', scale_i, scale_j, ni, nj);
257 p1b_tiles.emplace_back(35, 72, 'S', 'W', scale_i, scale_j, ni, nj);
258 return p1b_tiles;
259 }
260
generate_p1b_wr2_tiles()261 std::vector<volm_tile> volm_tile::generate_p1b_wr2_tiles()
262 {
263 std::vector<volm_tile> p1b_tiles;
264 unsigned ni = 3601;
265 unsigned nj = 3601;
266 float scale_i = 1;
267 float scale_j = 1;
268 p1b_tiles.emplace_back(12, 76, 'N', 'E', scale_i, scale_j, ni, nj);
269 p1b_tiles.emplace_back(12, 77, 'N', 'E', scale_i, scale_j, ni, nj);
270 p1b_tiles.emplace_back(13, 76, 'N', 'E', scale_i, scale_j, ni, nj);
271 p1b_tiles.emplace_back(13, 77, 'N', 'E', scale_i, scale_j, ni, nj);
272 return p1b_tiles;
273 }
274
generate_p1b_wr3_tiles()275 std::vector<volm_tile> volm_tile::generate_p1b_wr3_tiles()
276 {
277 std::vector<volm_tile> p1b_tiles;
278 unsigned ni = 3601;
279 unsigned nj = 3601;
280 float scale_i = 1;
281 float scale_j = 1;
282 p1b_tiles.emplace_back(30, 35, 'N', 'E', scale_i, scale_j, ni, nj);
283 p1b_tiles.emplace_back(31, 35, 'N', 'E', scale_i, scale_j, ni, nj);
284 p1b_tiles.emplace_back(31, 36, 'N', 'E', scale_i, scale_j, ni, nj);
285 p1b_tiles.emplace_back(32, 35, 'N', 'E', scale_i, scale_j, ni, nj);
286 p1b_tiles.emplace_back(32, 36, 'N', 'E', scale_i, scale_j, ni, nj);
287 return p1b_tiles;
288 }
289
generate_p1b_wr4_tiles()290 std::vector<volm_tile> volm_tile::generate_p1b_wr4_tiles()
291 {
292 std::vector<volm_tile> p1b_tiles;
293 unsigned ni = 3601;
294 unsigned nj = 3601;
295 float scale_i = 1;
296 float scale_j = 1;
297 p1b_tiles.emplace_back(5, 124, 'N', 'E', scale_i, scale_j, ni, nj);
298 p1b_tiles.emplace_back(5, 125, 'N', 'E', scale_i, scale_j, ni, nj);
299 p1b_tiles.emplace_back(6, 124, 'N', 'E', scale_i, scale_j, ni, nj);
300 p1b_tiles.emplace_back(6, 125, 'N', 'E', scale_i, scale_j, ni, nj);
301 p1b_tiles.emplace_back(6, 126, 'N', 'E', scale_i, scale_j, ni, nj);
302 p1b_tiles.emplace_back(7, 125, 'N', 'E', scale_i, scale_j, ni, nj);
303 p1b_tiles.emplace_back(7, 126, 'N', 'E', scale_i, scale_j, ni, nj);
304 return p1b_tiles;
305 }
306
generate_p1b_wr5_tiles()307 std::vector<volm_tile> volm_tile::generate_p1b_wr5_tiles()
308 {
309 std::vector<volm_tile> p1b_tiles;
310 unsigned ni = 3601;
311 unsigned nj = 3601;
312 float scale_i = 1;
313 float scale_j = 1;
314 p1b_tiles.emplace_back(22, 120, 'N', 'E', scale_i, scale_j, ni, nj);
315 p1b_tiles.emplace_back(23, 120, 'N', 'E', scale_i, scale_j, ni, nj);
316 p1b_tiles.emplace_back(24, 120, 'N', 'E', scale_i, scale_j, ni, nj);
317 p1b_tiles.emplace_back(24, 121, 'N', 'E', scale_i, scale_j, ni, nj);
318 p1b_tiles.emplace_back(25, 121, 'N', 'E', scale_i, scale_j, ni, nj);
319 return p1b_tiles;
320 }
321
generate_p1b_wr_tiles(int world_id)322 std::vector<volm_tile> volm_tile::generate_p1b_wr_tiles(int world_id)
323 {
324 std::vector<volm_tile> tiles;
325 switch (world_id) {
326 case 1: { tiles = volm_tile::generate_p1b_wr1_tiles(); break; }
327 case 2: { tiles = volm_tile::generate_p1b_wr2_tiles(); break; }
328 case 3: { tiles = volm_tile::generate_p1b_wr3_tiles(); break; }
329 case 4: { tiles = volm_tile::generate_p1b_wr4_tiles(); break; }
330 case 5: { tiles = volm_tile::generate_p1b_wr5_tiles(); break; }
331 default: {
332 std::cerr << "Unknown world id: " << world_id << std::endl;
333 } }
334 return tiles;
335 }
336
generate_p2a_wr8_tiles()337 std::vector<volm_tile> volm_tile::generate_p2a_wr8_tiles()
338 {
339 std::vector<volm_tile> tiles;
340 unsigned ni = 3601, nj = 3601;
341 float scale_i = 1.0f, scale_j = 1.0f;
342 tiles.emplace_back(33, 145, 'S', 'E', scale_i, scale_j, ni, nj);
343 tiles.emplace_back(33, 146, 'S', 'E', scale_i, scale_j, ni, nj);
344 tiles.emplace_back(33, 147, 'S', 'E', scale_i, scale_j, ni, nj);
345 tiles.emplace_back(33, 148, 'S', 'E', scale_i, scale_j, ni, nj);
346 tiles.emplace_back(33, 149, 'S', 'E', scale_i, scale_j, ni, nj);
347 tiles.emplace_back(33, 150, 'S', 'E', scale_i, scale_j, ni, nj);
348 tiles.emplace_back(33, 151, 'S', 'E', scale_i, scale_j, ni, nj);
349 tiles.emplace_back(33, 152, 'S', 'E', scale_i, scale_j, ni, nj);
350 tiles.emplace_back(34, 145, 'S', 'E', scale_i, scale_j, ni, nj);
351 tiles.emplace_back(34, 146, 'S', 'E', scale_i, scale_j, ni, nj);
352 tiles.emplace_back(34, 147, 'S', 'E', scale_i, scale_j, ni, nj);
353 tiles.emplace_back(34, 148, 'S', 'E', scale_i, scale_j, ni, nj);
354 tiles.emplace_back(34, 149, 'S', 'E', scale_i, scale_j, ni, nj);
355 tiles.emplace_back(34, 150, 'S', 'E', scale_i, scale_j, ni, nj);
356 tiles.emplace_back(34, 151, 'S', 'E', scale_i, scale_j, ni, nj);
357 tiles.emplace_back(35, 145, 'S', 'E', scale_i, scale_j, ni, nj);
358 tiles.emplace_back(35, 146, 'S', 'E', scale_i, scale_j, ni, nj);
359 tiles.emplace_back(35, 147, 'S', 'E', scale_i, scale_j, ni, nj);
360 tiles.emplace_back(35, 148, 'S', 'E', scale_i, scale_j, ni, nj);
361 tiles.emplace_back(35, 149, 'S', 'E', scale_i, scale_j, ni, nj);
362 tiles.emplace_back(35, 150, 'S', 'E', scale_i, scale_j, ni, nj);
363 tiles.emplace_back(35, 151, 'S', 'E', scale_i, scale_j, ni, nj);
364 tiles.emplace_back(36, 145, 'S', 'E', scale_i, scale_j, ni, nj);
365 tiles.emplace_back(36, 146, 'S', 'E', scale_i, scale_j, ni, nj);
366 tiles.emplace_back(36, 147, 'S', 'E', scale_i, scale_j, ni, nj);
367 tiles.emplace_back(36, 148, 'S', 'E', scale_i, scale_j, ni, nj);
368 tiles.emplace_back(36, 149, 'S', 'E', scale_i, scale_j, ni, nj);
369 tiles.emplace_back(36, 150, 'S', 'E', scale_i, scale_j, ni, nj);
370 tiles.emplace_back(37, 145, 'S', 'E', scale_i, scale_j, ni, nj);
371 tiles.emplace_back(37, 146, 'S', 'E', scale_i, scale_j, ni, nj);
372 tiles.emplace_back(37, 147, 'S', 'E', scale_i, scale_j, ni, nj);
373 tiles.emplace_back(37, 148, 'S', 'E', scale_i, scale_j, ni, nj);
374 tiles.emplace_back(37, 149, 'S', 'E', scale_i, scale_j, ni, nj);
375 tiles.emplace_back(38, 145, 'S', 'E', scale_i, scale_j, ni, nj);
376 tiles.emplace_back(38, 146, 'S', 'E', scale_i, scale_j, ni, nj);
377 tiles.emplace_back(38, 147, 'S', 'E', scale_i, scale_j, ni, nj);
378 tiles.emplace_back(38, 148, 'S', 'E', scale_i, scale_j, ni, nj);
379 tiles.emplace_back(38, 149, 'S', 'E', scale_i, scale_j, ni, nj);
380 return tiles;
381 }
382
generate_p2a_wr9_tiles()383 std::vector<volm_tile> volm_tile::generate_p2a_wr9_tiles()
384 {
385 std::vector<volm_tile> tiles;
386 unsigned ni = 3601, nj = 3601;
387 float scale_i = 1.0f, scale_j = 1.0f;
388 tiles.emplace_back(35, 147, 'S', 'E', scale_i, scale_j, ni, nj);
389 tiles.emplace_back(35, 148, 'S', 'E', scale_i, scale_j, ni, nj);
390 tiles.emplace_back(35, 149, 'S', 'E', scale_i, scale_j, ni, nj);
391 tiles.emplace_back(35, 150, 'S', 'E', scale_i, scale_j, ni, nj);
392 tiles.emplace_back(34, 147, 'S', 'E', scale_i, scale_j, ni, nj);
393 tiles.emplace_back(34, 148, 'S', 'E', scale_i, scale_j, ni, nj);
394 tiles.emplace_back(34, 149, 'S', 'E', scale_i, scale_j, ni, nj);
395 tiles.emplace_back(33, 147, 'S', 'E', scale_i, scale_j, ni, nj);
396 tiles.emplace_back(33, 148, 'S', 'E', scale_i, scale_j, ni, nj);
397 tiles.emplace_back(33, 149, 'S', 'E', scale_i, scale_j, ni, nj);
398 return tiles;
399 }
400
generate_p2a_wr10_tiles()401 std::vector<volm_tile> volm_tile::generate_p2a_wr10_tiles()
402 {
403 std::vector<volm_tile> tiles;
404 unsigned ni = 3601, nj = 3601;
405 float scale_i = 1.0f, scale_j = 1.0f;
406 tiles.emplace_back(35, 149, 'S', 'E', scale_i, scale_j, ni, nj);
407 tiles.emplace_back(35, 150, 'S', 'E', scale_i, scale_j, ni, nj);
408 tiles.emplace_back(35, 151, 'S', 'E', scale_i, scale_j, ni, nj);
409 tiles.emplace_back(34, 149, 'S', 'E', scale_i, scale_j, ni, nj);
410 tiles.emplace_back(34, 150, 'S', 'E', scale_i, scale_j, ni, nj);
411 tiles.emplace_back(34, 151, 'S', 'E', scale_i, scale_j, ni, nj);
412 tiles.emplace_back(33, 149, 'S', 'E', scale_i, scale_j, ni, nj);
413 tiles.emplace_back(33, 150, 'S', 'E', scale_i, scale_j, ni, nj);
414 tiles.emplace_back(33, 151, 'S', 'E', scale_i, scale_j, ni, nj);
415 tiles.emplace_back(33, 152, 'S', 'E', scale_i, scale_j, ni, nj);
416 return tiles;
417 }
418
generate_p2a_wr11_tiles()419 std::vector<volm_tile> volm_tile::generate_p2a_wr11_tiles()
420 {
421 std::vector<volm_tile> tiles;
422 unsigned ni = 3601, nj = 3601;
423 float scale_i = 1.0f, scale_j = 1.0f;
424 tiles.emplace_back(35, 145, 'S', 'E', scale_i, scale_j, ni, nj);
425 tiles.emplace_back(35, 146, 'S', 'E', scale_i, scale_j, ni, nj);
426 tiles.emplace_back(35, 147, 'S', 'E', scale_i, scale_j, ni, nj);
427 tiles.emplace_back(34, 145, 'S', 'E', scale_i, scale_j, ni, nj);
428 tiles.emplace_back(34, 146, 'S', 'E', scale_i, scale_j, ni, nj);
429 tiles.emplace_back(34, 147, 'S', 'E', scale_i, scale_j, ni, nj);
430 tiles.emplace_back(33, 145, 'S', 'E', scale_i, scale_j, ni, nj);
431 tiles.emplace_back(33, 146, 'S', 'E', scale_i, scale_j, ni, nj);
432 tiles.emplace_back(33, 147, 'S', 'E', scale_i, scale_j, ni, nj);
433 return tiles;
434 }
435
generate_p2a_wr12_tiles()436 std::vector<volm_tile> volm_tile::generate_p2a_wr12_tiles()
437 {
438 std::vector<volm_tile> tiles;
439 unsigned ni = 3601, nj = 3601;
440 float scale_i = 1.0f, scale_j = 1.0f;
441 tiles.emplace_back(38, 148, 'S', 'E', scale_i, scale_j, ni, nj);
442 tiles.emplace_back(38, 149, 'S', 'E', scale_i, scale_j, ni, nj);
443 tiles.emplace_back(37, 147, 'S', 'E', scale_i, scale_j, ni, nj);
444 tiles.emplace_back(37, 148, 'S', 'E', scale_i, scale_j, ni, nj);
445 tiles.emplace_back(37, 149, 'S', 'E', scale_i, scale_j, ni, nj);
446 tiles.emplace_back(36, 147, 'S', 'E', scale_i, scale_j, ni, nj);
447 tiles.emplace_back(36, 148, 'S', 'E', scale_i, scale_j, ni, nj);
448 tiles.emplace_back(36, 149, 'S', 'E', scale_i, scale_j, ni, nj);
449 tiles.emplace_back(36, 150, 'S', 'E', scale_i, scale_j, ni, nj);
450 tiles.emplace_back(35, 147, 'S', 'E', scale_i, scale_j, ni, nj);
451 tiles.emplace_back(35, 148, 'S', 'E', scale_i, scale_j, ni, nj);
452 tiles.emplace_back(35, 149, 'S', 'E', scale_i, scale_j, ni, nj);
453 tiles.emplace_back(35, 150, 'S', 'E', scale_i, scale_j, ni, nj);
454 return tiles;
455 }
456
generate_p2a_wr13_tiles()457 std::vector<volm_tile> volm_tile::generate_p2a_wr13_tiles()
458 {
459 std::vector<volm_tile> tiles;
460 unsigned ni = 3601, nj = 3601;
461 float scale_i = 1.0f, scale_j = 1.0f;
462 tiles.emplace_back(38, 145, 'S', 'E', scale_i, scale_j, ni, nj);
463 tiles.emplace_back(38, 146, 'S', 'E', scale_i, scale_j, ni, nj);
464 tiles.emplace_back(37, 145, 'S', 'E', scale_i, scale_j, ni, nj);
465 tiles.emplace_back(37, 146, 'S', 'E', scale_i, scale_j, ni, nj);
466 tiles.emplace_back(37, 147, 'S', 'E', scale_i, scale_j, ni, nj);
467 tiles.emplace_back(36, 145, 'S', 'E', scale_i, scale_j, ni, nj);
468 tiles.emplace_back(36, 146, 'S', 'E', scale_i, scale_j, ni, nj);
469 tiles.emplace_back(36, 147, 'S', 'E', scale_i, scale_j, ni, nj);
470 tiles.emplace_back(35, 145, 'S', 'E', scale_i, scale_j, ni, nj);
471 tiles.emplace_back(35, 146, 'S', 'E', scale_i, scale_j, ni, nj);
472 tiles.emplace_back(35, 147, 'S', 'E', scale_i, scale_j, ni, nj);
473 return tiles;
474 }
475
generate_p2a_wr_tiles(int world_id)476 std::vector<volm_tile> volm_tile::generate_p2a_wr_tiles(int world_id)
477 {
478 std::vector<volm_tile> tiles;
479 tiles.clear();
480 switch (world_id) {
481 case 8: { tiles = volm_tile::generate_p2a_wr8_tiles(); break; }
482 case 9: { tiles = volm_tile::generate_p2a_wr9_tiles(); break; }
483 case 10: { tiles = volm_tile::generate_p2a_wr10_tiles(); break; }
484 case 11: { tiles = volm_tile::generate_p2a_wr11_tiles(); break; }
485 case 12: { tiles = volm_tile::generate_p2a_wr12_tiles(); break; }
486 case 13: { tiles = volm_tile::generate_p2a_wr13_tiles(); break; }
487 default: { std::cerr << "in volm_tile::generate_p2a_wr_tiles: unknown world id: " << world_id << std::endl; }
488 }
489 return tiles;
490 }
491
generate_tiles(unsigned const & world_id,std::vector<volm_tile> & tiles)492 bool volm_tile::generate_tiles(unsigned const& world_id, std::vector<volm_tile>& tiles)
493 {
494 tiles.clear();
495 switch (world_id)
496 {
497 case 1: { tiles = volm_tile::generate_p1b_wr1_tiles(); break; }
498 case 2: { tiles = volm_tile::generate_p1b_wr2_tiles(); break; }
499 case 3: { tiles = volm_tile::generate_p1b_wr3_tiles(); break; }
500 case 4: { tiles = volm_tile::generate_p1b_wr4_tiles(); break; }
501 case 5: { tiles = volm_tile::generate_p1b_wr5_tiles(); break; }
502 case 6: { tiles = volm_tile::generate_p1_wr2_tiles(); break; }
503 case 7: { tiles = volm_tile::generate_p1_wr1_tiles(); break; }
504 case 8: { tiles = volm_tile::generate_p2a_wr8_tiles(); break; }
505 case 9: { tiles = volm_tile::generate_p2a_wr9_tiles(); break; }
506 case 10: { tiles = volm_tile::generate_p2a_wr10_tiles(); break; }
507 case 11: { tiles = volm_tile::generate_p2a_wr11_tiles(); break; }
508 case 12: { tiles = volm_tile::generate_p2a_wr12_tiles(); break; }
509 case 13: { tiles = volm_tile::generate_p2a_wr13_tiles(); break; }
510 default: { return false; }
511 }
512 return true;
513 }
514
get_string()515 std::string volm_tile::get_string()
516 {
517 std::stringstream str;
518 str << hemisphere_ << lat_ << direction_ << std::setfill('0') << std::setw(3) << lon_
519 << "_S" << scale_i_ << 'x' << scale_j_;
520 return str.str();
521 }
522
img_to_global(unsigned i,unsigned j,double & lon,double & lat)523 void volm_tile::img_to_global(unsigned i, unsigned j, double& lon, double& lat)
524 {
525 cam_.img_to_global(i, j, lon, lat);
526 if (direction_ == 'W')
527 lon = -lon;
528 if (hemisphere_ == 'S')
529 lat = -lat;
530 }
531
532 //: do the conversions in reporting lon,lat
lower_left_lon()533 double volm_tile::lower_left_lon()
534 {
535 if (direction_ == 'W')
536 return -lon_;
537 else
538 return lon_;
539 }
540
lower_left_lat()541 double volm_tile::lower_left_lat()
542 {
543 if (hemisphere_ == 'S')
544 return -lat_;
545 else
546 return lat_;
547 }
548
549
global_to_img(double lon,double lat,unsigned & i,unsigned & j)550 bool volm_tile::global_to_img(double lon, double lat, unsigned& i, unsigned& j)
551 {
552 double u,v; double dummy_elev = 0;
553
554 if (direction_ == 'W')
555 lon = -lon;
556 if (hemisphere_ == 'S')
557 lat = -lat;
558 cam_.global_to_img(lon, lat, dummy_elev, u, v);
559 if (u < 0 || v < 0 || u >= this->ni_ || v >= this->nj_)
560 return false;
561 i = (unsigned)std::floor( (int)(u*100+0.5)/100+0.5); // truncation up to 0.01 floating precion
562 j = (unsigned)std::floor( (int)(v*100+0.5)/100+0.5); // this may be ceil cause image direction is in reverse in latitude
563 if (j == this->nj_) j--; // v may be larger than nj+0.5 due to the floating point precision
564 return true;
565 }
566
567 //: calculate width of the tile
calculate_width()568 double volm_tile::calculate_width()
569 {
570 vpgl_utm utm;
571 double x, x2, y, y2; int utm_zone;
572 utm.transform(lat_, lon_, x, y, utm_zone);
573 utm.transform(lat_, lon_+scale_i_, x2, y2, utm_zone);
574 return std::abs(x-x2);
575 }
576
577 //: calculate width of the tile
calculate_height()578 double volm_tile::calculate_height()
579 {
580 vpgl_utm utm;
581 double x, x2, y, y2; int utm_zone;
582 utm.transform(lat_, lon_, x, y, utm_zone);
583 utm.transform(lat_+scale_j_, lon_, x2, y2, utm_zone);
584 return std::abs(y-y2);
585 }
586
get_uncertainty_region(float lambda_i,float lambda_j,float cutoff,vbl_array_2d<bool> & mask,vbl_array_2d<float> & kernel)587 void volm_tile::get_uncertainty_region(float lambda_i, float lambda_j, float cutoff, vbl_array_2d<bool>& mask, vbl_array_2d<float>& kernel)
588 {
589 #if 0
590 //brip_vil_float_ops::extrema_kernel_mask(lambda_i, lambda_j, 0.0, kernel, mask, cutoff);
591 //brip_vil_float_ops::gaussian_kernel_mask(lambda_i, lambda_j, 0.0, kernel, mask, cutoff, true);
592 #else
593 assert(lambda_i == lambda_j);
594 brip_vil_float_ops::gaussian_kernel_square_mask(lambda_i, kernel, mask, cutoff);
595 #endif
596 auto nrows = (unsigned)mask.rows();
597 auto ncols = (unsigned)mask.cols();
598 float kernel_max = kernel[nrows/2][ncols/2];
599 // normalize kernel
600 for (unsigned i = 0; i < ncols; ++i)
601 for (unsigned j = 0; j < nrows; ++j)
602 kernel[j][i] /= kernel_max;
603 }
604
605 // mark the uncertainty region around a given location using a gaussian mask, the center has the designated score, the rest diminishes with respect to a cutoff degree
mark_uncertainty_region(int i,int j,float score,vbl_array_2d<bool> & mask,vbl_array_2d<float> & kernel,vil_image_view<unsigned int> & img)606 void volm_tile::mark_uncertainty_region(int i, int j, float score, vbl_array_2d<bool>& mask, vbl_array_2d<float>& kernel, vil_image_view<unsigned int>& img)
607 {
608 auto nrows = (unsigned)mask.rows();
609 auto ncols = (unsigned)mask.cols();
610
611 int js = (int)std::floor(j - (float)nrows/2.0f + 0.5f);
612 int is = (int)std::floor(i - (float)ncols/2.0f + 0.5f);
613 int je = (int)std::floor(j + (float)nrows/2.0f + 0.5f);
614 int ie = (int)std::floor(i + (float)ncols/2.0f + 0.5f);
615
616 int ni = (int)img.ni();
617 int nj = (int)img.nj();
618 if (score > 0.0f) {
619 for (int ii = is; ii < ie; ++ii)
620 for (int jj = js; jj < je; ++jj) {
621 int mask_i = ii - is;
622 int mask_j = jj - js;
623 if (mask[mask_j][mask_i] && ii >= 0 && jj >= 0 && ii < ni && jj < nj) {
624 float val = score*kernel[mask_j][mask_i];
625 unsigned int pix_val = (unsigned int)(val*volm_io::SCALE_VALUE) + 1; // scale it
626 if (img(ii,jj) > 0)
627 img(ii,jj) = (img(ii,jj)+pix_val)/2; // overwrites whatever values was in the image
628 else
629 img(ii,jj) = pix_val;
630 }
631 }
632 }
633 }
634
635 // just overwrite whatever was in the image, cause it assumes that current value of the pixel is compared to this score before this function is called
mark_uncertainty_region(int i,int j,float score,vbl_array_2d<bool> & mask,vbl_array_2d<float> & kernel,vil_image_view<vxl_byte> & img)636 void volm_tile::mark_uncertainty_region(int i, int j, float score, vbl_array_2d<bool>& mask, vbl_array_2d<float>& kernel, vil_image_view<vxl_byte>& img)
637 {
638 auto nrows = (unsigned)mask.rows();
639 auto ncols = (unsigned)mask.cols();
640
641 int js = (int)std::floor(j - (float)nrows/2.0f + 0.5f);
642 int is = (int)std::floor(i - (float)ncols/2.0f + 0.5f);
643 int je = (int)std::floor(j + (float)nrows/2.0f + 0.5f);
644 int ie = (int)std::floor(i + (float)ncols/2.0f + 0.5f);
645
646 int ni = (int)img.ni();
647 int nj = (int)img.nj();
648 if (score > 0.0f) {
649
650 for (int ii = is; ii < ie; ++ii)
651 for (int jj = js; jj < je; ++jj) {
652 int mask_i = ii - is;
653 int mask_j = jj - js;
654 if (mask[mask_j][mask_i] && ii >= 0 && jj >= 0 && ii < ni && jj < nj) {
655 float val = score*kernel[mask_j][mask_i];
656 //unsigned int pix_val = (unsigned int)(val*volm_io::SCALE_VALUE) + 1; // scale it
657 unsigned char pix_val = (unsigned char)(val*volm_io::SCALE_VALUE) + 1;
658 #if 0
659 if (pix_val < volm_io::UNKNOWN) {
660 pix_val = volm_io::STRONG_NEGATIVE;
661 }
662 else if (pix_val == volm_io::UNKNOWN) {
663 pix_val = volm_io::UNKNOWN;
664 }
665 else {
666 pix_val = volm_io::STRONG_POSITIVE;
667 }
668 if (img(ii,jj) > 0)
669 img(ii,jj) = (img(ii,jj)+pix_val)/2;
670 else
671 img(ii,jj) = pix_val; // overwrites whatever value was in the image
672 #endif // 0
673 img(ii,jj) = pix_val;
674 }
675 }
676 }
677 }
678
679 // create a kml file of the tile as a box and with circular marks throughout at every n pixels in each direction
write_kml(const std::string & name,int n)680 void volm_tile::write_kml(const std::string& name, int n)
681 {
682 std::ofstream ofs(name.c_str());
683 bkml_write::open_document(ofs);
684
685 double lon, lat;
686 this->img_to_global(0,0,lon,lat);
687 vnl_double_2 ul; ul[0] = lat; ul[1] = lon;
688 this->img_to_global(0, this->nj_, lon, lat);
689 vnl_double_2 ll; ll[0] = lat; ll[1] = lon;
690 this->img_to_global(this->ni_, this->nj_, lon, lat);
691 vnl_double_2 lr; lr[0] = lat; lr[1] = lon;
692 this->img_to_global(this->ni_, 0, lon, lat);
693 vnl_double_2 ur; ur[0] = lat; ur[1] = lon;
694 bkml_write::write_box(ofs, this->get_string(), this->get_string(), ul, ur, ll, lr);
695 if (n != 0) {
696 double arcsecond = (n/2.0) * (1.0/3600.0);
697 for (unsigned i = 0; i < this->ni_; i+=n)
698 for (unsigned j = 0; j < this->nj_; j+=n) {
699 this->img_to_global(i,j,lon,lat);
700 ul[0] = lat; ul[1] = lon;
701 ll[0] = lat-arcsecond; ll[1] = lon;
702 lr[0] = lat-arcsecond; lr[1] = lon+arcsecond;
703 ur[0] = lat; ur[1] = lon+arcsecond;
704 bkml_write::write_box(ofs, this->get_string(), this->get_string(), ul, ur, ll, lr);
705 }
706 }
707 bkml_write::close_document(ofs);
708 ofs.close();
709 }
710
711 #if 0 //TODO
712 volm_tile volm_tile::parse_string(std::string& filename)
713 {
714 std::string name = vul_file::strip_directory(filename);
715 name = name.substr(name.find_first_of('_')+1, name.size());
716 if (name.find_first_of('N') != name.end())
717 std::string n = name.substr(name.find_first_of('N')+1, name.find_first_of('W'));
718 float lon, lat, scale;
719 std::stringstream str(n); str >> lat;
720 n = name.substr(name.find_first_of('W')+1, name.find_first_of('_'));
721 std::stringstream str2(n); str2 >> lon;
722 n = name.substr(name.find_first_of('x')+1, name.find_last_of('.'));
723 std::stringstream str3(n); str3 >> scale;
724 std::cout << " lat: " << lat << " lon: " << lon << " scale: " << scale << std::endl;
725 }
726 #endif
727
728
729 //: Binary save self to stream.
b_write(vsl_b_ostream & os) const730 void volm_tile::b_write(vsl_b_ostream &os) const
731 {
732 vsl_b_write(os, version());
733 vsl_b_write(os, lat_);
734 vsl_b_write(os, lon_);
735 vsl_b_write(os, hemisphere_);
736 vsl_b_write(os, direction_);
737 vsl_b_write(os, scale_i_);
738 vsl_b_write(os, scale_j_);
739 vsl_b_write(os, ni_);
740 vsl_b_write(os, nj_);
741 cam_.b_write(os);
742 }
743
744 //: Binary load self from stream.
b_read(vsl_b_istream & is)745 void volm_tile::b_read(vsl_b_istream &is)
746 {
747 if (!is) return;
748 short ver;
749 vsl_b_read(is, ver);
750 switch (ver)
751 {
752 case 1:
753 vsl_b_read(is, lat_);
754 vsl_b_read(is, lon_);
755 vsl_b_read(is, hemisphere_);
756 vsl_b_read(is, direction_);
757 vsl_b_read(is, scale_i_);
758 vsl_b_read(is, scale_j_);
759 vsl_b_read(is, ni_);
760 vsl_b_read(is, nj_);
761 cam_.b_read(is);
762 break;
763
764 default:
765 std::cerr << "I/O ERROR: vpgl_geo_camera::b_read(vsl_b_istream&)\n"
766 << " Unknown version number "<< ver << '\n';
767 is.is().clear(std::ios::badbit); // Set an unrecoverable IO error on stream
768 return;
769 }
770 }
771
bbox()772 vgl_box_2d<float> volm_tile::bbox()
773 {
774 vgl_point_2d<float> p1((float)this->lower_left_lon(), (float)this->lower_left_lat());
775 vgl_point_2d<float> p2(p1.x()+this->scale_i_, p1.y()+this->scale_j_);
776 vgl_box_2d<float> box(p1, p2);
777 return box;
778 }
779
bbox_double()780 vgl_box_2d<double> volm_tile::bbox_double()
781 {
782 vgl_point_2d<double> p1(this->lower_left_lon(), this->lower_left_lat());
783 vgl_point_2d<double> p2(p1.x()+this->scale_i_, p1.y()+this->scale_j_);
784 vgl_box_2d<double> box(p1, p2);
785 return box;
786 }
787