Skip to content

Commit

Permalink
Merge the branch for computing vertical distance into the simplificat…
Browse files Browse the repository at this point in the history
…ion branch
  • Loading branch information
ytSong74 committed Feb 27, 2024
2 parents cdbff63 + 27272f9 commit c83fc7e
Show file tree
Hide file tree
Showing 7 changed files with 312 additions and 3 deletions.
8 changes: 5 additions & 3 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,9 @@ target_link_libraries(test_edge_costs _terrain_trees_apps)
add_executable(test_smooth EXCLUDE_FROM_ALL sources/tests/test_smooth.cpp sources/utilities/utility_functions.h sources/utilities/utility_functions.cpp)
target_link_libraries(test_smooth _terrain_trees_apps)

add_executable(test_distance EXCLUDE_FROM_ALL sources/tests/test_vertical_distance.cpp sources/utilities/utility_functions.h sources/utilities/utility_functions.cpp)
target_link_libraries(test_distance _terrain_trees_apps)

# add_executable(test_morse EXCLUDE_FROM_ALL sources/tests/test_morse.cpp sources/utilities/utility_functions.h sources/utilities/utility_functions.cpp)
# target_link_libraries(test_morse _terrain_trees_apps)

Expand All @@ -80,10 +83,9 @@ target_link_libraries(test_smooth _terrain_trees_apps)
add_executable(test_morse_simplification EXCLUDE_FROM_ALL sources/tests/test_morse_simplification.cpp sources/utilities/utility_functions.h sources/utilities/utility_functions.cpp)
target_link_libraries(test_morse_simplification _terrain_trees_apps)


add_executable(test_calc_simp_stats EXCLUDE_FROM_ALL sources/tests/test_calc_node_stats_simplification.cpp sources/utilities/utility_functions.h sources/utilities/utility_functions.cpp)
target_link_libraries(test_calc_simp_stats _terrain_trees_apps)

#add_custom_target(tests DEPENDS test_index_generation test_curvature test_soup test_batched_queries test_terrain_feature_extraction test_spatial_queries test_edge_contraction test_parallel_edge_contraction test_morse test_morse_simplification test_multivariate_morse test_para_grad_edge_contraction)
add_custom_target(tests DEPENDS test_calc_simp_stats test_parallel_edge_contraction test_morse_simplification test_para_grad_edge_contraction test_smooth test_edge_costs)

add_custom_target(tests_simplification DEPENDS test_calc_simp_stats test_parallel_edge_contraction test_morse_simplification test_para_grad_edge_contraction test_smooth test_edge_costs test_distance)

Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
#include "distance_calculator.h"

Distance_Calculator::Distance_Calculator()
{
stat_max = 0;
stat_avg = 0;
stat_min = std::numeric_limits<double>::max();
count_external = 0;
}

Distance_Calculator::~Distance_Calculator()
{

}
147 changes: 147 additions & 0 deletions sources/core_library/sources/edge_contraction/distance_calculator.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,147 @@
#ifndef DISTANCE_CALCULATOR_H
#define DISTANCE_CALCULATOR_H

#pragma once

#include "terrain_trees/node_v.h"
#include "terrain_trees/node_t.h"
#include "terrain_trees/prt_tree.h"
#include "statistics/statistics.h"
#include "utilities/usage.h"
#include "utilities/cli_parameters.h"
#include "queries/topological_queries.h"
#include "simplification_aux_structures.h"
#include "utilities/timer.h"
#include <limits>

