ConcentricShells.h

Go to the documentation of this file.

#ifndef ConcentricShells_HEADER
#define ConcentricShells_HEADER

#ifdef HAVE_CONFIG_H
# include <config.h>
#endif

#include <boost/utility.hpp>
#include <utilities/SmallList.h>
#include <math.h> // sqrt
#include <utilities/foreach.h>
#include <utilities/bit-hacks.h>

#if 0
#define dprintf(...) printf(__VA_ARGS__)
#else
#define dprintf(...)
#endif

template <class Item, class Point>
class ConcentricShells : public boost::noncopyable {
  typedef hudson::SmallList<Item*> ItemList;

  struct Annulus {
    const int log_radius2;
    ItemList items;
    Annulus(int log_r2): log_radius2(log_r2) { }
    void add(Item *item) {
      items.push_front(item);
    }
  };

  typedef hudson::SmallList<Annulus> Annuli;

  Annuli annuli_;

  void newAnnulus(const typename Annuli::iterator& next, int log_r2, Item *item) {
    typename Annuli::iterator it = annuli_.insert(next, Annulus(log_r2));
    (*it).add(item);

    dprintf("Inserted annulus with radius %d; next has radius %d\n",
        log_r2,
        (++it) == annuli_.end()
          ? std::numeric_limits<int>::min()
          :  (*it).log_radius2
        );
  }

  public:

  /***************************************************************************
   * Add an item to the set.
   *
   * It will be stored in some annulus, as needed.
   * Takes time up to linear in the number of annuli.  For input vertices, this
   * is O(lg cfs(center)/lfs(center))).  For Steiner points, this is O(1).
   * If we cared about this term, we could make it be O(1) using a hash table
   * or an array like in BucketPQ.
   *
   * We don't store items at distance exactly zero.  TODO: instead, assert that
   * distance is non-zero, but first check for pointer equality that the center
   * is not the same as the item.
   ***************************************************************************/
  void add(const Point& center, Item *item) {

    double dist2 = item->toPoint().dist2(center);
    if (dist2 == 0.0) { return; }
    int log_d2 = bithacks::log2_norm(dist2);

    dprintf("Adding to center %s, item %s, annulus %d\n",
        center.toString().c_str(),
        item->toString().c_str(), log_d2);

    if (annuli_.empty()) {
      /* No annuli?  Better create one. */
      dprintf("Annuli empty\n");
      newAnnulus(annuli_.begin(), log_d2, item);
    }
    else if ( annuli_.front().log_radius2 < log_d2 ) {
      /* No annulus this big?  Create one. */
      dprintf("Annuli all small: %d\n", annuli_.front().log_radius2);
      newAnnulus(annuli_.begin(), log_d2, item);
    }
    else {
      /* Search for the proper annulus, or one a bit smaller. */
      typename Annuli::iterator it = annuli_.begin();
      typename Annuli::iterator end = annuli_.end();
      while (it != end && (*it).log_radius2 > log_d2) { ++it; }
      if (it != end && (*it).log_radius2 == log_d2) {
        // found!
        dprintf("Found correct annulus\n");
        (*it).add(item);
      } else {
        // not found; need a new annulus
        dprintf("No proper annulus\n");
        newAnnulus(it, log_d2, item);
      }
    }
  }

