diff options
author | Ulrich Bauer <mail@ulrich-bauer.org> | 2016-07-31 01:45:21 +0200 |
---|---|---|
committer | Ulrich Bauer <mail@ulrich-bauer.org> | 2016-07-31 01:45:21 +0200 |
commit | 351c09724300e8024decaff7fb9bf972b22c6107 (patch) | |
tree | 8023a373d40a1d7e6cf4feea7057b3f99bb55728 /ripser.cpp | |
parent | 945ee3a0707654cda7ad2a18fdce78254ce40dcf (diff) |
point cloud support
Diffstat (limited to 'ripser.cpp')
-rw-r--r-- | ripser.cpp | 271 |
1 files changed, 129 insertions, 142 deletions
@@ -15,7 +15,6 @@ You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. */ - //#define ASSEMBLE_REDUCTION_MATRIX //#define USE_COEFFICIENTS @@ -30,13 +29,13 @@ //#define USE_GOOGLE_HASHMAP #include <algorithm> -#include <iostream> -#include <fstream> -#include <sstream> #include <cassert> -#include <queue> #include <cmath> +#include <fstream> +#include <iostream> #include <numeric> +#include <queue> +#include <sstream> #include <unordered_map> #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<std::vector<value_t>> points; - - euclidean_distance_matrix(std::vector<std::vector<value_t>> _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>(), [] (value_t u, value_t v) { return (u-v)*(u-v); }); - return std::sqrt(result); - } - - size_t size() const { return points.size(); } + std::vector<std::vector<value_t>> points; + + euclidean_distance_matrix(std::vector<std::vector<value_t>> _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>(), + [](value_t u, value_t v) { return (u - v) * (u - v); })); + } + + size_t size() const { return points.size(); } }; template <> void compressed_distance_matrix<LOWER_TRIANGULAR>::init_rows() { @@ -341,7 +339,7 @@ template <> void compressed_distance_matrix<UPPER_TRIANGULAR>::init_rows() { template <> value_t compressed_distance_matrix<UPPER_TRIANGULAR>::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<UPPER_TRIANGULAR>::operator()(const index_t i template <> value_t compressed_distance_matrix<LOWER_TRIANGULAR>::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<LOWER_TRIANGULAR>::operator()(const index_t i return 0; } -typedef compressed_distance_matrix<LOWER_TRIANGULAR> - compressed_lower_distance_matrix; -typedef compressed_distance_matrix<UPPER_TRIANGULAR> - compressed_upper_distance_matrix; +typedef compressed_distance_matrix<LOWER_TRIANGULAR> compressed_lower_distance_matrix; +typedef compressed_distance_matrix<UPPER_TRIANGULAR> compressed_upper_distance_matrix; template <typename Heap> diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) { if (column.empty()) @@ -460,31 +456,31 @@ void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce hash_map<index_t, index_t>& 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<diameter_index_t>()); + + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), + greater_diameter_or_smaller_index<diameter_index_t>()); #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif } @@ -516,11 +512,13 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce, #ifdef ASSEMBLE_REDUCTION_MATRIX std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, - smaller_index<diameter_entry_t>> reduction_column; + smaller_index<diameter_entry_t>> + reduction_column; #endif std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, - greater_diameter_or_smaller_index<diameter_entry_t>> working_coboundary; + greater_diameter_or_smaller_index<diameter_entry_t>> + 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<std::vector<value_t>> points; - - std::string line; - value_t value; - while (std::getline(input_stream, line)) { - std::vector<value_t> 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<value_t> 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<std::vector<value_t>> points; + + std::string line; + value_t value; + while (std::getline(input_stream, line)) { + std::vector<value_t> 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<value_t> 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<value_t> 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<value_t> 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<value_t> 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<int64_t>(input_stream) != 8067171840) { - std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl; - exit(-1); - } - - if (read<int64_t>(input_stream) != 7) { - std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl; - exit(-1); - } - - index_t n = read<int64_t>(input_stream); - - std::vector<value_t> distances; - - for (int i = 0; i < n; ++i) - for (int j = 0; j < n; ++j) - if (i > j) distances.push_back(read<double>(input_stream)); else read<double>(input_stream); - - compressed_lower_distance_matrix dist(distances); - - std::cout << "distance matrix with " << n << " points" << std::endl; - -#endif - + + if (read<int64_t>(input_stream) != 8067171840) { + std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl; + exit(-1); + } + + if (read<int64_t>(input_stream) != 7) { + std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl; + exit(-1); + } + + index_t n = read<int64_t>(input_stream); + + std::vector<value_t> distances; + + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) + if (i > j) + distances.push_back(read<double>(input_stream)); + else + read<double>(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); |