Skip to content

Commit

Permalink
Merge branch 'dev' into alpha
Browse files Browse the repository at this point in the history
  • Loading branch information
@larc
larc committed Mar 21, 2024
2 parents 5ab234b + 93c3c8d commit 0f367f4
Show file tree
Hide file tree
Showing 35 changed files with 899 additions and 929 deletions.
40 changes: 21 additions & 19 deletions include/gproshan/geodesics/geodesics_ptp.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ namespace gproshan {
#ifdef __CUDACC__

__global__
void relax_ptp(const CHE * mesh, real_t * new_dist, real_t * old_dist, index_t * new_clusters, index_t * old_clusters, const index_t start, const index_t end, const index_t * sorted = nullptr);
void relax_ptp(const che * mesh, real_t * new_dist, real_t * old_dist, index_t * new_clusters, index_t * old_clusters, const index_t start, const index_t end, const index_t * sorted = nullptr);

__global__
void relative_error(real_t * error, const real_t * new_dist, const real_t * old_dist, const index_t start, const index_t end, const index_t * sorted = nullptr);
Expand Down Expand Up @@ -97,30 +97,30 @@ template<class T>
__forceinline__
#endif
__host_device__
real_t update_step(const CHE * mesh, const T * dist, const uvec3 & x)
real_t update_step(const che * mesh, const T * dist, const uvec3 & x)
{
const vec<T, 3> X[2] = {mesh->GT[x[0]] - mesh->GT[x[2]],
mesh->GT[x[1]] - mesh->GT[x[2]]
const vec<T, 3> X[2] = {mesh->point(x[0]) - mesh->point(x[2]),
mesh->point(x[1]) - mesh->point(x[2])
};

const vec<T, 2> & t = {dist[x[0]], dist[x[1]]};
const vec<T, 2> t = {dist[x[0]], dist[x[1]]};

mat<T, 2> q;
q[0][0] = dot(X[0], X[0]);
q[0][1] = dot(X[0], X[1]);
q[1][0] = dot(X[1], X[0]);
q[1][1] = dot(X[1], X[1]);

const T & det = q[0][0] * q[1][1] - q[0][1] * q[1][0];
const T det = q[0][0] * q[1][1] - q[0][1] * q[1][0];

mat<T, 2> Q;
Q[0][0] = q[1][1] / det;
Q[0][1] = -q[0][1] / det;
Q[1][0] = -q[1][0] / det;
Q[1][1] = q[0][0] / det;

const T & delta = t[0] * (Q[0][0] + Q[1][0]) + t[1] * (Q[0][1] + Q[1][1]);
const T & dis = delta * delta -
const T delta = t[0] * (Q[0][0] + Q[1][0]) + t[1] * (Q[0][1] + Q[1][1]);
const T dis = delta * delta -
(Q[0][0] + Q[0][1] + Q[1][0] + Q[1][1]) *
(t[0] * t[0] * Q[0][0] + t[0] * t[1] * (Q[1][0] + Q[0][1]) + t[1] * t[1] * Q[1][1] - 1);

Expand All @@ -130,13 +130,13 @@ real_t update_step(const CHE * mesh, const T * dist, const uvec3 & x)
T p = (delta + sqrt(dis)) / (Q[0][0] + Q[0][1] + Q[1][0] + Q[1][1]);
#endif

const vec<T, 2> & tp = t - p;
const vec<T, 3> & n = { tp[0] * (X[0][0]*Q[0][0] + X[1][0]*Q[1][0]) + tp[1] * (X[0][0]*Q[0][1] + X[1][0]*Q[1][1]),
const vec<T, 2> tp = t - p;
const vec<T, 3> n = { tp[0] * (X[0][0]*Q[0][0] + X[1][0]*Q[1][0]) + tp[1] * (X[0][0]*Q[0][1] + X[1][0]*Q[1][1]),
tp[0] * (X[0][1]*Q[0][0] + X[1][1]*Q[1][0]) + tp[1] * (X[0][1]*Q[0][1] + X[1][1]*Q[1][1]),
tp[0] * (X[0][2]*Q[0][0] + X[1][2]*Q[1][0]) + tp[1] * (X[0][2]*Q[0][1] + X[1][2]*Q[1][1])
};

const vec<T, 2> & c = Q * vec<T, 2>{dot(X[0], n), dot(X[1], n)};
const vec<T, 2> c = Q * vec<T, 2>{dot(X[0], n), dot(X[1], n)};

if(t[0] == INFINITY || t[1] == INFINITY || dis < 0 || c[0] >= 0 || c[1] >= 0)
{
Expand All @@ -153,35 +153,37 @@ template<class T>
__forceinline__
#endif
__host_device__
void relax_ptp(const CHE * mesh, T * new_dist, T * old_dist, index_t * new_clusters, index_t * old_clusters, const index_t v)
void relax_ptp(const che * mesh, T * new_dist, T * old_dist, index_t * new_clusters, index_t * old_clusters, const index_t v)
{
real_t & ndv = new_dist[v] = old_dist[v];
if(new_clusters) new_clusters[v] = old_clusters[v];

real_t d;
for_star(he, mesh, v)
for(const index_t he: mesh->star(v))
{
const uvec3 i = {mesh->VT[he_next(he)], mesh->VT[he_prev(he)], mesh->VT[he]};
const uvec3 i = { mesh->halfedge(he_next(he)),
mesh->halfedge(he_prev(he)),
mesh->halfedge(he)
};

d = update_step(mesh, old_dist, i);
real_t d = update_step(mesh, old_dist, i);

if(d < ndv)
{
ndv = d;
if(new_clusters)
new_clusters[v] = old_dist[i.y()] < old_dist[i.x()] ? old_clusters[i.y()] : old_clusters[i.x()];
new_clusters[v] = old_clusters[old_dist[i.y()] < old_dist[i.x()] ? i.y() : i.x()];
}
}
}


template<class T>
#ifdef __CUDACC__
index_t run_ptp(const CHE * mesh, const std::vector<index_t> & sources,
index_t run_ptp(const che * mesh, const std::vector<index_t> & sources,
const std::vector<index_t> & limits, T * error, T ** dist, index_t ** clusters,
const index_t * idx, index_t * sorted, const f_ptp<T> & fun = nullptr)
#else
index_t run_ptp(const CHE * mesh, const std::vector<index_t> & sources,
index_t run_ptp(const che * mesh, const std::vector<index_t> & sources,
const std::vector<index_t> & limits, T ** dist, index_t ** clusters,
const index_t * idx, index_t * sorted, const f_ptp<T> & fun = nullptr)
#endif
Expand Down
14 changes: 7 additions & 7 deletions include/gproshan/geodesics/heat_method.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@
#define HEAT_METHOD_H

#include <gproshan/mesh/che.h>
#include <gproshan/include_arma.h>

#include <armadillo>
#include <cholmod.h>


Expand All @@ -29,24 +29,24 @@ enum heat_method_opt {

double heat_method(real_t * dist, const che * mesh, const std::vector<index_t> & sources, const heat_method_opt & opt);

void compute_divergence(const che * mesh, const a_mat & u, a_mat & div);
arma::vec compute_divergence(const che * mesh, const arma::vec & u);

/// cholmod Keenan implementation
/// base on the code https://github.com/larc/dgpdec-course/tree/master/Geodesics
double solve_positive_definite(a_mat & x, const a_sp_mat & A, const a_mat & b, cholmod_common * context);
double solve_positive_definite(arma::mat & x, const arma::sp_mat & A, const arma::mat & b, cholmod_common * context);

cholmod_dense * arma_2_cholmod(const a_mat & m, cholmod_common * context);
cholmod_dense * arma_2_cholmod(const arma::mat & m, cholmod_common * context);

cholmod_sparse * arma_2_cholmod(const a_sp_mat & m, cholmod_common * context);
cholmod_sparse * arma_2_cholmod(const arma::sp_mat & m, cholmod_common * context);

#ifdef GPROSHAN_CUDA

///
double solve_positive_definite_gpu(a_mat & x, const a_sp_mat & A, const a_mat & b);
double solve_positive_definite_gpu(arma::mat & x, const arma::sp_mat & A, const arma::mat & b);

/// host and device support
/// https://docs.nvidia.com/cuda/cusolver/index.html#cusolver-lt-t-gt-csrlsvchol
double solve_positive_definite_cusolver(const int m, const int nnz, const real_t * hA_values, const int * hA_col_ptrs, const int * hA_row_indices, const real_t * hb, real_t * hx, const bool host = 0);
double solve_positive_definite_cusolver(const int m, const int nnz, const double * hA_values, const int * hA_col_ptrs, const int * hA_row_indices, const double * hb, double * hx, const bool host = 0);

#endif // GPROSHAN_CUDA

Expand Down
2 changes: 1 addition & 1 deletion include/gproshan/geometry/mat.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ class mat
}

__host_device__
const T & operator () (const index_t i, const index_t j) const
T operator () (const index_t i, const index_t j) const
{
return rows[i][j];
}
Expand Down
44 changes: 22 additions & 22 deletions include/gproshan/geometry/vec.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ class vec

public:
__host_device__
vec(const T & val = 0)
vec(const T val = 0)
{
for(T & v: values)
v = val;
Expand All @@ -34,12 +34,12 @@ class vec
vec(const std::initializer_list<T> & list)
{
int i = -1;
for(const T & v: list)
for(T v: list)
values[++i] = v;
}

__host_device__
vec(const vec<T, N - 1> & v, const T & val = 0)
vec(const vec<T, N - 1> & v, const T val = 0)
{
for(index_t i = 0; i < N - 1; ++i)
values[i] = v[i];
Expand All @@ -61,15 +61,15 @@ class vec
}

__host_device__
const T & operator [] (const index_t i) const
T operator [] (const index_t i) const
{
assert(i < N);
return values[i];
}

///< concatenate with comma operator
__host_device__
vec<T, N + 1> operator , (const T & a) const
vec<T, N + 1> operator , (const T a) const
{
return {*this, a};
}
Expand All @@ -82,7 +82,7 @@ class vec
}

__host_device__
const T & x() const
T x() const
{
assert(N > 0);
return values[0];
Expand All @@ -96,7 +96,7 @@ class vec
}

__host_device__
const T & y() const
T y() const
{
assert(N > 1);
return values[1];
Expand All @@ -110,7 +110,7 @@ class vec
}

__host_device__
const T & z() const
T z() const
{
assert(N > 2);
return values[2];
Expand All @@ -122,7 +122,7 @@ class vec
T norm() const
{
T res = 0;
for(const T & v: values)
for(T v: values)
res += v * v;
return std::sqrt(res);
}
Expand Down Expand Up @@ -258,15 +258,15 @@ class vec
///< scalar product
template<class T, size_t N>
__host_device__
vec<T, N> operator * (const T & a, const vec<T, N> & v)
vec<T, N> operator * (const T a, const vec<T, N> & v)
{
return v * a;
}

///< scalar product
template<class T, size_t N>
__host_device__
vec<T, N> operator / (const T & a, const vec<T, N> & v)
vec<T, N> operator / (const T a, const vec<T, N> & v)
{
vec<T, N> res;
for(index_t i = 0; i < N; ++i)
Expand All @@ -287,14 +287,10 @@ template<class T>
__host_device__
vec<T, 3> cross(const vec<T, 3> & u, const vec<T, 3> & v)
{
const T & ux = u[0];
const T & uy = u[1];
const T & uz = u[2];
const T & vx = v[0];
const T & vy = v[1];
const T & vz = v[2];

return {uy * vz - uz * vy, uz * vx - ux * vz, ux * vy - uy * vx};
return {u[1] * v[2] - u[2] * v[1],
u[2] * v[0] - u[0] * v[2],
u[0] * v[1] - u[1] * v[0]
};
}

///< dot product
Expand Down Expand Up @@ -350,9 +346,13 @@ std::istream & operator >> (std::istream & is, vec<T, N> & v)
}


using vec2 = vec<real_t, 2>;
using vec3 = vec<real_t, 3>;
using vec4 = vec<real_t, 4>;
using vec2 = vec<float, 2>;
using vec3 = vec<float, 3>;
using vec4 = vec<float, 4>;

using dvec2 = vec<double, 2>;
using dvec3 = vec<double, 3>;
using dvec4 = vec<double, 4>;

using uvec2 = vec<unsigned int, 2>;
using uvec3 = vec<unsigned int, 3>;
Expand Down
73 changes: 70 additions & 3 deletions include/gproshan/laplacian/laplacian.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,83 @@
#define LAPLACIAN_H

#include <gproshan/mesh/che.h>
#include <gproshan/include_arma.h>

#include <armadillo>


// geometry processing and shape analysis framework
namespace gproshan {


void laplacian(const che * mesh, a_sp_mat & L, a_sp_mat & A);
template<class T>
void laplacian(const che * mesh, arma::SpMat<T> & L, arma::SpMat<T> & A)
{
size_t n_edges = mesh->n_edges;
size_t n_vertices = mesh->n_vertices;

arma::umat DI(2, 2 * n_edges);
arma::Col<T> DV(2 * n_edges);

arma::umat SI(2, n_edges);
arma::Col<T> SV(n_edges);

#pragma omp parallel for
for(index_t e = 0; e < n_edges; ++e)
{
index_t i = e << 1;

DI(0, i) = e;
DI(1, i) = mesh->edge_u(e);
DV(i) = -1;

++i;

DI(0, i) = e;
DI(1, i) = mesh->edge_v(e);
DV(i) = 1;

SI(0, e) = SI(1, e) = e;
SV(e) = (mesh->cotan(mesh->edge_he_0(e)) + mesh->cotan(mesh->edge_he_1(e))) / 2;
}

arma::SpMat D(DI, DV, n_edges, n_vertices);
arma::SpMat S(SI, SV, n_edges, n_edges);

L = D.t() * S * D;

A.eye(n_vertices, n_vertices);

#pragma omp parallel for
for(index_t v = 0; v < n_vertices; ++v)
A(v, v) = mesh->area_vertex(v);
}

template<class T>
size_t eigs_laplacian(const che * mesh, arma::Col<T> & eigval, arma::Mat<T> & eigvec,
arma::SpMat<T> & L, arma::SpMat<T> & A, const size_t k)
{
laplacian(mesh, L, A);

std::string feigval = tmp_file_path(mesh->name_size() + ".eigval");
std::string feigvec = tmp_file_path(mesh->name_size() + ".eigvec");

if(!eigval.load(feigval) || !eigvec.load(feigvec) || eigval.n_elem < k)
{
if(!eigs_sym(eigval, eigvec, L, k, "sm"))
return 0;

eigval.save(feigval);
eigvec.save(feigvec);
}

if(k < eigval.n_elem)
{
eigval = eigval.head(k);
eigvec = eigvec.head_cols(k);
}

size_t eigs_laplacian(const che * mesh, a_vec & eigval, a_mat & eigvec, a_sp_mat & L, a_sp_mat & A, const size_t k);
return eigval.n_elem;
}


} // namespace gproshan
Expand Down
19 changes: 0 additions & 19 deletions include/gproshan/mesh/che.cuh
View file

This file was deleted.

Loading

0 comments on commit 0f367f4

Please sign in to comment.