From 351c09724300e8024decaff7fb9bf972b22c6107 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 31 Jul 2016 01:45:21 +0200 Subject: point cloud support --- ripser.cpp | 271 +++++++++++++++++++++++++++++-------------------------------- 1 file changed, 129 insertions(+), 142 deletions(-) (limited to 'ripser.cpp') diff --git a/ripser.cpp b/ripser.cpp index 78035ba..7d4c28d 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -15,7 +15,6 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . */ - //#define ASSEMBLE_REDUCTION_MATRIX //#define USE_COEFFICIENTS @@ -30,13 +29,13 @@ //#define USE_GOOGLE_HASHMAP #include -#include -#include -#include #include -#include #include +#include +#include #include +#include +#include #include #ifdef USE_GOOGLE_HASHMAP @@ -301,26 +300,25 @@ public: init_rows(); } - value_t operator()(const index_t i, const index_t j) const; + value_t operator()(const index_t i, const index_t j) const; size_t size() const { return rows.size(); } }; class euclidean_distance_matrix { public: - std::vector> points; - - euclidean_distance_matrix(std::vector> _points) - : points(std::move(_points)) { } - - value_t operator()(const index_t i, const index_t j) const { - value_t result; - result = std::inner_product((points[i]).begin(), (points[i]).end(), (points[j]).begin(), result, - std::plus(), [] (value_t u, value_t v) { return (u-v)*(u-v); }); - return std::sqrt(result); - } - - size_t size() const { return points.size(); } + std::vector> points; + + euclidean_distance_matrix(std::vector> _points) + : points(std::move(_points)) {} + + value_t operator()(const index_t i, const index_t j) const { + return std::sqrt(std::inner_product( + points[i].begin(), points[i].end(), points[j].begin(), value_t(), std::plus(), + [](value_t u, value_t v) { return (u - v) * (u - v); })); + } + + size_t size() const { return points.size(); } }; template <> void compressed_distance_matrix::init_rows() { @@ -341,7 +339,7 @@ template <> void compressed_distance_matrix::init_rows() { template <> value_t compressed_distance_matrix::operator()(const index_t i, - const index_t j) const { + const index_t j) const { if (i < j) return rows[i][j]; else if (i > j) @@ -352,7 +350,7 @@ value_t compressed_distance_matrix::operator()(const index_t i template <> value_t compressed_distance_matrix::operator()(const index_t i, - const index_t j) const { + const index_t j) const { if (i > j) return rows[i][j]; else if (i < j) @@ -361,10 +359,8 @@ value_t compressed_distance_matrix::operator()(const index_t i return 0; } -typedef compressed_distance_matrix - compressed_lower_distance_matrix; -typedef compressed_distance_matrix - compressed_upper_distance_matrix; +typedef compressed_distance_matrix compressed_lower_distance_matrix; +typedef compressed_distance_matrix compressed_upper_distance_matrix; template diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) { if (column.empty()) @@ -460,31 +456,31 @@ void assemble_columns_to_reduce(std::vector& columns_to_reduce hash_map& pivot_column_index, const Comparator& comp, index_t dim, index_t n, value_t threshold, const binomial_coeff_table& binomial_coeff) { - index_t num_simplices = binomial_coeff(n, dim + 2); - - columns_to_reduce.clear(); - + index_t num_simplices = binomial_coeff(n, dim + 2); + + columns_to_reduce.clear(); + #ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "assembling " << num_simplices << " columns" << std::flush << "\r"; -#endif - - for (index_t index = 0; index < num_simplices; ++index) { - if (pivot_column_index.find(index) == pivot_column_index.end()) { - value_t diameter = comp.diameter(index); - if (diameter <= threshold) columns_to_reduce.push_back(std::make_pair(diameter, index)); - } - } - + std::cout << "\033[K" + << "assembling " << num_simplices << " columns" << std::flush << "\r"; +#endif + + for (index_t index = 0; index < num_simplices; ++index) { + if (pivot_column_index.find(index) == pivot_column_index.end()) { + value_t diameter = comp.diameter(index); + if (diameter <= threshold) columns_to_reduce.push_back(std::make_pair(diameter, index)); + } + } + #ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "sorting " << num_simplices << " columns" << std::flush << "\r"; + std::cout << "\033[K" + << "sorting " << num_simplices << " columns" << std::flush << "\r"; #endif - - std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), - greater_diameter_or_smaller_index()); + + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), + greater_diameter_or_smaller_index()); #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif } @@ -516,11 +512,13 @@ void compute_pairs(std::vector& columns_to_reduce, #ifdef ASSEMBLE_REDUCTION_MATRIX std::priority_queue, - smaller_index> reduction_column; + smaller_index> + reduction_column; #endif std::priority_queue, - greater_diameter_or_smaller_index> working_coboundary; + greater_diameter_or_smaller_index> + working_coboundary; value_t diameter = get_diameter(column_to_reduce); @@ -762,111 +760,100 @@ int main(int argc, char** argv) { } } - std::ifstream input_stream(filename, std::ios_base::binary | std::ios_base::in); + std::ifstream input_stream(filename); if (input_stream.fail()) { std::cerr << "couldn't open file " << filename << std::endl; exit(-1); } - + #ifdef FILE_FORMAT_POINT_CLOUD - - std::vector> points; - - std::string line; - value_t value; - while (std::getline(input_stream, line)) { - std::vector point; - std::istringstream s(line); - while (s >> value) point.push_back(value); - points.push_back(point); - } - -// euclidean_distance_matrix dist(points); -// index_t n = dist.size(); - - euclidean_distance_matrix eucl_dist(points); - - index_t n = eucl_dist.size(); - - std::vector distances; - - for (int i = 0; i < n; ++i) - for (int j = 0; j < i; ++j) - if (i > j) distances.push_back(eucl_dist(i,j)); - - compressed_lower_distance_matrix dist(distances); - - std::cout << "point cloud with " << n << " points" << std::endl; - -#endif - - + + std::vector> points; + + std::string line; + value_t value; + while (std::getline(input_stream, line)) { + std::vector point; + std::istringstream s(line); + while (s >> value) point.push_back(value); + if (!point.empty()) points.push_back(point); + } + + euclidean_distance_matrix dist(points); + index_t n = dist.size(); + + std::cout << "point cloud with " << n << " points" << std::endl; + +#endif + #ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV - - std::vector distances; - value_t value; - while (input_stream >> value) { - distances.push_back(value); - input_stream.ignore(); - } - - compressed_lower_distance_matrix dist(distances); - - index_t n = dist.size(); - - std::cout << "distance matrix with " << n << " points" << std::endl; - -#endif - + + std::vector distances; + value_t value; + while (input_stream >> value) { + distances.push_back(value); + input_stream.ignore(); + } + + compressed_lower_distance_matrix dist(distances); + + index_t n = dist.size(); + + std::cout << "distance matrix with " << n << " points" << std::endl; + +#endif + #ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV - - std::vector distances; - value_t value; - while (input_stream >> value) { - distances.push_back(value); - input_stream.ignore(); - } - - compressed_upper_distance_matrix dist(distances); - - index_t n = dist.size(); - - std::cout << "distance matrix with " << n << " points" << std::endl; - -#endif - - + + std::vector distances; + value_t value; + while (input_stream >> value) { + distances.push_back(value); + input_stream.ignore(); + } + + compressed_upper_distance_matrix dist(distances); + + index_t n = dist.size(); + + std::cout << "distance matrix with " << n << " points" << std::endl; + +#endif + #ifdef FILE_FORMAT_DIPHA - - if (read(input_stream) != 8067171840) { - std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl; - exit(-1); - } - - if (read(input_stream) != 7) { - std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl; - exit(-1); - } - - index_t n = read(input_stream); - - std::vector distances; - - for (int i = 0; i < n; ++i) - for (int j = 0; j < n; ++j) - if (i > j) distances.push_back(read(input_stream)); else read(input_stream); - - compressed_lower_distance_matrix dist(distances); - - std::cout << "distance matrix with " << n << " points" << std::endl; - -#endif - + + if (read(input_stream) != 8067171840) { + std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl; + exit(-1); + } + + if (read(input_stream) != 7) { + std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl; + exit(-1); + } + + index_t n = read(input_stream); + + std::vector distances; + + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) + if (i > j) + distances.push_back(read(input_stream)); + else + read(input_stream); + + compressed_lower_distance_matrix dist(distances); + + std::cout << "distance matrix with " << n << " points" << std::endl; + +#endif + #ifndef FILE_FORMAT_POINT_CLOUD auto result = std::minmax_element(dist.distances.begin(), dist.distances.end()); - std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl; + std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl; #endif - + dim_max = std::min(dim_max, n - 2); binomial_coeff_table binomial_coeff(n, dim_max + 2); -- cgit v1.2.3