點雲的基本幾何計算

阿新 • • 發佈：2019-01-18

1.計算法向量

function normal = estimateNormal(data, tree, query, radius, min_neighbors)
% ESTIMATENORMAL Given a point cloud and query point, estimate the surface 
% normal by performing an eigendecomposition of the covariance matrix created 
% from the nearest neighbors of the query point for a fixed radius.
%
%  Example: 
%       data = randn(256,3);
%       tree = KDTreeSearcher(data);
%       query = [data(1,1) data(1,2) data(1,3)];
%       radius = 1.0;
%       min_neighbors = 8;
%       normal = estimateNormal(data, tree, query, radius, min_neighbors)
%
%  Copyright (c) 2014, William J. Beksi < 
[email protected]>
% 

% Find neighbors within the fixed radius 
idx = rangesearch(tree, query, radius);                                                   
idxs = idx{1};
neighbors = [data(idxs(:),1) data(idxs(:),2) data(idxs(:),3)];

if size(neighbors) < min_neighbors
    normal = {1};
    return;
end

% Compute the covariance matrix C
C = cov(neighbors);

% Compute the eigenvector of C
[v, lambda] = eig(C);

% Find the eigenvector corresponding to the minimum eigenvalue in C
[~, i] = min(diag(lambda));

% Normalize
normal = v(:,i) ./ norm(v(:,i));

end

2.計算曲率

　　曲線的曲率（curvature）就是針對曲線上某個點的切線方向角對弧長的轉動率，通過微分來定義，表明曲線偏離直線的程度。數學上表明曲線在某一點的彎曲程度的數值。曲率越大，表示曲線的彎曲程度越大。曲率的倒數就是曲率半徑。

　　平均曲率、高斯曲率、主曲率計算

　　賀美芳, et al. (2005.8). "散亂點雲資料的曲率估算及應用."

  1 ScalarType Neighbourhood::computeCurvature(unsigned neighbourIndex, CC_CURVATURE_TYPE cType)
  2 {
  3     switch 
 (cType)
  4     {
  5     case GAUSSIAN_CURV:
  6     case MEAN_CURV:
  7         {
  8             //we get 2D1/2 quadric parameters
  9             const PointCoordinateType* H = getQuadric();
 10             if (!H)
 11                 return NAN_VALUE;
 12 
 13             //compute centroid
 14             const CCVector3* G = getGravityCenter();
 15 
 16             //we compute curvature at the input neighbour position + we recenter it by the way
 17             CCVector3 Q = *m_associatedCloud->getPoint(neighbourIndex) - *G;
 18 
 19             const unsigned char X = m_quadricEquationDirections.x;
 20             const unsigned char Y = m_quadricEquationDirections.y;
 21 
 22             //z = a+b.x+c.y+d.x^2+e.x.y+f.y^2
 23             //const PointCoordinateType& a = H[0];
 24             const PointCoordinateType& b = H[1];
 25             const PointCoordinateType& c = H[2];
 26             const PointCoordinateType& d = H[3];
 27             const PointCoordinateType& e = H[4];
 28             const PointCoordinateType& f = H[5];
 29 
 30             //See "CURVATURE OF CURVES AND SURFACES – A PARABOLIC APPROACH" by ZVI HAR’EL
 31             const PointCoordinateType  fx    = b + (d*2) * Q.u[X] + (e  ) * Q.u[Y];    // b+2d*X+eY
 32             const PointCoordinateType  fy    = c + (e  ) * Q.u[X] + (f*2) * Q.u[Y];    // c+2f*Y+eX
 33             const PointCoordinateType  fxx    = d*2;                                    // 2d
 34             const PointCoordinateType  fyy    = f*2;                                    // 2f
 35             const PointCoordinateType& fxy    = e;                                    // e
 36 
 37             const PointCoordinateType fx2 = fx*fx;
 38             const PointCoordinateType fy2 = fy*fy;
 39             const PointCoordinateType q = (1 + fx2 + fy2);
 40 
 41             switch (cType)
 42             {
 43             case GAUSSIAN_CURV:
 44                 {
 45                     //to sign the curvature, we need a normal!
 46                     PointCoordinateType K = fabs( fxx*fyy - fxy*fxy ) / (q*q);
 47                     return static_cast<ScalarType>(K);
 48                 }
 49 
 50             case MEAN_CURV:
 51                 {
 52                     //to sign the curvature, we need a normal!
 53                     PointCoordinateType H = fabs( ((1+fx2)*fyy - 2*fx*fy*fxy + (1+fy2)*fxx) ) / (2*sqrt(q)*q);
 54                     return static_cast<ScalarType>(H);
 55                 }
 56 
 57             default:
 58                 assert(false);
 59             }
 60         }
 61         break;
 62 
 63     case NORMAL_CHANGE_RATE:
 64         {
 65             assert(m_associatedCloud);
 66             unsigned pointCount = (m_associatedCloud ? m_associatedCloud->size() : 0);
 67 
 68             //we need at least 4 points
 69             if (pointCount < 4)
 70             {
 71                 //not enough points!
 72                 return pointCount == 3 ? 0 : NAN_VALUE;
 73             }
 74 
 75             //we determine plane normal by computing the smallest eigen value of M = 1/n * S[(p-µ)*(p-µ)']
 76             CCLib::SquareMatrixd eigVectors;
 77             std::vector<double> eigValues;
 78             if (!Jacobi<double>::ComputeEigenValuesAndVectors(computeCovarianceMatrix(), eigVectors, eigValues))
 79             {
 80                 //failure
 81                 return NAN_VALUE;
 82             }
 83 
 84             //compute curvature as the rate of change of the surface
 85             double e0 = eigValues[0];
 86             double e1 = eigValues[1];
 87             double e2 = eigValues[2];
 88             double sum = fabs(e0+e1+e2);
 89             if (sum < ZERO_TOLERANCE)
 90             {
 91                 return NAN_VALUE;
 92             }
 93 
 94             double eMin = std::min(std::min(e0,e1),e2);
 95             return static_cast<ScalarType>(fabs(eMin) / sum);
 96         }
 97         break;
 98 
 99     default:
100         assert(false);
101     }
102 
103     return NAN_VALUE;
104 }

