python點云地面點濾波(Progressive Morphological Filter)算法介紹（PCL庫）

本篇博客參考Keqi Zhang的文章“A Progressive Morphological Filter for Removing Nonground Measurements From Airborne LIDAR Data”以去除大型建筑、樹木等常見地物。

不方便的小伙伴可以在此：資源鏈接下載。

1. 引言

機載LiDAR可以獲取快速、低成本地獲取大區(qū)域的高精度地形測量值。為了獲取高精度DTM/DEM需要區(qū)分測量點中的地面點（由地面直接返回）及非地面點（建筑、車、植被）。眾多學者采用了各種各樣的方法來進行"點云地面點濾波"。（此篇博客中也進行了相關介紹，不再驁述）

2. Morphological Filters（形態(tài)學濾波）

 2.1 膨脹/腐蝕

膨脹/腐蝕是其中的兩個基礎操作，通俗的說膨脹/腐蝕可以擴大/減小特征的尺寸，并以此組合為開/閉操作。針對LiDAR測量點p(x, y, z)，高程 z 值在(x, y)處對應的膨脹操作可以定義為：

式中：(xp, yp, zp) 代表p點的相鄰點，w為操作窗口(可以為一維“線”也可以為二維“矩形/圓/其他形狀等”)。膨脹操作完成后會輸出p點在窗口w內具有最大高程值的近鄰點。
與之類似的，腐蝕操作為在p點窗口w內找到具有最低高程值的近鄰點?？梢酝ㄟ^下式進行定義：

了解膨脹/腐蝕這兩個基礎操作之后，可以通過對其進行簡單組合來來形成開/閉操作，其中開操作為先進行腐蝕再進行膨脹操作，而閉操作為先膨脹再進行腐蝕。在LiDAR數(shù)據處理中使用了“開”操作，處理效果如下圖中所示：

可以從圖中得知：“虛線”是先進行“腐蝕”操作所形成的表面，這個表面剔除了“樹木”點，但是大型建筑物卻變得不完整?！皩嵕€”是對“腐蝕”操作結果進行“膨脹”操作所形成的表面，可以看出其又恢復了大型建筑的形狀?；诖耍覀兛梢灾?，“開操作”具備去除地面上的細小地物，保留大型地物的能力，這種能力對于后續(xù)處理是非常重要的。

2.2 形態(tài)學濾波

上述的“開操作”只是去除了細小地物，保留了大型地物，并沒有去除所有非地面點去除，而且僅僅通過一個“開操作”也不可能實現(xiàn)對復雜地表的提取。因此，怎么利用好“開操作”的能力進行下一步驟的提取是進一步提升的關鍵。
Kilian等人提出，可以在“開操作”處理后的結果中的每一個“窗口”內找到一個“最低點”，然后此窗口內的其他點若落在“最低點”的一個高程范圍內則認為是地面點。這個高程范圍通常根據機載LiDAR系統(tǒng)掃描的精度來定義，正常為20-30cm。
此方法中有兩個方面對最后的結果好壞非常重要：

1.濾波窗口的尺寸，如果窗口尺寸設置的比較小，則可以很好的保留地面起伏的細節(jié)，但是只能去除像汽車、樹木等細小地物，而對建筑物則去除效果較差（屋頂通常被認為是地面）。相反，若窗口尺寸設置的較大，則會過度去除一些“地面點”，例如，一條道路與一條小水溝相鄰，若窗口尺寸大于道路的寬度，則道路可能就會被認為是非地面點（因為小水溝中的點高程較低，會被認為是窗口內的最低點，而道路點較高，被判斷為非地面點）。此外，一些局部的小山丘、沙丘都極可能被“切除”。

2.建筑與樹木在特定/局部區(qū)域的分布。

注：一個最理想的情況是我們可以設置一個“窗口”，這個“窗口”的尺寸可以足夠小，能夠保留地面細節(jié)。同時，還可以足夠大，能夠去除建筑、汽車、樹木等地物。但是，這種理想情況在實際數(shù)據集中國并不存在。

為了解決這一問題，Kilian等人繼續(xù)提出了可以通過改變窗口大小來多次進行濾波。每個點都被賦予一個與窗口大小相關的權重，窗口尺寸越大，點的權重越高。這種方法雖然得到了更好一些的效果，但是沒有從"point level"進行區(qū)分地面點與非地面點。（"point level"區(qū)分的地面點與非地面點之后可以通過插值的方法使得DEM/DTM的生成效果更好。）

