包括兩種計算方法：精确計算和近似計算（思考：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點與最近鄰點之間距離的平方。