[CC]点云密度计算

本文转载自查看原文 2016-09-15 15:51 3424 PointCloud/ CloudCompare

包括两种计算方法：精确计算和近似计算（思考：local density=单位面积的点数 vs local density =1/单个点所占的面积）

每种方法可以实现三种模式的点云密度计算，CC里面的点云计算依赖于给定的近邻半径所对应的最佳八叉树层级（通过findBestLevelForAGivenNeighbourhoodSizeExtraction()方法实现）

在GeometricalAnalysisTools类文件中实现。

//volume of a unit sphere
static double s_UnitSphereVolume = 4.0 * M_PI / 3.0;

1.精确计算

 1 int GeometricalAnalysisTools::computeLocalDensity(    GenericIndexedCloudPersist* theCloud,  2  Density densityType,  3  PointCoordinateType kernelRadius,  4                                                     GenericProgressCallback* progressCb/*=0*/,  5                                                     DgmOctree* inputOctree/*=0*/)  6 {  7     if (!theCloud)  8         return -1;  9 
10     unsigned numberOfPoints = theCloud->size(); 11     if (numberOfPoints < 3) 12         return -2; 13 
14     //compute the right dimensional coef based on the expected output
15     double dimensionalCoef = 1.0; 16     switch (densityType) 17  { 18     case DENSITY_KNN: 19         dimensionalCoef = 1.0; 20         break; 21     case DENSITY_2D: 22         dimensionalCoef = M_PI * (static_cast<double>(kernelRadius) * kernelRadius); 23         break; 24     case DENSITY_3D: 25         dimensionalCoef = s_UnitSphereVolume * ((static_cast<double>(kernelRadius) * kernelRadius) * kernelRadius); 26         break; 27     default: 28         assert(false); 29         return -5; 30  } 31 
32     DgmOctree* theOctree = inputOctree; 33     if (!theOctree) 34  { 35         theOctree = new DgmOctree(theCloud); 36         if (theOctree->build(progressCb) < 1) 37  { 38             delete theOctree; 39             return -3; 40  } 41  } 42 
43     theCloud->enableScalarField(); 44 
45     //determine best octree level to perform the computation
46     unsigned char level = theOctree->findBestLevelForAGivenNeighbourhoodSizeExtraction(kernelRadius); 47 
48     //parameters
49     void* additionalParameters[] = {    static_cast<void*>(&kernelRadius), 50                                         static_cast<void*>(&dimensionalCoef) }; 51 
52     int result = 0; 53 
54     if (theOctree->executeFunctionForAllCellsAtLevel( level, 55                                                         &computePointsDensityInACellAtLevel, 56  additionalParameters, 57                                                         true, 58  progressCb, 59                                                         "Local Density Computation") == 0) 60  { 61         //something went wrong
62         result = -4; 63  } 64 
65     if (!inputOctree) 66         delete theOctree; 67 
68     return result; 69 }

GeometricalAnalysisTools::computeLocalDensity

 1 //"PER-CELL" METHOD: LOCAL DENSITY  2 //ADDITIONNAL PARAMETERS (2):  3 // [0] -> (PointCoordinateType*) kernelRadius : spherical neighborhood radius  4 // [1] -> (ScalarType*) sphereVolume : spherical neighborhood volume
 5 bool GeometricalAnalysisTools::computePointsDensityInACellAtLevel(    const DgmOctree::octreeCell& cell,  6                                                                     void** additionalParameters,  7                                                                     NormalizedProgress* nProgress/*=0*/)  8 {  9     //parameter(s)
10     PointCoordinateType radius = *static_cast<PointCoordinateType*>(additionalParameters[0]); 11     double dimensionalCoef = *static_cast<double*>(additionalParameters[1]); 12     
13     assert(dimensionalCoef > 0); 14 
15     //structure for nearest neighbors search
16  DgmOctree::NearestNeighboursSphericalSearchStruct nNSS; 17     nNSS.level = cell.level; 18     nNSS.prepare(radius,cell.parentOctree->getCellSize(nNSS.level)); 19     cell.parentOctree->getCellPos(cell.truncatedCode,cell.level,nNSS.cellPos,true); 20     cell.parentOctree->computeCellCenter(nNSS.cellPos,cell.level,nNSS.cellCenter); 21 
22     unsigned n = cell.points->size(); //number of points in the current cell 23     
24     //for each point in the cell
25     for (unsigned i=0; i<n; ++i) 26  { 27         cell.points->getPoint(i,nNSS.queryPoint); 28 
29         //look for neighbors inside a sphere 30         //warning: there may be more points at the end of nNSS.pointsInNeighbourhood than the actual nearest neighbors (neighborCount)!
31         unsigned neighborCount = cell.parentOctree->findNeighborsInASphereStartingFromCell(nNSS,radius,false); 32         //数目/体积
33         ScalarType density = static_cast<ScalarType>(neighborCount/dimensionalCoef); 34         cell.points->setPointScalarValue(i,density); 35 
36         if (nProgress && !nProgress->oneStep()) 37  { 38             return false; 39  } 40  } 41 
42     return true; 43 }

computePointsDensityInACellAtLevel