class Distance_Calculator
{
public:
Distance_Calculator();
~Distance_Calculator();
template<class T> void vertical_distance(T& tree, vector<Vertex> &points);
inline void print_stats(){
cout<<"avg: "<<stat_avg<<endl;
cout<<"max: "<<stat_max<<endl;
};
private:
template<class N> coord_type point_interpolate(N& n, Box& dom, int level, Point& p, Mesh& mesh, Spatial_Subdivision& division, bool& external);
template<class N> coord_type point_interpolate_leaf(N &n, Point &p, Mesh &mesh, bool& external);

inline coord_type interpolate_triangle(int t_id, Point& p, Mesh& mesh){
Triangle& t = mesh.get_triangle(t_id);
Vertex v1 = mesh.get_vertex(t.TV(0));
Vertex v2 = mesh.get_vertex(t.TV(1));
Vertex v3 = mesh.get_vertex(t.TV(2));
coord_type divider = (v2.get_y() - v3.get_y()) * (v1.get_x() - v3.get_x()) + (v3.get_x() - v2.get_x()) * (v1.get_y() - v3.get_y());
coord_type w1 = (v2.get_y() - v3.get_y()) * (p.get_x() - v3.get_x()) + (v3.get_x() - v2.get_x()) * (p.get_y() - v3.get_y());
w1 /= divider;
coord_type w2 = (v3.get_y() - v1.get_y()) * (p.get_x() - v3.get_x()) + (v1.get_x() - v3.get_x()) * (p.get_y() - v3.get_y());
w2 /= divider;
coord_type w3 = 1 - w1 - w2;
coord_type elevation = w1 * v1.get_z() + w2 * v2.get_z() + w3 * v3.get_z();
return elevation;
};

coord_type stat_max;
coord_type stat_avg;
coord_type stat_min;
int count_external;
coord_type elev_max;
coord_type elev_min;
};

template<class T> void Distance_Calculator::vertical_distance(T& tree, vector<Vertex> &points){
coord_type stat_sum = 0;
elev_max = points[0].get_z();
elev_min = points[0].get_z();
#pragma omp parallel for reduction(max:stat_max) reduction(+:stat_sum) reduction(+:count_external) reduction(min:stat_min) reduction(min:elev_min) reduction(max:elev_max)
for(int i = 0; i < points.size(); i++){
bool external = false;
elev_max = points[i].get_z() > elev_max ? points[i].get_z():elev_max;
elev_min = points[i].get_z() < elev_min ? points[i].get_z():elev_min;
coord_type z = this->point_interpolate(tree.get_root(),tree.get_mesh().get_domain(),0,points[i], tree.get_mesh(),tree.get_subdivision(), external);
if(external) {
count_external++;
// cout<<points[i].get_x()<<", "<<points[i].get_y()<<endl;
continue;
}
coord_type dist = abs(z - points[i].get_z());
dist = dist > Zero ? dist:0;
stat_max = dist > stat_max ? dist:stat_max;
stat_min = dist < stat_min ? dist:stat_min;
stat_sum += dist;

}

stat_avg = stat_sum/double(points.size() - count_external);
cout<<"Original mesh statistics (max, min): " <<endl;
cout<<elev_max<<" "<<elev_min<<endl;
cout<<"=============Vertical distance statistics==========="<<endl;
cout<<"avg: "<<stat_avg<<endl;
cout<<"max: "<<stat_max<<endl;
cout<<"min: "<<stat_min<<endl;
cout<<"count external: "<<count_external<<endl;
cout<<"diagonal: "<<tree.get_mesh().get_domain().get_diagonal()<<endl;
cout<<"max wrt. diagonal "<< stat_max/tree.get_mesh().get_domain().get_diagonal()<<endl;
}

template<class N> coord_type Distance_Calculator::point_interpolate(N &n, Box &dom, int level, Point& p, Mesh &mesh, Spatial_Subdivision &division, bool& external)
{
if (n.is_leaf())
{
return this->point_interpolate_leaf(n,p,mesh,external);
}
else
{
for (int i = 0; i < division.son_number(); i++)
{
Box son_dom = division.compute_domain(dom,level,i);
int son_level = level +1;
if(son_dom.contains(p,mesh.get_domain().get_max()))
{
return this->point_interpolate(*n.get_son(i),son_dom,son_level,p,mesh,division,external);
}
}
external = true;
return 0;
}
}

