点云密度计算方法matlab,点云中的几何计算及matlab源码

点云密度计算⽅法matlab,点云中的⼏何计算及matlab源码1、计算法向量2、计算曲率曲线的曲率(curvature)就是针对曲线上某个点的切线⽅向⾓对弧长的转动率，通过微分来定义，表明曲线偏离直线的程度。数学上表明曲线在某⼀点的弯曲程度的数值。曲率越⼤，表⽰曲线的弯曲程度越⼤。曲率的倒数就是曲率半径。平均曲率、⾼斯曲率、主曲率计算贺美芳,etal.(2005.8)."散乱点云数据的曲率估算及应⽤."ScalarTypeNeighbourhood::computeCurvature(unsignedneighbourIndex,CC_CURVATURE_TYPEcType){switch(cType){caseGAUSSIAN_CURV:caseMEAN_CURV:{//weget2D1/2quadricparametersconstPointCoordinateType*H=getQuadric();if(!H)returnNAN_VALUE;//computecentroidconstCCVector3*G=getGravityCenter();//wecomputecurvatureattheinputneighbourposition+werecenteritbythewayCCVector3Q=*m_associatedCloud->getPoint(neighbourIndex)-*G;constunsignedcharX=m_quadricEquationDirections.x;constunsignedcharY=m_quadricEquationDirections.y;//z=a+b.x+c.y+d.x^2+e.x.y+f.y^2//constPointCoordinateType&a=H[0];constPointCoordinateType&b=H[1];constPointCoordinateType&c=H[2];constPointCoordinateType&d=H[3];constPointCoordinateType&e=H[4];constPointCoordinateType&f=H[5];//See"CURVATUREOFCURVESANDSURFACES–APARABOLICAPPROACH"byZVIHAR’ELconstPointCoordinateTypefx=b+(d*2)*Q.u[X]+(e)*Q.u[Y];//b+2d*X+eYconstPointCoordinateTypefy=c+(e)*Q.u[X]+(f*2)*Q.u[Y];//c+2f*Y+eX

constPointCoordinateTypefxx=d*2;//2dconstPointCoordinateTypefyy=f*2;//2fconstPointCoordinateType&fxy=e;//econstPointCoordinateTypefx2=fx*fx;constPointCoordinateTypefy2=fy*fy;constPointCoordinateTypeq=(1+fx2+fy2);switch(cType){caseGAUSSIAN_CURV:{//tosignthecurvature,weneedanormal!PointCoordinateTypeK=fabs(fxx*fyy-fxy*fxy)/(q*q);returnstatic_cast(K);}caseMEAN_CURV:{//tosignthecurvature,weneedanormal!PointCoordinateTypeH=fabs(((1+fx2)*fyy-2*fx*fy*fxy+(1+fy2)*fxx))/(2*sqrt(q)*q);returnstatic_cast(H);}default:assert(false);}}break;caseNORMAL_CHANGE_RATE:{assert(m_associatedCloud);unsignedpointCount=(m_associatedCloud?m_associatedCloud->size():0);//weneedatleast4pointsif(pointCount<4){//notenoughpoints!returnpointCount==3?0:NAN_VALUE;

}//wedetermineplanenormalbycomputingthesmallesteigenvalueofM=1/n*S[(p-µ)*(p-µ)']CCLib::SquareMatrixdeigVectors;std::vectoreigValues;if(!Jacobi::ComputeEigenValuesAndVectors(computeCovarianceMatrix(),eigVectors,eigValues)){//failurereturnNAN_VALUE;}//computecurvatureastherateofchangeofthesurfacedoublee0=eigValues[0];doublee1=eigValues[1];doublee2=eigValues[2];doublesum=fabs(e0+e1+e2);if(sum<ZERO_TOLERANCE){returnNAN_VALUE;}doubleeMin=std::min(std::min(e0,e1),e2);returnstatic_cast(fabs(eMin)/sum);}break;default:assert(false);}returnNAN_VALUE;}Neighbourhood::computeCurvature3、计算点云密度4.计算点云表⾯粗糙度//"PER-CELL"METHOD:ROUGHNESSESTIMATION(LEASTSQUARESPLANEFIT)//ADDITIONNALPARAMETERS(1)://[0]->(PointCoordinateType*)kernelRadius:neighbourhoodradiusboolGeometricalAnalysisTools::computePointsRoughnessInACellAtLevel(constDgmOctree::octreeCell&cell,

void**additionalParameters,NormalizedProgress*nProgress/*=0*/){//parameter(s)PointCoordinateTyperadius=*static_cast(additionalParameters[0]);//structurefornearestneighborssearchDgmOctree::NearestNeighboursSphericalSearchStructnNSS;nNSS.level=cell.level;nNSS.prepare(radius,cell.parentOctree->getCellSize(nNSS.level));cell.parentOctree->getCellPos(cell.truncatedCode,cell.level,nNSS.cellPos,true);cell.parentOctree->computeCellCenter(nNSS.cellPos,cell.level,nNSS.cellCenter);unsignedn=cell.points->size();//numberofpointsinthecurrentcell//foreachpointinthecellfor(unsignedi=0;i{ScalarTyped=NAN_VALUE;cell.points->getPoint(i,nNSS.queryPoint);//lookforneighborsinsideasphere//warning:theremaybemorepointsattheendofnNSS.pointsInNeighbourhoodthantheactualnearestneighbors(=neighborCount)!unsignedneighborCount=cell.parentOctree->findNeighborsInASphereStartingFromCell(nNSS,radius,false);if(neighborCount>3){//findthequerypointinthenearestneighborssetandplaceitattheendconstunsignedglobalIndex=cell.points->getPointGlobalIndex(i);unsignedlocalIndex=0;while(localIndex<neighborCount&&nNSS.pointsInNeighbourhood[localIndex].pointIndex!=globalIndex)++localIndex;//thequerypointshouldbeinthenearestneighborsset!assert(localIndex<neighborCount);if(localIndex+1<neighborCount)//noneedtoswapwithanotherpointifit'salreadyattheend!{std::swap(nNSS.pointsInNeighbourhood[localIndex],nNSS.pointsInNeighbourhood[neighborCount-1]);}

DgmOctreeReferenceCloudneighboursCloud(&nNSS.pointsInNeighbourhood,neighborCount-1);//wedon'ttakethequerypointintoaccount!NeighbourhoodZ(&neighboursCloud);constPointCoordinateType*lsPlane=Z.getLSPlane();if(lsPlane)d=fabs(DistanceComputationTools::computePoint2PlaneDistance(&nNSS.queryPoint,lsPlane));//swapthepointsbacktotheiroriginalposition(DGM:notnecessary)//if(localIndex+1<neighborCount)//{//std::swap(nNSS.pointsInNeighbourhood[localIndex],nNSS.pointsInNeighbourhood[neighborCount-1]);//}}cell.points->setPointScalarValue(i,d);if(nProgress&&!nProgress->oneStep()){returnfalse;}}returntrue;}computePointsRoughnessInACellAtLevel地⾯粗糙度是指在⼀个特定的区域内，地球表⾯积与其投影⾯积之⽐。它也是反映地表形态的⼀个宏观指标。5.计算点云重⼼//计算重⼼CCVector3GeometricalAnalysisTools::computeGravityCenter(GenericCloud*theCloud){assert(theCloud);unsignedcount=theCloud->size();if(count==0)returnCCVector3();CCVector3dsum(0,0,0);theCloud->placeIteratorAtBegining();constCCVector3*P=0;while((P=theCloud->getNextPoint()))

{sum+=CCVector3d::fromArray(P->u);}sum/=static_cast(count);returnCCVector3::fromArray(sum.u);}computeGravityCenter6.计算点云权重重⼼//计算权重中⼼CCVector3GeometricalAnalysisTools::computeWeightedGravityCenter(GenericCloud*theCloud,ScalarField*weights){assert(theCloud&&weights);unsignedcount=theCloud->size();if(count==0||!weights||weights->currentSize()<count)returnCCVector3();CCVector3dsum(0,0,0);theCloud->placeIteratorAtBegining();doublewSum=0;for(unsignedi=0;i<count;++i){constCCVector3*P=theCloud->getNextPoint();ScalarTypew=weights->getValue(i);if(!ScalarField::ValidValue(w))continue;sum+=CCVector3d::fromArray(P->u)*fabs(w);wSum+=w;}if(wSum!=0)sum/=wSum;returnCCVector3::fromArray(sum.u);}computeWeightedGravityCenter7.计算点云协⽅差//计算协⽅差矩阵

CCLib::SquareMatrixdGeometricalAnalysisTools::computeCovarianceMatrix(GenericCloud*theCloud,constPointCoordinateType*_gravityCenter){assert(theCloud);unsignedn=(theCloud?theCloud->size():0);if(n==0)returnCCLib::SquareMatrixd();CCLib::SquareMatrixdcovMat(3);covMat.clear();//gravitycenterCCVector3G=(_gravityCenter?CCVector3(_gravityCenter):computeGravityCenter(theCloud));//crosssums(weusedoublestoavoidoverflow)doublemXX=0;doublemYY=0;doublemZZ=0;doublemXY=0;doublemXZ=0;doublemYZ=0;theCloud->placeIteratorAtBegining();for(unsignedi=0;i{constCCVector3*Q=theCloud->getNextPoint();CCVector3P=*Q-G;mXX+=static_cast(P.x*P.x);mYY+=static_cast(P.y*P.y);mZZ+=static_cast(P.z*P.z);mXY+=static_cast(P.x*P.y);mXZ+=static_cast(P.x*P.z);mYZ+=static_cast(P.y*P.z);}covMat.m_values[0][0]=mXX/static_cast(n);covMat.m_values[0][0]=mYY/static_cast(n);covMat.m_values[0][0]=mZZ/static_cast(n);covMat.m_values[1][0]=covMat.m_values[0][1]=mXY/static_cast(n);

covMat.m_values[2][0]=covMat.m_values[0][2]=mXZ/static_cast(n);covMat.m_values[2][1]=covMat.m_values[1][2]=mYZ/static_cast(n);returncovMat;}computeCovarianceMatrix8.计算点云互协⽅差//计算2个点云的互协⽅差CCLib::SquareMatrixdGeometricalAnalysisTools::computeCrossCovarianceMatrix(GenericCloud*P,GenericCloud*Q,constCCVector3&Gp,constCCVector3&Gq){assert(P&&Q);assert(Q->size()==P->size());//shortcutstooutputmatrixlinesCCLib::SquareMatrixdcovMat(3);double*l1=covMat.row(0);double*l2=covMat.row(1);double*l3=covMat.row(2);P->placeIteratorAtBegining();Q->placeIteratorAtBegining();//sumsunsignedcount=P->size();for(unsignedi=0;i{CCVector3Pt=*P->getNextPoint()-Gp;CCVector3Qt=*Q->getNextPoint()-Gq;l1[0]+=Pt.x*Qt.x;l1[1]+=Pt.x*Qt.y;l1[2]+=Pt.x*Qt.z;l2[0]+=Pt.y*Qt.x;l2[1]+=Pt.y*Qt.y;l2[2]+=Pt.y*Qt.z;l3[0]+=Pt.z*Qt.x;

l3[1]+=Pt.z*Qt.y;l3[2]+=Pt.z*Qt.z;}covMat.scale(1.0/static_cast(count));returncovMat;}computeCrossCovarianceMatrix//计算权重互协⽅差CCLib::SquareMatrixdGeometricalAnalysisTools::computeWeightedCrossCovarianceMatrix(GenericCloud*P,//dataGenericCloud*Q,//modelconstCCVector3&Gp,constCCVector3&Gq,ScalarField*coupleWeights/*=0*/){assert(P&&Q);assert(Q->size()==P->size());assert(coupleWeights);assert(coupleWeights->currentSize()==P->size());//shortcutstooutputmatrixlinesCCL