interval.go

Documentation: github.com/golang/geo/s1

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