template<class N> coord_type Distance_Calculator::point_interpolate_leaf(N &n, Point &p, Mesh &mesh, bool& external)
{
Box bb;
pair<itype,itype> run;

for(ivect_iter it=n.get_t_array_begin(); it!=n.get_t_array_end(); ++it)
{
if(n.get_run_bounding_box(it,bb,mesh,run))
{
if(bb.contains(p,mesh.get_domain().get_max()))
{
for(itype t_id=run.first; t_id<=run.second; t_id++)
{
if(Geometry_Wrapper::point_in_triangle(t_id,p,mesh))
{
coord_type z = interpolate_triangle(t_id, p, mesh);
return z;
}

}
}
}
else
{
if(Geometry_Wrapper::point_in_triangle(*it,p,mesh))
{
coord_type z = interpolate_triangle(*it, p, mesh);
return z;
}
}
}

external = true;
return 0;
}

#endif
37 changes: 37 additions & 0 deletions sources/core_library/sources/geometry/geometry.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -215,6 +215,43 @@ int Geometry::ClipTriangle2D_strict(double minX, double minY, double maxX, doubl
**
**/
int Geometry::PointInTriangle(double xp, double yp, double * v1, double * v2, double * v3)
{
// we check if the point is a vertex of the triangle
if (xp == v1[0] && yp == v1[1]) return 1;
if (xp == v2[0] && yp == v2[1]) return 1;
if (xp == v3[0] && yp == v3[1]) return 1;
v2[0] -= v1[0]; v3[0] -= v1[0];
v2[1] -= v1[1]; v3[1] -= v1[1];
xp -= v1[0]; yp -= v1[1];
v1[0] = 0; v1[1] = 0;
// static double * v[4] = {v1, v2, v3, v4};
int d1, d2, d3, orientation;
orientation = DetSign3D(v1[0], v1[1], 1,
v2[0], v2[1], 1,
v3[0], v3[1], 1);

d1 = DetSign3D(xp, yp, 1,
v2[0], v2[1], 1,
v3[0], v3[1], 1);

if (d1 != orientation && d1 != 0) return 0;

d2 = DetSign3D(v1[0], v1[1], 1,
xp, yp, 1,
v3[0], v3[1], 1);

if (d2 != orientation && d2 != 0) return 0;

d3 = DetSign3D(v1[0], v1[1], 1,
v2[0], v2[1], 1,
xp, yp, 1);

if (d3 != orientation && d3 != 0) return 0;

return 1;
}

int Geometry::PointInTriangle_old(double xp, double yp, double * v1, double * v2, double * v3)
{
// we check if the point is a vertex of the triangle
if (xp == v1[0] && yp == v1[1]) return 1;
Expand Down
1 change: 1 addition & 0 deletions sources/core_library/sources/geometry/geometry.h
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ class Geometry
Here the triangle is given as two arrays x[3], y[3] containing the coordinates of its three vertices
*/
static int PointInTriangle (double xp, double yp, double * v1, double * v2, double * v3);
static int PointInTriangle_old (double xp, double yp, double * v1, double * v2, double * v3);
static int PointInTriangle_strict (double xp, double yp, double * v1, double * v2, double * v3);


Expand Down
107 changes: 107 additions & 0 deletions sources/tests/test_vertical_distance.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
#include "utilities/utility_functions.h"

using namespace utility_functions;

void load_tree(PRT_Tree& tree, cli_parameters &cli, bool rebuild_tree);
void interpolate_vertices(PRT_Tree& tree, cli_parameters& cli, vector<Vertex>& input_vertices, vector<coord_type>& distances);

int main(int argc, char** argv)
{
cli_parameters cli;
cli.mesh_path = argv[1]; // simplified mesh
cli.debug_mode = false;
cerr<<"[OBJECTIVE] this unit-test reads the simplified mesh and the original mesh. For each vertex in the original mesh, it calculates its vertical distance to the simplified mesh."
<<"The output is the maximum of all vertice distances. "<<endl;

cli.division_type = QUAD;
cli.crit_type = "pr";
cli.v_per_leaf = atoi(argv[2]);
PRT_Tree ptree = PRT_Tree(cli.v_per_leaf,cli.division_type);
bool rebuild = false;
if(argc >= 5 && strcmp(argv[4], "-r") == 0)
rebuild = true;

load_tree(ptree, cli, rebuild);

ptree.init_leaves_list(ptree.get_root());

string orig_mesh_path = argv[3];
stringstream output_name;
Mesh orig_mesh;
if (!Reader::read_vertices(orig_mesh, orig_mesh_path))
{
cout << "[ERROR] Loading original mesh file. Execution Stopped." << endl;
return -1;
}
vector<coord_type> distances(orig_mesh.get_vertices_num(), 0);
vector<Vertex> vertices = orig_mesh.get_vertices_array();
Timer time;

time.start();
Distance_Calculator distance_calculator;
distance_calculator.vertical_distance(ptree, vertices);
time.stop();
time.print_elapsed_time("[TIME] Calculating vertical distance ");
cerr << "[MEMORY] peak for calculating vertical distance: " << to_string(MemoryUsage().get_Virtual_Memory_in_MB()) << " MBs" << std::endl;

return (EXIT_SUCCESS);
}

void load_tree(PRT_Tree& tree, cli_parameters &cli, bool rebuild_tree)
{
Timer time;
if (!Reader::read_mesh(tree.get_mesh(), cli.mesh_path))
{
cout << "[ERROR] Loading mesh file. Execution Stopped." << endl;
return;
}

stringstream base_info;
base_info << cli.v_per_leaf << " " << cli.t_per_leaf << " " << cli.crit_type << " ";
stringstream base;
base << get_path_without_file_extension(cli.mesh_path);
stringstream tree_info;
tree_info << base_info.str() << "[TIME] Building ";

stringstream out;
out << base.str() << "_" << SpatialDecType2string(cli.division_type) << "_" << cli.crit_type;
out << "_v_" << cli.v_per_leaf << "_.tree";

cli.tree_path=out.str();

if (rebuild_tree || !Reader::read_tree(tree, tree.get_root(), cli.tree_path))
{
cerr << "[ERROR] Loading .tree file." << endl;
cerr << "[GENERATION] tree from triangle mesh" << endl;
time.start();
tree.build_tree();
time.stop();
time.print_elapsed_time(tree_info.str());
Writer::write_tree(out.str(), tree.get_root(), tree.get_subdivision());
}
else
cout << "[NOTICE] Found corresponding .tree file. Loaded tree from file successfully"<<endl;


cerr << "[MEMORY] peak for encoding the Terrain tree: " << to_string(MemoryUsage().get_Virtual_Memory_in_MB()) << " MBs" << std::endl;


stringstream out2;
out2 << base.str();
out2 << "_" << SpatialDecType2string(cli.division_type) << "_" << cli.crit_type << "_v_" << cli.v_per_leaf << "_tree.vtk";



time.start();
Reindexer reindexer = Reindexer();
reindexer.reindex_tree_and_mesh(tree,false,cli.original_vertex_indices,
false,cli.original_triangle_indices);
time.stop();
time.print_elapsed_time("[TIME] Index and Mesh Reindexing ");
cerr << "[MEMORY] peak for Index and Mesh Reindexing: " <<
to_string(MemoryUsage().get_Virtual_Memory_in_MB()) << " MBs" << std::endl;
}




1 change: 1 addition & 0 deletions sources/utilities/utility_functions.h
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@
#include "edge_contraction/contraction_simplifier.h"
#include "edge_contraction/gradient_aware_contraction.h"
#include "edge_contraction/smooth.h"
#include "edge_contraction/distance_calculator.h"

#include "morse/forman_gradient.h"
#include "morse/forman_gradient_computation.h"
Expand Down

0 comments on commit c83fc7e

Please sign in to comment.