Merge branch 'dev' into alpha

.gitignore

@@ -13,3 +13,5 @@ tmp/*
 
 include/gproshan/config.h
 
+
+*.optixir

CMakeLists.txt

@@ -19,16 +19,16 @@ if(CUDAToolkit_FOUND)
 	enable_language(CUDA)
 	include_directories(SYSTEM ${CUDAToolkit_INCLUDE_DIRS})
 
-	set(CMAKE_CUDA_STANDARD 17)
+	set(CMAKE_CUDA_STANDARD 20)
 	set(CMAKE_CUDA_STANDARD_REQUIRED ON)
 
-	set(CMAKE_CUDA_FLAGS "-Xcompiler -fopenmp")
+	set(CMAKE_CUDA_FLAGS "-use_fast_math -Xcompiler -fopenmp")
 	set(CMAKE_CUDA_SEPARABLE_COMPILATION ON)
 	if(NOT DEFINED ENV{GITHUB_ACTIONS})
 		set(CMAKE_CUDA_ARCHITECTURES native)
 	endif()
 
-	find_package(OptiX 8)
+	find_package(OptiX 8.1)
 	if(OptiX_INCLUDE)
 		include_directories(SYSTEM ${OptiX_INCLUDE})
 	endif(OptiX_INCLUDE)
@@ -48,7 +48,7 @@ configure_file(	${gproshan_SOURCE_DIR}/include/gproshan/config.h.in
 				)
 
 
-find_package(embree 4 REQUIRED)
+find_package(embree 4.3 REQUIRED)
 find_package(OpenMP REQUIRED)
 find_package(OpenGL REQUIRED)
 find_package(GLEW REQUIRED)

cmake/FindOptiX.cmake

@@ -1,30 +1,31 @@
 #
-# Copyright (c) 2023, NVIDIA CORPORATION. All rights reserved.
-#
+# SPDX-FileCopyrightText: Copyright (c) 2010 - 2024 NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+# SPDX-License-Identifier: BSD-3-Clause
+# 
 # Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions
-# are met:
-#  * Redistributions of source code must retain the above copyright
-#    notice, this list of conditions and the following disclaimer.
-#  * Redistributions in binary form must reproduce the above copyright
-#    notice, this list of conditions and the following disclaimer in the
-#    documentation and/or other materials provided with the distribution.
-#  * Neither the name of NVIDIA CORPORATION nor the names of its
-#    contributors may be used to endorse or promote products derived
-#    from this software without specific prior written permission.
+# modification, are permitted provided that the following conditions are met:
 #
-# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY
-# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-# PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
-# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
-# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
-# OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
-# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+# 1. Redistributions of source code must retain the above copyright notice, this
+# list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+# this list of conditions and the following disclaimer in the documentation
+# and/or other materials provided with the distribution.
 #
+# 3. Neither the name of the copyright holder nor the names of its
+# contributors may be used to endorse or promote products derived from
+# this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 # Locate the OptiX distribution.  Search relative to the SDK first, then look in the system.
 

include/gproshan/geodesics/geodesics_ptp.h

@@ -13,11 +13,7 @@
 
 
 #ifdef __CUDACC__
-	#include <thrust/count.h>
-	#include <thrust/device_vector.h>
-	#include <thrust/execution_policy.h>
-
-	#define NT 64
+	#define NT 256
 	#define NB(x) (x + NT - 1) / NT
 #endif // __CUDACC__
 
@@ -34,18 +30,7 @@ __global__
 void relax_ptp(const che * mesh, float * new_dist, float * 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(float * error, const float * new_dist, const float * old_dist, const index_t start, const index_t end, const index_t * sorted = nullptr);
-
-struct is_ok
-{
-	const float * error = nullptr;
-
-	__host_device__
-	bool operator()(const float val) const;
-
-	__host_device__
-	bool operator()(const index_t val) const;
-};
+void relative_error(unsigned int * g_count, const float * new_dist, const float * old_dist, const index_t start, const index_t end, const index_t * sorted = nullptr);
 
 #endif // __CUDACC__
 
@@ -104,13 +89,9 @@ template<class T>
 __forceinline__
 #endif
 __host_device__
-float update_step(const che * mesh, const T * dist, const uvec3 & x)
+float update_step(const mat<T, 3> & points, const vec<T, 2> & t)
 {
-	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, 3> X[2] = {points[0] - points[2], points[1] - points[2]};
 
 	mat<T, 2> q;
 	q[0][0] = dot(X[0], X[0]);
@@ -143,7 +124,7 @@ float update_step(const che * mesh, const T * dist, const uvec3 & x)
 
 	if(t[0] == INFINITY || t[1] == INFINITY || dis < 0 || c[0] >= 0 || c[1] >= 0)
 	{
-		const vec<T, 2> & dp = {dist[x[0]] + norm(X[0]), dist[x[1]] + norm(X[1])};
+		const vec<T, 2> & dp = {t[0] + norm(X[0]), t[1] + norm(X[1])};
 		p = dp[dp[1] < dp[0]];
 	}
 
@@ -156,40 +137,45 @@ 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, const index_t * sorted, const index_t v, T * new_dist, T * old_dist, index_t * new_clusters, index_t * old_clusters)
 {
-	float & ndv = new_dist[v] = old_dist[v];
 	if(new_clusters) new_clusters[v] = old_clusters[v];
 
+	T ndv = old_dist[v];
+
+	mat<T, 3> X;
+	X[2] = mesh->point(v);
 	for(const index_t he: mesh->star(v))
 	{
-		const uvec3 i = {	mesh->halfedge(he_next(he)),
-							mesh->halfedge(he_prev(he)),
-							mesh->halfedge(he)
-							};
+		const uvec2 x = {mesh->halfedge(he_next(he)), mesh->halfedge(he_prev(he))};
+		const vec<T, 2> t = {old_dist[x[0]], old_dist[x[1]]};
 
-		float d = update_step(mesh, old_dist, i);
+		X[0] = mesh->point(x[0]);
+		X[1] = mesh->point(x[1]);
+
+		T d = update_step(X, t);
 
 		if(d < ndv)
 		{
 			ndv = d;
 			if(new_clusters)
-				new_clusters[v] = old_clusters[old_dist[i.y()] < old_dist[i.x()] ? i.y() : i.x()];
+				new_clusters[v] = old_clusters[x[t[1] < t[0]]];
 		}
 	}
+
+	new_dist[v] = ndv;
 }
 
 
-template<class T>
 #ifdef __CUDACC__
-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
+	__managed__ index_t count;
+#endif
+
+
+template<class T>
 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
 {
 #ifdef __CUDACC__
 	T * h_dist = dist[2];
@@ -225,10 +211,13 @@ index_t run_ptp(const che * mesh, const std::vector<index_t> & sources,
 	}
 #endif
 
+#ifndef __CUDACC__
+	index_t count = 0;
+#endif
+
 	const int max_iter = size(limits) << 1;
 
 	int iter = -1;
-	index_t count = 0;
 	index_t i = 1;
 	index_t j = 2;
 	while(i < j && ++iter < max_iter)
@@ -249,16 +238,12 @@ index_t run_ptp(const che * mesh, const std::vector<index_t> & sources,
 		relax_ptp<<< NB(end - start), NT >>>(mesh, new_dist, old_dist, new_cluster, old_cluster, start, end, sorted);
 		cudaDeviceSynchronize();
 
-		relative_error<<< NB(n_cond), NT >>>(error, new_dist, old_dist, start, start + n_cond, sorted);
+		relative_error<<< NB(n_cond), NT >>>(&count, new_dist, old_dist, start, start + n_cond, sorted);
 		cudaDeviceSynchronize();
-
-		count = sorted ? thrust::count_if(thrust::device, sorted + start, sorted + start + n_cond, is_ok{error})
-						: thrust::count_if(thrust::device, error + start, error + start + n_cond, is_ok{});
 	#else
 		#pragma omp parallel for
 		for(index_t v = start; v < end; ++v)
-			relax_ptp(mesh, new_dist, old_dist, new_cluster, old_cluster, sorted ? sorted[v] : v);
-
+			relax_ptp(mesh, sorted, sorted ? sorted[v] : v, new_dist, old_dist, new_cluster, old_cluster);
 
 		count = 0;
 		#pragma omp parallel for

include/gproshan/mesh/quaternion.h → include/gproshan/geometry/quaternion.h

@@ -10,55 +10,39 @@
 namespace gproshan {
 
 
-using vertex = vec3;
-
-
 class quaternion
 {
-	public:
-		float s;
-		vertex v;
+	private:
+		float s = 0;
+		vec3 v;
 
 	public:
-		quaternion(float s = 0, float vi = 0, float vj = 0, float vk = 0);
-		quaternion(float s, const vertex & v);
-		quaternion(const vertex & v);
+		quaternion(const vec3 & v = {});
+		quaternion(float s, const vec3 & v = {});
 
-		operator const vertex & () const;
-		const quaternion & operator = (float s);
-		const quaternion & operator = (const vertex & v);
+		operator const vec3 & () const;
 		float & operator [] (int index);
 		float operator [] (int index) const;
-		float & re(void);
-		float re(void) const;
-		vertex & im(void);
-		const vertex & im(void) const;
 
 		quaternion operator + (const quaternion & q) const;
 		quaternion operator - (const quaternion & q) const;
-		quaternion operator - (void) const;
+		quaternion operator - () const;
 		quaternion operator * (float c) const;
 		quaternion operator / (float c) const;
-		void operator += (const quaternion & q);
-		void operator += (float c);
-		void operator -= (const quaternion & q);
-		void operator -= (float c);
-		void operator *= (float c);
-		void operator /= (float c);
 		quaternion operator * (const quaternion & q) const;
-		void operator *= (const quaternion & q);
 
 		quaternion conj() const;
 		quaternion inv() const;
 		float norm() const;
 		float norm2() const;
-		quaternion unit() const;
-		void normalize();
 
 	friend std::ostream & operator << (std::ostream & os, const quaternion & q);
 	friend std::istream & operator >> (std::istream & is, quaternion & q);
 };
 
+
+float norm(const quaternion & q);
+quaternion normalize(const quaternion & q);
 quaternion operator * (float c, const quaternion & q);
 
 

include/gproshan/mesh/che.h

@@ -128,6 +128,13 @@ class che
 		return GT[v];
 	}
 
+	__host_device__
+	const vertex & vertex_he(const index_t he) const
+	{
+		assert(he < n_half_edges);
+		return GT[VT[he]];
+	}
+
 	__host_device__
 	const vertex & normal(const index_t v) const
 	{
@@ -266,7 +273,6 @@ class che
 	index_t edge_v(const index_t e) const;
 	index_t edge_he_0(const index_t e) const;
 	index_t edge_he_1(const index_t e) const;
-	const vertex & vertex_he(const index_t he) const;
 	const vertex & vertex_edge_u(const index_t e) const;
 	const vertex & vertex_edge_v(const index_t e) const;
 	index_t evt(const index_t v) const;

include/gproshan/pointcloud/knn.h

@@ -74,10 +74,12 @@ float mean_mean_knn_distant(const point * pc, const size_t n_points, const size_
 float mean_knn_area_radius(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
 float median_knn_area_radius(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
 
-
+float voronoi_radius(const point * pc, const int * id, const size_t n, const mat4 & model_mat);
 float median_pair_dist(const point * pc, const int * id, const size_t n, const mat4 & model_mat);
 float mean_knn(const point * pc, const int * id, const size_t n, const mat4 & model_mat);
 
+std::vector<float> anisotropic(const point * pc, const size_t n_points, const k3tree & nn, const int k);
+
 
 const char * radius_str(void *, int opt);
 

include/gproshan/raytracing/embree.h

@@ -22,15 +22,17 @@ class embree : public raytracing
 						MEAN,
 						MEDIAN,
 						AREA,
-						MEDIAN_PAIRS
+						MEDIAN_PAIRS,
+						VORONOI
 						};
 
 		struct pc_opts
 		{
+			knn_opt opt		= NONE;
 			bool enable		= false;
 			bool normals	= false;
+			bool anisotropy = false;
 			float radius	= 0.01;
-			knn_opt opt		= NONE;
 			float scale		= 1;
 			int knn			= 8;
 

include/gproshan/raytracing/optix.h

@@ -52,7 +52,7 @@ class optix : public raytracing
 		std::vector<unsigned char *> tex_data;
 
 	public:
-		optix(const std::string & ptx = "/src/optix.ptx");
+		optix(const std::string & program = "/src/optix.optixir");
 		optix(const std::vector<const che *> & meshes, const std::vector<mat4> & model_mats);
 		~optix();
 

include/gproshan/raytracing/utils.h

@@ -183,7 +183,7 @@ struct t_eval_hit
 	bool scatter_diffuse(vec<T, 3> & dir, random<T> & rnd)
 	{
 		// random unit sphere
-		const T theta = rnd() * 2.f * M_PI;
+		const T theta = rnd() * 2.f * 3.141592654f;
 		const T phi = acosf(2.f * rnd() - 1.f);
 		const T r = cbrtf(rnd());