1 #ifndef vsph_sph_box_2d_h_
2 #define vsph_sph_box_2d_h_
3 //:
4 // \file
5 // \brief An axis-aligned box on the unit sphere.
6 // \author J.L. Mundy
7 // \date February 1, 2013
8 // \verbatim
9 //  Modifications none
10 // \endverbatim
11 
12 #include <iostream>
13 #include <vector>
14 #include <vsl/vsl_binary_io.h>
15 #ifdef _MSC_VER
16 #  include <vcl_msvc_warnings.h>
17 #endif
18 #include <vsph/vsph_sph_point_2d.h>
19 #include <vgl/vgl_vector_3d.h>
20 
21 //: An axis-aligned box on the unit sphere.
22 // \theta is elevation, \phi is azimuth
23 //
24 // Note that the creation of intervals on the azimuth circle is not
25 // well-defined. In contrast to intervals on the real line,
26 // there is an ambiguity as to which portion of the circle is
27 // the bounded set of angles. That is two points divide the circle
28 // into two arcs and the bounded set could be either one. Thus it is
29 // necessary to have three points. Points a and b define the two arcs and
30 // point c determines which arc contains the bounded set.
31 //
32 // To infer the bounded set from a sequence of points it is necessary to assume
33 // the initial three points are "compact," that is that they all lie in the
34 // smaller of the two arcs. That is |b - a|<180 and a<c<b . In this case,
35 // the predicate < means counter clockwise, i.e., a<b indicates that
36 // rotation to go from a to b is less than 180 and is counter-clockwise.
37 //
38 class vsph_sph_box_2d
39 {
40  public:
41   //: Default constructor
42   vsph_sph_box_2d();
43   //: Specify units
44   vsph_sph_box_2d(bool in_radians);
45   //: copy constructor
46   vsph_sph_box_2d(const vsph_sph_box_2d& sbox);
47   //: Constructor from three points.
48   // Three points, pa, pb and pc, are needed to define an
49   // unambiguous order on azimuth having a cut at +-180.
50   // Two of the points, pa and pb, bound the interval and the third
51   // point, c, defines which complementary arc of the circle is "in" the box.
52   // pa < pb by definition. pc lies on the  smaller or larger arc to
53   // break the ambiguity
54   vsph_sph_box_2d(vsph_sph_point_2d const& pa, vsph_sph_point_2d const& pb,
55                   vsph_sph_point_2d const& pc);
56 
57   ~vsph_sph_box_2d() = default;
58 
59   void set(double min_theta, double max_theta, double a_phi, double b_phi,
60            double c_phi, bool in_radians = true);
61 
62   vsph_sph_box_2d& operator= (const vsph_sph_box_2d & rhs);
63 
64   //: the angle units (radians or degrees) maintained by *this box
in_radians()65   bool in_radians() const {return in_radians_;}
66 
67   //: bounds on azimuth and elevation, if in_radians == false.
68   // the return is in degrees.
69   double min_phi(bool in_radians = true) const;
70   double min_theta(bool in_radians = true) const;
71   double max_phi(bool in_radians = true) const;
72   double max_theta(bool in_radians = true) const;
73 
74   double a_phi(bool in_radians = true) const;
75   double b_phi(bool in_radians = true) const;
76   double c_phi(bool in_radians = true) const;
77 
78   vsph_sph_point_2d min_point(bool in_radians = true) const;
79   vsph_sph_point_2d max_point(bool in_radians = true) const;
80 
81   //: is box empty, i.e. no points have been added
82   bool is_empty() const;
83 
84   //: add point to update box bounds.
85   // assumes the box already has the three point basis
86   // or the added point sequence is from a "compact" set on the circle
87   void add( double theta, double phi, bool in_radians = true);
add(vsph_sph_point_2d const & pt)88   void add( vsph_sph_point_2d const& pt) { add(pt.theta_, pt.phi_, pt.in_radians_); }
89   //: are two boxes the same box (same spherical domain)
90   bool operator==(const vsph_sph_box_2d& other) const;
91 
92   //: does the box have enough added points to uniquely define interval bounds
93   bool defined() const;
94 
95   //: is an azimuth angle contained in the current bounded azimuth interval
96   bool in_interval(double phi, bool in_radians=true) const;
97   //: does the box contain the specified spherical point ?
98   bool contains(double const& theta, double const& phi, bool in_radians = true) const;
99   bool contains(vsph_sph_point_2d const& p) const;
100 
101   //: does *this  box fully contain box b
102   bool contains(vsph_sph_box_2d const& b) const;
103 
104   //: area on the surface of unit sphere
105   double area() const;
106 
107   //: the spherical point corresponding to the center of the box
108   vsph_sph_point_2d center(bool in_radians = true) const;
109 
110   //: transform the box on the unit sphere
111   vsph_sph_box_2d transform(double t_theta, double t_phi, double scale,
112                             bool in_radians) const;
113   vsph_sph_box_2d transform(double t_theta,
114                             double t_phi, double scale,
115                             double theta_c,double phi_c,
116                             bool in_radians) const;
117 
118   //:subdivide box into potentially a 2x2 set of sub-boxes. Don't subdivide
119   // along a given axis if resulting angular range is less than min_ang
120   bool sub_divide(std::vector<vsph_sph_box_2d>& sub_boxes,
121                   double min_ang = 0.035) const;
122 
123   //: decompose box into approximately planar quadrilaterals
124   void planar_quads(std::vector<vgl_vector_3d<double> >& verts,
125                     std::vector<std::vector<int> >& quads,
126                     double tol = 0.01) const;
127 
128   //: display the box as a set of planar quadrilaterals in vrml
129   void display_box(std::ostream& os, float r, float g, float b,
130                    double tol = 0.01,
131                    double factor =1.0) const;
132 
133   //: display a set of boxes
134   static void display_boxes(std::string const& path,
135                             std::vector<vsph_sph_box_2d> const& boxes,
136                             std::vector<std::vector<float> > colors,
137                             double tol = 0.01,
138                             double factor =1.0);
139 
140   //: support for binary I/O
141   void print(std::ostream& os, bool in_radians = false) const;
142 
143   void b_read(vsl_b_istream& is);
144 
145   void b_write(vsl_b_ostream& os);
146 
version()147   short version() const {return 1;}
148  private:
149   //: the value of pi in the units of *this box
150   double pye() const;
151 
152   //: reduce phi to the range +-180 (assumed in *this angle units)
153   double reduce_phi(double phi) const;
154 
155   //: bounds of ccw traversal of phi interval in *this angle units
156   void phi_bounds(double& phi_start, double& phi_end) const;
157 
158   //: assign phi_a, phi_b so that b is ccw of a
159   void b_ccw_a();
160   //: assign comparisons
161   void set_comparisons();
162   //: update the current theta bounds
163   void update_theta(double th);
164   //: the azimuth angle ph is outside the current interval so extend it
165   bool extend_interval(double ph);
166   //: the angle units of *this box
167   bool in_radians_;
168   //: the three azimuth angles that define an unambigous bounded set
169   // a and b define two circular intervals and c determines the bounded interval
170   double a_phi_, b_phi_, c_phi_;
171   double min_phi_, max_phi_;
172   //: the bounds on elevation
173   double min_th_, max_th_;
174   bool a_ge_zero_, b_le_zero_, b_lt_zero_;
175   bool c_gt_a_, c_lt_b_;
176   bool c_le_pi_, c_ge_mpi_;
177   bool c_ge_zero_;
178 };
179 
180 //: return a box that represents the intersection of two boxes (could be empty)
181 // note that it is possible to have two disjoint intervals in phi (2 boxes)
182 bool intersection(vsph_sph_box_2d const& b1, vsph_sph_box_2d const& b2,
183                   std::vector<vsph_sph_box_2d>& boxes);
184 
185 //: return the area of the intersection
186 double intersection_area(vsph_sph_box_2d const& b1, vsph_sph_box_2d const& b2);
187 
188 std::ostream& operator<<(std::ostream& os, vsph_sph_box_2d const& p);
189 
190 #endif //vsph_sph_box_2d_h_
191