2. 近似计算

 1 int GeometricalAnalysisTools::computeLocalDensityApprox(GenericIndexedCloudPersist* theCloud,
 2                                                         Density densityType,
 3                                                         GenericProgressCallback* progressCb/*=0*/,
 4                                                         DgmOctree* inputOctree/*=0*/)
 5 {
 6     if (!theCloud)
 7         return -1;
 8 
 9     unsigned numberOfPoints = theCloud->size();
10     if (numberOfPoints < 3)
11         return -2;
12 
13     DgmOctree* theOctree = inputOctree;
14     if (!theOctree)
15     {
16         theOctree = new DgmOctree(theCloud);
17         if (theOctree->build(progressCb) < 1)
18         {
19             delete theOctree;
20             return -3;
21         }
22     }
23 
24     theCloud->enableScalarField();
25 
26     //determine best octree level to perform the computation
27     unsigned char level = theOctree->findBestLevelForAGivenPopulationPerCell(3);
28 
29     //parameters
30     void* additionalParameters[] = { static_cast<void*>(&densityType) };
31 
32     int result = 0;
33 
34     if (theOctree->executeFunctionForAllCellsAtLevel(    level,
35                                                         &computeApproxPointsDensityInACellAtLevel,
36                                                         additionalParameters,
37                                                         true,
38                                                         progressCb,
39                                                         "Approximate Local Density Computation") == 0)
40     {
41         //something went wrong
42         result = -4;
43     }
44 
45     if (!inputOctree)
46         delete theOctree;
47 
48     return result;
49 }

View Code

 1 //"PER-CELL" METHOD: APPROXIMATE LOCAL DENSITY
 2 //ADDITIONAL PARAMETERS (0): NONE
 3 bool GeometricalAnalysisTools::computeApproxPointsDensityInACellAtLevel(const DgmOctree::octreeCell& cell,
 4                                                                         void** additionalParameters,
 5                                                                         NormalizedProgress* nProgress/*=0*/)
 6 {
 7     //extract additional parameter(s)
 8     Density densityType = *static_cast<Density*>(additionalParameters[0]);
 9     
10     DgmOctree::NearestNeighboursSearchStruct nNSS;
11     nNSS.level                                = cell.level;
12     nNSS.alreadyVisitedNeighbourhoodSize    = 0;
13     nNSS.minNumberOfNeighbors                = 2;
14     cell.parentOctree->getCellPos(cell.truncatedCode,cell.level,nNSS.cellPos,true);
15     cell.parentOctree->computeCellCenter(nNSS.cellPos,cell.level,nNSS.cellCenter);
16 
17     unsigned n = cell.points->size();
18     for (unsigned i=0; i<n; ++i)
19     {
20         cell.points->getPoint(i,nNSS.queryPoint);
21 
22         //the first point is always the point itself!
23         if (cell.parentOctree->findNearestNeighborsStartingFromCell(nNSS) > 1)
24         {
25             double R2 = nNSS.pointsInNeighbourhood[1].squareDistd;
26 
27             ScalarType density = NAN_VALUE;
28             if (R2 > ZERO_TOLERANCE)
29             {
30                 switch (densityType)
31                 {
32                 case DENSITY_KNN:
33                     {
34                         //we return in fact the (inverse) distance to the nearest neighbor
35                         density = static_cast<ScalarType>(1.0 / sqrt(R2));
36                     }
37                     break;
38                 case DENSITY_2D:
39                     {
40                         //circle area (2D approximation)
41                         double circleArea = M_PI * R2;
42                         density = static_cast<ScalarType>(1.0 / circleArea);
43                     }
44                     break;
45                 case DENSITY_3D:
46                     {
47                         //sphere area
48                         double sphereArea =  s_UnitSphereVolume * R2 * sqrt(R2);
49                         density = static_cast<ScalarType>(1.0 / sphereArea);
50                     }
51                     break;
52                 default:
53                     assert(false);
54                     break;
55                 }
56             }
57             cell.points->setPointScalarValue(i,density);
58         }
59         else
60         {
61             //shouldn't happen! Apart if the cloud has only one point...
62             cell.points->setPointScalarValue(i,NAN_VALUE);
63         }
64 
65         if (nProgress && !nProgress->oneStep())
66         {
67             return false;
68         }
69     }
70 
71     return true;
72 }

computeApproxPointsDensityInACellAtLevel

double R2 = nNSS.pointsInNeighbourhood[1].squareDistd; //索引为1的点，表示最近邻点。pi点与最近邻点之间距离的平方。

免责声明！

本站转载的文章为个人学习借鉴使用，本站对版权不负任何法律责任。如果侵犯了您的隐私权益，请联系本站邮箱yoyou2525@163.com删除。

猜您在找 matlab——功率谱密度计算 2020年四月点密度、线密度与核密度 PCL的PNG文件和计算点云重心点云/网格模型的体积计算 LiDAR点云利用Arcgis基于DEM法计算土方量 CloudCompare点云体积计算方法 MATLAB实例：散点密度图 PPI是什么？pixels per inch像素密度是什么？PPI如何计算？沟谷网络的提取及沟壑密度的计算（ArcPy实现）点云滤波简介