---

more/num/point_map.h

Go to the documentation of this file.

//  Copyright (C) 2000--2001  Petter Urkedal (petter.urkedal@matfys.lth.se)

//  This file is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.

//  This file is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

//  As a special exception, you may use this file as part of a free
//  software library without restriction.  Specifically, if other files
//  instantiate templates or use macros or inline functions from this
//  file, or you compile this file and link it with other files to
//  produce an executable, this file does not by itself cause the
//  resulting executable to be covered by the GNU General Public
//  License.  This exception does not however invalidate any other
//  reasons why the executable file might be covered by the GNU General
//  Public License.

//  $Id: point_map.h,v 1.1 2002/05/30 18:01:40 petter_urkedal Exp $

//  Status: experimental


#include <limits>
#include <iterator>
#include <map>
#include <more/gen/iterator.h>
#include <more/gen/view.h>
#include <more/gen/utility.h>
#include <more/math/math.h>
#include <more/num/vectin.h>

#ifndef MORE_POINT_MAP_H
#define MORE_POINT_MAP_H


namespace more {
namespace num {

#ifdef CXX_HAVE_LONG_LONG
  typedef unsigned long long pindex_type;
#else
  typedef unsigned long pindex_type;
#endif

  template<typename T, int D>
    inline pindex_type
    fold(more::vectin<T, D> const& x) {
        const double scale = ldexp(.49, 8*sizeof(pindex_type)/D);
        // since std::numeric_limits<pindex_type>::max() may not work
        // for long long.
        pindex_type res = 0;
        for (int i = 0; i < D; ++i)
            res |= more::bitpow((pindex_type)((1.0+x[i])*scale), D) << i;
        return res;
    }

  template< typename Position, typename Data,
            typename Metric = euclidian_metric<Position> >
    struct point_map {
        // XXX todo: if trunc_cnt/c.size() becomes non-negligible,
        // reconstruct the container with a smaller scale.
        typedef Position position_type;
        typedef Position key_type;
        typedef Metric metric_function;
        typedef typename metric_function::result_type metric_type;
        static const int ndim = Position::CTC_size;
        typedef more::pindex_type pindex_type;
        typedef Data data_type;
        typedef std::pair<Position, Data> value_type;
#if !defined(MORE_BACKWARD) || MORE_BACKWARD >= 20000601
      private:
#endif
        typedef std::map<pindex_type, value_type> container;
        typedef typename container::iterator full_iterator;
        typedef typename container::const_iterator full_const_iterator;
        typedef typename container::value_type full_value_type;
      public:

        struct vicinity_filter
            : std::unary_function<full_value_type, bool> {
            vicinity_filter() {}
            vicinity_filter(position_type const& x_, metric_type d_,
                            metric_function const& met)
                : x(x_), d(d_), metric(met) {}
            vicinity_filter(vicinity_filter const& flt)
                : x(flt.x), d(flt.d), metric(flt.metric) {}

            bool operator()(full_value_type const& v) {
                return metric(v.second.first, x) <= d;
            }

          private:
            position_type x;
            metric_type d;
            metric_function metric;
        };

        typedef gen::transform_if_view< full_iterator, vicinity_filter,
                                        gen::select2nd<full_value_type> >
                vicinity_view;
        typedef gen::transforming_iterator< full_iterator,
                                            gen::select2nd<full_value_type> >
                iterator;
        typedef gen::transforming_iterator
                  < full_const_iterator,
                    gen::select2nd<const full_value_type> >
                const_iterator;

        point_map()
            : scale(1.0), trunc_cnt(0) {}
        point_map(metric_type radius_max)
            : scale(1.0/radius_max), trunc_cnt(0) {}

#if defined (MORE_BACKWARD) && MORE_BACKWARD < 20000601
        void insert(position_type const& x, data_type const& y) {
            position_type z = scale_and_truncate(x);
            c.insert(full_value_type(fold(z), value_type(x, y)));
        }
#endif
        void insert(value_type const& x) {
            position_type z = scale_and_truncate(x.first);
            c.insert(full_value_type(fold(z), x));
        }
#if defined (MORE_BACKWARD) && MORE_BACKWARD < 20000601
        void insert(full_iterator hint,
                    position_type const& x, data_type const& y) {
            position_type z = scale_and_truncate(x);
            c.insert(hint, full_value_type(fold(z), value_type(x, y)));
        }
#endif
        void insert(iterator hint, value_type const& x) {
            position_type z = scale_and_truncate(x);
            c.insert(hint.sub_iterator(), full_value_type(fold(z), x));
        }

