diff --git a/include/gproshan/pointcloud/knn.h b/include/gproshan/pointcloud/knn.h
index 31fe8af8..e9d4972f 100644
--- a/include/gproshan/pointcloud/knn.h
+++ b/include/gproshan/pointcloud/knn.h
@@ -56,6 +56,7 @@ class k3tree
 
 	public:
 		k3tree(const point * pc, const size_t n_points, const size_t k = 8, const std::vector<point> & query = {});
+		k3tree(const point * pc, const size_t n_points, const point * query, const size_t n_query, const size_t k = 8);
 		~k3tree();
 
 		const int * operator [] (const index_t i) const;
@@ -64,7 +65,23 @@ class k3tree
 };
 
 
-real_t pc_median_pairwise_distant(const point * pc, const size_t n_points, const mat4 & model_mat = mat4::identity());
+real_t mean_knn_distant(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
+real_t median_knn_distant(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
+real_t median_median_knn_distant(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
+real_t mean_median_knn_distant(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
+real_t median_mean_knn_distant(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
+real_t mean_mean_knn_distant(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
+real_t mean_knn_area_radius(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
+real_t median_knn_area_radius(const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
+
+
+real_t median_pair_dist(const point * pc, const int * id, const size_t n, const mat4 & model_mat);
+real_t mean_knn(const point * pc, const int * id, const size_t n, const mat4 & model_mat);
+
+
+const char * radius_str(void *, int opt);
+
+real_t radius(const int opt, const point * pc, const size_t n_points, const size_t k = 8, const mat4 & model_mat = mat4::identity());
 
 
 } // namespace gproshan
diff --git a/include/gproshan/raytracing/embree.h b/include/gproshan/raytracing/embree.h
index 8b01b240..a7010a2e 100644
--- a/include/gproshan/raytracing/embree.h
+++ b/include/gproshan/raytracing/embree.h
@@ -17,12 +17,22 @@ namespace gproshan::rt {
 class embree : public raytracing
 {
 	public:
+		enum knn_opt {	NONE,
+						MAX,
+						MEAN,
+						MEDIAN,
+						AREA,
+						MEDIAN_PAIRS
+						};
 		struct pc_opts
 		{
 			bool enable		= false;
 			bool normals	= false;
 			float radius	= 0.01;
-			int knn			= 0;
+			knn_opt opt		= NONE;
+			float scale		= 1;
+			int knn			= 8;
+
 
 			pc_opts() {};
 		};
diff --git a/src/gproshan/pointcloud/knn.cpp b/src/gproshan/pointcloud/knn.cpp
index 2ff08e62..56bcc6a5 100644
--- a/src/gproshan/pointcloud/knn.cpp
+++ b/src/gproshan/pointcloud/knn.cpp
@@ -71,8 +71,11 @@ std::vector<index_t> grid::operator () (const point & p, int k) const
 }
 
 
-///< Implementation using flann, by default compute all knn
 k3tree::k3tree(const point * pc, const size_t n_points, const size_t k, const std::vector<point> & query)
+		: k3tree(pc, n_points, size(query) ? query.data() : nullptr, size(query), k) {}
+
+///< Implementation using flann, by default compute all knn
+k3tree::k3tree(const point * pc, const size_t n_points, const point * query, const size_t n_query, const size_t k)
 {
 	double time_build, time_query;
 
@@ -84,8 +87,8 @@ k3tree::k3tree(const point * pc, const size_t n_points, const size_t k, const st
 	gproshan_log_var(time_build);
 
 	TIC(time_query);
-		const point * q = size(query) ? query.data() : pc;
-		const size_t n_results = size(query) ? size(query) : n_points;
+		const point * q = query && n_query ? query : pc;
+		const size_t n_results = query && n_query ? n_query : n_points;
 
 		flann::Matrix<real_t> mq((real_t *) q, n_results, 3);
 
@@ -93,7 +96,7 @@ k3tree::k3tree(const point * pc, const size_t n_points, const size_t k, const st
 		flann::Matrix dists(new real_t[n_results * k], n_results, k);
 
 		flann::SearchParams params;
-		params.cores = 16;
+		params.cores = 0;
 		index.knnSearch(mq, indices, dists, k, params);
 	TOC(time_query);
 	gproshan_log_var(time_query);
@@ -124,21 +127,197 @@ int k3tree::operator () (const index_t i, const index_t j) const
 }
 
 
-real_t pc_median_pairwise_distant(const point * pc, const size_t n_points, const mat4 & model_mat)
+real_t mean_knn_distant(const point * pc, const size_t n_points, const size_t k, const mat4 & model_mat)
+{
+	k3tree p2nn(pc, n_points, k + 1);
+
+	real_t mean = 0;
+
+	#pragma omp parallel for
+	for(index_t i = 0; i < n_points; ++i)
+	{
+		#pragma omp atomic
+		mean += length(model_mat * (pc[i] - pc[p2nn(i, k)], 0));
+	}
+
+	return mean / n_points;
+}
+
+real_t median_knn_distant(const point * pc, const size_t n_points, const size_t k, const mat4 & model_mat)
+{
+	k3tree nn(pc, n_points, k + 1);
+
+	std::vector<real_t> dist(n_points);
+
+	#pragma omp parallel for
+	for(index_t i = 0; i < n_points; ++i)
+		dist[i] = length(model_mat * (pc[i] - pc[nn(i, k)], 0));
+
+	std::ranges::sort(dist);
+
+	return dist[size(dist) >> 1];
+}
+
+real_t median_median_knn_distant(const point * pc, const size_t n_points, const size_t k, const mat4 & model_mat)
 {
-	k3tree p2nn(pc, n_points, 2);
+	k3tree nn(pc, n_points, k + 1);
 
 	std::vector<real_t> dist(n_points);
 
 	#pragma omp parallel for
 	for(index_t i = 0; i < n_points; ++i)
-		dist[i] = length(model_mat * (pc[i] - pc[p2nn(i, 1)], 0));
+		dist[i] = length(model_mat * (pc[i] - pc[nn(i, (k >> 1) + 1)], 0));
+
+	std::ranges::sort(dist);
+
+	return dist[size(dist) >> 1];
+}
+
+
+real_t mean_median_knn_distant(const point * pc, const size_t n_points, const size_t k, const mat4 & model_mat)
+{
+	k3tree nn(pc, n_points, k + 1);
+
+	real_t mean = 0;
+
+	#pragma omp parallel for
+	for(index_t i = 0; i < n_points; ++i)
+	{
+		#pragma omp atomic
+		mean += length(model_mat * (pc[i] - pc[nn(i, (k >> 1) + 1)], 0));
+	}
+
+	return mean / n_points;
+}
+
+real_t median_mean_knn_distant(const point * pc, const size_t n_points, const size_t k, const mat4 & model_mat)
+{
+	k3tree nn(pc, n_points, k + 1);
+
+	std::vector<real_t> dist(n_points);
+
+	#pragma omp parallel for
+	for(index_t i = 0; i < n_points; ++i)
+		dist[i] = mean_knn(pc, nn(i), k, model_mat);
+
+	std::ranges::sort(dist);
+
+	return dist[size(dist) >> 1];
+}
+
+real_t mean_mean_knn_distant(const point * pc, const size_t n_points, const size_t k, const mat4 & model_mat)
+{
+	k3tree nn(pc, n_points, k + 1);
+
+	real_t mean = 0;
+
+	#pragma omp parallel for
+	for(index_t i = 0; i < n_points; ++i)
+	{
+		#pragma omp atomic
+		mean += mean_knn(pc, nn(i), k, model_mat);
+	}
+
+	return mean / n_points;
+}
+
+real_t mean_knn_area_radius(const point * pc, const size_t n_points, const size_t k, const mat4 & model_mat)
+{
+	k3tree nn(pc, n_points, k + 1);
+
+	real_t mean_r = 0;
+
+	#pragma omp parallel for
+	for(index_t i = 0; i < n_points; ++i)
+	{
+		real_t r = length(model_mat * (pc[i] - pc[nn(i, k)], 0));
+
+		#pragma omp atomic
+		mean_r += sqrt(r * r / k);
+	}
+
+	gproshan_log_var(mean_r);
+
+	return mean_r / n_points;
+}
+
+real_t median_knn_area_radius(const point * pc, const size_t n_points, const size_t k, const mat4 & model_mat)
+{
+	k3tree nn(pc, n_points, k + 1);
+
+	std::vector<real_t> radius(n_points);
+
+	#pragma omp parallel for
+	for(index_t i = 0; i < n_points; ++i)
+	{
+		real_t r = length(model_mat * (pc[i] - pc[nn(i, k)], 0));
+
+		radius[i] = sqrt(r * r / k);
+	}
+
+	std::ranges::sort(radius);
+
+	return radius[size(radius) >> 1];
+}
+
+
+real_t median_pair_dist(const point * pc, const int * id, const size_t n, const mat4 & model_mat)
+{
+	std::vector<real_t> dist;
+	dist.reserve((n * (n - 1)) >> 1);
+
+	for(index_t i = 0; i < n; ++i)
+	for(index_t j = i + 1; j < n; ++j)
+		dist.push_back(length(model_mat * (pc[id[i]] - pc[id[j]], 0)));
 
 	std::ranges::sort(dist);
 
 	return dist[size(dist) >> 1];
 }
 
+real_t mean_knn(const point * pc, const int * id, const size_t n, const mat4 & model_mat)
+{
+	real_t mean = 0;
+	for(index_t i = 1; i <= n; ++i)
+		mean += length(model_mat * (pc[id[0]] - pc[id[i]], 0));
+
+	return mean / n;
+}
+
+
+const char * radius_str(void *, int opt)
+{
+	static const char * str[] = {	"none",
+									"mean_knn_distant",
+									"median_knn_distant",
+									"median_median_knn_distant",
+									"mean_median_knn_distant",
+									"median_mean_knn_distant",
+									"mean_mean_knn_distant",
+									"mean_knn_area_radius",
+									"median_knn_area_radius"
+									};
+
+	return str[opt];
+}
+
+real_t radius(const int opt, const point * pc, const size_t n_points, const size_t k, const mat4 & model_mat)
+{
+	using fknn = real_t (*) (const point *, const size_t, const size_t, const mat4 &);
+	static const fknn f[] = {	nullptr,
+								mean_knn_distant,
+								median_knn_distant,
+								median_median_knn_distant,
+								mean_median_knn_distant,
+								median_mean_knn_distant,
+								mean_mean_knn_distant,
+								mean_knn_area_radius,
+								median_knn_area_radius
+								};
+
+	return opt ? f[opt](pc, n_points, k, model_mat) : 0;
+}
+
 
 } // namespace gproshan
 
diff --git a/src/gproshan/raytracing/embree.cpp b/src/gproshan/raytracing/embree.cpp
index d68b9ec5..a63606d0 100644
--- a/src/gproshan/raytracing/embree.cpp
+++ b/src/gproshan/raytracing/embree.cpp
@@ -216,19 +216,48 @@ index_t embree::add_pointcloud(const che * mesh, const mat4 & model_mat, const p
 																			);
 
 
-	knn::k3tree * nn = !pc.knn	? nullptr
-								: new knn::k3tree(&mesh->point(0), mesh->n_vertices, pc.knn + 1);
+	knn::k3tree * nn = pc.opt == NONE ? nullptr : new knn::k3tree(&mesh->point(0), mesh->n_vertices, pc.knn + 1);
 
 
 	#pragma omp parallel for
 	for(index_t i = 0; i < mesh->n_vertices; ++i)
 		pxyzr[i] = model_mat * (mesh->point(i), 1);
 
-	const real_t r = pc.radius < 0 && !nn ? knn::pc_median_pairwise_distant(&mesh->point(0), mesh->n_vertices, model_mat) : pc.radius;
-
 	#pragma omp parallel for
 	for(index_t i = 0; i < mesh->n_vertices; ++i)
-		pxyzr[i][3] = nn ? length(model_mat * (mesh->point(i) - mesh->point((*nn)(i, pc.knn)), 0)) : r;
+	{
+		real_t r = 0;
+
+		switch(pc.opt)
+		{
+			case NONE:
+					r = pc.radius;
+					break;
+
+			case MAX:
+					r = length(model_mat * (mesh->point(i) - mesh->point((*nn)(i, pc.knn)), 0));
+					break;
+
+			case MEAN:
+					r = knn::mean_knn(&mesh->point(0), (*nn)(i), pc.knn, model_mat);
+					break;
+
+			case MEDIAN:
+					r = length(model_mat * (mesh->point(i) - mesh->point((*nn)(i, pc.knn >> 1)), 0));
+					break;
+
+			case AREA:
+					r = length(model_mat * (mesh->point(i) - mesh->point((*nn)(i, pc.knn)), 0));
+					r = sqrt(r * r / pc.knn);
+					break;
+
+			case MEDIAN_PAIRS:
+					r = knn::median_pair_dist(&mesh->point(0), (*nn)(i), pc.knn, model_mat);
+					break;
+		};
+
+		pxyzr[i][3] = pc.scale * r;
+	}
 
 
 	delete nn;
diff --git a/src/gproshan/viewer/viewer.cpp b/src/gproshan/viewer/viewer.cpp
index 1f5511f5..33bea2a3 100644
--- a/src/gproshan/viewer/viewer.cpp
+++ b/src/gproshan/viewer/viewer.cpp
@@ -14,6 +14,7 @@
 #include <gproshan/mesh/che_xyz.h>
 #include <gproshan/mesh/che_pts.h>
 #include <gproshan/mesh/che_pcd.h>
+#include <gproshan/pointcloud/knn.h>
 #include <gproshan/viewer/scene_viewer.h>
 #include <gproshan/viewer/glfw_keys.h>
 
@@ -978,6 +979,7 @@ bool viewer::m_setup_raytracing(viewer * view)
 	static int rt = 0;
 	static double time = 0;
 	static rt::embree::pc_opts pc;
+	static int opt_radius = 0;
 
 	ImGui::SliderInt("depth", (int *) &view->render_params.depth, 1, 1 << 5);
 	ImGui::SliderInt("n_samples", (int *) &view->render_params.n_samples, 1, 1 << 5);
@@ -986,8 +988,15 @@ bool viewer::m_setup_raytracing(viewer * view)
 	if(rt == R_EMBREE && (mesh.render_pointcloud || mesh->is_pointcloud()))
 	{
 		ImGui::Indent();
-		ImGui::SliderInt("pc.knn", &pc.knn, 0, 1 << 5);
-		ImGui::SliderFloat("pc.radius", &pc.radius, 0, 1);
+		ImGui::Combo("pc.opt", (int *) &pc.opt, "NONE\0MAX\0MEAN\0MEDIAN\0AREA\0MEDIAN_PAIRS\0\0");
+		if(pc.opt == rt::embree::NONE)
+		{
+			ImGui::SliderFloat("pc.radius", &pc.radius, 0, 1);
+			if(ImGui::Combo("knn_radius", &opt_radius, knn::radius_str, nullptr, 5, 10))
+				pc.radius = knn::radius(opt_radius, &mesh->point(0), mesh->n_vertices, pc.knn, mesh.model_mat);
+		}
+		ImGui::SliderInt("pc.knn", &pc.knn, 0, 1 << 6);
+		ImGui::SliderFloat("pc.scale", &pc.scale, 0, 10);
 		ImGui::Checkbox("pc.normals", &pc.normals);
 		ImGui::Unindent();
 	}