summaryrefslogtreecommitdiff
path: root/GudhUI/utils/Furthest_point_epsilon_net.h
blob: 98346daa67db85be57666f4b73b98fc836a875a7 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
/* This file is part of the Gudhi Library. The Gudhi library 
 *    (Geometric Understanding in Higher Dimensions) is a generic C++ 
 *    library for computational topology.
 *
 *    Author(s):       David Salinas
 *
 *    Copyright (C) 2014  INRIA Sophia Antipolis-Mediterranee (France)
 *
 *    This program 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 3 of the License, or
 *    (at your option) any later version.
 *
 *    This program 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, see <http://www.gnu.org/licenses/>.
 * 
 */

#ifndef UTILS_FURTHEST_POINT_EPSILON_NET_H_
#define UTILS_FURTHEST_POINT_EPSILON_NET_H_

#include <vector>
#include <algorithm>  // for sort

#include "utils/UI_utils.h"

/**
 * Computes an epsilon net with furthest point strategy.
 */
template<typename SkBlComplex> class Furthest_point_epsilon_net {
 private:
  SkBlComplex& complex_;
  typedef typename SkBlComplex::Vertex_handle Vertex_handle;
  typedef typename SkBlComplex::Edge_handle Edge_handle;

  /**
   * Let V be the set of vertices.
   * Initially v0 is one arbitrarly vertex and the set V0 is {v0}.
   * Then Vk is computed as follows.
   * First we compute the vertex pk that is the furthest from Vk
   * then Vk = Vk \cup pk.
   * The radius of pk is its distance to Vk and its meeting vertex
   * is the vertex of Vk for which this distance is achieved.
   */
  struct Net_filtration_vertex {
    Vertex_handle vertex_handle;
    Vertex_handle meeting_vertex;
    double radius;

    Net_filtration_vertex(Vertex_handle vertex_handle_,
                          Vertex_handle meeting_vertex_,
                          double radius_) :
        vertex_handle(vertex_handle_), meeting_vertex(meeting_vertex_), radius(radius_) { }

    bool operator<(const Net_filtration_vertex& other) const {
      return radius < other.radius;
    }
  };

 public:
  std::vector<Net_filtration_vertex> net_filtration_;

  /**
   * @brief Modify complex to be the expansion of the k-nearest neighbor
   * symetric graph.
   */
  Furthest_point_epsilon_net(SkBlComplex& complex) :
      complex_(complex) {
    if (!complex.empty()) {
      init_filtration();
      for (std::size_t k = 2; k < net_filtration_.size(); ++k) {
        update_radius_value(k);
      }
    }
  }

  // xxx does not work if complex not full

  double radius(Vertex_handle v) {
    return net_filtration_[v.vertex].radius;
  }

 private:
  void init_filtration() {
    Vertex_handle v0 = *(complex_.vertex_range().begin());
    net_filtration_.reserve(complex_.num_vertices());
    for (auto v : complex_.vertex_range()) {
      if (v != v0)
        net_filtration_.push_back(Net_filtration_vertex(v,
                                                        Vertex_handle(-1),
                                                        squared_eucl_distance(v, v0)));
    }
    net_filtration_.push_back(Net_filtration_vertex(v0, Vertex_handle(-1), 1e10));
    auto n = net_filtration_.size();
    sort_filtration(n - 1);
  }

  void update_radius_value(int k) {
    int n = net_filtration_.size();
    int index_to_update = n - k;
    for (int i = 0; i < index_to_update; ++i) {
      net_filtration_[i].radius = (std::min)(net_filtration_[i].radius,
                                             squared_eucl_distance(net_filtration_[i].vertex_handle,
                                                                   net_filtration_[index_to_update].vertex_handle));
    }
    sort_filtration(n - k);
  }

  /**
   * sort all i first elements.
   */
  void sort_filtration(int i) {
    std::sort(net_filtration_.begin(), net_filtration_.begin() + i);
  }

  double squared_eucl_distance(Vertex_handle v1, Vertex_handle v2) const {
    return std::sqrt(Geometry_trait::Squared_distance_d()(complex_.point(v1), complex_.point(v2)));
  }

  void print_filtration() const {
    for (auto v : net_filtration_) {
      std::cerr << "v=" << v.vertex_handle << "-> d=" << v.radius << std::endl;
    }
  }
};

#endif  // UTILS_FURTHEST_POINT_EPSILON_NET_H_