3. Progressive Morphological Filters

由上述2.2節(jié)中的分析可知，傳統(tǒng)的形態(tài)學濾波很難通過一個固定大小的窗口去檢測出各種尺寸變化的不同地物。這一問題可以通過逐漸改變窗口大小來解決。
如下圖中所示，首先使用一個尺寸為l1的窗口來對原始數(shù)據進行開操作。由圖中的“虛線”可以看出樹木等尺寸小于l1的地物被去除，且地形特征中小于l1的部分被“切除”（山丘頂部高程被替換成了l1中的最小值），但是，尺寸大于l1的建筑物被保留了下來。接著，進行下一次迭代，窗口尺寸變?yōu)閘2，對上一次的處理結果進行開操作處理，結果從“實線”中可以看出，l2大于建筑的尺寸，所以建筑也被去除，但同時山丘頂部被“切除”的范圍更大。

需要注意的是： 通過逐漸增加窗口尺寸解決了去除不同大小地物的問題，但是又引入了"山丘"頂部等小于窗口尺寸的地形特征部分被“切除”的問題。

為了解決這一新出現(xiàn)的問題，可以通過引入一個高度差閾值來解決。建筑屋頂和地面點之間的高程差是“突變”（abrupt change），而地面高程是“漸變”（gradual change）。因此，二者之間高程變化中的明顯區(qū)別可以幫助我們進行區(qū)分。假設dhp,1代表原始LiDAR測量值與在任意給定p點處第一次迭代表面之間的高程差，dhT,1代表高程差閾值，則如果dhp,1 ≤ dhT,1點p就被認為是地面點，反之如果dhp,1 > dhT,1就認為點p是一個非地面點。此后，令dhmax(t),1為當前迭代中初始地面點與濾波表面之間差值的最大值，則如果選取的dhT,1 > dhmax(t),1則所有的測量值都會保留。
在第二次迭代中假設上一次濾波表面和本次濾波表面的最大高程差為dhmax(t),2，則如果dhmax(t),2 dhT,2,則高程差值在dhmax(t),2范圍內的地面點都會被保留。類似的，假設在上次迭代和本次迭代之間建筑高程差值最小為dhmin(b),2（通常近似為建筑的高度），如果dhmin(b),2 > dhT,2，則建筑就會被移除。
通常設置dhT,k為研究區(qū)域第k次迭代中建筑物的最矮高程值。以dhT,k作為閾值，對于第k次迭代中的任意點p如果dhp,k dhT,k則將其標記為地面點，反之為非地面點。通過這種方式，不同尺寸的建筑物（樹）可以隨著迭代窗口尺寸的增加逐步被去除。
綜上所述，Progressive Morphological Filters的詳細流程如下圖所示：

我們可以對上圖總結以下四個步驟：

1. 加載不規(guī)則點云，劃分為規(guī)則網格，在每個網格中選取高程最低點（如果網格中沒有點則根據最近鄰點進行插值），構建最小高程表面。
2. 使用輸入的初始濾波窗口尺寸、1）中獲取的最小高程表面作為第一次迭代的參數(shù)進行第一次迭代。隨后，在后續(xù)的迭代中，以前一次獲取的濾波表面及3）中計算的濾波窗口尺寸作為輸入。每次迭代的輸出有兩部分：a) 形態(tài)學濾波后得到的更加光滑的表面；b) 基于不同閾值所檢測出來的非地面點。
3. 計算新的濾波窗口尺寸及不同的高程插值閾值。重復步驟2)、步驟3)直到窗口尺寸大于預設的最大窗口。
4. 基于去除非地面測量值的數(shù)據進行生成DEM/DTM。

注意：每一次迭代中的“開操作”實際都是作用在步驟1）所劃分網格中的點，所以Progressive Morphological
Filters是"point level"來對LiDAR測量值進行濾波處理的。

3.1 參數(shù)計算(窗口尺寸/高程差閾值)

在上述步驟3）中我們要變化窗口尺寸 wk和高程差閾值 dhT,k兩個參數(shù)的值，以進行下一次迭代，那么這兩個值是怎么計算的呢？