Neighbourhood::computeCurvature

3.計算點雲密度

4.計算點雲表面粗糙度

　　表面粗糙度

 1 //"PER-CELL" METHOD: ROUGHNESS ESTIMATION (LEAST SQUARES PLANE FIT)
 2 //ADDITIONNAL PARAMETERS (1):
 3 // [0] -> (PointCoordinateType*) kernelRadius : neighbourhood radius
 4 bool GeometricalAnalysisTools::computePointsRoughnessInACellAtLevel(const DgmOctree::octreeCell& cell, 
 5                                                                     void** additionalParameters,
 6                                                                     NormalizedProgress* nProgress/*=0*/)
 7 {
 8     //parameter(s)
 9     PointCoordinateType radius = *static_cast<PointCoordinateType*>(additionalParameters[0]);
10 
11     //structure for nearest neighbors search
12     DgmOctree::NearestNeighboursSphericalSearchStruct nNSS;
13     nNSS.level = cell.level;
14     nNSS.prepare(radius,cell.parentOctree->getCellSize(nNSS.level));
15     cell.parentOctree->getCellPos(cell.truncatedCode,cell.level,nNSS.cellPos,true);
16     cell.parentOctree->computeCellCenter(nNSS.cellPos,cell.level,nNSS.cellCenter);
17 
18     unsigned n = cell.points->size(); //number of points in the current cell
19     
20     //for each point in the cell
21     for (unsigned i=0; i<n; ++i)
22     {
23         ScalarType d = NAN_VALUE;
24         cell.points->getPoint(i,nNSS.queryPoint);
25 
26         //look for neighbors inside a sphere
27         //warning: there may be more points at the end of nNSS.pointsInNeighbourhood than the actual nearest neighbors (= neighborCount)!
28         unsigned neighborCount = cell.parentOctree->findNeighborsInASphereStartingFromCell(nNSS,radius,false);
29         if (neighborCount > 3)
30         {
31             //find the query point in the nearest neighbors set and place it at the end
32             const unsigned globalIndex = cell.points->getPointGlobalIndex(i);
33             unsigned localIndex = 0;
34             while (localIndex < neighborCount && nNSS.pointsInNeighbourhood[localIndex].pointIndex != globalIndex)
35                 ++localIndex;
36             //the query point should be in the nearest neighbors set!
37             assert(localIndex < neighborCount);
38             if (localIndex+1 < neighborCount) //no need to swap with another point if it's already at the end!
39             {
40                 std::swap(nNSS.pointsInNeighbourhood[localIndex],nNSS.pointsInNeighbourhood[neighborCount-1]);
41             }
42 
43             DgmOctreeReferenceCloud neighboursCloud(&nNSS.pointsInNeighbourhood,neighborCount-1); //we don't take the query point into account!
44             Neighbourhood Z(&neighboursCloud);
45 
46             const PointCoordinateType* lsPlane = Z.getLSPlane();
47             if (lsPlane)
48                 d = fabs(DistanceComputationTools::computePoint2PlaneDistance(&nNSS.queryPoint,lsPlane));
49 
50             //swap the points back to their original position (DGM: not necessary)
51             //if (localIndex+1 < neighborCount)
52             //{
53             //    std::swap(nNSS.pointsInNeighbourhood[localIndex],nNSS.pointsInNeighbourhood[neighborCount-1]);
54             //}
55         }
56 
57         cell.points->setPointScalarValue(i,d);
58 
59         if (nProgress && !nProgress->oneStep())
60         {
61             return false;
62         }
63     }
64 
65     return true;
66 }

