diff options
Diffstat (limited to 'src/Nerve_GIC/example/density_smoother.cpp')
-rw-r--r-- | src/Nerve_GIC/example/density_smoother.cpp | 155 |
1 files changed, 155 insertions, 0 deletions
diff --git a/src/Nerve_GIC/example/density_smoother.cpp b/src/Nerve_GIC/example/density_smoother.cpp new file mode 100644 index 00000000..4f3c4137 --- /dev/null +++ b/src/Nerve_GIC/example/density_smoother.cpp @@ -0,0 +1,155 @@ +#include <cstdlib> +#include <string.h> +#include <stdio.h> +#include <iostream> +#include <sstream> +#include <fstream> +#include <vector> +#include <algorithm> +#include <set> +#include <map> +#include <limits> +#include <assert.h> +#include <math.h> +#include <iterator> +#include <omp.h> +#include <list> +#include <random> + +using namespace std; + +double EuclideanDistance(const vector<double> & V1, const vector<double> & V2){ + assert(V1.size() == V2.size()); + double res = 0; + for(int i = 0; i < V1.size(); i++) res += pow(V1[i]-V2[i],2); + return sqrt(res); +} + +int main(int argc, char** argv){ + + char* const cloud = argv[1]; + char* const method = argv[2]; + + if(argc <= 4){cout << "Usage: <cloud_file> <DTM/GDE> <k/h> <threshold>" << endl; return 0;} + if(argc >= 6){cout << "Too many arguments !" << endl; return 0;} + + ifstream Cloud(cloud); + vector<vector<double> > C, distances; C.clear(); double d; + string line; + + int nb = 0; cout << "Reading cloud..." << endl; + while(getline(Cloud,line)){ + if(nb%100000 == 0) cout << " " << nb << endl; + vector<double> P; + stringstream stream(line); + while(stream >> d) P.push_back(d); + C.push_back(P); nb++; + } + cout << endl; + + for(int i = 0; i < nb; i++){vector<double> vect(nb); distances.push_back(vect);} + + char name[100]; sprintf(name,"%s_dist",cloud); + ifstream input(name, ios::out | ios::binary); + + if(input.good()){ + + cout << " Reading distances..." << endl; + + for(int i = 0; i < nb; i++){ + for (int j = 0; j < nb; j++){ + input.read((char*) &d,8); distances[i][j] = d; + } + } + + input.close(); + cout << " Done." << endl; + + } + + else{ + + cout << " Computing distances..." << endl; + input.close(); ofstream output(name, ios::out | ios::binary); + + for(int i = 0; i < nb-1; i++){ + if( (int) floor( 100*((double) i)/((double) nb)+1 ) %10 == 0 ){cout << "\r" << floor( 100*((double) i)/((double) nb) +1) << "%" << flush;} + for (int j = i+1; j < nb; j++){ + d = EuclideanDistance(C[i],C[j]); distances[i][j] = d; distances[j][i] = d; + output.write((char*) &d,8); + } + } + + output.close(); + cout << endl << " Done." << endl; + + } + + if(strcmp(method,"DTM") == 0){ + + int k = atoi(argv[3]); + double thresh = atof(argv[4]); + + char nameo[100]; sprintf(nameo,"%s_DTM_%d_%g",cloud,k,thresh); + ofstream output(nameo); vector<double> densities; double M = 0; + + cout << "Computing DTM..." << endl; + for(int i = 0; i < nb; i++){ + if( (int) floor( 100*((double) i)/((double) nb)+1 ) %10 == 0 ){cout << "\r" << floor( 100*((double) i)/((double) nb) +1) << "%" << flush;} + double density = 0; + sort(distances[i].begin(),distances[i].end()); + for(int j = 0; j < k; j++) + density += pow(distances[i][j],2); + density = sqrt(k)/sqrt(density); + if(density >= M) M = density; + densities.push_back(density); + + } + cout << endl << " Done." << endl; + + for(int i = 0; i < nb; i++){ + if(densities[i] >= thresh*M){ + for(int j = 0; j < C[i].size(); j++) output << C[i][j] << " "; + output << endl; + } + } + + output.close(); + + } + + if(strcmp(method,"GaussianDE") == 0){ + + double bandwidth = atof(argv[3]); + double thresh = atof(argv[4]); + + char nameo[100]; sprintf(nameo,"%s_GDE_%g_%g",cloud,bandwidth,thresh); + ofstream output(nameo); vector<double> densities; double M = 0; + + cout << "Computing GDE..." << endl; + for(int i = 0; i < nb; i++){ + if( (int) floor( 100*((double) i)/((double) nb)+1 ) %10 == 0 ){cout << "\r" << floor( 100*((double) i)/((double) nb) +1) << "%" << flush;} + double density = 0; + for(int j = 0; j < nb; j++) + density += exp(-pow(distances[i][j],2)/(2*pow(bandwidth,2))); + density /= sqrt(2*3.1415)*nb*bandwidth; + if(density >= M) M = density; + densities.push_back(density); + + } + cout << endl << " Done." << endl; + + for(int i = 0; i < nb; i++){ + if(densities[i] >= thresh*M){ + for(int j = 0; j < C[i].size(); j++) output << C[i][j] << " "; + output << endl; + } + } + + output.close(); + + } + + return 0; + +} |