3.1.1 窗口尺寸

首先是窗口尺寸 wk有兩種計算方式，第一種是：

式中，k為迭代次數(shù)，b是初始窗口大小（由用戶進行輸入），最后+1是為了保證 wk為一個奇數(shù)，窗口對稱。但是，如果一個研究區(qū)域具有非常大的地物，這種增加窗口尺寸速度太慢則會耗費較多時間。因此，可以使用第二種方式，通過指數(shù)增長來改變窗口大小，計算如下式：

同樣的，式中k為迭代次數(shù)，b是初始窗口大小（由用戶進行輸入），這種方式的增長速度較第一種方式快很多。

3.1.2 高程差閾值

高程差閾值與研究區(qū)域的地形坡度密不可分，因此可以通過下式進行計算：

式中，dh0為初始高程差閾值，s為坡度，c為格網大小，dhmax為最大高程差閾值。

在城市區(qū)域，樹木、汽車相對于建筑的尺寸小很多，所以通常是最后濾除建筑，最大高程差閾值dhmax可以設置為一個固定值（如最矮建筑物高度）。而在山區(qū)，主要的非地面點為植被，所以并沒有必要設置固定的最大高程差閾值dhmax，于是dhmax通常被設置為測區(qū)內的最大高程差。

此外，坡度s通過第k次迭代的最大高程差dhmax(t),k，以及窗口尺寸wk進行計算，如下式所示：

3.2 參數(shù)輸入/輸出

3.2.1 參數(shù)輸入

•  原始LiDAR點云數(shù)據，每個點都由(x, y, z)進行表示
• 劃分格網大小c 參數(shù)b（計算窗口尺寸）
• 最大窗口尺寸（判斷是否停止迭代）
• 地形坡度s
• 初始高程差閾值dh0
• 最大高程差值dhmax

3.2.1 參數(shù)輸出

• 地面點
• 非地面點

 3.3 PCL官方示例代碼

#include iostream>
#include pcl/io/pcd_io.h>
#include pcl/point_types.h>
#include pcl/filters/extract_indices.h>
#include pcl/segmentation/progressive_morphological_filter.h>

int
main (int argc, char** argv)
{
  pcl::PointCloudpcl::PointXYZ>::Ptr cloud (new pcl::PointCloudpcl::PointXYZ>);
  pcl::PointCloudpcl::PointXYZ>::Ptr cloud_filtered (new pcl::PointCloudpcl::PointXYZ>);
  pcl::PointIndicesPtr ground (new pcl::PointIndices);

  // Fill in the cloud data
  pcl::PCDReader reader;
  // Replace the path below with the path where you saved your file
  reader.readpcl::PointXYZ> ("samp11-utm.pcd", *cloud);

  std::cerr  "Cloud before filtering: "  std::endl;
  std::cerr  *cloud  std::endl;

  // Create the filtering object
  pcl::ProgressiveMorphologicalFilterpcl::PointXYZ> pmf;
  pmf.setInputCloud (cloud);
  pmf.setMaxWindowSize (20);
  pmf.setSlope (1.0f);
  pmf.setInitialDistance (0.5f);
  pmf.setMaxDistance (3.0f);
  pmf.extract (ground->indices);

  // Create the filtering object
  pcl::ExtractIndicespcl::PointXYZ> extract;
  extract.setInputCloud (cloud);
  extract.setIndices (ground);
  extract.filter (*cloud_filtered);

  std::cerr  "Ground cloud after filtering: "  std::endl;
  std::cerr  *cloud_filtered  std::endl;

  pcl::PCDWriter writer;
  writer.writepcl::PointXYZ> ("samp11-utm_ground.pcd", *cloud_filtered, false);

  // Extract non-ground returns
  extract.setNegative (true);
  extract.filter (*cloud_filtered);

  std::cerr  "Object cloud after filtering: "  std::endl;
  std::cerr  *cloud_filtered  std::endl;

  writer.writepcl::PointXYZ> ("samp11-utm_object.pcd", *cloud_filtered, false);

  return (0);
}

到此這篇關于python點云地面點濾波(Progressive Morphological Filter)算法介紹（PCL庫）的文章就介紹到這了,更多相關python點云地面點濾波PCL庫內容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關文章希望大家以后多多支持腳本之家！