1 // This is bbas/bpgl/bpgl_geotif_camera.h 2 #ifndef bpgl_geotif_camera_h_ 3 #define bpgl_geotif_camera_h_ 4 //: 5 // \file 6 // \brief Converts a vpgl_camera to a geographic camera based on geotiff header information 7 // \author Joseph Mundy 8 // \date November 1, 2019 9 // 10 // The motivation is that 3-d data can be represented as a geotif image, or a point cloud derived from 11 // the geotif image. It is convenient to be able to project a pixel from the geotif image into another image 12 // based on a general camera viewpoint with respect to the geotif data. Similarly, a 3-d point in a point cloud 13 // derived from the geotif data, or in some cases used to create the geotif data, can be projected with respect to a general 14 // camera view. 15 // 16 // the geotif header defines the mapping from geotif image coordinates to 2-d geographic coordinates, 17 // and the image values define the global elevation at each 2-d location, i.e., z(x,y) 18 // 19 // for the point cloud, a lvcs is required to relate the point coordinates to geographic coordinates, not necessarily 20 // in the same as the geographic coordinate system as the geotif image. 21 // 22 #include <iostream> 23 #include <string> 24 #include <memory> 25 #ifdef _MSC_VER 26 # include <vcl_msvc_warnings.h> 27 #endif 28 #include <vnl/vnl_matrix.h> 29 #include <vpgl/vpgl_camera.h> 30 #include <vpgl/file_formats/vpgl_geo_camera.h> 31 #include <vil/vil_image_resource.h> 32 #include <vgl/vgl_box_2d.h> 33 #include <vgl/vgl_polygon.h> 34 // templated since may want to project <float> points 35 template <class T> 36 class bpgl_geotif_camera : vpgl_camera<T> 37 { 38 public: 39 //: default constructor - may want a container of cameras 40 // the default member values represent the most common case, e.g. a local_rational_camera bpgl_geotif_camera()41 bpgl_geotif_camera():has_lvcs_(true), gcam_has_wgs84_cs_(true), 42 elev_org_at_zero_(true), is_utm_(false), project_local_points_(true),scale_defined_(false), projection_enabled_(false){} 43 44 virtual ~bpgl_geotif_camera() = default; 45 46 //: factory methods to allow for failure conditions 47 48 // if the lvcs is not defined (null) then the general camera must be a rational camera with WGS84 CS, 49 // or a local rational camera with an internal lvcs that defines the local Cartesian CS 50 // Notes: 1) input points may have been generated with local CS elevation origin = 0. In this case 51 // the local rational camera lvcs elevation origin is correct, and is used by default. 52 // If the elevation origin is global due to post processing, elev_org_at_zero = false and 53 // the necessary adjustment to elevation values is made. 54 // 2) the lvcs CS may not match the GEOTIFF header CS so extra internal conversion may be necessary, 55 // e.g., the lvcs CS is WGS84 and the GEOTIFF header CS is UTM 56 // 3) The geographic bounds of a geotif camera are determined from: a) the image bounds or b)from the DSM bounds or 57 // c) from RPC scale and offset bounds, in this order, when defined 58 // 59 bool construct_from_geotif(vpgl_camera<T> const& general_cam, vil_image_resource_sptr resc, 60 vgl_box_2d<T> const& image_bounds = vgl_box_2d<T>(), bool elev_org_at_zero = false, 61 vpgl_lvcs_sptr lvcs_ptr = nullptr); 62 63 //: construct from a 4x4 transform matrix. 64 // if the lvcs is not defined then the camera must be a RPC camera in the WGS84 CS 65 // the UTM CS case is indicated by hemisphere flag being either 0 - Northern Hemisphere or 1 - Southern Hemisphere 66 // when the lvcs_ptr is null, UTM coordinates must be converted to WGS84 before projecting through the camera 67 // The geographic bounds of a geotif camera are determined from: a) the image bounds or b)from the DSM bounds or 68 // c) from RPC scale and offset bounds, in this order, when defined 69 bool construct_from_matrix(vpgl_camera<T> const& general_cam, vnl_matrix<T> const& geo_transform_matrix, 70 vgl_box_2d<T> const& image_bounds = vgl_box_2d<T>(), bool elev_org_at_zero = false, 71 vpgl_lvcs_sptr lvcs_ptr = nullptr,int hemisphere_flag = -1, int zone = 0); 72 73 //: construct without a camera to provide geographic transforms 74 // when projection functions are not needed 75 bool construct_geo_data_only(vil_image_resource_sptr resc); 76 77 //: are the input points in a local CS ? project_local_points()78 bool project_local_points() const {return project_local_points_;} 79 80 //: project from local or global 3-d coordinates, to an image location (u, v) 81 // if project_local_points() == true, coordinates are in a local CS otherwise in a global CS 82 virtual void project(const T x, const T y, const T z, T& u, T& v) const; 83 84 //: project from an image location in the GEOTIFF image to an image location (u, v), 85 // as projected by the general camera. tifz is the height value of the GEOTIFF image at (tifu, tifv) 86 void project_gtif_to_image(const T tifu, const T tifv, const T tifz , T& u, T& v) const; 87 88 //: legal C++ because the return type is covariant with vpgl_camera<T>* clone(void)89 virtual bpgl_geotif_camera<T>* clone(void) const {return new bpgl_geotif_camera<T>(*this);} 90 91 //: accessors 92 // in some versions of DSM algorithms the elevations are relative to 93 // the local CS, i.e. the origin for z is 0; elevation_origin_at_zero()94 bool elevation_origin_at_zero() const {return elev_org_at_zero_;} 95 96 //: is the cs of the DSM UTM? is_utm()97 bool is_utm() const {return is_utm_;} 98 99 //: does the general camera used for projection have a lvcs? has_lvcs()100 bool has_lvcs() const {return has_lvcs_;} 101 102 //: is the general camera coordinate system wgs84? general_cam_has_wgs84_cs()103 bool general_cam_has_wgs84_cs() const {return gcam_has_wgs84_cs_;} 104 105 //: return a pointer to the general camera, e.g. a local rational camera general_camera()106 std::shared_ptr<vpgl_camera<T> > general_camera() const {return general_cam_;} 107 108 //: return the lvcs pointer if defined, otherwise null lvcs_ptr()109 vpgl_lvcs_sptr lvcs_ptr() {return lvcs_ptr_;} 110 111 //: the geographic matrix extracted from the geotiff header matrix()112 vnl_matrix<T> matrix() const {return matrix_;} 113 114 //: spacing between dsm samples in meters dsm_spacing()115 T dsm_spacing() const{return dsm_spacing_;} 116 117 //: utm CS information utm_zone()118 int utm_zone() const {return utm_zone_;} hemisphere_flag()119 int hemisphere_flag() const {return hemisphere_flag_;} 120 121 //: construct a lvcs at the lower_left corner of the geotiff image array 122 // can be either wgs84 or utm depending on CS of geotiff image 123 vpgl_lvcs_sptr lower_left_lvcs(T elev_ll = T(0)) const; 124 125 //: if not empty, defines the geographic region of validity of the camera 126 // both a polygonal region and an enclosing bounding box are provided geo_bb()127 vgl_box_2d<T> geo_bb() const {return geo_bb_;} geo_boundary()128 vgl_polygon<T> geo_boundary() const {return geo_boundary_;} 129 130 //: =====Transforms between DSM image and global geo coordinates==== 131 // [e.g., (u,v)->(lon, lat) or (lon, lat, elev)->(u,v)] 132 133 //: map dsm image location to global geo X-Y coordinates (wgs84 or UTM) 134 void dsm_to_global(T i, T j, T& gx, T& gy) const; 135 136 //: map global geo X-Y (wgs84 or UTM) to dsm u,v (uses GEOTIFF matrix) 137 void global_to_dsm(T gx, T gy, T& i, T& j) const; 138 139 //: intialize geographic info from the geotiff header 140 // camera projection not enabled but the above functions are available 141 bool init_from_geotif(vil_image_resource_sptr const& resc); 142 143 //===================================================================== 144 protected: 145 //: extract the geographic matrix from the geotiff header 146 bool construct_matrix(T sx, T sy, T sz, std::vector<std::vector<T> > tiepoints); 147 148 //: map local coordinate to global coordinates 149 bool local_to_global(T lx, T ly, T lz, T& gx, T& gy, T& gz) const; 150 //: the opposite direction 151 bool global_to_local(T gx, T gy, T gz, T& lx, T& ly, T& lz) const; 152 153 //: the elevation value at the CS origin 154 T elevation_origin() const; 155 156 //: set the DSM spacing using the wgs CS - convert deg./pix to meters/pix 157 bool set_spacing_from_wgs_matrix(); 158 159 // set the geographic bounds given a rational camera 160 static bool geo_bounds_from_rational_cam(vpgl_camera<T>* rat_cam_ptr, vgl_box_2d<T> const& image_bounds, 161 vgl_box_2d<T>& geo_bb, vgl_polygon<T>& geo_boundary); 162 163 // set the geographic bounds when the general camera is a local camera (e.g. affine) 164 bool geo_bounds_from_local_cam(std::shared_ptr<vpgl_camera<T> >const& lcam_ptr); 165 166 //:a general camera is defined, e.g. RPC 167 bool projection_enabled_; 168 169 //: the dsm grid spacing in meters 170 T dsm_spacing_; 171 172 //: if true, the input points are in a local CS, or if false in a global CS 173 bool project_local_points_; 174 175 //: the DSM is constructed with a zero elevation origin, i.e. local Cartesian elevation 176 bool elev_org_at_zero_; 177 178 //: the points to be projected are in a local CS 179 bool has_lvcs_; 180 181 //: the general camera to be used to project 182 std::shared_ptr<vpgl_camera<T> > general_cam_; 183 184 //: the general camera projects using WGS84 global coordinates 185 bool gcam_has_wgs84_cs_; 186 187 //: the local vertical CS to convert local to global 188 vpgl_lvcs_sptr lvcs_ptr_; 189 190 vnl_matrix<T> matrix_; 191 bool scale_defined_; 192 bool is_utm_; 193 int utm_zone_; 194 int hemisphere_flag_; //0 North, 1 South 195 vgl_box_2d<T> image_bounds_; 196 vgl_box_2d<T> geo_bb_; 197 vgl_polygon<T> geo_boundary_; 198 }; 199 200 #define BPGL_GEOTIF_CAMERA_INSTANTIATE(T) extern "please include vgl/bpgl_geotif_camera.txx first" 201 202 203 #endif // bpgl_geotif_camera_h_ 204