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

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

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

1. 引言

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

2. Morphological Filters（形態學濾波）

 2.1 膨脹/腐蝕

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

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

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

可以從圖中得知：“虛線”是先進行“腐蝕”操作所形成的表面，這個表面剔除瞭“樹木”點，但是大型建築物卻變得不完整。“實線”是對“腐蝕”操作結果進行“膨脹”操作所形成的表面，可以看出其又恢復瞭大型建築的形狀。基於此，我們可以知道，“開操作”具備去除地面上的細小地物，保留大型地物的能力，這種能力對於後續處理是非常重要的。

2.2 形態學濾波

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

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

2.建築與樹木在特定/局部區域的分佈。

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

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

3. Progressive Morphological Filters

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

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

為瞭解決這一新出現的問題，可以通過引入一個高度差閾值來解決。建築屋頂和地面點之間的高程差是“突變”（abrupt change），而地面高程是“漸變”（gradual change）。因此，二者之間高程變化中的明顯區別可以幫助我們進行區分。假設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為研究區域第k次迭代中建築物的最矮高程值。以dhT,k作為閾值，對於第k次迭代中的任意點p如果dhp,k < dhT,k則將其標記為地面點，反之為非地面點。通過這種方式，不同尺寸的建築物（樹）可以隨著迭代窗口尺寸的增加逐步被去除。
綜上所述，Progressive Morphological Filters的詳細流程如下圖所示：

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

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

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

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

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

3.1.1 窗口尺寸

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

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

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

3.1.2 高程差閾值

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

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

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

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

3.2 參數輸入/輸出

3.2.1 參數輸入

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

3.2.1 參數輸出

• 地面點
• 非地面點

 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::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_filtered (new pcl::PointCloud<pcl::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.read<pcl::PointXYZ> ("samp11-utm.pcd", *cloud);

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

  // Create the filtering object
  pcl::ProgressiveMorphologicalFilter<pcl::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::ExtractIndices<pcl::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.write<pcl::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.write<pcl::PointXYZ> ("samp11-utm_object.pcd", *cloud_filtered, false);

  return (0);
}

到此這篇關於python點雲地面點濾波(Progressive Morphological Filter)算法介紹（PCL庫）的文章就介紹到這瞭,更多相關python點雲地面點濾波PCL庫內容請搜索WalkonNet以前的文章或繼續瀏覽下面的相關文章希望大傢以後多多支持WalkonNet！

推薦閱讀：

• 點雲地面點濾波(Cloth Simulation Filter, CSF)"佈料"濾波算法介紹
• python 點雲地面點濾波-progressive TIN densification(PTD)算法介紹
• python實現CSF地面點濾波算法原理解析
• C++輸入流和輸出流 超級詳細
• C語言lidar_align雷達裡程計校準功能詳解