簡介

上篇博客介紹了特征提取的原理以及坐標轉化和插值的源碼理解，這篇將介紹特征提取的后續算法模塊。

源碼

1、新數據進來進行坐標轉換和插補等工作

見【LeGO-LOAM論文閱讀（二）–特征提取（一）】

2、進行光滑度計算

光滑度計算公式：

S取10，代表參與計算點數，計算點左右各有5個，r表示深度距離，計算點云開始5個點后到點云結束5個點前的點云平滑度，計算方式與論文中給的并不一樣。

// 計算光滑性，這里的計算沒有完全按照公式進行，// 缺少除以總點數i和r[i]void calculateSmoothness(){int cloudSize = segmentedCloud->points.size();for (int i = 5; i < cloudSize - 5; i++) {float diffRange = segInfo.segmentedCloudRange[i-5] + segInfo.segmentedCloudRange[i-4]+ segInfo.segmentedCloudRange[i-3] + segInfo.segmentedCloudRange[i-2]+ segInfo.segmentedCloudRange[i-1] - segInfo.segmentedCloudRange[i] * 10+ segInfo.segmentedCloudRange[i+1] + segInfo.segmentedCloudRange[i+2]+ segInfo.segmentedCloudRange[i+3] + segInfo.segmentedCloudRange[i+4]+ segInfo.segmentedCloudRange[i+5]; cloudCurvature[i] = diffRange*diffRange;

初始化點云標簽，初始化為0，surfPointsFlat標記為-1，surfPointsLessFlatScan為不大于0的標簽，cornerPointsSharp標記為2，cornerPointsLessSharp標記為1。

// 在markOccludedPoints()函數中對該參數進行重新修改cloudNeighborPicked[i] = 0;// 在extractFeatures()函數中會對標簽進行修改，// 初始化為0，surfPointsFlat標記為-1，surfPointsLessFlatScan為不大于0的標簽// cornerPointsSharp標記為2，cornerPointsLessSharp標記為1cloudLabel[i] = 0;cloudSmoothness[i].value = cloudCurvature[i];cloudSmoothness[i].ind = i;}}

3、標記阻塞點。

阻塞點指點云之間相互遮擋而且距離比較近的點。將被遮擋的點云連續五個標志位cloudNeighborPicked標記為1。depth大的那個就是被遮擋的那個。具體參數設置為什么是10和0.3我也不清楚，應該是是試驗過的，各位可以考慮話說弄成別的試試

