1 // This is brl/bseg/boxm2/pro/processes/boxm2_roi_init_geotiff_process.cxx
2 #include <cmath>
3 #include <fstream>
4 #include <ios>
5 #include <iostream>
6 #include <algorithm>
7 #include <bprb/bprb_func_process.h>
8 //:
9 // \file
10 // \brief  A process to retrieve geocam from a geotiff image and lvcs from the scene
11 // Also creates a generic camera
12 // and crops the portion of the scene for the given scene and return it
13 //
14 // \author Ozge C. Ozcanli
15 // \date May 07, 2012
16 
17 #include <boxm2/boxm2_scene.h>
18 
19 #include <vpgl/file_formats/vpgl_geo_camera.h>
20 #include <vpgl/algo/vpgl_camera_convert.h>
21 #include "vil/vil_image_view.h"
22 #include "vil/vil_image_resource.h"
23 #include <vil/file_formats/vil_geotiff_header.h>
24 #include <vil/file_formats/vil_tiff.h>
25 #include "vil/vil_resample_bilin.h"
26 #include "vil/vil_load.h"
27 #ifdef _MSC_VER
28 #  include "vcl_msvc_warnings.h"
29 #endif
30 
31 namespace boxm2_roi_init_geotiff_process_globals
32 {
33   constexpr unsigned n_inputs_ = 4;
34   constexpr unsigned n_outputs_ = 2;
35 }
36 
boxm2_roi_init_geotiff_process_cons(bprb_func_process & pro)37 bool boxm2_roi_init_geotiff_process_cons(bprb_func_process& pro)
38 {
39   using namespace boxm2_roi_init_geotiff_process_globals;
40 
41   //process takes 4 inputs
42   std::vector<std::string> input_types_(n_inputs_);
43   input_types_[0] = "boxm2_scene_sptr";
44   input_types_[1] = "vpgl_camera_double_sptr";  // geocam
45   input_types_[2] = "vcl_string";  // tiff image
46   input_types_[3] = "unsigned";  // resolution level, if 0 orig img res, 1 => downsample by 2^1, 2 => by 2^2
47 
48   // process has 2 outputs:
49   std::vector<std::string>  output_types_(n_outputs_);
50   output_types_[0] = "vpgl_camera_double_sptr";   // output generic cam
51   output_types_[1] = "vil_image_view_base_sptr";
52 
53   return pro.set_input_types(input_types_) && pro.set_output_types(output_types_);
54 }
55 
boxm2_roi_init_geotiff_process(bprb_func_process & pro)56 bool boxm2_roi_init_geotiff_process(bprb_func_process& pro)
57 {
58   using namespace boxm2_roi_init_geotiff_process_globals;
59 
60   if ( pro.n_inputs() < n_inputs_ ) {
61     std::cout << pro.name() << ": The input number should be " << n_inputs_<< std::endl;
62     return false;
63   }
64   //get the inputs
65   boxm2_scene_sptr scene = pro.get_input<boxm2_scene_sptr>(0);
66   vpgl_camera_double_sptr cam = pro.get_input<vpgl_camera_double_sptr>(1);
67   std::string geotiff_fname = pro.get_input<std::string>(2);
68   auto level = pro.get_input<unsigned>(3);
69 
70   vil_image_resource_sptr img_res = vil_load_image_resource(geotiff_fname.c_str());
71   vpgl_geo_camera* geocam = nullptr;
72   if (cam) {
73     std::cout << "Using the provided loaded camera.\n";
74     geocam = dynamic_cast<vpgl_geo_camera*> (cam.ptr());
75   } else {
76     std::cout << "Using the camera initialized from the image.\n";
77     vpgl_lvcs_sptr lvcs = new vpgl_lvcs(scene->lvcs());
78     vpgl_geo_camera::init_geo_camera(img_res, lvcs, geocam);
79   }
80 
81   if (!geocam) {
82     std::cerr << "In boxm2_roi_init_geotiff_process() - the geocam could not be initialized!\n";
83     return false;
84   }
85 
86   double e1,n1,e2,n2,elev,lat,lon;
87   scene->lvcs().get_origin(lat, lon, elev);
88   geocam->img_four_corners_in_utm(img_res->ni(), img_res->nj(), elev, e1, n1, e2, n2);
89   std::cout << "img res ni: " << img_res->ni() << " nj: " << img_res->nj() << std::endl;
90 
91   vgl_box_3d<double> scene_bbox = scene->bounding_box();
92   std::cout.setf(std::ios::fixed, std::ios::floatfield);
93 
94   // crop the image
95   vgl_box_2d<double> proj_bbox;
96   double u,v;
97   geocam->project(scene_bbox.min_x(), scene_bbox.min_y(), scene_bbox.min_z(), u, v);
98   proj_bbox.add(vgl_point_2d<double>(u,v));
99   geocam->project(scene_bbox.max_x(), scene_bbox.max_y(), scene_bbox.max_z(), u, v);
100   proj_bbox.add(vgl_point_2d<double>(u,v));
101   int min_i = int(std::max(0.0, std::floor(proj_bbox.min_x())));
102   int min_j = int(std::max(0.0, std::floor(proj_bbox.min_y())));
103   int max_i = int(std::min(img_res->ni()-1.0, std::ceil(proj_bbox.max_x())));
104   int max_j = int(std::min(img_res->nj()-1.0, std::ceil(proj_bbox.max_y())));
105   int crop_ni = max_i - min_i + 1;
106   int crop_nj = max_j - min_j + 1;
107   if (crop_ni < 0 || crop_nj < 0) {
108     std::cerr << "Error: boxm2_roi_init_geotiff_process: the image does not contain scene!\n"
109              << "scene bbox:\n " << scene_bbox << '\n';
110     double n,e; int utm_zone;
111     std::cerr.setf(std::ios::fixed, std::ios::floatfield); std::cerr.precision(2);
112     geocam->local_to_utm(scene_bbox.min_x(), scene_bbox.min_y(), scene_bbox.min_z(), e, n, utm_zone);
113     std::cerr << "in UTM bbox transforms to " << n << " N " << e << " E and ";
114     geocam->local_to_utm(scene_bbox.max_x(), scene_bbox.max_y(), scene_bbox.max_z(), e, n, utm_zone);
115     std::cerr << n << " N " << e << " E\n";
116     return false;
117   }
118 
119   vil_image_view_base_sptr crop_base = img_res->get_view((unsigned int)min_i, crop_ni, (unsigned int)min_j, crop_nj);
120   if (!crop_base) {
121     std::cerr << "Error: boxm2_roi_init_geotiff_process: could not crop the scene in this image!\n";
122     return false;
123   }
124 
125   // downsample the image based on the level
126   if (level < 32) {
127     crop_ni /= (1L<<level);
128     crop_nj /= (1L<<level);
129   }
130   else {
131     crop_ni = (unsigned int) (crop_ni*std::pow(0.5,static_cast<double>(level)));
132     crop_nj = (unsigned int) (crop_nj*std::pow(0.5,static_cast<double>(level)));
133   }
134 
135   // now use this camera to generate generic cam
136   vpgl_generic_camera<double> gcam;
137   // make the geocam local
138   geocam->translate(min_i, min_j, 0);
139 
140   vpgl_generic_camera_convert::convert(*geocam, crop_ni, crop_nj, scene_bbox.max_z()+1.0, gcam, level);
141 
142   vpgl_camera_double_sptr out = new vpgl_generic_camera<double>(gcam);
143   pro.set_output_val<vpgl_camera_double_sptr>(0, out);
144 
145   auto* crop_view = dynamic_cast<vil_image_view<vxl_byte>*>(crop_base.ptr());
146   if (!crop_view) {
147     auto* crop_view_f = dynamic_cast<vil_image_view<float>*>(crop_base.ptr());
148     if (!crop_view_f) {
149       std::cerr << "Error: boxm2_roi_init_geotiff_process: could not cast first return image to a vil_image_view<float>\n";
150       return false;
151     }
152     auto* out_img = new vil_image_view<float>(crop_ni, crop_nj, crop_view_f->nplanes());
153     vil_resample_bilin(*crop_view_f, *out_img, crop_ni, crop_nj);
154     pro.set_output_val<vil_image_view_base_sptr>(1,new vil_image_view<float>(*out_img));
155     return true;
156   }
157 
158   //: downsample the image based on the level
159   auto* out_img = new vil_image_view<vxl_byte>(crop_ni, crop_nj, crop_view->nplanes());
160   vil_resample_bilin(*crop_view, *out_img, crop_ni, crop_nj);
161   pro.set_output_val<vil_image_view_base_sptr>(1,new vil_image_view<vxl_byte>(*out_img));
162 
163   return true;
164 }
165