How to compute the density of a 3D point cloud?

Question

Ali 2020-7-12

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/563603-how-to-compute-the-density-of-a-3d-point-cloud

评论： Ali ，2020-7-12

3_buildings.zip

I'm trying to write a program that operates on 3D point clouds (.ply, .pcd), I need to know the density of the given point cloud file and compare it with a threshold or certain percentage to decide for some operations. How can I do that?

I have found this for 3D point clouds on the web:

Two methods can be used to compute the density:

either 'Precise': the density is estimated by counting for each point the number of neighbors N (inside a sphere of radius R)
or 'Approximate': the density is simply estimated by determining the distance to the nearest neighbor (which is generally much faster). This distance is considered as being equivalent to the above spherical neighborhood radius R (and N = 1).

However I'm not sure how to put this into matlab code, I guess I should use pcfitplane for sphere fitting and findNearestNeighbors. Would appreciate a step by step example.

0 个评论
显示 -1更早的评论隐藏 -1更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Thiago Henrique Gomes Lobato 2020-7-12

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/563603-how-to-compute-the-density-of-a-3d-point-cloud#answer_464621

编辑：Thiago Henrique Gomes Lobato 2020-7-12

You can just directly apply the definition you gave considering that the density is N/Volume. The easiest is the second one, which would be:

AllPoints = % your points
K = 1;
for idx=1:length(AllPoints)
[~,r] = findNearestNeighbors(AllPoints,AllPoints(idx,:),K);
density(idx) = 1/(4*pi*r.^3/3);
end

The first one is a little more complicated but at the same time not so much:

AllPoints = % your points
R = 1; % depends on your data
for idx=1:length(AllPoints)
Distances = sqrt( sum( (AllPoints-AllPoints(idx,:)).^2 ,2) );
Ninside   = length( find(Distances<=R) );
density(idx) = Ninside/(4*pi*R.^3/3);
end

6 个评论
显示 5更早的评论隐藏 5更早的评论

Ali 2020-7-12

编辑：Ali 2020-7-12

Thank you for your quick reply @Thiago I managed to get the following piece of code working thanks to your help, the 'Approximate' method (your first block of code) produces an 1xN array (N being the total point count in the cloud) with all elements of "Inf", however the 'Precise' method produces the same array size with elements of different numbers ranging from 0 to 4 with some trailing fractions.

This may sound noob, but if I want to get a single number representing the whole point cloud's density what should I do with this array? Should I divide the sum of all values by the number of total points existing in the cloud or should I calculate the mean of this array?

Also the radius variable 'R' being equal to 1 in our example, is there a way to compute it for different given point clouds before running the for loop like finding the center of the point cloud and calculating the radius to the largest possible direction so that it covers the whole point cloud?

clc
clear 
close all
% input point cloud file gui
[FileName,PathName] = uigetfile({'*.ply';'*.pcd'},'Input point cloud');
ptCloud=pcread([PathName,FileName]);
AllPoints = ptCloud.Location;
% Approximate method
K = 1;
for idx=1:length(AllPoints)
[~,r] = findNearestNeighbors(ptCloud,AllPoints(idx,:),K);
densityApprox(idx) = 1/(4*pi*r.^3/3);
end
fprintf('\n Density (Approx.): %d \n', densityApprox);
% Precise method
R = 1; % depends on your data
for idx=1:length(AllPoints)
Distances = sqrt( sum( (AllPoints-AllPoints(idx,:)).^2 ,2) );
Ninside   = length( find(Distances<=R) );
densityPrecise(idx) = Ninside/(4*pi*R.^3/3);
end
fprintf('\n Density (Precise): %d \n', densityPrecise);