computePointsRoughnessInACellAtLevel

　　地面粗糙度是指在一個特定的區域內，地球表面積與其投影面積之比。它也是反映地表形態的一個巨集觀指標。

5.計算點雲重心

 1 //計算重心
 2 CCVector3 GeometricalAnalysisTools::computeGravityCenter(GenericCloud* theCloud)
 3 {
 4     assert(theCloud);
 5     
 6     unsigned count = theCloud->size();
 7     if (count == 0)
 8         return CCVector3();
 9 
10     CCVector3d sum(0,0,0);
11 
12     theCloud->placeIteratorAtBegining();
13     const CCVector3 *P = 0;
14     while ((P = theCloud->getNextPoint()))
15     {
16         sum += CCVector3d::fromArray(P->u);
17     }
18 
19     sum /= static_cast<double>(count);
20     return CCVector3::fromArray(sum.u);
21 }

computeGravityCenter

6.計算點雲權重重心

 1 //計算權重中心
 2 CCVector3 GeometricalAnalysisTools::computeWeightedGravityCenter(GenericCloud* theCloud, ScalarField* weights)
 3 {
 4     assert(theCloud && weights);
 5 
 6     unsigned count = theCloud->size();
 7     if (count == 0 || !weights || weights->currentSize() < count)
 8         return CCVector3();
 9 
10     CCVector3d sum(0, 0, 0);
11 
12     theCloud->placeIteratorAtBegining();
13     double wSum = 0;
14     for (unsigned i = 0; i < count; ++i)
15     {
16         const CCVector3* P = theCloud->getNextPoint();
17         ScalarType w = weights->getValue(i);
18         if (!ScalarField::ValidValue(w))
19             continue;
20         sum += CCVector3d::fromArray(P->u) * fabs(w);
21         wSum += w;
22     }
23 
24     if (wSum != 0)
25         sum /= wSum;
26 
27     return CCVector3::fromArray(sum.u);
28 }

computeWeightedGravityCenter