        vicinity_view vicinity(position_type const& x, metric_type d);
        iterator find_closest(position_type const& x);

        const_iterator begin() const { return const_iterator(c.begin()); }
        const_iterator end() const { return const_iterator(c.end()); }
        iterator begin() { return iterator(c.begin()); }
        iterator end() { return iterator(c.end()); }

        void clear() { c.clear(); trunc_cnt = 0; scale = 1.0; }

      private:
        position_type
          scale_and_truncate(position_type x) {
              for (int i = 0; i < ndim; ++i) {
                  x[i] *= scale;
                  if (x[i] < -1.0) { ++trunc_cnt; x[i] = -1.0; }
                  else if (x[i] > 1.0) { ++trunc_cnt; x[i] = 1.0; }
              }
              return x;
          }
        void corner_indices(position_type const&, metric_type,
                            pindex_type&, pindex_type&);

        metric_type scale;
        int trunc_cnt;
        container c;
        metric_function metric;
    };


  template<typename T, typename U, typename Metric>
    inline void
    point_map<T, U, Metric>
    ::corner_indices(position_type const& x, metric_type d,
                     pindex_type& idx_inf, pindex_type& idx_sup) {
        position_type x_inf = x, x_sup = x;
        for (int i = 0; i < ndim; ++i) {
            x_inf[i] -= d;
            x_sup[i] += d;
        }
        position_type z_inf = scale*x_inf;
        position_type z_sup = scale*x_sup;
        for (int i = 0; i < ndim; ++i) {
            if (z_inf[i] < -1.0) {
                z_inf[i] = -1.0;
                if (z_sup[i] < -1.0)  z_sup[i] = -1.0;
            }
            if (z_sup[i] > 1.0) {
                z_sup[i] = 1.0;
                if (z_inf[i] > 1.0)  z_inf[i] = 1.0;
            }
        }
        idx_inf = more::max(fold(z_inf), c.begin()->first);
        full_iterator last = c.end(); --last;
        idx_sup = more::min(fold(z_sup), last->first);
    }

  template<typename T, typename U, typename Metric>
    inline typename point_map<T, U, Metric>::vicinity_view
    point_map<T, U, Metric>
    ::vicinity(position_type const& x, metric_type d) {
        pindex_type idx_inf, idx_sup;
        corner_indices(x, d, idx_inf, idx_sup);
        full_iterator it_inf = c.lower_bound(idx_inf);
        if (it_inf != c.begin()) --it_inf;
        return vicinity_view(it_inf, c.upper_bound(idx_sup),  
                             vicinity_filter(x, d, metric));
    }

  template<typename T, typename U, typename Metric>
    inline typename point_map<T, U, Metric>::iterator
    point_map<T, U, Metric>
    ::find_closest(position_type const& x) {
        if (c.size() < 2) {
            if (c.size() == 1) return iterator(c.begin());
            else throw std::logic_error
                   ("point_map::find_closest in empty map.");
        }
        pindex_type idx = fold(scale_and_truncate(x));
        full_iterator it_inc = c.lower_bound(idx);
        full_iterator it_dec = it_inc;
        full_iterator it_min = it_inc;
        metric_type d_min = metric(it_min->second.first, x);
        pindex_type idx_inf, idx_sup;
        corner_indices(x, d_min, idx_inf, idx_sup);

        bool changed;
        do {
            changed = 0;
            if (it_dec->first > idx_inf) {
                changed = true;
                --it_dec;
                metric_type d = metric(it_dec->second.first, x);
                if (d < d_min) {
                    d_min = d;
                    it_min = it_dec;
                    corner_indices(x, d_min, idx_inf, idx_sup);
                }
            }
            if (it_inc->first < idx_sup) {
                changed = true;
                ++it_inc;
                metric_type d = metric(it_inc->second.first, x);
                if (d < d_min) {
                    d_min = d;
                    it_min = it_inc;
                    corner_indices(x, d_min, idx_inf, idx_sup);
                }
            }
        }
        while (changed);

        return iterator(it_min);
    }
}}
#endif

---