1 // This is bbas/bpgl/bpgl_geotif_camera.hxx
2 #ifndef bpgl_geotif_camera_hxx_
3 #define bpgl_geotif_camera_hxx_
4 //:
5 // \file
6 #include <iostream>
7 #include <cmath>
8 #include <vector>
9 #include <fstream>
10 #include "bpgl_geotif_camera.h"
11 #include <vpgl/vpgl_utm.h>
12 #include <vpgl/vpgl_rational_camera.h>
13 #include <vpgl/vpgl_local_rational_camera.h>
14 #include <vpgl/algo/vpgl_backproject.h>
15 #include <vil/file_formats/vil_tiff.h>
16 #include <vgl/vgl_point_2d.h>
17 #include "bpgl_lon_lat_camera.h"
18 #include <vnl/vnl_inverse.h>
19
20 #ifdef _MSC_VER
21 # include <vcl_msvc_warnings.h>
22 #endif
23
24 template <class T>
construct_matrix(T sx,T sy,T sz,std::vector<std::vector<T>> tiepoints)25 bool bpgl_geotif_camera<T>::construct_matrix(T sx, T sy, T sz, std::vector<std::vector<T> > tiepoints){
26 matrix_.set_size(4, 4);
27 T fv = T(0);
28 matrix_.fill(fv);
29 // use tiepoints and scale values to create a transformation matrix
30 // for now use the first tiepoint if there are more than one
31 bool valid_matrix = (tiepoints.size() > 0) && (tiepoints[0].size() == 6);
32 if (!valid_matrix) {
33 std::cerr << "GEOTIFF matrix not valid" << std::endl;
34 return false;
35 }
36 T I = tiepoints[0][0];
37 T J = tiepoints[0][1];
38 T K = tiepoints[0][2];
39 T X = tiepoints[0][3];
40 T Y = tiepoints[0][4];
41 T Z = tiepoints[0][5];
42
43 // Define a transformation matrix as follows:
44 //
45 // |- -|
46 // | Sx 0.0 0.0 Tx |
47 // | | Tx = X - I*Sx
48 // | 0.0 -Sy 0.0 Ty | Ty = Y + J*Sy
49 // | | Tz = Z - K*Sz
50 // | 0.0 0.0 Sz Tz |
51 // | |
52 // | 0.0 0.0 0.0 1.0 |
53 // |- -|
54
55 T Tx = X - I*sx;
56 T Ty = Y + J*sy;
57 T Tz = Z - K*sz;
58
59 matrix_[0][0] = sx;
60 matrix_[1][1] = -1*sy;
61 matrix_[2][2] = sz;
62 matrix_[3][3] = 1.0;
63 matrix_[0][3] = Tx;
64 matrix_[1][3] = Ty;
65 matrix_[2][3] = Tz;
66 return true;
67 }
68
69 template <class T>
init_from_geotif(vil_image_resource_sptr const & resc)70 bool bpgl_geotif_camera<T>::init_from_geotif(vil_image_resource_sptr const& resc){
71 // check if the image is GEOTIFF
72 auto* geotiff_tiff = dynamic_cast<vil_tiff_image*> (resc.ptr());
73 if (!geotiff_tiff) {
74 std::cerr << "Can't cast resource to TIFF" << std::endl;
75 return false;
76 }
77
78 // check if the tiff file is geotiff
79 if (!geotiff_tiff->is_GEOTIFF()) {
80 std::cerr << "TIFF image is not GEOTIFF!\n";
81 return false;
82 }
83
84 // retrieve header
85 vil_geotiff_header* gtif = geotiff_tiff->get_geotiff_header();
86 if (!gtif) {
87 std::cerr << "no geotiff header! "<< std::endl;
88 return false;
89 }
90
91 vil_geotiff_header::GTIF_HEMISPH h;
92
93 std::vector < std::vector<double> > dtiepoints;
94 gtif->gtif_tiepoints(dtiepoints);
95
96 // create a transformation matrix
97 // if there is a transformation matrix in GEOTIFF, use that
98 double* matrix_dvalues;
99 double dsx = 1.0, dsy = 1.0, dsz = 1.0;
100 scale_defined_ = false;
101 if (gtif->gtif_trans_matrix(matrix_dvalues)) {
102 std::cout << "Transfer matrix is given, using that...." << std::endl;
103 T vals[16];
104 for (size_t i = 0; i < 16; ++i)
105 vals[i] = static_cast<T>(matrix_dvalues[i]);
106 matrix_.copy_in(vals);
107 }else if (gtif->gtif_pixelscale(dsx, dsy, dsz)) {
108 scale_defined_ = true;
109 T sx = static_cast<T>(dsx);
110 T sy = static_cast<T>(dsy);
111 T sz = static_cast<T>(dsz);
112 if(sz == T(0)) sz = T(1);
113 if (dtiepoints.size() == 0) {
114 std::cerr << "null tiepoint array - fatal" << std::endl;
115 return false;
116 }
117 std::vector<std::vector<T> > tiepoints(dtiepoints.size(), std::vector<T>(dtiepoints[0].size()));
118 for (size_t r = 0; r < dtiepoints.size(); ++r)
119 for (size_t c = 0; c < dtiepoints[r].size(); ++c)
120 tiepoints[r][c] = static_cast<T>(dtiepoints[r][c]);
121 if(!this->construct_matrix(sx, sy, sz, tiepoints))
122 return false;
123 }
124 else {
125 std::cout << "Transform matrix cannot be formed..\n";
126 return false;
127 }
128
129 // check if the model type is geographic and also the units
130 if (gtif->GCS_WGS84_MET_DEG()) {
131 is_utm_ = false;
132 if(!this->set_spacing_from_wgs_matrix())
133 return false;
134 return true;
135 }
136
137 // otherwise check if it is projected to UTM and figure out the zone
138 if (gtif->PCS_WGS84_UTM_zone(utm_zone_, h) || gtif->PCS_NAD83_UTM_zone(utm_zone_, h))
139 {
140 hemisphere_flag_ = static_cast<int>(h);
141 is_utm_ = true;
142 dsm_spacing_ = matrix_[0][0];
143 return true;
144 }
145 std::cout << "Only PCS_WGS84_UTM , PCS_NAD83_UTM, and GCS_WGS_84 with linear units in meters, "
146 << "angular units in degrees are supported" << std::endl;
147 return false;
148 }
149 template <class T>
set_spacing_from_wgs_matrix()150 bool bpgl_geotif_camera<T>::set_spacing_from_wgs_matrix(){
151 if(is_utm_){
152 std::cerr << "attempting to set dsm_spacing from a utm matrix" << std::endl;
153 return false;
154 }
155 vpgl_lvcs lvcs(matrix_[1][3], matrix_[0][3], double(0),
156 vpgl_lvcs::wgs84, vpgl_lvcs::DEG, vpgl_lvcs::METERS);
157 T lon0, lat0, lon1, lat1;
158 this->dsm_to_global(T(0), T(0), lon0, lat0);
159 this->dsm_to_global(T(100000), T(0), lon1, lat1);
160 double dlx0, dlx1, dly0, dly1, dlz ;
161 lvcs.global_to_local(double(lon0), double(lat0), double(0),
162 vpgl_lvcs::wgs84, dlx0, dly0, dlz);
163 lvcs.global_to_local(double(lon1), double(lat1), double(0),
164 vpgl_lvcs::wgs84, dlx1, dly1, dlz);
165 dsm_spacing_ = static_cast<T>(sqrt((dlx1-dlx0)*(dlx1-dlx0) +
166 (dly1-dly0) * (dly1-dly0)))/T(100000);
167 return true;
168 }
169 // determine the geographic bounds of the camera region of interest
170 template <class T>
geo_bounds_from_rational_cam(vpgl_camera<T> * cam_ptr,vgl_box_2d<T> const & image_bounds,vgl_box_2d<T> & geo_bb,vgl_polygon<T> & geo_boundary)171 bool bpgl_geotif_camera<T>::geo_bounds_from_rational_cam(vpgl_camera<T>* cam_ptr, vgl_box_2d<T> const& image_bounds,
172 vgl_box_2d<T>& geo_bb, vgl_polygon<T>& geo_boundary){
173 vpgl_rational_camera<T>* rat_cam_ptr = dynamic_cast<vpgl_rational_camera<T>*>(cam_ptr);
174 if(!rat_cam_ptr){
175 std::cerr << "null rational camera pointer - fatal" << std::endl;
176 return false;
177 }
178 //cast camera<T> to camera<double> for backproject=======
179 std::vector<vpgl_scale_offset<T> > scl_off = rat_cam_ptr->scale_offsets();
180 std::vector<std::vector<T> > coef = rat_cam_ptr->coefficients();
181 std::vector<vpgl_scale_offset<double> > scale_d;
182 std::vector<std::vector<double> > coef_d;
183 for(typename std::vector<vpgl_scale_offset<T> >::iterator sit = scl_off.begin();
184 sit != scl_off.end(); ++sit){
185 vpgl_scale_offset<T>& so = *sit;
186 vpgl_scale_offset<double> sod(static_cast<double>(so.scale()), static_cast<double>(so.offset()));
187 scale_d.push_back(sod);
188 }
189 size_t np = coef.size();
190 for(size_t i = 0; i<np; ++i){
191 std::vector<T>& pcoef = coef[i];
192 size_t nc = pcoef.size();
193 std::vector<double> pcoef_d;
194 for(size_t j = 0; j<nc; ++j)
195 pcoef_d.push_back(static_cast<double>(pcoef[j]));
196 coef_d.push_back(pcoef_d);
197 }
198 vpgl_rational_camera<double> rat_cam_d(coef_d, scale_d);
199 //=======================================================
200
201 // define ground plane as bottom of valid 3-d region
202 // also define middle of ground plane as the initial guess for
203 // back projection
204 double sz = rat_cam_d.scale(vpgl_rational_camera<double>::Z_INDX);
205 double oz = rat_cam_d.offset(vpgl_rational_camera<double>::Z_INDX);
206 double ox = rat_cam_d.offset(vpgl_rational_camera<double>::X_INDX);
207 double oy = rat_cam_d.offset(vpgl_rational_camera<double>::Y_INDX);
208 double z0 = oz-sz;//min elevation
209 // default image bounds if input bounds are empty
210 vgl_box_2d<T> img_bb = image_bounds;
211 if(img_bb.is_empty()){
212 double su = rat_cam_d.scale(vpgl_rational_camera<double>::U_INDX);
213 double ou = rat_cam_d.offset(vpgl_rational_camera<double>::U_INDX);
214 double sv = rat_cam_d.scale(vpgl_rational_camera<double>::V_INDX);
215 double ov = rat_cam_d.offset(vpgl_rational_camera<double>::V_INDX);
216 double u0 = ou-su, v0 = ov-sv, u1 = ou+su, v1 = v0;
217 double u2 = u1, v2 = ov+sv,u3 = u0, v3 = v2;
218 vgl_point_2d<T> b0(u0, v0), b1(u1,v1), b2(u2, v2), b3(u3,v3);
219 img_bb.add(b0);img_bb.add(b1);img_bb.add(b2);img_bb.add(b3);
220 }
221 double minx = img_bb.min_x(), miny = img_bb.min_y();
222 double maxx = img_bb.max_x(), maxy = img_bb.max_y();
223 vgl_point_2d<double> p0(minx, miny), p1(maxx, miny), p2(maxx, maxy), p3(minx, maxy);
224 vgl_point_3d<double> P0, P1, P2, P3, Pi(ox, oy, z0);
225 bool good = true;
226 vgl_plane_3d<double> pl(T(0), T(0), T(1), -z0);
227 good = good && vpgl_backproject::bproj_plane(rat_cam_d, p0, pl, Pi, P0);
228 good = good && vpgl_backproject::bproj_plane(rat_cam_d, p1, pl, Pi, P1);
229 good = good && vpgl_backproject::bproj_plane(rat_cam_d, p2, pl, Pi, P2);
230 good = good && vpgl_backproject::bproj_plane(rat_cam_d, p3, pl, Pi, P3);
231 if(!good){
232 std::cerr << "Backprojection failed in geotif bounds" << std::endl;
233 return false;
234 }
235 std::vector<vgl_point_2d<T> > verts;
236 verts.emplace_back(P0.x(), P0.y()); verts.emplace_back(P1.x(), P1.y());
237 verts.emplace_back(P2.x(), P2.y()); verts.emplace_back(P3.x(), P3.y());
238 geo_boundary.push_back(verts);
239 geo_bb.add(verts[0]); geo_bb.add(verts[1]);
240 geo_bb.add(verts[2]); geo_bb.add(verts[3]);
241 return true;
242 }
243 // determine the geographic bounds of the camera region of interest
244 template <class T>
geo_bounds_from_local_cam(std::shared_ptr<vpgl_camera<T>> const & lcam_ptr)245 bool bpgl_geotif_camera<T>::geo_bounds_from_local_cam(std::shared_ptr<vpgl_camera<T> >const& lcam_ptr){
246 if(image_bounds_.is_empty()){
247 std::cerr << "can't define geo bounds for a local camera without image bounds" << std::endl;
248 return false;
249 }
250 if(!lvcs_ptr_){
251 std::cout << "can't define geo bounds for a local camers if lvcs_ptr is null" << std::endl;
252 return false;
253 }
254 // cast the local camera to vpgl_camera<double>
255 vpgl_camera<double>* cam_d;
256 vnl_matrix_fixed<double, 3, 4> md;
257 if(lcam_ptr->type_name() == "vpgl_proj_camera" || lcam_ptr->type_name() == "vpgl_affine_camera"){
258 vpgl_proj_camera<T>* pcam = dynamic_cast<vpgl_proj_camera<T>*>(lcam_ptr.get());
259 const vnl_matrix_fixed<T, 3, 4>& m = pcam->get_matrix();
260 for(size_t r = 0; r<3; ++r)
261 for(size_t c = 0; c<4; ++r)
262 md[r][c] = static_cast<double>(m[r][c]);
263 cam_d = new vpgl_proj_camera<double>(md);
264 }
265 if(!cam_d){
266 std::cerr << "local camera is not projective or a subclass of projective" << std::endl;
267 return false;
268 }
269 // the ground plane for a local camera is defined at z = 0
270 vgl_plane_3d<double> gpl(0.0, 0.0, 1.0, 0.0);
271 double minx = image_bounds_.min_x(), miny = image_bounds_.min_y();
272 double maxx = image_bounds_.max_x(), maxy = image_bounds_.max_y();
273 vgl_point_2d<double> p0(minx, miny), p1(maxx, miny), p2(maxx, maxy), p3(minx, maxy);
274 vgl_point_3d<double> P0, P1, P2, P3, Pi(0.0, 0.0, 0.0);
275 bool good = true;
276 good = good && vpgl_backproject::bproj_plane(cam_d, p0, gpl, Pi, P0);
277 good = good && vpgl_backproject::bproj_plane(cam_d, p1, gpl, Pi, P1);
278 good = good && vpgl_backproject::bproj_plane(cam_d, p2, gpl, Pi, P2);
279 good = good && vpgl_backproject::bproj_plane(cam_d, p3, gpl, Pi, P3);
280 if(!good){
281 std::cerr << "Backprojection failed in geotif bounds" << std::endl;
282 delete cam_d;
283 return false;
284 }
285 //convert local intersection points to global pts using the lvcs
286 vgl_point_2d<T> gp0, gp1, gp2, gp3;
287 T z0;
288 this->local_to_global(P0.x(), P0.y(), 0.0, gp0.x(), gp0.y(), z0);
289 this->local_to_global(P1.x(), P1.y(), 0.0, gp1.x(), gp1.y(), z0);
290 this->local_to_global(P2.x(), P2.y(), 0.0, gp2.x(), gp2.y(), z0);
291 this->local_to_global(P3.x(), P3.y(), 0.0, gp3.x(), gp3.y(), z0);
292 std::vector<vgl_point_2d<T> > verts;
293 verts.emplace_back(gp0.x(), gp0.y()); verts.emplace_back(gp1.x(), gp1.y());
294 verts.emplace_back(gp2.x(), gp2.y()); verts.emplace_back(gp3.x(), gp3.y());
295 geo_boundary_.push_back(verts);
296 geo_bb_.add(verts[0]); geo_bb_.add(verts[1]);
297 geo_bb_.add(verts[2]); geo_bb_.add(verts[3]);
298 delete cam_d;
299 return true;
300 }
301 //--------------------------------------
302 // factory constructors
303 //
304 template <class T>
construct_from_geotif(vpgl_camera<T> const & general_cam,vil_image_resource_sptr resc,vgl_box_2d<T> const & image_bounds,bool elev_org_at_zero,vpgl_lvcs_sptr lvcs_ptr)305 bool bpgl_geotif_camera<T>::construct_from_geotif(vpgl_camera<T> const& general_cam, vil_image_resource_sptr resc,
306 vgl_box_2d<T> const& image_bounds, bool elev_org_at_zero,
307 vpgl_lvcs_sptr lvcs_ptr) {
308 if(resc == nullptr){
309 std::cerr << "null geotiff resource - can't proceed" << std::endl;
310 return false;
311 }
312 projection_enabled_ = true;
313 elev_org_at_zero_ = elev_org_at_zero;
314 has_lvcs_ = static_cast<bool>(lvcs_ptr);
315 dsm_spacing_ = T(1);
316 if(image_bounds_.is_empty()){
317 T ni = static_cast<T>(resc->ni()), nj = static_cast<T>(resc->ni());
318 vgl_point_2d<T> min_p(T(0), T(0)), max_p(ni, nj);
319 image_bounds_.add(min_p); image_bounds_.add(max_p);
320 }
321 //case I - a camera which projects global geo coordinates - no local CS
322 // TO DO - there might exist a UTM camera in the future
323 if(general_cam.type_name() == "vpgl_rational_camera"||general_cam.type_name() == "bpgl_lon_lat_camera"&& !lvcs_ptr){
324 general_cam_ = std::shared_ptr<vpgl_camera<T> >(general_cam.clone());
325 gcam_has_wgs84_cs_= true;
326 project_local_points_ = false;
327 bpgl_lon_lat_camera<T>* lon_lat_ptr = dynamic_cast<bpgl_lon_lat_camera<T>*>(general_cam_.get());
328 if (lon_lat_ptr) {
329 geo_bb_ = lon_lat_ptr->geo_bb();
330 geo_boundary_ = lon_lat_ptr->geo_boundary();
331 }
332 vpgl_rational_camera<T>* rat_cam_ptr = dynamic_cast<vpgl_rational_camera<T>*>(general_cam_.get());
333 // backproject image bounds onto a horizontal plane to get bounds
334 if(rat_cam_ptr)
335 if(!geo_bounds_from_rational_cam(rat_cam_ptr, image_bounds_, geo_bb_, geo_boundary_)){
336 std::cerr << "failed to get geo bounds from rational camera" << std::endl;
337 return false;
338 }
339 this->init_from_geotif(resc);
340 return true;
341 }
342 // Case II - camera is local rational camera and points are in its local CS
343 // possibly with a elevation offset required if points have a z=0 reference elevation
344 // rather than global z
345 if(general_cam.is_a() == "vpgl_local_rational_camera"&&!lvcs_ptr){
346 const vpgl_local_rational_camera<T>& lrcam = dynamic_cast<const vpgl_local_rational_camera<T>&>(general_cam);
347 gcam_has_wgs84_cs_= true;
348 lvcs_ptr_ = new vpgl_lvcs(lrcam.lvcs());
349 has_lvcs_ = true;
350 project_local_points_ = true;
351 const vpgl_rational_camera<T>& rcam = dynamic_cast<const vpgl_rational_camera<T>&>(lrcam);
352 vpgl_rational_camera<T>* rcam_ptr = new vpgl_rational_camera<T>(rcam);
353 // backproject image bounds onto a horizontal plane to get bounds
354 if (!geo_bounds_from_rational_cam(rcam_ptr, image_bounds_, geo_bb_, geo_boundary_)) {
355 std::cerr << "failed to get geo bounds from rational camera" << std::endl;
356 return false;
357 }
358 delete rcam_ptr;
359 general_cam_ = std::shared_ptr<vpgl_camera<T> >(new vpgl_rational_camera<T>(rcam));
360 this->init_from_geotif(resc);
361 return true;
362 }
363 // Case III - camera is local and lvcs is specified, input points are in the global CS of the lvcs
364 // TODO define geo bounds
365 if(lvcs_ptr){
366 has_lvcs_ = true;
367 lvcs_ptr_ = lvcs_ptr;
368 gcam_has_wgs84_cs_= false;
369 project_local_points_ = false;
370 general_cam_ = std::shared_ptr<vpgl_camera<T> >(general_cam.clone());
371 if (!geo_bounds_from_local_cam(general_cam_)) {
372 std::cerr << "failed to get geo bounds from local camera" << std::endl;
373 return false;
374 }
375 this->init_from_geotif(resc);
376 return true;
377 }
378
379 std::cout << "camera is local but no lvcs specified - fatal" << std::endl;
380 return false;
381 }
382 template <class T>
construct_from_matrix(vpgl_camera<T> const & general_cam,vnl_matrix<T> const & geo_transform_matrix,vgl_box_2d<T> const & image_bounds,bool elev_org_at_zero,vpgl_lvcs_sptr lvcs_ptr,int hemisphere_flag,int zone)383 bool bpgl_geotif_camera<T>::construct_from_matrix(vpgl_camera<T> const& general_cam, vnl_matrix<T> const& geo_transform_matrix,
384 vgl_box_2d<T> const& image_bounds, bool elev_org_at_zero,
385 vpgl_lvcs_sptr lvcs_ptr, int hemisphere_flag , int zone)
386 {
387 projection_enabled_ = true;
388 scale_defined_ = true;
389 elev_org_at_zero_ = elev_org_at_zero;
390 lvcs_ptr_ = lvcs_ptr;
391 if(general_cam.is_a() == "vpgl_rational_camera"&& !lvcs_ptr){
392 project_local_points_ = false;
393 general_cam_ = std::shared_ptr<vpgl_camera<T> >(general_cam.clone());
394 gcam_has_wgs84_cs_= true;
395 matrix_ = geo_transform_matrix;
396 if(hemisphere_flag>=0){
397 is_utm_ = true;
398 utm_zone_ = zone;
399 hemisphere_flag_ = hemisphere_flag;
400 dsm_spacing_ = geo_transform_matrix[0][0];
401 return true;
402 }
403 return set_spacing_from_wgs_matrix();
404 }
405 // Case II - camera is local rational camera and points are in its local CS
406 // possibly with a elevation offset required if points have a z=0 reference elevation
407 // rather than global z
408 if(general_cam.is_a() == "vpgl_local_rational_camera"&&!lvcs_ptr){
409 const vpgl_local_rational_camera<T>& lrcam = dynamic_cast<const vpgl_local_rational_camera<T>&>(general_cam);
410 vpgl_lvcs_sptr lr_lvcs_ptr = new vpgl_lvcs(lrcam.lvcs());
411 const vpgl_rational_camera<T>& rcam = dynamic_cast<const vpgl_rational_camera<T>&>(lrcam);
412 project_local_points_ = true;
413 general_cam_ = std::shared_ptr<vpgl_camera<T> >(new vpgl_rational_camera<T>(rcam));
414 has_lvcs_ = true;
415 gcam_has_wgs84_cs_= true;
416 lvcs_ptr_ = lr_lvcs_ptr;
417 matrix_ = geo_transform_matrix;
418 if(hemisphere_flag >=0){
419 is_utm_ = true;
420 utm_zone_ = zone;
421 hemisphere_flag_ = hemisphere_flag;
422 dsm_spacing_ = geo_transform_matrix[0][0];
423 return true;
424 }
425 return set_spacing_from_wgs_matrix();
426 }
427 // Case III - camera is local and lvcs is specified - input points are in a global geographic CS
428 if(lvcs_ptr){
429 lvcs_ptr_ = lvcs_ptr;
430 has_lvcs_ = true;
431 gcam_has_wgs84_cs_= false;
432 general_cam_ = std::shared_ptr<vpgl_camera<T> >(general_cam.clone());
433 matrix_ = geo_transform_matrix;
434 project_local_points_ = false;
435 if(hemisphere_flag >=0){
436 is_utm_ = true;
437 utm_zone_ = zone;
438 hemisphere_flag_ = hemisphere_flag;
439 dsm_spacing_ = geo_transform_matrix[0][0];
440 return true;
441 }
442 return set_spacing_from_wgs_matrix();
443 }
444 std::cout << "camera is local but no lvcs specified - fatal" << std::endl;
445 return false;
446 }
447 template <class T>
construct_geo_data_only(vil_image_resource_sptr resc)448 bool bpgl_geotif_camera<T>::construct_geo_data_only(vil_image_resource_sptr resc){
449 if(!resc){
450 std::cerr << "null DSM resource - can't proceed" << std::endl;
451 return false;
452 }
453 projection_enabled_ = false;
454 T tni = static_cast<T>(resc->ni()), tnj = static_cast<T>(resc->nj());
455 vgl_point_2d<T> minp(T(0), T(0));
456 vgl_point_2d<T> maxp(tni, tnj);
457 image_bounds_.add(minp);image_bounds_.add(maxp);
458 return this->init_from_geotif(resc);
459 }
460
461 template <class T>
local_to_global(T lx,T ly,T lz,T & gx,T & gy,T & gz) const462 bool bpgl_geotif_camera<T>::local_to_global(T lx, T ly, T lz, T& gx, T& gy, T& gz) const{
463 double dlx = static_cast<double>(lx);
464 double dly = static_cast<double>(ly);
465 double dlz = static_cast<double>(lz);
466 double dgx, dgy, dgz;
467 if (lvcs_ptr_) {
468 if (lvcs_ptr_->get_cs_name() == vpgl_lvcs::utm) {
469 if (is_utm_) { // this CS is utm so keep using utm
470 lvcs_ptr_->local_to_global(dlx, dly, dlz, vpgl_lvcs::utm, dgx, dgy, dgz);
471 }else { // this CS is not UTM, convert to wgs84 as global
472 lvcs_ptr_->local_to_global(dlx, dly, dlz, vpgl_lvcs::wgs84, dgx, dgy,dgz);
473 }
474 }else {
475 lvcs_ptr_->local_to_global(dlx, dly, dlz, vpgl_lvcs::wgs84, dgx, dgy, dgz);
476 }
477 gx = static_cast<T>(dgx);
478 gy = static_cast<T>(dgy);
479 gz = static_cast<T>(dgz);
480 return true;
481 }
482 std::cerr << "global coordinate system not defined" << std::endl;
483 return false;
484 }
485 template <class T>
global_to_local(T gx,T gy,T gz,T & lx,T & ly,T & lz) const486 bool bpgl_geotif_camera<T>::global_to_local(T gx, T gy, T gz, T& lx, T& ly, T& lz) const {
487 if(!lvcs_ptr_){
488 std::cout << "No local vertical CS defined - can't map global to local" << std::endl;
489 return false;
490 }
491 double dgx = static_cast<double>(gx);
492 double dgy = static_cast<double>(gy);
493 double dgz = static_cast<double>(gz);
494 double dlx, dly, dlz;
495 bool good = false;
496 if(lvcs_ptr_->get_cs_name() == vpgl_lvcs::utm){
497 lvcs_ptr_->global_to_local(dgx, dgy, dgz, vpgl_lvcs::utm, dlx, dly, dlz);
498 good = true;
499 } else if(lvcs_ptr_->get_cs_name() == vpgl_lvcs::wgs84){
500 lvcs_ptr_->global_to_local(dgx, dgy, dgz, vpgl_lvcs::wgs84, dlx, dly, dlz);
501 good = true;
502 }
503 if (good) {
504 lx = static_cast<T>(dlx);
505 ly = static_cast<T>(dly);
506 lz = static_cast<T>(dlz);
507 return true;
508 }
509 std::cout << "lvcs name " << lvcs_ptr_->get_cs_name() << " not handled " << std::endl;
510 return false;
511 }
512 template <class T>
project(const T x,const T y,const T z,T & u,T & v) const513 void bpgl_geotif_camera<T>::project(const T x, const T y, const T z, T& u, T& v) const{
514 if(!projection_enabled_){
515 std::cerr << "Warning trying to project without a general camera defined" << std::endl;
516 return;
517 }
518 // case I - The camera and the input points have the same CS
519 if(!has_lvcs_){
520 general_cam_->project(x, y, z, u, v);
521 return;
522 }
523 // case II - the camera is global WGS 84,
524 // and input points are in a local CS
525 if (has_lvcs_&& gcam_has_wgs84_cs_) {
526 T gx, gy, gz;
527 T zadj = z;
528 if(!elev_org_at_zero_){
529 T elev_org = this->elevation_origin();
530 zadj = static_cast<T>(zadj - elev_org);
531 }
532 this->local_to_global(x, y, zadj, gx, gy, gz);
533 general_cam_->project(T(gx), T(gy), T(gz), u, v);
534 return;
535 }
536 // case III - the camera is global
537 // and input points are in a global CS
538 if (has_lvcs_) {
539 T lx, ly, lz;
540 if(!this->global_to_local(x, y, z, lx, ly, lz)){
541 u = T(-1); v = T(-1);
542 return;
543 }
544 general_cam_->project(T(lx), T(ly), T(lz), u, v);
545 return;
546 }
547 std::cerr << "invalid projection" << std::endl;
548 u = T(-1); v = T(-1);
549 }
550 //
551 // use the 4x4 GEOTIFF matrix to transform DSM image coordinates to global coordinates
552 // either UTM or WGS84
553 template <class T>
dsm_to_global(T i,T j,T & gx,T & gy) const554 void bpgl_geotif_camera<T>::dsm_to_global(T i, T j, T& gx, T& gy) const{
555 vnl_vector<T> v(4), res(4);
556 if (scale_defined_) {
557 v[0] = matrix_[0][3] + i*matrix_[0][0];
558 v[1] = matrix_[1][3] + j*matrix_[1][1];
559 }
560 else {
561 v[0] = matrix_[0][3] + i;
562 v[1] = matrix_[1][3] - j;
563 }
564 gx = v[0]; gy = v[1];
565 }
566 //
567 // use the 4x4 GEOTIFF matrix to transform global coordinates to DSM image coordinates
568 // either UTM or WGS84
569 template <class T>
global_to_dsm(T gx,T gy,T & i,T & j) const570 void bpgl_geotif_camera<T>::global_to_dsm(T gx, T gy, T& i, T& j) const{
571 vnl_vector<T> vec(4,T(0)), res(4);
572 T X=gx, Y=gy;
573 if (is_utm_) {
574 double deast = static_cast<double>(X);
575 double dnorth = static_cast<double>(Y);
576 double dlat, dlon, delev;
577 bool south_flag = hemisphere_flag_ > 0;
578 vpgl_utm utm;
579 utm.transform(utm_zone_, deast, dnorth, 0.0, dlat, dlon, delev, south_flag);
580 X = T(dlon); Y = T(dlat);
581 }
582 vec[0] = X;
583 vec[1] = Y;
584 vec[3] = 1;
585
586 if (scale_defined_) {
587 i = (vec[0] - matrix_[0][3])/matrix_[0][0];
588 j = (vec[1] - matrix_[1][3])/matrix_[1][1];
589 }
590 else {//TO DO check this case
591 vnl_matrix<T> matrix_inv = vnl_inverse(matrix_);
592 res = matrix_inv*vec;
593 i = res[0];
594 j = res[1];
595 }
596 }
597
598 // get the global elevation origin to handle the case where the input points
599 // were generated with a local z elevation origin of zero. This reference elevation
600 // is correct if the camera is a local camera by which the input points were
601 // generated. E.g. a local affine camera. Otherwise, the global offset of the lvcs
602 // has to be added to the input point elevations before projecting.
603 template <class T>
elevation_origin() const604 T bpgl_geotif_camera<T>::elevation_origin() const{
605 if(!lvcs_ptr_)
606 return T(0);
607 double ox, oy, oz;
608 if(is_utm_){
609 int temp;
610 lvcs_ptr_->get_utm_origin(ox, oy, oz, temp);
611 return oz;
612 }
613 lvcs_ptr_->get_origin(oy, ox, oz);
614 return oz;
615 }
616
617 template <class T>
project_gtif_to_image(const T tifu,const T tifv,const T tifz,T & u,T & v) const618 void bpgl_geotif_camera<T>::project_gtif_to_image(const T tifu, const T tifv, const T tifz, T& u, T& v) const {
619 if(!projection_enabled_){
620 std::cerr << "Warning trying to project without a general camera defined" << std::endl;
621 return;
622 }
623 // normal situation, wgs84 cs and lvcs
624 if (has_lvcs_ && gcam_has_wgs84_cs_ && !is_utm_) {
625 T lat = 0.0, lon = 0.0;
626 // map dsm image coordinates to lon lat
627 this->dsm_to_global(tifu, tifv, lon, lat);
628 T zadj = tifz;
629 if(elev_org_at_zero_){
630 T elev_org = this->elevation_origin();
631 zadj = static_cast<T>(zadj + elev_org);
632 }
633 general_cam_->project(T(lon), T(lat), zadj, u, v);
634 return;
635 }
636 // wgs84 cs and no lvcs, i.e. a global camera
637 if (!has_lvcs_ && gcam_has_wgs84_cs_ && !is_utm_) {
638 T lat = 0.0, lon = 0.0;
639 // map dsm image coordinates to lon lat
640 this->dsm_to_global(tifu, tifv, lon, lat);
641 T zadj = tifz;
642 if(elev_org_at_zero_){
643 T elev_org = this->elevation_origin();
644 zadj = static_cast<T>(zadj + elev_org);
645 }
646 general_cam_->project(T(lon), T(lat), zadj, u, v);
647 return;
648 }
649 // geotiff is in UTM so conversion is required before projecting
650 // through the WGS84 local camera
651 if (has_lvcs_ && gcam_has_wgs84_cs_ && is_utm_) {
652 T east = 0.0, north = 0.0;
653 this->dsm_to_global(tifu, tifv, east, north);
654 bool south_flag = hemisphere_flag_ > 0;
655 T zadj = tifz;
656 if (elev_org_at_zero_) {
657 T elev_org = this->elevation_origin();
658 zadj = static_cast<T>(zadj + elev_org);
659 }
660 T lat = 0.0, lon = 0.0, elev = 0.0;
661 double deast = static_cast<double>(east);
662 double dnorth = static_cast<double>(north);
663 double dzadj = static_cast<double>(zadj);
664 double dlat, dlon, delev;
665 vpgl_utm utm;
666 utm.transform(utm_zone_, deast, dnorth, dzadj, dlat, dlon, delev, south_flag);
667 general_cam_->project(T(dlon), T(dlat), T(delev), u, v);
668 return;
669 }// TO DO geotiff in UTM but wgs84 global camera with no lvcs
670 std::cout << "improper camera or geographic CS to project DSM points" << std::endl;
671 }
672 template <class T>
lower_left_lvcs(T elev_ll) const673 vpgl_lvcs_sptr bpgl_geotif_camera<T>::lower_left_lvcs(T elev_ll) const{
674 vpgl_lvcs_sptr ret;
675 if(image_bounds_.is_empty()){
676 std::cerr << "can't define lower left origin without image bounds" << std::endl;
677 return ret;
678 }
679 T j_ll = image_bounds_.max_y()-T(1);
680 T i_ll = T(0);
681 T gx, gy;
682 dsm_to_global(i_ll, j_ll, gx, gy);
683 if(!is_utm_){
684 ret = new vpgl_lvcs(gy, gx, elev_ll, vpgl_lvcs::wgs84, vpgl_lvcs::DEG, vpgl_lvcs::METERS);
685 return ret;
686 }
687 double lon, lat;
688 vpgl_utm utm;
689 utm.transform(utm_zone_, gx, gy, lat, lon, hemisphere_flag_>0);
690 ret = new vpgl_lvcs(lat, lon, elev_ll, vpgl_lvcs::utm, vpgl_lvcs::DEG, vpgl_lvcs::METERS);
691 return ret;
692 }
693 // Code for easy instantiation.
694 #undef BPGL_GEOTIF_CAMERA_INSTANTIATE
695 #define BPGL_GEOTIF_CAMERA_INSTANTIATE(T) \
696 template class bpgl_geotif_camera<T >
697
698
699
700 #endif // bpgl_geotif_camera_hxx_
701