7.計算點雲協方差

 1 //計算協方差矩陣
 2 CCLib::SquareMatrixd GeometricalAnalysisTools::computeCovarianceMatrix(GenericCloud* theCloud, const PointCoordinateType* _gravityCenter)
 3 {
 4     assert(theCloud);
 5     unsigned n = (theCloud ? theCloud->size() : 0);
 6     if (n==0)
 7         return CCLib::SquareMatrixd();
 8 
 9     CCLib::SquareMatrixd covMat(3);
10     covMat.clear();
11 
12     //gravity center
13     CCVector3 G = (_gravityCenter ?  CCVector3(_gravityCenter) : computeGravityCenter(theCloud));
14 
15     //cross sums (we use doubles to avoid overflow)
16     double mXX = 0;
17     double mYY = 0;
18     double mZZ = 0;
19     double mXY = 0;
20     double mXZ = 0;
21     double mYZ = 0;
22 
23     theCloud->placeIteratorAtBegining();
24     for (unsigned i=0;i<n;++i)
25     {
26         const CCVector3* Q = theCloud->getNextPoint();
27 
28         CCVector3 P = *Q-G;
29         mXX += static_cast<double>(P.x*P.x);
30         mYY += static_cast<double>(P.y*P.y);
31         mZZ += static_cast<double>(P.z*P.z);
32         mXY += static_cast<double>(P.x*P.y);
33         mXZ += static_cast<double>(P.x*P.z);
34         mYZ += static_cast<double>(P.y*P.z);
35     }
36 
37     covMat.m_values[0][0] = mXX/static_cast<double>(n);
38     covMat.m_values[0][0] = mYY/static_cast<double>(n);
39     covMat.m_values[0][0] = mZZ/static_cast<double>(n);
40     covMat.m_values[1][0] = covMat.m_values[0][1] = mXY/static_cast<double>(n);
41     covMat.m_values[2][0] = covMat.m_values[0][2] = mXZ/static_cast<double>(n);
42     covMat.m_values[2][1] = covMat.m_values[1][2] = mYZ/static_cast<double>(n);
43 
44     return covMat;
45 }

computeCovarianceMatrix

8.計算點雲互協方差

 1 //計算2個點雲的互協方差
 2 CCLib::SquareMatrixd GeometricalAnalysisTools::computeCrossCovarianceMatrix(GenericCloud* P,
 3                                                                             GenericCloud* Q,
 4                                                                             const CCVector3& Gp,
 5                                                                             const CCVector3& Gq)
 6 {
 7     assert(P && Q);
 8     assert(Q->size() == P->size());
 9 
10     //shortcuts to output matrix lines
11     CCLib::SquareMatrixd covMat(3);
12     double* l1 = covMat.row(0);
13     double* l2 = covMat.row(1);
14     double* l3 = covMat.row(2);
15 
16     P->placeIteratorAtBegining();
17     Q->placeIteratorAtBegining();
18 
19     //sums
20     unsigned count = P->size();
21     for (unsigned i=0; i<count; i++)
22     {
23         CCVector3 Pt = *P->getNextPoint() - Gp;
24         CCVector3 Qt = *Q->getNextPoint() - Gq;
25 
26         l1[0] += Pt.x*Qt.x;
27         l1[1] += Pt.x*Qt.y;
28         l1[2] += Pt.x*Qt.z;
29         l2[0] += Pt.y*Qt.x;
30         l2[1] += Pt.y*Qt.y;
31         l2[2] += Pt.y*Qt.z;
32         l3[0] += Pt.z*Qt.x;
33         l3[1] += Pt.z*Qt.y;
34         l3[2] += Pt.z*Qt.z;
35     }
36 
37     covMat.scale(1.0/static_cast<double>(count));
38 
39     return covMat;
40 }

