IntervalList.h
// IntervalList.h
//
// Author David Barrett-Lennard
// (C)opyright Cedanet Pty Ltd 2007
#pragma once
#ifndef Ceda_cxUtils_IntervalList_H
#define Ceda_cxUtils_IntervalList_H
#include "HalfOpenInterval.h"
#include "CedaAssert.h"
#include "xvector.h"
#include "xostream.h"
#define INTERVAL_USE_BINARY_SEARCH
namespace ceda
{
/*
An IntervalList<T> represents an ordered list of half open intervals in coordinate type T
[----------) [---) [---------------------) [-----)
t1 t2 | |
gaps are
nonzero
It is assumed that no intervals are empty (i.e. t1 < t2), no pair of intervals
overlap or are adjacent and they are sorted in increasing order
Note that this implies a uniqueness of representation for a set in type T. Therefore
IntervalLists compare equal as vectors of (t1,t2) pairs if and only if they compare equal as
sets.
*/
template <typename T>
class IntervalList : public xvector<HalfOpenInterval<T> >
{
private:
typedef xvector<HalfOpenInterval<T> > Base;
public:
using Base::begin;
using Base::end;
using Base::size;
using Base::empty;
using Base::front;
using Base::back;
using Base::clear;
using Base::insert;
using Base::erase;
using Base::push_back;
using value_type = T;
using Interval = HalfOpenInterval<T>;
using const_iterator = typename xvector<HalfOpenInterval<T> >::const_iterator;
using iterator = typename xvector<HalfOpenInterval<T> >::iterator;
void Write(xostream& os) const
{
os << '{';
for (const_iterator i = begin() ; i != end() ; ++i)
{
os << *i;
}
os << '}';
}
void AssertValid() const
{
#ifdef CEDA_CHECK_ASSERTIONS
// Ensure that all intervals are non-empty. i.e. t1 < t2
for (const_iterator i = begin() ; i != Base::end() ; ++i)
{
cxAssert(!i->empty());
}
// Ensure that intervals are sorted correctly, and there are nonzero gaps
// between the intervals
ssize_t n = Base::size();
for (ssize_t i=0 ; i < n-1 ; ++i)
{
cxAssert((*this)[i].t2 < (*this)[i+1].t1);
}
#endif // CEDA_CHECK_ASSERTIONS
}
T Extent() const
{
T dt = 0;
for (const_iterator i = Base::begin() ; i != Base::end() ; ++i) dt += i->size();
return dt;
}
void ScalarMultiply(T f)
{
for (iterator i = Base::begin() ; i != Base::end() ; ++i)
{
i->ScalarMultiply(f);
}
}
// Returns extent of the prefix (i.e. of (*this) intersect [-infinity,t)).
T PrefixExtent(T t) const
{
T dt = 0;
const_iterator w = Base::begin();
for(;;)
{
if (w == Base::end()) return dt;
if (t < w->t2)
{
if (w->t1 < t)
{
dt += (t - w->t1);
}
return dt;
}
dt += w->size();
++w;
}
}
/*
A safe version of push_back in DEBUG mode
push_back(h) should only be called when
1) h is nonempty; and
2) h is to the right (with a nonzero gap) of the last interval if any.
*/
void safe_push_back(const Interval& h)
{
#ifdef CEDA_CHECK_ASSERTIONS
cxAssert(!h.empty());
ssize_t n = Base::size();
if (n > 0)
{
cxAssert((*this)[n-1].t2 < h.t1);
}
#endif
push_back(h);
}
/*
AddBack is similar to push_back.
push_back(h) should only be called when
1) h is nonempty; and
2) h is to the right (with a nonzero gap) of the last interval if any.
AddBack is slightly less restrictive. It automatically discards empty
intervals, and doesn't require a nonzero gap from the last interval.
If h is adjacent to the last interval then the last interval grows to
contain h, rather than adding another interval to the list.
*/
void AddBack(const Interval& h)
{
if (h.t1 < h.t2)
{
ssize_t n = Base::size();
if (n > 0)
{
cxAssert((*this)[n-1].t2 <= h.t1);
if ((*this)[n-1].t2 == h.t1)
{
(*this)[n-1].t2 = h.t2;
return;
}
}
push_back(h);
}
}
// Returns the index position of the interval of this interval set that contains t,
// or else returns -1 if no interval contains t.
// Uses binary search.
ssize_t Find_linear(T t) const
{
for (ssize_t i=0 ; i < Base::size() ; ++i)
{
if ((*this)[i].Contains(t)) return i;
}
return -1;
}
// Returns the index position of the interval of this interval set that contains t,
// or else returns -1 if no interval contains t.
// Uses binary search.
ssize_t Find(T t) const
{
if (Base::empty() || t < (*this)[0].t1) return -1;
ssize_t i2 = Base::size()-1;
if (t >= (*this)[i2].t2) return -1;
ssize_t i1 = 0;
while(i1 <= i2)
{
ssize_t i = (i1+i2)/2;
if (t < (*this)[i].t1)
{
// Recurse into left side
i2 = i-1;
}
else if (t >= (*this)[i].t2)
{
// Recurse into right side
i1 = i+1;
}
else
{
// t is inside ith interval
return i;
}
}
return -1;
}
/*
t <= t2: F F T T T T
[----) [-----) [--------------) [--) [-----) [---)
| t1 t2
t
Find the index of the left most interval satisfying t <= t2, or -1 if there is no such
interval.
*/
ssize_t FindLUB_le2(T t) const
{
if (Base::empty()) return -1;
if (t <= (*this)[0].t2) return 0;
ssize_t i2 = Base::size()-1;
if (t > (*this)[i2].t2) return -1;
cxAssert(Base::size() >= 2);
ssize_t i1 = 1;
// Loop invariant:
// 1) 1 <= i1 <= i2 < Base::size(); and
// 2) t <= (*this)[i].t2 is
// a) false for i < i1
// b) unknown for i in [i1,i2)
// c) true for i >= i2
while(i1 < i2)
{
cxAssert(1 <= i1);
cxAssert(i1 <= i2);
cxAssert(i2 < Base::size());
cxAssert(t <= (*this)[i2].t2);
cxAssert(t > (*this)[i1-1].t2);
ssize_t i = (i1+i2)/2;
if (t <= (*this)[i].t2)
{
// Recurse into left side
i2 = i;
}
else
{
// Recurse into right side
i1 = i+1;
}
}
return i1;
}
/*
t < t2: F F T T T T
[----) [-----) [--------------) [--) [-----) [---)
| t1 t2
t
Find the index of the left most interval satisfying t < t2, or -1 if there is no such
interval.
*/
ssize_t FindLUB_l2(T t) const
{
if (Base::empty()) return -1;
if (t < (*this)[0].t2) return 0;
ssize_t i2 = Base::size()-1;
if (t >= (*this)[i2].t2) return -1;
cxAssert(Base::size() >= 2);
ssize_t i1 = 1;
// Loop invariant:
// 1) 1 <= i1 <= i2 < Base::size(); and
// 2) t < (*this)[i].t2 is
// a) false for i < i1
// b) unknown for i in [i1,i2)
// c) true for i >= i2
while(i1 < i2)
{
cxAssert(1 <= i1);
cxAssert(i1 <= i2);
cxAssert(i2 < Base::size());
cxAssert(t < (*this)[i2].t2);
cxAssert(t >= (*this)[i1-1].t2);
ssize_t i = (i1+i2)/2;
if (t < (*this)[i].t2)
{
// Recurse into left side
i2 = i;
}
else
{
// Recurse into right side
i1 = i+1;
}
}
return i1;
}
/*
t < t2: F F T T T T
[----) [-----) [--------------) [--) [-----) [---)
| t1 t2
t
i1 i2-1
Find the index of the left most interval in [i1,i2) satisfying t < t2, or -1 if there is no
such interval.
*/
ssize_t FindLUB_l2(ssize_t i1,ssize_t i2, T t) const
{
cxAssert(0 <= i1 && i1 <= i2 && i2 <= Base::size());
if (i1 == i2) return -1;
if (t < (*this)[i1].t2) return i1;
--i2;
if (t >= (*this)[i2].t2) return -1;
++i1;
cxAssert(i1 <= i2);
// Loop invariant:
// 1) 1 <= i1 <= i2 < Base::size(); and
// 2) t < (*this)[i].t2 is
// a) false for i < i1
// b) unknown for i in [i1,i2)
// c) true for i >= i2
while(i1 < i2)
{
cxAssert(1 <= i1);
cxAssert(i1 <= i2);
cxAssert(i2 < Base::size());
cxAssert(t < (*this)[i2].t2);
cxAssert(t >= (*this)[i1-1].t2);
ssize_t i = (i1+i2)/2;
if (t < (*this)[i].t2)
{
// Recurse into left side
i2 = i;
}
else
{
// Recurse into right side
i1 = i+1;
}
}
return i1;
}
// Returns the index position of the interval of this interval set that contains h which
// must be nonempty, or else returns -1 if no interval contains h.
// This simple implementation uses a linear scan
ssize_t Find_linear(const Interval& h) const
{
cxAssert(!h.empty());
for (ssize_t i=0 ; i < Base::size() ; ++i)
{
if ((*this)[i].Contains(h)) return i;
}
return -1;
}
// Returns the index position of the interval of this interval set that contains h which
// must be nonempty, or else returns -1 if no interval contains h.
// This implementation uses a binary search
ssize_t Find(const Interval& h) const
{
cxAssert(!h.empty());
if (Base::empty() || h.t1 < (*this)[0].t1) return -1;
ssize_t i2 = Base::size()-1;
if (h.t2 > (*this)[i2].t2) return -1;
ssize_t i1 = 0;
while(i1 <= i2)
{
ssize_t i = (i1+i2)/2;
// There are 4 cases
// h is below and not adjacent to i
// h is above and not adjacent to i
// h is contained inside i
// h partially intersects or is adjacent to i
if (h.t2 < (*this)[i].t1)
{
// h is below and not adjacent to i
// Recurse into left side
i2 = i-1;
}
else if (h.t1 > (*this)[i].t2)
{
// h is above and not adjacent to i
// Recurse into right side
i1 = i+1;
}
else if ((*this)[i].Contains(h))
{
// h is contained inside i
return i;
}
else
{
// h partially intersects or is adjacent to i
return -1;
}
}
return -1;
}
// Returns true if this interval set contains the given value
// Linear implementation
bool Contains_linear(T t) const
{
return Find_linear(t) != -1;
}
// Like Contains_linear except uses binary search for better performance
// on large lists.
bool Contains(T t) const
{
return Find(t) != -1;
}
// Returns true iff this interval set contains h
// Simple implementation using a linear scan
bool Contains_linear(const Interval& h) const
{
if (h.empty()) return true;
return Find_linear(h) != -1;
}
// Returns true iff this interval set contains h
// Implementation uses an efficient binary search
bool Contains(const Interval& h) const
{
if (h.empty()) return true;
return Find(h) != -1;
}
bool Contains(const IntervalList<T>& rhs) const
{
// Iterate through rhs, checking that each interval is contained
// in this Interval set
//
// todo:
// This is efficient if the size of rhs is small. Otherwise it can
// be better to step through both interval sets from left to right
// in a similar manner to merge sort.
for (const_iterator i = rhs.begin() ; i != rhs.end() ; ++i)
{
if (!Contains(*i)) return false;
}
return true;
}
// Return the range of this interval set, assuming it is non-empty
// [--------) [---) [----------------) [-----)
// --> range [------------------------------------------------)
Interval GetRange() const
{
cxAssert(!Base::empty());
return Interval( front().t1, back().t2 );
}
// Returns true if this IntervalList has a nonempty intersection
// with the given interval
// This implementation uses a simple linear scan.
bool DoesIntersect_linear(const Interval& h) const
{
if (h.empty()) return false;
for (const_iterator i = begin() ; i != end() ; ++i)
{
if (h.DoesIntersect(*i)) return true;
}
// No intersection
return false;
}
// Version that uses binary search
bool DoesIntersect(const Interval& h) const
{
if (h.empty() || Base::empty() || h.t2 <= (*this)[0].t1) return false;
ssize_t i2 = Base::size()-1;
if (h.t1 >= (*this)[i2].t2) return false;
ssize_t i1 = 0;
while(i1 <= i2)
{
ssize_t i = (i1+i2)/2;
if (h.t2 <= (*this)[i].t1)
{
// h is below i
// Recurse into left side
i2 = i-1;
}
else if (h.t1 >= (*this)[i].t2)
{
// h is above i
// Recurse into right side
i1 = i+1;
}
else
{
// h intersects i
return true;
}
}
// No intersection
return false;
}
bool DoesIntersect(const IntervalList<T>& rhs) const
{
// Iterate through rhs, testing whether any interval intersects
// this Interval set
//
// todo:
// This is efficient if the size of rhs is small. Otherwise it can
// be better to step through both interval sets from left to right
// in a similar manner to merge sort.
for (const_iterator i = rhs.begin() ; i != rhs.end() ; ++i)
{
if (DoesIntersect(*i)) return true;
}
return false;
}
// Apply a positive shift to all intervals in this interval set
void operator+=(T offset)
{
for (iterator i = begin() ; i < end() ; ++i)
{
*i += offset;
}
}
// Apply a negative shift to all intervals in this interval set
void operator-=(T offset)
{
for (iterator i = begin() ; i < end() ; ++i)
{
*i -= offset;
}
}
// Take the union of this set of intervals with the given interval h
void UnionWith(const Interval& h)
{
if (!h.empty())
{
ssize_t n = Base::size();
// Find first interval that goes past beginning of h
// i.e. find index i of the left most interval satisfying h.t1 <= (*this)[i].t2
ssize_t i = FindLUB_le2(h.t1);
if (i == -1)
{
// h is completely above all the intervals - therefore append to the
// end of the list
push_back(h);
}
else if (h.t2 < (*this)[i].t1)
{
// Interval i has skipped completely past h, so therefore h needs to be
// inserted into the set
insert(begin()+i, h);
}
else
{
// h and interval i overlap or touch so we adjust interval i to include h
// If h starts before ith interval then extend ith interval to the left
if (h.t1 < (*this)[i].t1)
{
(*this)[i].t1 = h.t1;
}
if (h.t2 > (*this)[i].t2)
{
// (*this)[i] needs to be extended to the right
(*this)[i].t2 = h.t2;
// Search ahead for an interval that goes past the end of h
// i.e. find left most index j in [i+1,n) satisfying h.t2 < (*this)[j].t2
ssize_t j = FindLUB_l2(i+1,n,h.t2);
if (j == -1) j = n;
if (j < n)
{
if ((*this)[j].t1 <= h.t2)
{
// Interval j can be merged with interval i
(*this)[i].t2 = (*this)[j].t2;
++j;
}
}
// Erase all intervals in range [i+1,j)
erase(begin()+i+1, begin()+j);
}
}
}
}
// Take the difference between this interval set and the given interval h
void Subtract(const Interval& h)
{
if (!h.empty())
{
ssize_t n = Base::size();
// Find first interval that goes past the beginning of h
// i.e. find left most index i satisfying h.t1 < (*this)[i].t2
ssize_t i = FindLUB_l2(h.t1);
if (i != -1)
{
cxAssert(i < n);
if ((*this)[i].t1 < h.t1)
{
if ((*this)[i].t2 > h.t2)
{
// Interval i is broken up into two pieces by h
Interval v(h.t2, (*this)[i].t2);
(*this)[i].t2 = h.t1;
insert(begin()+i+1, v);
return;
}
// The top part of interval i needs to be chopped away
(*this)[i].t2 = h.t1;
++i;
}
if (i < n)
{
// Find first interval that goes past the end of h
// i.e. find left most index j in [i,n) satisfying h.t2 < (*this)[j].t2
ssize_t j = FindLUB_l2(i,n,h.t2);
if (j == -1)
{
j = n;
}
else
{
if ((*this)[j].t1 < h.t2)
{
// The bottom part of interval j needs to be chopped away
(*this)[j].t1 = h.t2;
}
}
// Erase all intervals in range [i,j)
erase(begin()+i, begin()+j);
}
}
}
}
void UnionWith(const IntervalList<T>& rhs)
{
if (&rhs != this)
{
/*
todo: inefficient
*/
for (const_iterator i = rhs.begin() ; i != rhs.end() ; ++i)
{
UnionWith(*i);
}
}
}
void Subtract(const IntervalList<T>& rhs)
{
if (&rhs != this)
{
/*
todo: inefficient
*/
for (const_iterator i = rhs.begin() ; i != rhs.end() ; ++i)
{
Subtract(*i);
}
}
else
{
clear();
}
}
void IntersectWith(const IntervalList<T>& rhs)
{
/*
todo: inefficient
*/
IntervalList<T> W;
Intersect(W, *this, rhs);
*this = W;
}
// Apply shift to all intervals in [w1,w2)
static void Shift(iterator w1, iterator w2, T offset)
{
while(w1 != w2)
{
*w1 += offset;
++w1;
}
}
// Insert interval [t,t+n) into this interval set, causing intervals to the right of
// position t to be shifted to the right by n.
// 'filled' indicates whether 1's or 0's are inserted in the interval [t,n)
void Insert(T t, T n, bool filled)
{
cxAssert(n >= 0);
if (n > 0)
{
iterator w = begin();
if (filled)
{
// Find the left most interval *w = [t1,t2) satisfying t <= t2
#ifdef INTERVAL_USE_BINARY_SEARCH
ssize_t i = FindLUB_le2(t);
if (i == -1)
{
/*
nonzero gap
| |
[---------) [---)
^
t
*/
safe_push_back(Interval(t,t+n));
return;
}
w += i;
#else
for(;;)
{
if (w == end())
{
/*
nonzero gap
| |
[---------) [---)
^
t
*/
safe_push_back(Interval(t,t+n));
return;
}
if (t <= w->t2) break;
// w is to left of t with a nonzero gap. Therefore can skip w
++w;
}
#endif
cxAssert(w != end());
cxAssert(t <= w->t2);
if (t < w->t1)
{
/*
Both gaps nonzero
(note that gap1 only applicable if w != begin())
|gap1| gap2 |
| | |
[---------) [---) [--w-)
^
t
*/
w = insert(w, Interval(t,t+n)); // w now points at inserted element
}
else
{
/*
w->t1 <= t <= w->t2
The insertion coalesces with w
w
[---------) [---) [------------) [----)
^
t
*/
w->t2 += n;
}
++w;
}
else
{
// Find the left most interval *w = [t1,t2) satisfying t < t2
#ifdef INTERVAL_USE_BINARY_SEARCH
ssize_t i = FindLUB_l2(t);
if (i == -1)
{
/*
0 <= gap
| |
[---------) [---)
^
t
*/
return;
}
w += i;
#else
for(;;)
{
if (w == end())
{
/*
0 <= gap
| |
[---------) [---)
^
t
*/
return;
}
if (t < w->t2) break;
// w is to left of t. Therefore can skip w
++w;
}
#endif
cxAssert(w != end());
cxAssert(t < w->t2);
if (w->t1 < t)
{
/*
w->t1 < t < w->t2
The insertion splits w into two pieces
w
[---------) [---) [------------) [----)
^
t
*/
w = insert(w, Interval(w->t1,t)); // w now points at inserted element
++w;
cxAssert(w != end());
w->t1 = t;
}
}
Shift(w,end(),n); // Shift all remaining intervals if any
}
}
/*
Erase [t,t+n) from this interval set. Intervals after t+n are shifted by -n.
W <-- W intersect [-infinity,t)
union
shift( W intersect [t+n, infinity), -n )
|<------- erase ------->|<------ shift by -n ---------->
| |
[-----) [--------) [--) [-------------) [--------) [---)
| |
t t+n
*/
void Erase(T t, T n)
{
cxAssert(0 <= n);
if (n > 0)
{
/*
Skip intervals that are to the left of t, with a nonzero gap. These intervals
cannot be affected by this erase so we can skip past them and never consider them
again.
nonzero gap
| |
[----) [--) [--------) |
|
t
(we require a nonzero gap because an erase immediately to the right of an interval
could eat up a gap causing it to coalesce with another interval)
*/
#ifdef INTERVAL_USE_BINARY_SEARCH
// Using binary search, find the index of the left most interval satisfying
// t <= t2, or -1 if there is no such interval.
ssize_t i = FindLUB_le2(t);
if (i == -1)
{
/*
gap > 0
| |
[---------) [---) |
|
t
*/
return; // The erase has no effect!
}
iterator w = begin() + i;
#else
iterator w = begin();
for(;;)
{
if (w == end())
{
/*
gap > 0
| |
[---------) [---) |
|
t
*/
return; // The erase has no effect!
}
if (t <= w->t2) break;
// w is to left of t. Therefore can skip w
++w;
}
#endif
cxAssert(w != end());
cxAssert(t <= w->t2);
/*
If t > w->t1 then we need to keep a prefix [w->t1,t) of w. This might end
up being all of w. We cannot add this to the result yet because the prefix
may coalesce with something after the erase. We choose to simply record the
size of the prefix so we can account for it later.
prefix
| |
[--------w-------)
|
t
*/
const T prefix = t - w->t1;
/*
Scan through the intervals to be erased. An interval is deleted if gap >= 0
gap >= 0
| |
[-w--) [----) [--) [----) [---w2---)
|
t+n
*/
const T t2 = t+n;
// Find w2, the left most interval in [w,end()) satisfying t+n < w2.t2
#ifdef INTERVAL_USE_BINARY_SEARCH
ssize_t i2 = FindLUB_l2(i,Base::size(), t2);
if (i2 == -1)
{
if (prefix > 0)
{
w->t2 = t; // Erase suffix of w
++w; // Keep w
}
erase(w,end());
return;
}
cxAssert(i <= i2 && i2 < Base::size());
iterator w2 = begin() + i2;
#else
iterator w2 = w;
while(w2->t2 <= t2)
{
if (++w2 == end())
{
if (prefix > 0)
{
w->t2 = t; // Erase suffix of w
++w; // Keep w
}
erase(w,end());
return;
}
}
#endif
cxAssert(w2 != end());
cxAssert(t+n < w2->t2);
if (w2->t1 <= t2)
{
/*
[-----w2------)
|
t+n
*/
cxAssert(w2->t1 <= t+n && t+n < w2->t2);
w2->t1 = t+n; // erase prefix of w2
if (prefix > 0) w2->t1 -= prefix; // Coalesce w2 with prefix of w.
}
else
{
/*
[-----w2------)
|
t+n
*/
if (prefix > 0)
{
cxAssert(w != w2);
// Prefix of w cannot coalesce with w2
w->t2 = t; // Erase suffix of w
++w; // Keep w
}
}
Shift(w2,end(),-n); // Shift all remaining intervals if any
erase(w,w2); // Erase [w,w2) which may be nothing if w == w2
}
}
};
// W = W1 intersect h
template <typename T>
void Intersect(IntervalList<T>& W, const IntervalList<T>& W1, const HalfOpenInterval<T>& h)
{
// In-place not supported
cxAssert(&W != &W1);
W.clear();
if (!h.empty() && !W1.empty())
{
// Skip intervals before h
// i.e. find left most interval *w1 = [t1,t2) satisfying h.t1 < t2
#ifdef INTERVAL_USE_BINARY_SEARCH
ssize_t i = W1.FindLUB_l2(h.t1);
if (i == -1)
{
return;
}
typename IntervalList<T>::const_iterator w1 = W1.begin() + i;
#else
typename IntervalList<T>::const_iterator w1 = W1.begin();
while(w1->t2 <= h.t1)
{
if (++w1 == W1.end()) return;
}
#endif
#if 1
// Optimised version with linear scan. Note that there is no point performing
// binary search because we need to add these to W.
if (w1->t1 < h.t2)
{
W.safe_push_back(Intersect(*w1,h));
if (++w1 == W1.end()) return;
while(w1->t1 < h.t2)
{
if (h.t2 <= w1->t2)
{
W.safe_push_back(HalfOpenInterval<T>(w1->t1, h.t2));
return;
}
else
{
W.safe_push_back(*w1);
if (++w1 == W1.end()) return;
}
}
}
#else
/*
The code below is sub optimal because
a) There is only a need to compute max(w1->t1,h.t1) initially. After that
point it can be assumed max(w1->t1,h.t1) = w1->t1.
b) As soon as it is found that min(w1->t2,h.t2) = h.t2, it is possible to return
without incrementing w1 or testing w1->t1 < h.t2 again.
*/
// Process intervals not after h
while(w1->t1 < h.t2)
{
// These intervals are not before and not after h. Therefore there is a non-empty
// intersection which can be added to the result.
W.safe_push_back(Intersect(*w1,h));
if (++w1 == W1.end()) return;
}
#endif
}
}
} // namespace ceda
#endif // include guard