1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#include "s1interval.h"
4
5#include "base/logging.h"
6
7S1Interval S1Interval::FromPoint(double p) {
8 if (p == -M_PI) p = M_PI;
9 return S1Interval(p, p, ARGS_CHECKED);
10}
11
12double S1Interval::GetCenter() const {
13 double center = 0.5 * (lo() + hi());
14 if (!is_inverted()) return center;
15 // Return the center in the range (-Pi, Pi].
16 return (center <= 0) ? (center + M_PI) : (center - M_PI);
17}
18
19double S1Interval::GetLength() const {
20 double length = hi() - lo();
21 if (length >= 0) return length;
22 length += 2 * M_PI;
23 // Empty intervals have a negative length.
24 return (length > 0) ? length : -1;
25}
26
27S1Interval S1Interval::Complement() const {
28 if (lo() == hi()) return Full(); // Singleton.
29 return S1Interval(hi(), lo(), ARGS_CHECKED); // Handles empty and full.
30}
31
32double S1Interval::GetComplementCenter() const {
33 if (lo() != hi()) {
34 return Complement().GetCenter();
35 } else { // Singleton.
36 return (hi() <= 0) ? (hi() + M_PI) : (hi() - M_PI);
37 }
38}
39
40bool S1Interval::FastContains(double p) const {
41 if (is_inverted()) {
42 return (p >= lo() || p <= hi()) && !is_empty();
43 } else {
44 return p >= lo() && p <= hi();
45 }
46}
47
48bool S1Interval::Contains(double p) const {
49 // Works for empty, full, and singleton intervals.
50 DCHECK_LE(fabs(p), M_PI);
51 if (p == -M_PI) p = M_PI;
52 return FastContains(p);
53}
54
55bool S1Interval::InteriorContains(double p) const {
56 // Works for empty, full, and singleton intervals.
57 DCHECK_LE(fabs(p), M_PI);
58 if (p == -M_PI) p = M_PI;
59
60 if (is_inverted()) {
61 return p > lo() || p < hi();
62 } else {
63 return (p > lo() && p < hi()) || is_full();
64 }
65}
66
67bool S1Interval::Contains(S1Interval const& y) const {
68 // It might be helpful to compare the structure of these tests to
69 // the simpler Contains(double) method above.
70
71 if (is_inverted()) {
72 if (y.is_inverted()) return y.lo() >= lo() && y.hi() <= hi();
73 return (y.lo() >= lo() || y.hi() <= hi()) && !is_empty();
74 } else {
75 if (y.is_inverted()) return is_full() || y.is_empty();
76 return y.lo() >= lo() && y.hi() <= hi();
77 }
78}
79
80bool S1Interval::InteriorContains(S1Interval const& y) const {
81 if (is_inverted()) {
82 if (!y.is_inverted()) return y.lo() > lo() || y.hi() < hi();
83 return (y.lo() > lo() && y.hi() < hi()) || y.is_empty();
84 } else {
85 if (y.is_inverted()) return is_full() || y.is_empty();
86 return (y.lo() > lo() && y.hi() < hi()) || is_full();
87 }
88}
89
90bool S1Interval::Intersects(S1Interval const& y) const {
91 if (is_empty() || y.is_empty()) return false;
92 if (is_inverted()) {
93 // Every non-empty inverted interval contains Pi.
94 return y.is_inverted() || y.lo() <= hi() || y.hi() >= lo();
95 } else {
96 if (y.is_inverted()) return y.lo() <= hi() || y.hi() >= lo();
97 return y.lo() <= hi() && y.hi() >= lo();
98 }
99}
100
101bool S1Interval::InteriorIntersects(S1Interval const& y) const {
102 if (is_empty() || y.is_empty() || lo() == hi()) return false;
103 if (is_inverted()) {
104 return y.is_inverted() || y.lo() < hi() || y.hi() > lo();
105 } else {
106 if (y.is_inverted()) return y.lo() < hi() || y.hi() > lo();
107 return (y.lo() < hi() && y.hi() > lo()) || is_full();
108 }
109}
110
111inline static double PositiveDistance(double a, double b) {
112 // Compute the distance from "a" to "b" in the range [0, 2*Pi).
113 // This is equivalent to (drem(b - a - M_PI, 2 * M_PI) + M_PI),
114 // except that it is more numerically stable (it does not lose
115 // precision for very small positive distances).
116 double d = b - a;
117 if (d >= 0) return d;
118 // We want to ensure that if b == Pi and a == (-Pi + eps),
119 // the return result is approximately 2*Pi and not zero.
120 return (b + M_PI) - (a - M_PI);
121}
122
123double S1Interval::GetDirectedHausdorffDistance(S1Interval const& y) const {
124 if (y.Contains(*this)) return 0.0; // this includes the case *this is empty
125 if (y.is_empty()) return M_PI; // maximum possible distance on S1
126
127 double y_complement_center = y.GetComplementCenter();
128 if (Contains(y_complement_center)) {
129 return PositiveDistance(y.hi(), y_complement_center);
130 } else {
131 // The Hausdorff distance is realized by either two hi() endpoints or two
132 // lo() endpoints, whichever is farther apart.
133 double hi_hi = S1Interval(y.hi(), y_complement_center).Contains(hi()) ?
134 PositiveDistance(y.hi(), hi()) : 0;
135 double lo_lo = S1Interval(y_complement_center, y.lo()).Contains(lo()) ?
136 PositiveDistance(lo(), y.lo()) : 0;
137 DCHECK(hi_hi > 0 || lo_lo > 0);
138 return max(hi_hi, lo_lo);
139 }
140}
141
142void S1Interval::AddPoint(double p) {
143 DCHECK_LE(fabs(p), M_PI);
144 if (p == -M_PI) p = M_PI;
145
146 if (FastContains(p)) return;
147 if (is_empty()) {
148 set_hi(p);
149 set_lo(p);
150 } else {
151 // Compute distance from p to each endpoint.
152 double dlo = PositiveDistance(p, lo());
153 double dhi = PositiveDistance(hi(), p);
154 if (dlo < dhi) {
155 set_lo(p);
156 } else {
157 set_hi(p);
158 }
159 // Adding a point can never turn a non-full interval into a full one.
160 }
161}
162
163S1Interval S1Interval::FromPointPair(double p1, double p2) {
164 DCHECK_LE(fabs(p1), M_PI);
165 DCHECK_LE(fabs(p2), M_PI);
166 if (p1 == -M_PI) p1 = M_PI;
167 if (p2 == -M_PI) p2 = M_PI;
168 if (PositiveDistance(p1, p2) <= M_PI) {
169 return S1Interval(p1, p2, ARGS_CHECKED);
170 } else {
171 return S1Interval(p2, p1, ARGS_CHECKED);
172 }
173}
174
175S1Interval S1Interval::Expanded(double radius) const {
176 DCHECK_GE(radius, 0);
177 if (is_empty()) return *this;
178
179 // Check whether this interval will be full after expansion, allowing
180 // for a 1-bit rounding error when computing each endpoint.
181 if (GetLength() + 2 * radius >= 2 * M_PI - 1e-15) return Full();
182
183 S1Interval result(drem(lo() - radius, 2*M_PI), drem(hi() + radius, 2*M_PI));
184 if (result.lo() <= -M_PI) result.set_lo(M_PI);
185 return result;
186}
187
188S1Interval S1Interval::Union(S1Interval const& y) const {
189 // The y.is_full() case is handled correctly in all cases by the code
190 // below, but can follow three separate code paths depending on whether
191 // this interval is inverted, is non-inverted but contains Pi, or neither.
192
193 if (y.is_empty()) return *this;
194 if (FastContains(y.lo())) {
195 if (FastContains(y.hi())) {
196 // Either this interval contains y, or the union of the two
197 // intervals is the Full() interval.
198 if (Contains(y)) return *this; // is_full() code path
199 return Full();
200 }
201 return S1Interval(lo(), y.hi(), ARGS_CHECKED);
202 }
203 if (FastContains(y.hi())) return S1Interval(y.lo(), hi(), ARGS_CHECKED);
204
205 // This interval contains neither endpoint of y. This means that either y
206 // contains all of this interval, or the two intervals are disjoint.
207 if (is_empty() || y.FastContains(lo())) return y;
208
209 // Check which pair of endpoints are closer together.
210 double dlo = PositiveDistance(y.hi(), lo());
211 double dhi = PositiveDistance(hi(), y.lo());
212 if (dlo < dhi) {
213 return S1Interval(y.lo(), hi(), ARGS_CHECKED);
214 } else {
215 return S1Interval(lo(), y.hi(), ARGS_CHECKED);
216 }
217}
218
219S1Interval S1Interval::Intersection(S1Interval const& y) const {
220 // The y.is_full() case is handled correctly in all cases by the code
221 // below, but can follow three separate code paths depending on whether
222 // this interval is inverted, is non-inverted but contains Pi, or neither.
223
224 if (y.is_empty()) return Empty();
225 if (FastContains(y.lo())) {
226 if (FastContains(y.hi())) {
227 // Either this interval contains y, or the region of intersection
228 // consists of two disjoint subintervals. In either case, we want
229 // to return the shorter of the two original intervals.
230 if (y.GetLength() < GetLength()) return y; // is_full() code path
231 return *this;
232 }
233 return S1Interval(y.lo(), hi(), ARGS_CHECKED);
234 }
235 if (FastContains(y.hi())) return S1Interval(lo(), y.hi(), ARGS_CHECKED);
236
237 // This interval contains neither endpoint of y. This means that either y
238 // contains all of this interval, or the two intervals are disjoint.
239
240 if (y.FastContains(lo())) return *this; // is_empty() okay here
241 DCHECK(!Intersects(y));
242 return Empty();
243}
244
245bool S1Interval::ApproxEquals(S1Interval const& y, double max_error) const {
246 if (is_empty()) return y.GetLength() <= max_error;
247 if (y.is_empty()) return GetLength() <= max_error;
248 return (fabs(drem(y.lo() - lo(), 2 * M_PI)) +
249 fabs(drem(y.hi() - hi(), 2 * M_PI))) <= max_error;
250}
251