computeCrossCovarianceMatrix

 1 //計算權重互協方差
 2 CCLib::SquareMatrixd GeometricalAnalysisTools::computeWeightedCrossCovarianceMatrix(GenericCloud* P, //data
 3                                                                                     GenericCloud* Q, //model
 4                                                                                     const CCVector3& Gp,
 5                                                                                     const CCVector3& Gq,
 6                                                                                     ScalarField* coupleWeights/*=0*/)
 7 {
 8     assert(P && Q);
 9     assert(Q->size() == P->size());
10     assert(coupleWeights);
11     assert(coupleWeights->currentSize() == P->size());
12 
13     //shortcuts to output matrix lines
14     CCLib::SquareMatrixd covMat(3);
15     double* r1 = covMat.row(0);
16     double* r2 = covMat.row(1);
17     double* r3 = covMat.row(2);
18 
19     P->placeIteratorAtBegining();
20     Q->placeIteratorAtBegining();
21 
22     //sums
23     unsigned count = P->size();
24     double wSum = 0.0; //we will normalize by the sum
25     for (unsigned i = 0; i<count; i++)
26     {
27         CCVector3d Pt = CCVector3d::fromArray((*P->getNextPoint() - Gp).u);
28         CCVector3 Qt = *Q->getNextPoint() - Gq;
29 
30         //Weighting scheme for cross-covariance is inspired from
31         //https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Weighted_sample_covariance
32         double wi = 1.0;
33         if (coupleWeights)
34         {
35             ScalarType w = coupleWeights->getValue(i);
36             if (!ScalarField::ValidValue(w))
37                 continue;
38             wi = fabs(w);
39         }
40 
41         //DGM: we virtually make the P (data) point nearer if it has a lower weight
42         Pt *= wi;
43         wSum += wi;
44 
45         //1st row
46         r1[0] += Pt.x * Qt.x;
47         r1[1] += Pt.x * Qt.y;
48         r1[2] += Pt.x * Qt.z;
49         //2nd row
50         r2[0] += Pt.y * Qt.x;
51         r2[1] += Pt.y * Qt.y;
52         r2[2] += Pt.y * Qt.z;
53         //3rd row
54         r3[0] += Pt.z * Qt.x;
55         r3[1] += Pt.z * Qt.y;
56         r3[2] += Pt.z * Qt.z;
57     }
58 
59     if (wSum != 0.0)
60         covMat.scale(1.0/wSum);
61 
62     return covMat;
63 }

computeWeightedCrossCovarianceMatrix

點雲的基本幾何計算

1.計算法向量原檔案 function normal = estimateNormal(data, tree, query, radius, min_neighbors) % ESTIMATENORMAL Given a point cloud and query point, estimate

[PCL]2 點雲法向量計算NormalEstimation

參考：http://www.cnblogs.com/yhlx125/p/5137850.html 從GitHub的程式碼版本庫下載原始碼https://github.com/PointCloudLibrary/pcl，用CMake生成VS專案，檢視PCL的原始碼位於pcl_featu

點雲模型的局部曲面幾何資訊的提取（點雲法向量的計算）

估算法向量法向量是三維點雲資料具有的一個很重要的區域性特性，在點雲的許多處理中是必不可少的資訊。三維掃描獲取的初始取樣點集只記錄了各取樣點的空間三維座標，而不存在任何連線關係，求解法向量是處理點雲資料的第一步。在三維軟體中求解取樣點的法向量，一般要求建立點雲資料對應的網路模型

點雲/網格模型的體積計算

pro nis top 方便 vol ref org splay 會有點雲體積計算　　有時用激光掃描設備掃描零件或者用無人機進行測量後會想知道它們的體積。比如下面的土堆：　　如果掃描得到的數據是一系列三維點雲，那麽體積就比較難求，因為如何定義物體的邊界比較

由深度圖計算點雲的原理

本文整理了兩種計算點雲資料的方法，其核心都是根據針孔成像原理而來。一、ROS中的方法首先，要了解下世界座標到影象的對映過程，考慮世界座標點M(Xw,Yw,Zw)對映到影象點m(u,v)的過程，如下圖所示：參考針孔成像原理其中u,v為影象座標系下的任意座標點。

極簡多檢視幾何（3）基於基本矩陣計算攝像機矩陣

MVG: p256 假設F為基本矩陣，S為任意反對稱矩陣，則可以將一對攝像機矩陣定義為：其中e’為極點，因為只有在如下條件滿足時，P’的秩為3，因此S可以表示為極點的反對稱矩陣的形式進而得到如下求解公式：

點雲密度計算

