1 // This is brl/bbas/volm/volm_geo_index2.cxx
2 #include "volm_geo_index2.h"
3 //:
4 // \file
5 #include <bkml/bkml_write.h>
6 #include "vgl/vgl_intersection.h"
7 #include "vpgl/vpgl_utm.h"
8 #include "volm_io.h"
9 #include "vgl/vgl_area.h"
10
11 #include <utility>
12
13 // function to check whether the given box intersect with a line defined by a vector of points
is_intersect(vgl_box_2d<double> const & box,std::vector<vgl_point_2d<double>> const & line)14 static bool is_intersect(vgl_box_2d<double> const& box, std::vector<vgl_point_2d<double> > const& line)
15 {
16 if (box.min_x() == 0.375 && box.min_y() == 0.5)
17 unsigned i = 1;
18 unsigned n_line_segs = line.size()-1;
19 for (unsigned i = 0; i < n_line_segs; i++) {
20 vgl_point_2d<double> s = line[i]; vgl_point_2d<double> e = line[i+1];
21 vgl_line_segment_2d<double> line_segment(s,e);
22 vgl_point_2d<double> p0, p1;
23 if (vgl_intersection(box, line_segment, p0, p1) != 0)
24 return true;
25 if (box.contains(s) || box.contains(e))
26 return true;
27 }
28 return false;
29 }
30
31 // write node to kml
write_to_kml_node(std::ofstream & ofs,const volm_geo_index2_node_sptr & n,unsigned current_depth,unsigned depth,std::string explanation)32 void volm_geo_index2::write_to_kml_node(std::ofstream& ofs, const volm_geo_index2_node_sptr& n, unsigned current_depth, unsigned depth, std::string explanation)
33 {
34 if (!n)
35 return;
36 if (current_depth == depth) {
37 // note that in the extent bounding box, x -- lon, y -- lat (kml format: lat lon)
38 vnl_double_2 ul, ll, lr, ur;
39 ll[0] = n->extent_.min_point().y(); ll[1] = n->extent_.min_point().x();
40 ul[0] = n->extent_.max_point().y(); ul[1] = n->extent_.min_point().x();
41 lr[0] = n->extent_.min_point().y(); lr[1] = n->extent_.max_point().x();
42 ur[0] = n->extent_.max_point().y(); ur[1] = n->extent_.max_point().x();
43 bkml_write::write_box(ofs, " ", std::move(explanation), ul, ur,ll,lr);
44 }
45 else {
46 for (const auto & c_idx : n->children_)
47 write_to_kml_node(ofs, c_idx, current_depth+1, depth);
48 }
49 }
50
prune_tree(const volm_geo_index2_node_sptr & root,vgl_polygon<double> const & poly)51 bool volm_geo_index2::prune_tree(const volm_geo_index2_node_sptr& root, vgl_polygon<double> const& poly)
52 {
53 // note that the tree will not be pruned if the root bounding box does not intersect with the polygon
54 if (!vgl_intersection(root->extent_, poly))
55 return false;
56
57 for (auto & i : root->children_) {
58 if (!i)
59 continue;
60 if (!prune_tree(i, poly)) // the child does not intersect with the polygon
61 i = nullptr; // sptr de-allocates this child
62 }
63 return true;
64 }
65
prune_tree(const volm_geo_index2_node_sptr & root,vgl_polygon<float> const & poly)66 bool volm_geo_index2::prune_tree(const volm_geo_index2_node_sptr& root, vgl_polygon<float> const& poly)
67 {
68 vgl_polygon<double> poly_double;
69 volm_io::convert_polygons(poly, poly_double);
70 return volm_geo_index2::prune_tree(root, poly_double);
71 }
72
prune_by_zone(const volm_geo_index2_node_sptr & root,unsigned utm_zone)73 bool volm_geo_index2::prune_by_zone(const volm_geo_index2_node_sptr& root, unsigned utm_zone)
74 {
75 // note the tree will not be pruned if the root is outside the utm_zone;
76 vpgl_utm u;
77 int zone1, zone2; double x, y;
78 u.transform(root->extent_.min_point().y(), root->extent_.min_point().x(), x, y, zone1);
79 u.transform(root->extent_.max_point().y(), root->extent_.max_point().x(), x, y, zone2);
80 if (zone1 != (int)utm_zone && zone2 != (int)utm_zone) // the whole quadtree is outside utm_zone
81 return false;
82
83 for (auto & i : root->children_) {
84 if (!i)
85 continue;
86 if (!prune_by_zone(i, utm_zone)) // the child i is not in the utm_zone
87 i = nullptr; // sptr deallocates this child
88 }
89 return true;
90 }
91
write_to_kml(const volm_geo_index2_node_sptr & root,unsigned depth,std::string const & file_name)92 void volm_geo_index2::write_to_kml(const volm_geo_index2_node_sptr& root, unsigned depth, std::string const& file_name)
93 {
94 std::ofstream ofs(file_name.c_str());
95 bkml_write::open_document(ofs);
96 write_to_kml_node(ofs, root, 0, depth);
97 bkml_write::close_document(ofs);
98 }
99
depth(const volm_geo_index2_node_sptr & node)100 unsigned volm_geo_index2::depth(const volm_geo_index2_node_sptr& node)
101 {
102 if (node->children_.empty()) // alreay at leaf level
103 return 0;
104 unsigned d = 0;
105 for (auto & i : node->children_) {
106 if (!i)
107 continue;
108 unsigned dd = depth(i);
109 if (dd > d)
110 d = dd;
111 }
112 return d+1;
113 }
114
write_to_text(std::ofstream & ofs,const volm_geo_index2_node_sptr & n)115 void write_to_text(std::ofstream& ofs, const volm_geo_index2_node_sptr& n)
116 {
117 ofs << std::setprecision(6) << std::fixed << n->extent_.min_point().x() << ' '
118 << std::setprecision(6) << std::fixed << n->extent_.min_point().y() << ' '
119 << std::setprecision(6) << std::fixed << n->extent_.max_point().x() << ' '
120 << std::setprecision(6) << std::fixed << n->extent_.max_point().y() << '\n'
121 << n->children_.size() << '\n';
122 for (auto & i : n->children_) {
123 if (!i) ofs << " 0";
124 else ofs << " 1";
125 }
126 ofs << '\n';
127 for (auto & i : n->children_) {
128 if (i)
129 write_to_text(ofs, i);
130 }
131 }
132
write(const volm_geo_index2_node_sptr & root,std::string const & file_name,double const & min_size)133 void volm_geo_index2::write(const volm_geo_index2_node_sptr& root, std::string const& file_name, double const& min_size)
134 {
135 std::ofstream ofs(file_name.c_str());
136 ofs << min_size << '\n';
137 write_to_text(ofs, root);
138 }
139
get_leaves(const volm_geo_index2_node_sptr & root,std::vector<volm_geo_index2_node_sptr> & leaves)140 void volm_geo_index2::get_leaves(const volm_geo_index2_node_sptr& root, std::vector<volm_geo_index2_node_sptr>& leaves)
141 {
142 if (!root)
143 return;
144 if (!root->children_.size())
145 leaves.push_back(root);
146 else {
147 bool at_least_one_child = false;
148 for (auto & i : root->children_) {
149 if (!i)
150 continue;
151 else {
152 get_leaves(i, leaves);
153 at_least_one_child = true;
154 }
155 }
156 if (!at_least_one_child)
157 leaves.push_back(root);
158 }
159 }
160
get_leaves(const volm_geo_index2_node_sptr & root,std::vector<volm_geo_index2_node_sptr> & leaves,vgl_box_2d<double> const & area)161 void volm_geo_index2::get_leaves(const volm_geo_index2_node_sptr& root, std::vector<volm_geo_index2_node_sptr>& leaves, vgl_box_2d<double> const& area)
162 {
163 if (!root) // the node is empty
164 return;
165 if (vgl_area(vgl_intersection(root->extent_,area)) == 0.0f) // the node doesn't intersect with box
166 return;
167
168 if (!root->children_.size()) // the node has no children and it intersects with box
169 leaves.push_back(root);
170 else {
171 bool at_least_one_child = false;
172 for (auto & i : root->children_) {
173 if (!i)
174 continue;
175 else {
176 get_leaves(i, leaves, area);
177 at_least_one_child = true;
178 }
179 }
180 if (!at_least_one_child) // the node has children but all children are empty
181 leaves.push_back(root);
182 }
183 }
184
get_leaves(const volm_geo_index2_node_sptr & root,std::vector<volm_geo_index2_node_sptr> & leaves,vgl_box_2d<float> const & area)185 void volm_geo_index2::get_leaves(const volm_geo_index2_node_sptr& root, std::vector<volm_geo_index2_node_sptr>& leaves, vgl_box_2d<float> const& area)
186 {
187 vgl_box_2d<double> area_double( (double)(area.min_point().x()), (double)(area.max_point().x()),
188 (double)(area.min_point().y()), (double)(area.max_point().y()) );
189 get_leaves(root, leaves, area_double);
190 }
191
get_leaves(const volm_geo_index2_node_sptr & root,std::vector<volm_geo_index2_node_sptr> & leaves,vgl_polygon<double> const & poly)192 void volm_geo_index2::get_leaves(const volm_geo_index2_node_sptr& root, std::vector<volm_geo_index2_node_sptr>& leaves, vgl_polygon<double> const& poly)
193 {
194 if (!root) // the node is empty
195 return;
196
197 // check whether the polygon intersects with current root
198 if (!vgl_intersection(root->extent_, poly)) // the node does not intersect with given polygon
199 return;
200
201 // the node intersects with the polygon
202 if (!root->children_.size()) { // the node intersects with the polygon and has no child
203 leaves.push_back(root);
204 }
205 else { // the node has children, go inside to its children
206 bool at_least_one_child = false;
207 for (auto & i : root->children_) {
208 if (!i) // the node has children but child i is empty
209 continue;
210 else {
211 get_leaves(i, leaves, poly); // check the intersection of child i and its following children with poly
212 at_least_one_child = true;
213 }
214 }
215 if (!at_least_one_child) // the node has children but all children are empty
216 leaves.push_back(root);
217 }
218 }
219
get_leaves(const volm_geo_index2_node_sptr & root,std::vector<volm_geo_index2_node_sptr> & leaves,vgl_polygon<float> const & poly)220 void volm_geo_index2::get_leaves(const volm_geo_index2_node_sptr& root, std::vector<volm_geo_index2_node_sptr>& leaves, vgl_polygon<float> const& poly)
221 {
222 vgl_polygon<double> poly_double;
223 volm_io::convert_polygons(poly, poly_double);
224 get_leaves(root, leaves, poly_double);
225 }
226
get_leaves(const volm_geo_index2_node_sptr & root,std::vector<volm_geo_index2_node_sptr> & leaves,std::vector<vgl_point_2d<double>> const & line)227 void volm_geo_index2::get_leaves(const volm_geo_index2_node_sptr& root, std::vector<volm_geo_index2_node_sptr>& leaves, std::vector<vgl_point_2d<double> > const& line)
228 {
229 if (!root) // the tree is empty
230 return;
231 if (!is_intersect(root->extent_, line)) // the root does not intersect with the line
232 return;
233
234 // the node intersects with the line
235 if (!root->children_.size()) { // the node intersects with the line and has no child
236 leaves.push_back(root);
237 }
238 else { // the node has children and check the intersection recursively
239 bool at_least_one_child = false;
240 for (auto & i : root->children_) {
241 if (!i) // the node has children but child i is empty
242 continue;
243 else {
244 get_leaves(i, leaves, line); // check the intersection of child i and its following children
245 at_least_one_child = true;
246 }
247 }
248 if (!at_least_one_child) // the node has children but all children are empty
249 leaves.push_back(root);
250 }
251 }
252
get_leaves(const volm_geo_index2_node_sptr & root,std::vector<volm_geo_index2_node_sptr> & leaves,std::vector<vgl_point_2d<float>> const & line)253 void volm_geo_index2::get_leaves(const volm_geo_index2_node_sptr& root, std::vector<volm_geo_index2_node_sptr>& leaves, std::vector<vgl_point_2d<float> > const& line)
254 {
255 // transfer double to float since our tree bounding box is float
256 std::vector<vgl_point_2d<double> > line_double;
257 auto num_pts = (unsigned)line.size();
258 for (unsigned i = 0; i < num_pts; i++)
259 line_double.emplace_back((double)line[i].x(), (double)line[i].y());
260 get_leaves(root, leaves, line_double);
261 }
262
get_leaf(const volm_geo_index2_node_sptr & root,volm_geo_index2_node_sptr & leaf,vgl_point_2d<double> const & point)263 void volm_geo_index2::get_leaf(const volm_geo_index2_node_sptr& root, volm_geo_index2_node_sptr& leaf, vgl_point_2d<double> const& point)
264 {
265 if (!root)
266 return;
267 if (!root->extent_.contains(point))
268 return;
269
270 if (!root->children_.size()) {
271 leaf = root;
272 }
273 else {
274 bool at_least_one_child = false;
275 for (auto & i : root->children_) {
276 if (!i)
277 continue;
278 else {
279 get_leaf(i, leaf, point);
280 at_least_one_child = true;
281 }
282 }
283 if (!at_least_one_child)
284 leaf = root;
285 }
286 }
287
get_leaf(const volm_geo_index2_node_sptr & root,volm_geo_index2_node_sptr & leaf,vgl_point_2d<float> const & point)288 void volm_geo_index2::get_leaf(const volm_geo_index2_node_sptr& root, volm_geo_index2_node_sptr& leaf, vgl_point_2d<float> const& point)
289 {
290 vgl_point_2d<double> pt_double((double)point.x(), (double)point.y());
291 get_leaf(root, leaf, pt_double);
292 }
293