void markOccludedPoints(){int cloudSize = segmentedCloud->points.size();for (int i = 5; i < cloudSize - 6; ++i){float depth1 = segInfo.segmentedCloudRange[i];float depth2 = segInfo.segmentedCloudRange[i+1];int columnDiff = std::abs(int(segInfo.segmentedCloudColInd[i+1] - segInfo.segmentedCloudColInd[i]));if (columnDiff < 10){// 選擇距離較遠的那些點，并將他們標記為1if (depth1 - depth2 > 0.3){cloudNeighborPicked[i - 5] = 1;cloudNeighborPicked[i - 4] = 1;cloudNeighborPicked[i - 3] = 1;cloudNeighborPicked[i - 2] = 1;cloudNeighborPicked[i - 1] = 1;cloudNeighborPicked[i] = 1;}else if (depth2 - depth1 > 0.3){cloudNeighborPicked[i + 1] = 1;cloudNeighborPicked[i + 2] = 1;cloudNeighborPicked[i + 3] = 1;cloudNeighborPicked[i + 4] = 1;cloudNeighborPicked[i + 5] = 1;cloudNeighborPicked[i + 6] = 1;}}

接著排除共面的情況，舍棄距離變換比較大的點

float diff1 = std::abs(segInfo.segmentedCloudRange[i-1] - segInfo.segmentedCloudRange[i]);float diff2 = std::abs(segInfo.segmentedCloudRange[i+1] - segInfo.segmentedCloudRange[i]);// 選擇距離變化較大的點，并將他們標記為1if (diff1 > 0.02 * segInfo.segmentedCloudRange[i] && diff2 > 0.02 * segInfo.segmentedCloudRange[i])cloudNeighborPicked[i] = 1;}}

4、特征抽取。

這部分對計算好平滑度的點云將進行特征提取，提取出線特征和面特征。由于點云過于龐大，后面處理的時候對點云進行了下采樣來提高計算效率。
清空上一幀特征：

void extractFeatures(){cornerPointsSharp->clear();cornerPointsLessSharp->clear();surfPointsFlat->clear();surfPointsLessFlat->clear();

將每一幀點云數據按水平角度分為6份，每份是60°的激光數據，sp儲存的是每份數據的起點，ep儲存的是每份數據的終點。

for (int i = 0; i < N_SCAN; i++) {surfPointsLessFlatScan->clear();for (int j = 0; j < 6; j++) {// sp和ep的含義是什么???startPointer,endPointer?int sp = (segInfo.startRingIndex[i] * (6 - j) + segInfo.endRingIndex[i] * j) / 6;int ep = (segInfo.startRingIndex[i] * (5 - j) + segInfo.endRingIndex[i] * (j + 1)) / 6 - 1;if (sp >= ep)continue;

將每份點云按照平滑度進行排序（由小到大）：

std::sort(cloudSmoothness.begin()+sp, cloudSmoothness.begin()+ep, by_value());

提取線特征的條件：不是被淘汰的點（cloudNeighborPicked=0）并且平滑度大于設定的閾值還有不能是地面點：

int largestPickedNum = 0;for (int k = ep; k >= sp; k--) {// 每次ind的值就是等于k??? 有什么意義?// 因為上面對cloudSmoothness進行了一次從小到大排序，所以ind不一定等于k了int ind = cloudSmoothness[k].ind;if (cloudNeighborPicked[ind] == 0 &&cloudCurvature[ind] > edgeThreshold &&segInfo.segmentedCloudGroundFlag[ind] == false) {

在滿足條件的點云中，將最大平滑度的兩個點云加入到cornerPointsSharp中，并將標志cloudLabel置為2

largestPickedNum++;if (largestPickedNum <= 2) {// 論文中nFe=2,cloudSmoothness已經按照從小到大的順序排列，// 所以這邊只要選擇最后兩個放進隊列即可// cornerPointsSharp標記為2cloudLabel[ind] = 2;cornerPointsSharp->push_back(segmentedCloud->points[ind]);cornerPointsLessSharp->push_back(segmentedCloud->points[ind]);

將平滑度前20的點云加入到cornerPointsLessSharp中，并將cloudLabel置為1

} else if (largestPickedNum <= 20) {// 塞20個點到cornerPointsLessSharp中去// cornerPointsLessSharp標記為1cloudLabel[ind] = 1;cornerPointsLessSharp->push_back(segmentedCloud->points[ind]);} else {break;}

將已經加入到cornerPointsSharp和cornerPointsLessSharp中的點淘汰，后面不再用：

cloudNeighborPicked[ind] = 1;

計算已經被選為線特征點云前后5個點之間的距離差值，若距離較近
同樣淘汰（避免兩個特征點表示的是同一個線特征）：

for (int l = 1; l <= 5; l++) {// 從ind+l開始后面5個點，每個點index之間的差值，// 確保columnDiff<=10,然后標記為我們需要的點int columnDiff = std::abs(int(segInfo.segmentedCloudColInd[ind + l] - segInfo.segmentedCloudColInd[ind + l - 1]));if (columnDiff > 10)break;cloudNeighborPicked[ind + l] = 1;}for (int l = -1; l >= -5; l--) {// 從ind+l開始前面五個點，計算差值然后標記int columnDiff = std::abs(int(segInfo.segmentedCloudColInd[ind + l] - segInfo.segmentedCloudColInd[ind + l + 1]));if (columnDiff > 10)break;cloudNeighborPicked[ind + l] = 1;}}}

提取面特征的條件，不能是淘汰點，平滑度小于設定閾值并且只能是地面點：

int smallestPickedNum = 0;for (int k = sp; k <= ep; k++) {int ind = cloudSmoothness[k].ind;// 平面點只從地面點中進行選擇? 為什么要這樣做?if (cloudNeighborPicked[ind] == 0 &&cloudCurvature[ind] < surfThreshold &&segInfo.segmentedCloudGroundFlag[ind] == true) {

將滿足條件的點加入surfPointsFlat同時標志cloudLabel設為-1：

cloudLabel[ind] = -1;surfPointsFlat->push_back(segmentedCloud->points[ind]);

將最小平滑度的四個點放入surfPointsFlat，同時設置為淘汰點。

// 論文中nFp=4，將4個最平的平面點放入隊列中smallestPickedNum++;if (smallestPickedNum >= 4) {break;}cloudNeighborPicked[ind] = 1;

同時判斷面特征點前后5個點并標記（和線特征一樣）：

for (int l = 1; l <= 5; l++) {// 從前面往后判斷是否是需要的鄰接點，是的話就進行標記int columnDiff = std::abs(int(segInfo.segmentedCloudColInd[ind + l] - segInfo.segmentedCloudColInd[ind + l - 1]));if (columnDiff > 10)break;cloudNeighborPicked[ind + l] = 1;}for (int l = -1; l >= -5; l--) {// 從后往前開始標記int columnDiff = std::abs(int(segInfo.segmentedCloudColInd[ind + l] - segInfo.segmentedCloudColInd[ind + l + 1]));if (columnDiff > 10)break;cloudNeighborPicked[ind + l] = 1;}}}

將cloudLabe<= 0的點先都加入 surfPointsLessFlatScan中：

for (int k = sp; k <= ep; k++) {if (cloudLabel[k] <= 0) {surfPointsLessFlatScan->push_back(segmentedCloud->points[k]);}}}

surfPointsLessFlatScan中的點云會非常的多，為了提高計算效率，對點云進行下采樣：

// surfPointsLessFlatScan中有過多的點云，如果點云太多，計算量太大// 進行下采樣，可以大大減少計算量surfPointsLessFlatScanDS->clear();downSizeFilter.setInputCloud(surfPointsLessFlatScan);downSizeFilter.filter(*surfPointsLessFlatScanDS);*surfPointsLessFlat += *surfPointsLessFlatScanDS;}}

5、發布四種點云信息。

分別發布平滑度大的點云、平滑度較大的點云、平滑度小的點云、平滑度較小的點云：

void publishCloud(){sensor_msgs::PointCloud2 laserCloudOutMsg;if (pubCornerPointsSharp.getNumSubscribers() != 0){pcl::toROSMsg(*cornerPointsSharp, laserCloudOutMsg);laserCloudOutMsg.header.stamp = cloudHeader.stamp;laserCloudOutMsg.header.frame_id = "/camera";pubCornerPointsSharp.publish(laserCloudOutMsg);}if (pubCornerPointsLessSharp.getNumSubscribers() != 0){pcl::toROSMsg(*cornerPointsLessSharp, laserCloudOutMsg);laserCloudOutMsg.header.stamp = cloudHeader.stamp;laserCloudOutMsg.header.frame_id = "/camera";pubCornerPointsLessSharp.publish(laserCloudOutMsg);}if (pubSurfPointsFlat.getNumSubscribers() != 0){pcl::toROSMsg(*surfPointsFlat, laserCloudOutMsg);laserCloudOutMsg.header.stamp = cloudHeader.stamp;laserCloudOutMsg.header.frame_id = "/camera";pubSurfPointsFlat.publish(laserCloudOutMsg);}if (pubSurfPointsLessFlat.getNumSubscribers() != 0){pcl::toROSMsg(*surfPointsLessFlat, laserCloudOutMsg);laserCloudOutMsg.header.stamp = cloudHeader.stamp;laserCloudOutMsg.header.frame_id = "/camera";pubSurfPointsLessFlat.publish(laserCloudOutMsg);}}

回到runFeatureAssociation函數里，判斷有沒有初始化·，即判斷是不是第一幀數據，systemInitedLM初始為0，：

if (!systemInitedLM) {checkSystemInitialization();return;}

將第一幀的曲率較大的點保存為上一幀用于匹配的線特征：

void checkSystemInitialization(){// 交換cornerPointsLessSharp和laserCloudCornerLastpcl::PointCloud<PointType>::Ptr laserCloudTemp = cornerPointsLessSharp;cornerPointsLessSharp = laserCloudCornerLast;laserCloudCornerLast = laserCloudTemp;

將第一幀的曲率較小的點保存為上一幀用于匹配的面特征

// 交換surfPointsLessFlat和laserCloudSurfLastlaserCloudTemp = surfPointsLessFlat;surfPointsLessFlat = laserCloudSurfLast;laserCloudSurfLast = laserCloudTemp;

將特征點保存到KD樹中：

kdtreeCornerLast->setInputCloud(laserCloudCornerLast);kdtreeSurfLast->setInputCloud(laserCloudSurfLast);

記錄特征點數

laserCloudCornerLastNum = laserCloudCornerLast->points.size();laserCloudSurfLastNum = laserCloudSurfLast->points.size();

發布特征點云消息：

sensor_msgs::PointCloud2 laserCloudCornerLast2;pcl::toROSMsg(*laserCloudCornerLast, laserCloudCornerLast2);laserCloudCornerLast2.header.stamp = cloudHeader.stamp;laserCloudCornerLast2.header.frame_id = "/camera";pubLaserCloudCornerLast.publish(laserCloudCornerLast2);sensor_msgs::PointCloud2 laserCloudSurfLast2;pcl::toROSMsg(*laserCloudSurfLast, laserCloudSurfLast2);laserCloudSurfLast2.header.stamp = cloudHeader.stamp;laserCloudSurfLast2.header.frame_id = "/camera";pubLaserCloudSurfLast.publish(laserCloudSurfLast2);

保存累積旋轉量：

transformSum[0] += imuPitchStart;transformSum[2] += imuRollStart;

checkSystemInitialization判斷完第一幀數據將systemInitedLM設為1

systemInitedLM = true; }

6、更新位姿

這里主要是保存當前點云中最后一個點的旋轉角、最后一個點相對于第一個點的位移以及速度：

void updateInitialGuess(){imuPitchLast = imuPitchCur;imuYawLast = imuYawCur;imuRollLast = imuRollCur;imuShiftFromStartX = imuShiftFromStartXCur;imuShiftFromStartY = imuShiftFromStartYCur;imuShiftFromStartZ = imuShiftFromStartZCur;imuVeloFromStartX = imuVeloFromStartXCur;imuVeloFromStartY = imuVeloFromStartYCur;imuVeloFromStartZ = imuVeloFromStartZCur;

但是我發現最后一個點的位移變換函數根本沒有被調用，我覺得這個函數是要在插值和坐標轉那被調用，可以并沒有，有沒有大佬能解釋下。

void ShiftToStartIMU(float pointTime){// 下面三個量表示的是世界坐標系下，從start到cur的坐標的漂移imuShiftFromStartXCur = imuShiftXCur - imuShiftXStart - imuVeloXStart * pointTime;imuShiftFromStartYCur = imuShiftYCur - imuShiftYStart - imuVeloYStart * pointTime;imuShiftFromStartZCur = imuShiftZCur - imuShiftZStart - imuVeloZStart * pointTime;// 從世界坐標系變換到start坐標系float x1 = cosImuYawStart * imuShiftFromStartXCur - sinImuYawStart * imuShiftFromStartZCur;float y1 = imuShiftFromStartYCur;float z1 = sinImuYawStart * imuShiftFromStartXCur + cosImuYawStart * imuShiftFromStartZCur;float x2 = x1;float y2 = cosImuPitchStart * y1 + sinImuPitchStart * z1;float z2 = -sinImuPitchStart * y1 + cosImuPitchStart * z1;imuShiftFromStartXCur = cosImuRollStart * x2 + sinImuRollStart * y2;imuShiftFromStartYCur = -sinImuRollStart * x2 + cosImuRollStart * y2;imuShiftFromStartZCur = z2;}

更新旋轉量（相鄰兩幀）：

// 關于下面負號的說明：// transformCur是在Cur坐標系下的 p_start=R*p_cur+t// R和t是在Cur坐標系下的// 而imuAngularFromStart是在start坐標系下的，所以需要加負號if (imuAngularFromStartX != 0 || imuAngularFromStartY != 0 || imuAngularFromStartZ != 0){transformCur[0] = - imuAngularFromStartY;transformCur[1] = - imuAngularFromStartZ;transformCur[2] = - imuAngularFromStartX;

更新位移量（兩幀之間）：

// 速度乘以時間，當前變換中的位移if (imuVeloFromStartX != 0 || imuVeloFromStartY != 0 || imuVeloFromStartZ != 0){transformCur[3] -= imuVeloFromStartX * scanPeriod;transformCur[4] -= imuVeloFromStartY * scanPeriod;transformCur[5] -= imuVeloFromStartZ * scanPeriod;}}

7、更新變換矩陣。

如果線特征或者面特征太少，直接退出：

void updateTransformation(){if (laserCloudCornerLastNum < 10 || laserCloudSurfLastNum < 100)return;

尋找面特征對應點findCorrespondingSurfFeatures， iterCount迭代次數？：

for (int iterCount1 = 0; iterCount1 < 25; iterCount1++) {laserCloudOri->clear();coeffSel->clear();// 找到對應的特征平面// 然后計算協方差矩陣，保存在coeffSel隊列中// laserCloudOri中保存的是對應于coeffSel的未轉換到開始時刻的原始點云數據findCorrespondingSurfFeatures(iterCount1);

面特征點數賦值：

void findCorrespondingSurfFeatures(int iterCount){int surfPointsFlatNum = surfPointsFlat->points.size();

將特征點坐標變換到開始時刻 TransformToStart：

for (int i = 0; i < surfPointsFlatNum; i++) {// 坐標變換到開始時刻，參數0是輸入，參數1是輸出TransformToStart(&surfPointsFlat->points[i], &pointSel);

注意假設了勻速運動模型，也就是說位姿的變化是線性的，而代碼中提到的s就是隨時間線性變化的因子。回顧之前點云強度的保存的數據的計算，point.intensity = int(segmentedCloud->points[i].intensity) + scanPeriod * relTime;，再結合s的定義， float s = 10 * (pi->intensity - int(pi->intensity));，發現s就是relTime，為當前點云時間減去初始點云相對于幀點云時間差的比值。

void TransformToStart(PointType const * const pi, PointType * const po){// intensity代表的是：整數部分ring序號，小數部分是當前點在這一圈中所花的時間// 關于intensity， 參考 adjustDistortion() 函數中的定義// s代表的其實是一個比例，s的計算方法應該如下：// s=(pi->intensity - int(pi->intensity))/SCAN_PERIOD// ===> SCAN_PERIOD=0.1(雷達頻率為10hz)// 以上理解感謝github用戶StefanGlaser// https://github.com/laboshinl/loam_velodyne/issues/29float s = 10 * (pi->intensity - int(pi->intensity));

計算當前點云和開始時刻的旋轉角度以及位移，這里用了勻速運動模型：

float rx = s * transformCur[0];float ry = s * transformCur[1];float rz = s * transformCur[2];float tx = s * transformCur[3];float ty = s * transformCur[4];float tz = s * transformCur[5];

計算特征點對應于開始時刻坐標系的坐標，先平移在旋轉：

float x1 = cos(rz) * (pi->x - tx) + sin(rz) * (pi->y - ty);float y1 = -sin(rz) * (pi->x - tx) + cos(rz) * (pi->y - ty);float z1 = (pi->z - tz);float x2 = x1;float y2 = cos(rx) * y1 + sin(rx) * z1;float z2 = -sin(rx) * y1 + cos(rx) * z1;po->x = cos(ry) * x2 - sin(ry) * z2;po->y = y2;po->z = sin(ry) * x2 + cos(ry) * z2;po->intensity = pi->intensity;}

每5此迭代重新匹配特征點：

if (iterCount % 5 == 0) {

在上一幀曲率較小的點（平面特征點）組成的KD樹中尋找一個最近鄰點：

kdtreeSurfLast->nearestKSearch(pointSel, 1, pointSearchInd, pointSearchSqDis);int closestPointInd = -1, minPointInd2 = -1, minPointInd3 = -1;// sq:平方，距離的平方值// 如果nearestKSearch找到的1(k=1)個鄰近點滿足條件if (pointSearchSqDis[0] < nearestFeatureSearchSqDist) {closestPointInd = pointSearchInd[0];

找到該最近鄰點的線束（第幾行第幾列）point.intensity保存里點云的行列值：

// point.intensity 保存的是什么值? 第幾次scan?// thisPoint.intensity = (float)rowIdn + (float)columnIdn / 10000.0;// fullInfoCloud->points[index].intensity = range;// 在imageProjection.cpp文件中有上述兩行代碼，兩種類型的值，應該存的是上面一個int closestPointScan = int(laserCloudSurfLast->points[closestPointInd].intensity);

在最近鄰點i 所在線束尋找再尋找另外一個最近點，以及在最近鄰點i ii所在線束的相鄰線束中尋找第三個最近點，組成匹配平面：

// 主要功能是找到2個scan之內的最近點，并將找到的最近點及其序號保存// 之前掃描的保存到minPointSqDis2，之后的保存到minPointSqDis2float pointSqDis, minPointSqDis2 = nearestFeatureSearchSqDist, minPointSqDis3 = nearestFeatureSearchSqDist;

首先在雷達點云存儲的正方向(ScanID增大的方向和雷達旋轉正方向)尋找其他兩個點，為什么加上2.5，因為上下兩個scan相差2°：

for (int j = closestPointInd + 1; j < surfPointsFlatNum; j++) {// int類型的值加上2.5? 為什么不直接加上2?// 四舍五入if (int(laserCloudSurfLast->points[j].intensity) > closestPointScan + 2.5) {break;}

計算距離，找到最近鄰的點，將之前掃描到的點保存在minPointSqDis2中，之后掃描到的保存在minPointSqDis3中：

pointSqDis = (laserCloudSurfLast->points[j].x - pointSel.x) * (laserCloudSurfLast->points[j].x - pointSel.x) + (laserCloudSurfLast->points[j].y - pointSel.y) * (laserCloudSurfLast->points[j].y - pointSel.y) + (laserCloudSurfLast->points[j].z - pointSel.z) * (laserCloudSurfLast->points[j].z - pointSel.z);if (int(laserCloudSurfLast->points[j].intensity) <= closestPointScan) {if (pointSqDis < minPointSqDis2) {minPointSqDis2 = pointSqDis;minPointInd2 = j;}} else {if (pointSqDis < minPointSqDis3) {minPointSqDis3 = pointSqDis;minPointInd3 = j;}}

反方向(ScanID減小的方向和雷達旋轉反方向尋找點，并確定總的最近鄰點：

// 往前找for (int j = closestPointInd - 1; j >= 0; j--) {if (int(laserCloudSurfLast->points[j].intensity) < closestPointScan - 2.5) {break;}pointSqDis = (laserCloudSurfLast->points[j].x - pointSel.x) * (laserCloudSurfLast->points[j].x - pointSel.x) + (laserCloudSurfLast->points[j].y - pointSel.y) * (laserCloudSurfLast->points[j].y - pointSel.y) + (laserCloudSurfLast->points[j].z - pointSel.z) * (laserCloudSurfLast->points[j].z - pointSel.z);if (int(laserCloudSurfLast->points[j].intensity) >= closestPointScan) {if (pointSqDis < minPointSqDis2) {minPointSqDis2 = pointSqDis;minPointInd2 = j;}} else {if (pointSqDis < minPointSqDis3) {minPointSqDis3 = pointSqDis;minPointInd3 = j;}}}}

找到平面匹配的三個點：

pointSearchSurfInd1[i] = closestPointInd;pointSearchSurfInd2[i] = minPointInd2;pointSearchSurfInd3[i] = minPointInd3;}

找到特征點的對應點以后的工作就是計算距離了，在這里還計算了平面法向量（就是點到平面的直線的方向向量）在各個軸的方向分解，在雅可比的計算中使用：

// 前后都能找到對應的最近點在給定范圍之內// 那么就開始計算距離// [pa,pb,pc]是tripod1，tripod2，tripod3這3個點構成的一個平面的方向量，// ps是模長，它是三角形面積的2倍if (pointSearchSurfInd2[i] >= 0 && pointSearchSurfInd3[i] >= 0) {tripod1 = laserCloudSurfLast->points[pointSearchSurfInd1[i]];tripod2 = laserCloudSurfLast->points[pointSearchSurfInd2[i]];tripod3 = laserCloudSurfLast->points[pointSearchSurfInd3[i]];float pa = (tripod2.y - tripod1.y) * (tripod3.z - tripod1.z) - (tripod3.y - tripod1.y) * (tripod2.z - tripod1.z);float pb = (tripod2.z - tripod1.z) * (tripod3.x - tripod1.x) - (tripod3.z - tripod1.z) * (tripod2.x - tripod1.x);float pc = (tripod2.x - tripod1.x) * (tripod3.y - tripod1.y) - (tripod3.x - tripod1.x) * (tripod2.y - tripod1.y);float pd = -(pa * tripod1.x + pb * tripod1.y + pc * tripod1.z);float ps = sqrt(pa * pa + pb * pb + pc * pc);pa /= ps;pb /= ps;pc /= ps;pd /= ps;// 距離沒有取絕對值// 兩個向量的點乘，分母除以ps中已經除掉了，// 加pd原因:pointSel與tripod1構成的線段需要相減float pd2 = pa * pointSel.x + pb * pointSel.y + pc * pointSel.z + pd;float s = 1;if (iterCount >= 5) {// /加上影響因子s = 1 - 1.8 * fabs(pd2) / sqrt(sqrt(pointSel.x * pointSel.x+ pointSel.y * pointSel.y + pointSel.z * pointSel.z));}if (s > 0.1 && pd2 != 0) {// [x,y,z]是整個平面的單位法量// intensity是平面外一點到該平面的距離coeff.x = s * pa;coeff.y = s * pb;coeff.z = s * pc;coeff.intensity = s * pd2;// 未經變換的點放入laserCloudOri隊列，距離，法向量值放入coeffSellaserCloudOri->push_back(surfPointsFlat->points[i]);coeffSel->push_back(coeff);}}}}

面特征優化（這塊我看暈了，以下參考LeGO-LOAM源碼解析5： featureAssociation（三））：
定義各個矩陣，matA表示J，matAtA表示H，matAtB表示b：

int pointSelNum = laserCloudOri->points.size();cv::Mat matA(pointSelNum, 3, CV_32F, cv::Scalar::all(0));cv::Mat matAt(3, pointSelNum, CV_32F, cv::Scalar::all(0));cv::Mat matAtA(3, 3, CV_32F, cv::Scalar::all(0));cv::Mat matB(pointSelNum, 1, CV_32F, cv::Scalar::all(0));cv::Mat matAtB(3, 1, CV_32F, cv::Scalar::all(0));cv::Mat matX(3, 1, CV_32F, cv::Scalar::all(0));

對每個線特征點求關于對應三個變量的解雅可比矩陣A和損失向量B，拼裝成為一個大的求解矩陣：

float srx = sin(transformCur[0]);float crx = cos(transformCur[0]);float sry = sin(transformCur[1]);float cry = cos(transformCur[1]);float srz = sin(transformCur[2]);float crz = cos(transformCur[2]);float tx = transformCur[3];float ty = transformCur[4];float tz = transformCur[5];float a1 = crx*sry*srz; float a2 = crx*crz*sry; float a3 = srx*sry; float a4 = tx*a1 - ty*a2 - tz*a3;float a5 = srx*srz; float a6 = crz*srx; float a7 = ty*a6 - tz*crx - tx*a5;float a8 = crx*cry*srz; float a9 = crx*cry*crz; float a10 = cry*srx; float a11 = tz*a10 + ty*a9 - tx*a8;float b1 = -crz*sry - cry*srx*srz; float b2 = cry*crz*srx - sry*srz;float b5 = cry*crz - srx*sry*srz; float b6 = cry*srz + crz*srx*sry;float c1 = -b6; float c2 = b5; float c3 = tx*b6 - ty*b5; float c4 = -crx*crz; float c5 = crx*srz; float c6 = ty*c5 + tx*-c4;float c7 = b2; float c8 = -b1; float c9 = tx*-b2 - ty*-b1;// 構建雅可比矩陣，求解for (int i = 0; i < pointSelNum; i++) {pointOri = laserCloudOri->points[i];coeff = coeffSel->points[i];float arx = (-a1*pointOri.x + a2*pointOri.y + a3*pointOri.z + a4) * coeff.x+ (a5*pointOri.x - a6*pointOri.y + crx*pointOri.z + a7) * coeff.y+ (a8*pointOri.x - a9*pointOri.y - a10*pointOri.z + a11) * coeff.z;float arz = (c1*pointOri.x + c2*pointOri.y + c3) * coeff.x+ (c4*pointOri.x - c5*pointOri.y + c6) * coeff.y+ (c7*pointOri.x + c8*pointOri.y + c9) * coeff.z;float aty = -b6 * coeff.x + c4 * coeff.y + b2 * coeff.z;float d2 = coeff.intensity;matA.at<float>(i, 0) = arx;matA.at<float>(i, 1) = arz;matA.at<float>(i, 2) = aty;matB.at<float>(i, 0) = -0.05 * d2;}

使用OpenCV函數進行求解，利用QR分解加速：

cv::transpose(matA, matAt);matAtA = matAt * matA;matAtB = matAt * matB;cv::solve(matAtA, matAtB, matX, cv::DECOMP_QR);

關于退化問題的討論，與LOAM中一致，請參考loam源碼解析5 ： laserOdometry（三）中的退化問題分析：

if (iterCount == 0) {cv::Mat matE(1, 3, CV_32F, cv::Scalar::all(0));cv::Mat matV(3, 3, CV_32F, cv::Scalar::all(0));cv::Mat matV2(3, 3, CV_32F, cv::Scalar::all(0));cv::eigen(matAtA, matE, matV);matV.copyTo(matV2);isDegenerate = false;float eignThre[3] = {10, 10, 10};for (int i = 2; i >= 0; i--) {if (matE.at<float>(0, i) < eignThre[i]) {for (int j = 0; j < 3; j++) {matV2.at<float>(i, j) = 0;}isDegenerate = true;} else {break;}}matP = matV.inv() * matV2;}if (isDegenerate) {cv::Mat matX2(3, 1, CV_32F, cv::Scalar::all(0));matX.copyTo(matX2);matX = matP * matX2;}

使用迭代計算的增量進行姿態更新：

transformCur[1] += matX.at<float>(0, 0);transformCur[3] += matX.at<float>(1, 0);transformCur[5] += matX.at<float>(2, 0);for(int i=0; i<6; i++){if(isnan(transformCur[i]))transformCur[i]=0;}

計算姿態是否合法，以及是否滿足迭代條件：

float deltaR = sqrt(pow(rad2deg(matX.at<float>(0, 0)), 2) +pow(rad2deg(matX.at<float>(1, 0)), 2));float deltaT = sqrt(pow(matX.at<float>(2, 0) * 100, 2));if (deltaR < 0.1 && deltaT < 0.1) {return false;}return true;}

線特征匹配與優化：

for (int iterCount2 = 0; iterCount2 < 25; iterCount2++) {laserCloudOri->clear();coeffSel->clear();// 找到對應的特征邊/角點// 尋找邊特征的方法和尋找平面特征的很類似，過程可以參照尋找平面特征的注釋findCorrespondingCornerFeatures(iterCount2);if (laserCloudOri->points.size() < 10)continue;// 通過角/邊特征的匹配，計算變換矩陣if (calculateTransformationCorner(iterCount2) == false)break;}}

這塊和面特征匹配優化的思想類似：
將曲率大的點（待匹配的特征點）的坐標全部轉換到當前幀的初始時刻：

void findCorrespondingCornerFeatures(int iterCount){int cornerPointsSharpNum = cornerPointsSharp->points.size();for (int i = 0; i < cornerPointsSharpNum; i++) {TransformToStart(&cornerPointsSharp->points[i], &pointSel);

每五次迭代重新匹配特征點：

if (iterCount % 5 == 0) {

在上一幀曲率較大的點組成的KD樹中尋找一個最近鄰點，找到該最近鄰點的線束，

kdtreeCornerLast->nearestKSearch(pointSel, 1, pointSearchInd, pointSearchSqDis);int closestPointInd = -1, minPointInd2 = -1;if (pointSearchSqDis[0] < nearestFeatureSearchSqDist) {closestPointInd = pointSearchInd[0];int closestPointScan = int(laserCloudCornerLast->points[closestPointInd].intensity);

在最近鄰點i ii所在線束的相鄰線束中尋找第二個最近點，組成匹配直線：
首先在雷達點云存儲的正方向(ScanID增大的方向和雷達旋轉正方向)尋找另外一個點：
然后在雷達點云存儲的反方向(ScanID減小的方向和雷達旋轉反方向)尋找其他兩個點：

float pointSqDis, minPointSqDis2 = nearestFeatureSearchSqDist;for (int j = closestPointInd + 1; j < cornerPointsSharpNum; j++) {if (int(laserCloudCornerLast->points[j].intensity) > closestPointScan + 2.5) {break;}pointSqDis = (laserCloudCornerLast->points[j].x - pointSel.x) * (laserCloudCornerLast->points[j].x - pointSel.x) + (laserCloudCornerLast->points[j].y - pointSel.y) * (laserCloudCornerLast->points[j].y - pointSel.y) + (laserCloudCornerLast->points[j].z - pointSel.z) * (laserCloudCornerLast->points[j].z - pointSel.z);if (int(laserCloudCornerLast->points[j].intensity) > closestPointScan) {if (pointSqDis < minPointSqDis2) {minPointSqDis2 = pointSqDis;minPointInd2 = j;}}}for (int j = closestPointInd - 1; j >= 0; j--) {if (int(laserCloudCornerLast->points[j].intensity) < closestPointScan - 2.5) {break;}pointSqDis = (laserCloudCornerLast->points[j].x - pointSel.x) * (laserCloudCornerLast->points[j].x - pointSel.x) + (laserCloudCornerLast->points[j].y - pointSel.y) * (laserCloudCornerLast->points[j].y - pointSel.y) + (laserCloudCornerLast->points[j].z - pointSel.z) * (laserCloudCornerLast->points[j].z - pointSel.z);if (int(laserCloudCornerLast->points[j].intensity) < closestPointScan) {if (pointSqDis < minPointSqDis2) {minPointSqDis2 = pointSqDis;minPointInd2 = j;}}}}

找到用于匹配的兩個點（形成邊緣線）：

pointSearchCornerInd1[i] = closestPointInd;pointSearchCornerInd2[i] = minPointInd2;}

找到特征點的對應點以后的工作就是計算距離了，在這里還計算了就是點到直線的方向向量在各個軸的方向分解，在雅可比的計算中使用：

if (pointSearchCornerInd2[i] >= 0) {tripod1 = laserCloudCornerLast->points[pointSearchCornerInd1[i]];tripod2 = laserCloudCornerLast->points[pointSearchCornerInd2[i]];float x0 = pointSel.x;float y0 = pointSel.y;float z0 = pointSel.z;float x1 = tripod1.x;float y1 = tripod1.y;float z1 = tripod1.z;float x2 = tripod2.x;float y2 = tripod2.y;float z2 = tripod2.z;float m11 = ((x0 - x1)*(y0 - y2) - (x0 - x2)*(y0 - y1));float m22 = ((x0 - x1)*(z0 - z2) - (x0 - x2)*(z0 - z1));float m33 = ((y0 - y1)*(z0 - z2) - (y0 - y2)*(z0 - z1));float a012 = sqrt(m11 * m11 + m22 * m22 + m33 * m33);float l12 = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) + (z1 - z2)*(z1 - z2));float la = ((y1 - y2)*m11 + (z1 - z2)*m22) / a012 / l12;float lb = -((x1 - x2)*m11 - (z1 - z2)*m33) / a012 / l12;float lc = -((x1 - x2)*m22 + (y1 - y2)*m33) / a012 / l12;float ld2 = a012 / l12;float s = 1;if (iterCount >= 5) {s = 1 - 1.8 * fabs(ld2);}if (s > 0.1 && ld2 != 0) {coeff.x = s * la; coeff.y = s * lb;coeff.z = s * lc;coeff.intensity = s * ld2;laserCloudOri->push_back(cornerPointsSharp->points[i]);coeffSel->push_back(coeff);}}}}

線特征優化（我已經麻了）：

bool calculateTransformationCorner(int iterCount){int pointSelNum = laserCloudOri->points.size();cv::Mat matA(pointSelNum, 3, CV_32F, cv::Scalar::all(0));cv::Mat matAt(3, pointSelNum, CV_32F, cv::Scalar::all(0));cv::Mat matAtA(3, 3, CV_32F, cv::Scalar::all(0));cv::Mat matB(pointSelNum, 1, CV_32F, cv::Scalar::all(0));cv::Mat matAtB(3, 1, CV_32F, cv::Scalar::all(0));cv::Mat matX(3, 1, CV_32F, cv::Scalar::all(0));// 以下為開始計算A,A=[J的偏導],J的偏導的計算公式是什么?float srx = sin(transformCur[0]);float crx = cos(transformCur[0]);float sry = sin(transformCur[1]);float cry = cos(transformCur[1]);float srz = sin(transformCur[2]);float crz = cos(transformCur[2]);float tx = transformCur[3];float ty = transformCur[4];float tz = transformCur[5];float b1 = -crz*sry - cry*srx*srz; float b2 = cry*crz*srx - sry*srz; float b3 = crx*cry; float b4 = tx*-b1 + ty*-b2 + tz*b3;float b5 = cry*crz - srx*sry*srz; float b6 = cry*srz + crz*srx*sry; float b7 = crx*sry; float b8 = tz*b7 - ty*b6 - tx*b5;float c5 = crx*srz;for (int i = 0; i < pointSelNum; i++) {pointOri = laserCloudOri->points[i];coeff = coeffSel->points[i];float ary = (b1*pointOri.x + b2*pointOri.y - b3*pointOri.z + b4) * coeff.x+ (b5*pointOri.x + b6*pointOri.y - b7*pointOri.z + b8) * coeff.z;float atx = -b5 * coeff.x + c5 * coeff.y + b1 * coeff.z;float atz = b7 * coeff.x - srx * coeff.y - b3 * coeff.z;float d2 = coeff.intensity;// A=[J的偏導]; B=[權重系數*(點到直線的距離)] 求解公式: AX=B// 為了讓左邊滿秩，同乘At-> At*A*X = At*BmatA.at<float>(i, 0) = ary;matA.at<float>(i, 1) = atx;matA.at<float>(i, 2) = atz;matB.at<float>(i, 0) = -0.05 * d2;}// transpose函數求得matA的轉置matAtcv::transpose(matA, matAt);matAtA = matAt * matA;matAtB = matAt * matB;// 通過QR分解的方法，求解方程AtA*X=AtB，得到Xcv::solve(matAtA, matAtB, matX, cv::DECOMP_QR);if (iterCount == 0) {cv::Mat matE(1, 3, CV_32F, cv::Scalar::all(0));cv::Mat matV(3, 3, CV_32F, cv::Scalar::all(0));cv::Mat matV2(3, 3, CV_32F, cv::Scalar::all(0));// 計算At*A的特征值和特征向量// 特征值存放在matE，特征向量matVcv::eigen(matAtA, matE, matV);matV.copyTo(matV2);// 退化的具體表現是指什么？isDegenerate = false;float eignThre[3] = {10, 10, 10};for (int i = 2; i >= 0; i--) {if (matE.at<float>(0, i) < eignThre[i]) {for (int j = 0; j < 3; j++) {matV2.at<float>(i, j) = 0;}// 存在比10小的特征值則出現退化isDegenerate = true;} else {break;}}matP = matV.inv() * matV2;}if (isDegenerate) {cv::Mat matX2(3, 1, CV_32F, cv::Scalar::all(0));matX.copyTo(matX2);matX = matP * matX2;}transformCur[1] += matX.at<float>(0, 0);transformCur[3] += matX.at<float>(1, 0);transformCur[5] += matX.at<float>(2, 0);for(int i=0; i<6; i++){if(isnan(transformCur[i]))transformCur[i]=0;}float deltaR = sqrt(pow(rad2deg(matX.at<float>(0, 0)), 2));float deltaT = sqrt(pow(matX.at<float>(1, 0) * 100, 2) +pow(matX.at<float>(2, 0) * 100, 2));if (deltaR < 0.1 && deltaT < 0.1) {return false;}return true;}

以上部分我不完全理解就不贅述了，總之經過上面的騷操作我們得到了當前幀的角度旋轉量和位移量。

8、積分總變換。

然后就算旋轉變化：

// 旋轉角的累計變化量void integrateTransformation(){float rx, ry, rz, tx, ty, tz; // AccumulateRotation作用// 將計算的兩幀之間的位姿“累加”起來，獲得相對于第一幀的旋轉矩陣// transformSum + (-transformCur) =(rx,ry,rz)AccumulateRotation(transformSum[0], transformSum[1], transformSum[2], -transformCur[0], -transformCur[1], -transformCur[2], rx, ry, rz);

其中 AccumulateRotation功能就是計算旋轉變化的：

void AccumulateRotation(float cx, float cy, float cz, float lx, float ly, float lz, float &ox, float &oy, float &oz){// 參考：https://www.cnblogs.com/ReedLW/p/9005621.html// 0--->(cx,cy,cz)--->(lx,ly,lz)// 從0時刻到(cx,cy,cz),然后在(cx,cy,cz)的基礎上又旋轉(lx,ly,lz)// 求最后總的位姿結果是什么？// R*p_cur = R_c*R_l*p_cur ==> R=R_c* R_l//// |cly*clz+sly*slx*slz clz*sly*slx-cly*slz clx*sly|// R_l=| clx*slz clx*clz -slx|// |cly*slx*slz-clz*sly cly*clz*slx+sly*slz cly*clx|// R_c=...// -srx=(ccx*scy,-scx,cly*clx)*(clx*slz,clx*clz,-slx)// ...// 然后根據R再來求(ox,oy,oz)float srx = cos(lx)*cos(cx)*sin(ly)*sin(cz) - cos(cx)*cos(cz)*sin(lx) - cos(lx)*cos(ly)*sin(cx);ox = -asin(srx);float srycrx = sin(lx)*(cos(cy)*sin(cz) - cos(cz)*sin(cx)*sin(cy)) + cos(lx)*sin(ly)*(cos(cy)*cos(cz) + sin(cx)*sin(cy)*sin(cz)) + cos(lx)*cos(ly)*cos(cx)*sin(cy);float crycrx = cos(lx)*cos(ly)*cos(cx)*cos(cy) - cos(lx)*sin(ly)*(cos(cz)*sin(cy) - cos(cy)*sin(cx)*sin(cz)) - sin(lx)*(sin(cy)*sin(cz) + cos(cy)*cos(cz)*sin(cx));oy = atan2(srycrx / cos(ox), crycrx / cos(ox));float srzcrx = sin(cx)*(cos(lz)*sin(ly) - cos(ly)*sin(lx)*sin(lz)) + cos(cx)*sin(cz)*(cos(ly)*cos(lz) + sin(lx)*sin(ly)*sin(lz)) + cos(lx)*cos(cx)*cos(cz)*sin(lz);float crzcrx = cos(lx)*cos(lz)*cos(cx)*cos(cz) - cos(cx)*sin(cz)*(cos(ly)*sin(lz) - cos(lz)*sin(lx)*sin(ly)) - sin(cx)*(sin(ly)*sin(lz) + cos(ly)*cos(lz)*sin(lx));oz = atan2(srzcrx / cos(ox), crzcrx / cos(ox));}

接著計算當前幀與初始幀的位移變化：

// 進行平移分量的更新float x1 = cos(rz) * (transformCur[3] - imuShiftFromStartX) - sin(rz) * (transformCur[4] - imuShiftFromStartY);float y1 = sin(rz) * (transformCur[3] - imuShiftFromStartX) + cos(rz) * (transformCur[4] - imuShiftFromStartY);float z1 = transformCur[5] - imuShiftFromStartZ;float x2 = x1;float y2 = cos(rx) * y1 - sin(rx) * z1;float z2 = sin(rx) * y1 + cos(rx) * z1;tx = transformSum[3] - (cos(ry) * x2 + sin(ry) * z2);ty = transformSum[4] - y2;tz = transformSum[5] - (-sin(ry) * x2 + cos(ry) * z2);

修正累計位姿變換（公式又讓我麻了，后面要是有心情在先詳細解釋這些數學原理）：

// 與accumulateRotatio聯合起來更新transformSum的rotation部分的工作// 可視為transformToEnd的下部分的逆過程PluginIMURotation(rx, ry, rz, imuPitchStart, imuYawStart, imuRollStart, imuPitchLast, imuYawLast, imuRollLast, rx, ry, rz); void PluginIMURotation(float bcx, float bcy, float bcz, float blx, float bly, float blz, float alx, float aly, float alz, float &acx, float &acy, float &acz){// 參考：https://www.cnblogs.com/ReedLW/p/9005621.html// -imuStart imuEnd 0// transformSum.rot= R * R * R// YXZ ZXY k+1// bcx,bcy,bcz (rx, ry, rz)構成了 R_(k+1)^(0)// blx,bly,blz（imuPitchStart, imuYawStart, imuRollStart） 構成了 R_(YXZ)^(-imuStart)// alx,aly,alz（imuPitchLast, imuYawLast, imuRollLast）構成了 R_(ZXY)^(imuEnd)float sbcx = sin(bcx);float cbcx = cos(bcx);float sbcy = sin(bcy);float cbcy = cos(bcy);float sbcz = sin(bcz);float cbcz = cos(bcz);float sblx = sin(blx);float cblx = cos(blx);float sbly = sin(bly);float cbly = cos(bly);float sblz = sin(blz);float cblz = cos(blz);float salx = sin(alx);float calx = cos(alx);float saly = sin(aly);float caly = cos(aly);float salz = sin(alz);float calz = cos(alz);float srx = -sbcx*(salx*sblx + calx*caly*cblx*cbly + calx*cblx*saly*sbly) - cbcx*cbcz*(calx*saly*(cbly*sblz - cblz*sblx*sbly) - calx*caly*(sbly*sblz + cbly*cblz*sblx) + cblx*cblz*salx) - cbcx*sbcz*(calx*caly*(cblz*sbly - cbly*sblx*sblz) - calx*saly*(cbly*cblz + sblx*sbly*sblz) + cblx*salx*sblz);acx = -asin(srx);float srycrx = (cbcy*sbcz - cbcz*sbcx*sbcy)*(calx*saly*(cbly*sblz - cblz*sblx*sbly) - calx*caly*(sbly*sblz + cbly*cblz*sblx) + cblx*cblz*salx) - (cbcy*cbcz + sbcx*sbcy*sbcz)*(calx*caly*(cblz*sbly - cbly*sblx*sblz) - calx*saly*(cbly*cblz + sblx*sbly*sblz) + cblx*salx*sblz) + cbcx*sbcy*(salx*sblx + calx*caly*cblx*cbly + calx*cblx*saly*sbly);float crycrx = (cbcz*sbcy - cbcy*sbcx*sbcz)*(calx*caly*(cblz*sbly - cbly*sblx*sblz) - calx*saly*(cbly*cblz + sblx*sbly*sblz) + cblx*salx*sblz) - (sbcy*sbcz + cbcy*cbcz*sbcx)*(calx*saly*(cbly*sblz - cblz*sblx*sbly) - calx*caly*(sbly*sblz + cbly*cblz*sblx) + cblx*cblz*salx) + cbcx*cbcy*(salx*sblx + calx*caly*cblx*cbly + calx*cblx*saly*sbly);acy = atan2(srycrx / cos(acx), crycrx / cos(acx));float srzcrx = sbcx*(cblx*cbly*(calz*saly - caly*salx*salz) - cblx*sbly*(caly*calz + salx*saly*salz) + calx*salz*sblx) - cbcx*cbcz*((caly*calz + salx*saly*salz)*(cbly*sblz - cblz*sblx*sbly) + (calz*saly - caly*salx*salz)*(sbly*sblz + cbly*cblz*sblx) - calx*cblx*cblz*salz) + cbcx*sbcz*((caly*calz + salx*saly*salz)*(cbly*cblz + sblx*sbly*sblz) + (calz*saly - caly*salx*salz)*(cblz*sbly - cbly*sblx*sblz) + calx*cblx*salz*sblz);float crzcrx = sbcx*(cblx*sbly*(caly*salz - calz*salx*saly) - cblx*cbly*(saly*salz + caly*calz*salx) + calx*calz*sblx) + cbcx*cbcz*((saly*salz + caly*calz*salx)*(sbly*sblz + cbly*cblz*sblx) + (caly*salz - calz*salx*saly)*(cbly*sblz - cblz*sblx*sbly) + calx*calz*cblx*cblz) - cbcx*sbcz*((saly*salz + caly*calz*salx)*(cblz*sbly - cbly*sblx*sblz) + (caly*salz - calz*salx*saly)*(cbly*cblz + sblx*sbly*sblz) - calx*calz*cblx*sblz);acz = atan2(srzcrx / cos(acx), crzcrx / cos(acx));}

最終位姿變換：

transformSum[0] = rx;transformSum[1] = ry;transformSum[2] = rz;transformSum[3] = tx;transformSum[4] = ty;transformSum[5] = tz;}

9、發布里程計及上一幀點云信息

發布里程計：

void publishOdometry(){geometry_msgs::Quaternion geoQuat = tf::createQuaternionMsgFromRollPitchYaw(transformSum[2], -transformSum[0], -transformSum[1]);// rx,ry,rz轉化為四元數發布laserOdometry.header.stamp = cloudHeader.stamp;laserOdometry.pose.pose.orientation.x = -geoQuat.y;laserOdometry.pose.pose.orientation.y = -geoQuat.z;laserOdometry.pose.pose.orientation.z = geoQuat.x;laserOdometry.pose.pose.orientation.w = geoQuat.w;laserOdometry.pose.pose.position.x = transformSum[3];laserOdometry.pose.pose.position.y = transformSum[4];laserOdometry.pose.pose.position.z = transformSum[5];pubLaserOdometry.publish(laserOdometry);// laserOdometryTrans 是用于tf廣播laserOdometryTrans.stamp_ = cloudHeader.stamp;laserOdometryTrans.setRotation(tf::Quaternion(-geoQuat.y, -geoQuat.z, geoQuat.x, geoQuat.w));laserOdometryTrans.setOrigin(tf::Vector3(transformSum[3], transformSum[4], transformSum[5]));tfBroadcaster.sendTransform(laserOdometryTrans);}

發布點云：
把特征點云投影到每幀的結尾時刻：

void publishCloudsLast(){updateImuRollPitchYawStartSinCos();int cornerPointsLessSharpNum = cornerPointsLessSharp->points.size();for (int i = 0; i < cornerPointsLessSharpNum; i++) {// TransformToEnd的作用是將k+1時刻的less特征點轉移至k+1時刻的sweep的結束位置處的雷達坐標系下TransformToEnd(&cornerPointsLessSharp->points[i], &cornerPointsLessSharp->points[i]);}int surfPointsLessFlatNum = surfPointsLessFlat->points.size();for (int i = 0; i < surfPointsLessFlatNum; i++) {TransformToEnd(&surfPointsLessFlat->points[i], &surfPointsLessFlat->points[i]);}

用KD樹存儲當前幀點云：

pcl::PointCloud<PointType>::Ptr laserCloudTemp = cornerPointsLessSharp;cornerPointsLessSharp = laserCloudCornerLast;laserCloudCornerLast = laserCloudTemp;laserCloudTemp = surfPointsLessFlat;surfPointsLessFlat = laserCloudSurfLast;laserCloudSurfLast = laserCloudTemp;laserCloudCornerLastNum = laserCloudCornerLast->points.size();laserCloudSurfLastNum = laserCloudSurfLast->points.size();if (laserCloudCornerLastNum > 10 && laserCloudSurfLastNum > 100) {kdtreeCornerLast->setInputCloud(laserCloudCornerLast);kdtreeSurfLast->setInputCloud(laserCloudSurfLast);}frameCount++;

發布各類點云（發布界外點云、線特征點云、面特征點云）：

if (frameCount >= skipFrameNum + 1) {frameCount = 0;// 調整坐標系，x=y,y=z,z=xadjustOutlierCloud();sensor_msgs::PointCloud2 outlierCloudLast2;pcl::toROSMsg(*outlierCloud, outlierCloudLast2);outlierCloudLast2.header.stamp = cloudHeader.stamp;outlierCloudLast2.header.frame_id = "/camera";pubOutlierCloudLast.publish(outlierCloudLast2);sensor_msgs::PointCloud2 laserCloudCornerLast2;pcl::toROSMsg(*laserCloudCornerLast, laserCloudCornerLast2);laserCloudCornerLast2.header.stamp = cloudHeader.stamp;laserCloudCornerLast2.header.frame_id = "/camera";pubLaserCloudCornerLast.publish(laserCloudCornerLast2);sensor_msgs::PointCloud2 laserCloudSurfLast2;pcl::toROSMsg(*laserCloudSurfLast, laserCloudSurfLast2);laserCloudSurfLast2.header.stamp = cloudHeader.stamp;laserCloudSurfLast2.header.frame_id = "/camera";pubLaserCloudSurfLast.publish(laserCloudSurfLast2);}}

參考：
LeGO-LOAM源碼解析5： featureAssociation（三）
坐標轉換與imu融合

總結

以上是生活随笔為你收集整理的【LeGO-LOAM论文阅读（二）--特征提取（二）】的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。