  /************************************************************
    * Reassign items *from* this set *to* the other set.
    * We look only at annuli that are at least in part closer
    * to the other center than to our own center.
   ************************************************************/
  void reassignTo(const Point& center, ConcentricShells& other, const Point& othercenter) {
    // find the midpoint between center and other.
    double pq2 = center.dist2(othercenter);
    assert(pq2 > 0.0);
    double half2 = pq2 / 4.0; // half squared, hence 1/4
    int last_annulus = bithacks::log2_norm(half2);

    /* We can't use BOOST_FOREACH since we modify the list in-place. */
    typename Annuli::iterator annulus_it = annuli_.begin();

    /* Keep going until all points in the annulus are closer
     * to our center than to the other center. */
    while (annulus_it != annuli_.end()
        && (*annulus_it).log_radius2 >= last_annulus) {
      ItemList& l = (*annulus_it).items;
      typename ItemList::iterator item_it = l.begin();
      while (item_it != l.end() /* end() must be recomputed: we modify the list */) {
        /* Reassign u if it's closer to newv than to this. */
        Item *item = *item_it;
        double up2 = center.dist2(item->toPoint());
        double uq2 = othercenter.dist2(item->toPoint());
        if (uq2 < up2) {
          item_it = l.erase(item_it);
          other.add(othercenter, item);
        } else {
          ++item_it;
        }
      }

      /* Erase the annulus if it's now empty. */
      if (l.empty()) {
        annulus_it = annuli_.erase(annulus_it);
      } else {
        ++annulus_it;
      }
    }
  }

  private:
  static int computeSmallestAnnulus(double pq2, double r2) {
    if (pq2 <= r2) {
      // B includes p, so every annulus intersects.
      return std::numeric_limits<int>::min();
    } else {
      double pB2 = pq2 + r2 - 2.0 * sqrt(pq2 * r2);
      return bithacks::log2_norm(pB2);
    }
  }

  /***************************************************************************
   * Find the closest item to q that is also within distance sqrt(r2) of q.
   *
   * Time is up to linear in the number of items, but it often much less
   * because we can avoid looking at any annuli that do not overlap B(q, r).
   ***************************************************************************/
  public:
  std::pair<Item*, double> findClosest(const Point& center, const Point& q, double r2) const {
    double pq2 = center.dist2(q);
    int smallest_annulus = computeSmallestAnnulus(pq2, r2);

    Item *best = NULL;
    // r2 stores the distance from q to best

    /* Keep going until all points in the annulus are outside B(q,r) */
    BOOST_FOREACH(const Annulus& annulus, annuli_) {
      /* Break if the annulus is too small to matter. */
      if (annulus.log_radius2 < smallest_annulus) { break; }

      /* Otherwise, iterate over all the items and keep track of the closest.
       * Update r to shrink B(q,r) as we find items. */
      BOOST_FOREACH(Item *item, annulus.items) {
        double iq2 = item->toPoint().dist2(q);
        if (iq2 < r2) {
          best = item;
          r2 = iq2;
          smallest_annulus = computeSmallestAnnulus(pq2, r2);
        }
      }
    }
    return make_pair(best, r2);
  }

  /***************************************************************************
   * Find an item that is within distance sqrt(r2) of q.
   *
   * Same as above, but short-cut.
   ***************************************************************************/
  public:
  std::pair<Item*, double> findPoint(const Point& center, const Point& q, double r2) const {
    double pq2 = center.dist2(q);
    int smallest_annulus = computeSmallestAnnulus(pq2, r2);

    // Loop over the annuli in order of size.
    BOOST_FOREACH(const Annulus& annulus, annuli_) {
      // Stop if the remaining annuli don't intersect B.
      if (annulus.log_radius2 < smallest_annulus) { break; }

      // Stop if we find an item in range.
      const ItemList& l = annulus.items;
      BOOST_FOREACH(Item *item, l) {
        double d2 = item->toPoint().dist2(q);
        if (d2 < r2) {
          return std::make_pair(item, d2);
        }
      }
    }
    return std::make_pair((Item*)0, r2);
  }

  /***************************************************************************
   * Return whether there are any items in this set.
   *
   * O(1) time.
   ***************************************************************************/
  bool empty() const {
    // in add we only create annuli when they have members; in reassign we
    // delete empty ones.  So the existence of an annulus is sufficient
    // evidence for there being an item.
    return annuli_.empty();
  }

  /***************************************************************************
   * Return an item in the outermost shell; it is within a factor of two of
   * the farthest item.
   *
   * O(1) time.
   ***************************************************************************/
  Item *getFarItem() const {
    assert(!empty());
    return annuli_.front().items.front();
  }
};

#undef dprintf

#endif

---