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 output15     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 computation46     unsigned char level = theOctree->findBestLevelForAGivenNeighbourhoodSizeExtraction(kernelRadius);47 48     //parameters49     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 wrong62         result = -4;63     }64 65     if (!inputOctree)66         delete theOctree;67 68     return result;69 }

 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 search16     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 cell23     24     //for each point in the cell25     for (unsigned i=0; i<n; ++i)26     {27         cell.points->getPoint(i,nNSS.queryPoint);28 29         //look for neighbors inside a sphere30         //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 }

 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 computation27     unsigned char level = theOctree->findBestLevelForAGivenPopulationPerCell(3);28 29     //parameters30     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 wrong42         result = -4;43     }44 45     if (!inputOctree)46         delete theOctree;47 48     return result;49 }

 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 neighbor35                         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 area48                         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         else60         {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 }