interval.go (14140B)
1 // Copyright 2014 Google Inc. All rights reserved. 2 // 3 // Licensed under the Apache License, Version 2.0 (the "License"); 4 // you may not use this file except in compliance with the License. 5 // You may obtain a copy of the License at 6 // 7 // http://www.apache.org/licenses/LICENSE-2.0 8 // 9 // Unless required by applicable law or agreed to in writing, software 10 // distributed under the License is distributed on an "AS IS" BASIS, 11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 // See the License for the specific language governing permissions and 13 // limitations under the License. 14 15 package s1 16 17 import ( 18 "math" 19 "strconv" 20 ) 21 22 // An Interval represents a closed interval on a unit circle (also known 23 // as a 1-dimensional sphere). It is capable of representing the empty 24 // interval (containing no points), the full interval (containing all 25 // points), and zero-length intervals (containing a single point). 26 // 27 // Points are represented by the angle they make with the positive x-axis in 28 // the range [-π, π]. An interval is represented by its lower and upper 29 // bounds (both inclusive, since the interval is closed). The lower bound may 30 // be greater than the upper bound, in which case the interval is "inverted" 31 // (i.e. it passes through the point (-1, 0)). 32 // 33 // The point (-1, 0) has two valid representations, π and -π. The 34 // normalized representation of this point is π, so that endpoints 35 // of normal intervals are in the range (-π, π]. We normalize the latter to 36 // the former in IntervalFromEndpoints. However, we take advantage of the point 37 // -π to construct two special intervals: 38 // The full interval is [-π, π] 39 // The empty interval is [π, -π]. 40 // 41 // Treat the exported fields as read-only. 42 type Interval struct { 43 Lo, Hi float64 44 } 45 46 // IntervalFromEndpoints constructs a new interval from endpoints. 47 // Both arguments must be in the range [-π,π]. This function allows inverted intervals 48 // to be created. 49 func IntervalFromEndpoints(lo, hi float64) Interval { 50 i := Interval{lo, hi} 51 if lo == -math.Pi && hi != math.Pi { 52 i.Lo = math.Pi 53 } 54 if hi == -math.Pi && lo != math.Pi { 55 i.Hi = math.Pi 56 } 57 return i 58 } 59 60 // IntervalFromPointPair returns the minimal interval containing the two given points. 61 // Both arguments must be in [-π,π]. 62 func IntervalFromPointPair(a, b float64) Interval { 63 if a == -math.Pi { 64 a = math.Pi 65 } 66 if b == -math.Pi { 67 b = math.Pi 68 } 69 if positiveDistance(a, b) <= math.Pi { 70 return Interval{a, b} 71 } 72 return Interval{b, a} 73 } 74 75 // EmptyInterval returns an empty interval. 76 func EmptyInterval() Interval { return Interval{math.Pi, -math.Pi} } 77 78 // FullInterval returns a full interval. 79 func FullInterval() Interval { return Interval{-math.Pi, math.Pi} } 80 81 // IsValid reports whether the interval is valid. 82 func (i Interval) IsValid() bool { 83 return (math.Abs(i.Lo) <= math.Pi && math.Abs(i.Hi) <= math.Pi && 84 !(i.Lo == -math.Pi && i.Hi != math.Pi) && 85 !(i.Hi == -math.Pi && i.Lo != math.Pi)) 86 } 87 88 // IsFull reports whether the interval is full. 89 func (i Interval) IsFull() bool { return i.Lo == -math.Pi && i.Hi == math.Pi } 90 91 // IsEmpty reports whether the interval is empty. 92 func (i Interval) IsEmpty() bool { return i.Lo == math.Pi && i.Hi == -math.Pi } 93 94 // IsInverted reports whether the interval is inverted; that is, whether Lo > Hi. 95 func (i Interval) IsInverted() bool { return i.Lo > i.Hi } 96 97 // Invert returns the interval with endpoints swapped. 98 func (i Interval) Invert() Interval { 99 return Interval{i.Hi, i.Lo} 100 } 101 102 // Center returns the midpoint of the interval. 103 // It is undefined for full and empty intervals. 104 func (i Interval) Center() float64 { 105 c := 0.5 * (i.Lo + i.Hi) 106 if !i.IsInverted() { 107 return c 108 } 109 if c <= 0 { 110 return c + math.Pi 111 } 112 return c - math.Pi 113 } 114 115 // Length returns the length of the interval. 116 // The length of an empty interval is negative. 117 func (i Interval) Length() float64 { 118 l := i.Hi - i.Lo 119 if l >= 0 { 120 return l 121 } 122 l += 2 * math.Pi 123 if l > 0 { 124 return l 125 } 126 return -1 127 } 128 129 // Assumes p ∈ (-π,π]. 130 func (i Interval) fastContains(p float64) bool { 131 if i.IsInverted() { 132 return (p >= i.Lo || p <= i.Hi) && !i.IsEmpty() 133 } 134 return p >= i.Lo && p <= i.Hi 135 } 136 137 // Contains returns true iff the interval contains p. 138 // Assumes p ∈ [-π,π]. 139 func (i Interval) Contains(p float64) bool { 140 if p == -math.Pi { 141 p = math.Pi 142 } 143 return i.fastContains(p) 144 } 145 146 // ContainsInterval returns true iff the interval contains oi. 147 func (i Interval) ContainsInterval(oi Interval) bool { 148 if i.IsInverted() { 149 if oi.IsInverted() { 150 return oi.Lo >= i.Lo && oi.Hi <= i.Hi 151 } 152 return (oi.Lo >= i.Lo || oi.Hi <= i.Hi) && !i.IsEmpty() 153 } 154 if oi.IsInverted() { 155 return i.IsFull() || oi.IsEmpty() 156 } 157 return oi.Lo >= i.Lo && oi.Hi <= i.Hi 158 } 159 160 // InteriorContains returns true iff the interior of the interval contains p. 161 // Assumes p ∈ [-π,π]. 162 func (i Interval) InteriorContains(p float64) bool { 163 if p == -math.Pi { 164 p = math.Pi 165 } 166 if i.IsInverted() { 167 return p > i.Lo || p < i.Hi 168 } 169 return (p > i.Lo && p < i.Hi) || i.IsFull() 170 } 171 172 // InteriorContainsInterval returns true iff the interior of the interval contains oi. 173 func (i Interval) InteriorContainsInterval(oi Interval) bool { 174 if i.IsInverted() { 175 if oi.IsInverted() { 176 return (oi.Lo > i.Lo && oi.Hi < i.Hi) || oi.IsEmpty() 177 } 178 return oi.Lo > i.Lo || oi.Hi < i.Hi 179 } 180 if oi.IsInverted() { 181 return i.IsFull() || oi.IsEmpty() 182 } 183 return (oi.Lo > i.Lo && oi.Hi < i.Hi) || i.IsFull() 184 } 185 186 // Intersects returns true iff the interval contains any points in common with oi. 187 func (i Interval) Intersects(oi Interval) bool { 188 if i.IsEmpty() || oi.IsEmpty() { 189 return false 190 } 191 if i.IsInverted() { 192 return oi.IsInverted() || oi.Lo <= i.Hi || oi.Hi >= i.Lo 193 } 194 if oi.IsInverted() { 195 return oi.Lo <= i.Hi || oi.Hi >= i.Lo 196 } 197 return oi.Lo <= i.Hi && oi.Hi >= i.Lo 198 } 199 200 // InteriorIntersects returns true iff the interior of the interval contains any points in common with oi, including the latter's boundary. 201 func (i Interval) InteriorIntersects(oi Interval) bool { 202 if i.IsEmpty() || oi.IsEmpty() || i.Lo == i.Hi { 203 return false 204 } 205 if i.IsInverted() { 206 return oi.IsInverted() || oi.Lo < i.Hi || oi.Hi > i.Lo 207 } 208 if oi.IsInverted() { 209 return oi.Lo < i.Hi || oi.Hi > i.Lo 210 } 211 return (oi.Lo < i.Hi && oi.Hi > i.Lo) || i.IsFull() 212 } 213 214 // Compute distance from a to b in [0,2π], in a numerically stable way. 215 func positiveDistance(a, b float64) float64 { 216 d := b - a 217 if d >= 0 { 218 return d 219 } 220 return (b + math.Pi) - (a - math.Pi) 221 } 222 223 // Union returns the smallest interval that contains both the interval and oi. 224 func (i Interval) Union(oi Interval) Interval { 225 if oi.IsEmpty() { 226 return i 227 } 228 if i.fastContains(oi.Lo) { 229 if i.fastContains(oi.Hi) { 230 // Either oi ⊂ i, or i ∪ oi is the full interval. 231 if i.ContainsInterval(oi) { 232 return i 233 } 234 return FullInterval() 235 } 236 return Interval{i.Lo, oi.Hi} 237 } 238 if i.fastContains(oi.Hi) { 239 return Interval{oi.Lo, i.Hi} 240 } 241 242 // Neither endpoint of oi is in i. Either i ⊂ oi, or i and oi are disjoint. 243 if i.IsEmpty() || oi.fastContains(i.Lo) { 244 return oi 245 } 246 247 // This is the only hard case where we need to find the closest pair of endpoints. 248 if positiveDistance(oi.Hi, i.Lo) < positiveDistance(i.Hi, oi.Lo) { 249 return Interval{oi.Lo, i.Hi} 250 } 251 return Interval{i.Lo, oi.Hi} 252 } 253 254 // Intersection returns the smallest interval that contains the intersection of the interval and oi. 255 func (i Interval) Intersection(oi Interval) Interval { 256 if oi.IsEmpty() { 257 return EmptyInterval() 258 } 259 if i.fastContains(oi.Lo) { 260 if i.fastContains(oi.Hi) { 261 // Either oi ⊂ i, or i and oi intersect twice. Neither are empty. 262 // In the first case we want to return i (which is shorter than oi). 263 // In the second case one of them is inverted, and the smallest interval 264 // that covers the two disjoint pieces is the shorter of i and oi. 265 // We thus want to pick the shorter of i and oi in both cases. 266 if oi.Length() < i.Length() { 267 return oi 268 } 269 return i 270 } 271 return Interval{oi.Lo, i.Hi} 272 } 273 if i.fastContains(oi.Hi) { 274 return Interval{i.Lo, oi.Hi} 275 } 276 277 // Neither endpoint of oi is in i. Either i ⊂ oi, or i and oi are disjoint. 278 if oi.fastContains(i.Lo) { 279 return i 280 } 281 return EmptyInterval() 282 } 283 284 // AddPoint returns the interval expanded by the minimum amount necessary such 285 // that it contains the given point "p" (an angle in the range [-π, π]). 286 func (i Interval) AddPoint(p float64) Interval { 287 if math.Abs(p) > math.Pi { 288 return i 289 } 290 if p == -math.Pi { 291 p = math.Pi 292 } 293 if i.fastContains(p) { 294 return i 295 } 296 if i.IsEmpty() { 297 return Interval{p, p} 298 } 299 if positiveDistance(p, i.Lo) < positiveDistance(i.Hi, p) { 300 return Interval{p, i.Hi} 301 } 302 return Interval{i.Lo, p} 303 } 304 305 // Define the maximum rounding error for arithmetic operations. Depending on the 306 // platform the mantissa precision may be different than others, so we choose to 307 // use specific values to be consistent across all. 308 // The values come from the C++ implementation. 309 var ( 310 // epsilon is a small number that represents a reasonable level of noise between two 311 // values that can be considered to be equal. 312 epsilon = 1e-15 313 // dblEpsilon is a smaller number for values that require more precision. 314 dblEpsilon = 2.220446049e-16 315 ) 316 317 // Expanded returns an interval that has been expanded on each side by margin. 318 // If margin is negative, then the function shrinks the interval on 319 // each side by margin instead. The resulting interval may be empty or 320 // full. Any expansion (positive or negative) of a full interval remains 321 // full, and any expansion of an empty interval remains empty. 322 func (i Interval) Expanded(margin float64) Interval { 323 if margin >= 0 { 324 if i.IsEmpty() { 325 return i 326 } 327 // Check whether this interval will be full after expansion, allowing 328 // for a rounding error when computing each endpoint. 329 if i.Length()+2*margin+2*dblEpsilon >= 2*math.Pi { 330 return FullInterval() 331 } 332 } else { 333 if i.IsFull() { 334 return i 335 } 336 // Check whether this interval will be empty after expansion, allowing 337 // for a rounding error when computing each endpoint. 338 if i.Length()+2*margin-2*dblEpsilon <= 0 { 339 return EmptyInterval() 340 } 341 } 342 result := IntervalFromEndpoints( 343 math.Remainder(i.Lo-margin, 2*math.Pi), 344 math.Remainder(i.Hi+margin, 2*math.Pi), 345 ) 346 if result.Lo <= -math.Pi { 347 result.Lo = math.Pi 348 } 349 return result 350 } 351 352 // ApproxEqual reports whether this interval can be transformed into the given 353 // interval by moving each endpoint by at most ε, without the 354 // endpoints crossing (which would invert the interval). Empty and full 355 // intervals are considered to start at an arbitrary point on the unit circle, 356 // so any interval with (length <= 2*ε) matches the empty interval, and 357 // any interval with (length >= 2*π - 2*ε) matches the full interval. 358 func (i Interval) ApproxEqual(other Interval) bool { 359 // Full and empty intervals require special cases because the endpoints 360 // are considered to be positioned arbitrarily. 361 if i.IsEmpty() { 362 return other.Length() <= 2*epsilon 363 } 364 if other.IsEmpty() { 365 return i.Length() <= 2*epsilon 366 } 367 if i.IsFull() { 368 return other.Length() >= 2*(math.Pi-epsilon) 369 } 370 if other.IsFull() { 371 return i.Length() >= 2*(math.Pi-epsilon) 372 } 373 374 // The purpose of the last test below is to verify that moving the endpoints 375 // does not invert the interval, e.g. [-1e20, 1e20] vs. [1e20, -1e20]. 376 return (math.Abs(math.Remainder(other.Lo-i.Lo, 2*math.Pi)) <= epsilon && 377 math.Abs(math.Remainder(other.Hi-i.Hi, 2*math.Pi)) <= epsilon && 378 math.Abs(i.Length()-other.Length()) <= 2*epsilon) 379 380 } 381 382 func (i Interval) String() string { 383 // like "[%.7f, %.7f]" 384 return "[" + strconv.FormatFloat(i.Lo, 'f', 7, 64) + ", " + strconv.FormatFloat(i.Hi, 'f', 7, 64) + "]" 385 } 386 387 // Complement returns the complement of the interior of the interval. An interval and 388 // its complement have the same boundary but do not share any interior 389 // values. The complement operator is not a bijection, since the complement 390 // of a singleton interval (containing a single value) is the same as the 391 // complement of an empty interval. 392 func (i Interval) Complement() Interval { 393 if i.Lo == i.Hi { 394 // Singleton. The interval just contains a single point. 395 return FullInterval() 396 } 397 // Handles empty and full. 398 return Interval{i.Hi, i.Lo} 399 } 400 401 // ComplementCenter returns the midpoint of the complement of the interval. For full and empty 402 // intervals, the result is arbitrary. For a singleton interval (containing a 403 // single point), the result is its antipodal point on S1. 404 func (i Interval) ComplementCenter() float64 { 405 if i.Lo != i.Hi { 406 return i.Complement().Center() 407 } 408 // Singleton. The interval just contains a single point. 409 if i.Hi <= 0 { 410 return i.Hi + math.Pi 411 } 412 return i.Hi - math.Pi 413 } 414 415 // DirectedHausdorffDistance returns the Hausdorff distance to the given interval. 416 // For two intervals i and y, this distance is defined by 417 // h(i, y) = max_{p in i} min_{q in y} d(p, q), 418 // where d(.,.) is measured along S1. 419 func (i Interval) DirectedHausdorffDistance(y Interval) Angle { 420 if y.ContainsInterval(i) { 421 return 0 // This includes the case i is empty. 422 } 423 if y.IsEmpty() { 424 return Angle(math.Pi) // maximum possible distance on s1. 425 } 426 yComplementCenter := y.ComplementCenter() 427 if i.Contains(yComplementCenter) { 428 return Angle(positiveDistance(y.Hi, yComplementCenter)) 429 } 430 431 // The Hausdorff distance is realized by either two i.Hi endpoints or two 432 // i.Lo endpoints, whichever is farther apart. 433 hiHi := 0.0 434 if IntervalFromEndpoints(y.Hi, yComplementCenter).Contains(i.Hi) { 435 hiHi = positiveDistance(y.Hi, i.Hi) 436 } 437 438 loLo := 0.0 439 if IntervalFromEndpoints(yComplementCenter, y.Lo).Contains(i.Lo) { 440 loLo = positiveDistance(i.Lo, y.Lo) 441 } 442 443 return Angle(math.Max(hiHi, loLo)) 444 } 445 446 // Project returns the closest point in the interval to the given point p. 447 // The interval must be non-empty. 448 func (i Interval) Project(p float64) float64 { 449 if p == -math.Pi { 450 p = math.Pi 451 } 452 if i.fastContains(p) { 453 return p 454 } 455 // Compute distance from p to each endpoint. 456 dlo := positiveDistance(p, i.Lo) 457 dhi := positiveDistance(i.Hi, p) 458 if dlo < dhi { 459 return i.Lo 460 } 461 return i.Hi 462 }