double computeCloudResolution (const pcl::PointCloud<PointType>::ConstPtr &cloud) { double res = 0.0; int n_points = 0; int nres;

點雲的平滑與法線計算

需要平滑的情況： 1、用RGB-D鐳射掃描器等裝置掃描物體，尤其是比較小的物體時，往往會有測量誤差。這些誤差所造成的不規則資料如果直接拿來曲面重建的話，會使得重建的曲面不光滑或者有漏洞，而且這種不規則資料很難用前面我們提到過的統計分析等濾波方法消除； 2、後處理過程中，對同一個物體從不同方向

[CC]點雲密度計算

包括兩種計算方法：精確計算和近似計算（思考：local density=單位面積的點數 vs local density =1/單個點所佔的面積）每種方法可以實現三種模式的點雲密度計算，CC裡面的點雲端計算依賴於給定的近鄰半徑所對應的最佳八叉樹層級（通過findBestLevelForAGiv

計算點雲的最小BBOX

參考1 參考2 #include <iostream> #include <pcl/io/pcd_io.h> #include <pcl/point_types.h> #include <pcl/visualizat

VTK修煉之道57：圖形基本操作進階_點雲配準技術（LandMark標記點演算法和座標系顯示方法）

#include <vtkAutoInit.h> VTK_MODULE_INIT(vtkRenderingOpenGL); VTK_MODULE_INIT(vtkRenderingFreeType); VTK_MODULE_INIT(vtkInteractionStyle); #include

PCL中計算點雲的法向量並顯示

// NormalEstimation.cpp : Defines the entry point for the console application. // #include "stdafx.h" #include <pcl/point_t

orientation_points_xld (Operator)----計算點雲輪廓或點雲多邊形的方向

計算點雲輪廓或點雲多邊形的方向（該點雲的順序不會考慮其中）①如果該XLD與自身相交或者如果自身相交能用一條線把其首尾相連，使用orientation_points_xld 運算元，可以使用 test_self_intersection_xld運算元測試該XLD能否自身相交②如

通過Kinect的深度影象資料計算三維點雲

在可以透過 OpenNI 讀取到 Kinect 的深度、色彩資訊之後，其實就可以試著用這些資訊，來重建 3D 的環境做顯示了～不過實際上，在前面的範例中所讀到的深度資訊，都算是原始資料，而且座標軸也都是感應器二維影像的座標系統，如果要重建 3D 場景的話，這些資訊都還是需要換算的；所幸，OpenNI 在 D

點分治基本膜版

!= 方式 turn n) space ring 命運 code cstring 點分治基本的幾個步驟吧找根（保證平衡度最優）從找到的根處理起，每次刪去已有的根，分治其每一個子樹非暴力算法解決樹形問題的有效方式 1 #include<

百度地圖API一：根據標註點坐標範圍計算顯示縮放級別zoom自適應顯示地圖

var spa get bsp pan nts viewport 百度 getview 百度地圖中根據頁面中的point，自動設置縮放級別和視圖中心，將所有的point在視圖範圍內展示。 var points = [point1, point2,point3];

1.ArcGis幾何圖形之幾何計算

nal contain pre true 包含 geo summary relation sta 1 /// <summary> 2 /// 檢測幾何圖形A是否包含幾何圖形B 3 /// </s

VisualSFM+PMVS生成稠密點雲

我的電腦調用 span visual bit .cn 電腦菜單 32bit 利用相機拍攝一個場景不同角度的圖片，使用VisualSFM能夠得到稀疏點雲，如果想要得到稠密點雲，可以在VisualSFM中加入PMVS的應用程序，PMVS會作為一個插件運行將稀疏點雲插成稠密的

數據埋點最基本需要獲取的用戶數據有哪些？

根據可能 -a 統計代碼自己一個著作權 mage 由於數據埋點最基本需要獲取的用戶數據有哪些？原文出處https://www.zhihu.com/question/20412632 作者：劉大大鏈接：https://www.zhihu.com/question/

PCL—點雲分割（最小割算法）

number 作用早就有效好的介紹不同的優勢 bsp 1.點雲分割的精度　　在之前的兩個章節裏介紹了基於采樣一致的點雲分割和基於臨近搜索的點雲分割算法。基於采樣一致的點雲分割算法顯然是意識流的，它只能割出大概的點雲（可能是杯子的一部分，但杯把兒肯定沒分割出來）

點雲的基本幾何計算

相關推薦