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

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

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

1. 引言

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

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

 2.1 膨脹/腐蝕

膨脹/腐蝕是其中的兩個(gè)基礎(chǔ)操作，通俗的說(shuō)膨脹/腐蝕可以擴(kuò)大/減小特征的尺寸，并以此組合為開(kāi)/閉操作。針對(duì)LiDAR測(cè)量點(diǎn)p(x, y, z)，高程 z 值在(x, y)處對(duì)應(yīng)的膨脹操作可以定義為：

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

了解膨脹/腐蝕這兩個(gè)基礎(chǔ)操作之后，可以通過(guò)對(duì)其進(jìn)行簡(jiǎn)單組合來(lái)來(lái)形成開(kāi)/閉操作，其中開(kāi)操作為先進(jìn)行腐蝕再進(jìn)行膨脹操作，而閉操作為先膨脹再進(jìn)行腐蝕。在LiDAR數(shù)據(jù)處理中使用了“開(kāi)”操作，處理效果如下圖中所示：

可以從圖中得知：“虛線”是先進(jìn)行“腐蝕”操作所形成的表面，這個(gè)表面剔除了“樹(shù)木”點(diǎn)，但是大型建筑物卻變得不完整?！皩?shí)線”是對(duì)“腐蝕”操作結(jié)果進(jìn)行“膨脹”操作所形成的表面，可以看出其又恢復(fù)了大型建筑的形狀?；诖?，我們可以知道，“開(kāi)操作”具備去除地面上的細(xì)小地物，保留大型地物的能力，這種能力對(duì)于后續(xù)處理是非常重要的。

2.2 形態(tài)學(xué)濾波

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

1.濾波窗口的尺寸，如果窗口尺寸設(shè)置的比較小，則可以很好的保留地面起伏的細(xì)節(jié)，但是只能去除像汽車(chē)、樹(shù)木等細(xì)小地物，而對(duì)建筑物則去除效果較差（屋頂通常被認(rèn)為是地面）。相反，若窗口尺寸設(shè)置的較大，則會(huì)過(guò)度去除一些“地面點(diǎn)”，例如，一條道路與一條小水溝相鄰，若窗口尺寸大于道路的寬度，則道路可能就會(huì)被認(rèn)為是非地面點(diǎn)（因?yàn)樾∷疁现械狞c(diǎn)高程較低，會(huì)被認(rèn)為是窗口內(nèi)的最低點(diǎn)，而道路點(diǎn)較高，被判斷為非地面點(diǎn)）。此外，一些局部的小山丘、沙丘都極可能被“切除”。

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

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

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

3. Progressive Morphological Filters

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

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

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

我們可以對(duì)上圖總結(jié)以下四個(gè)步驟：

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

注意：每一次迭代中的“開(kāi)操作”實(shí)際都是作用在步驟1）所劃分網(wǎng)格中的點(diǎn)，所以Progressive Morphological
Filters是"point level"來(lái)對(duì)LiDAR測(cè)量值進(jìn)行濾波處理的。

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

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

3.1.1 窗口尺寸

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

式中，k為迭代次數(shù)，b是初始窗口大?。ㄓ捎脩暨M(jìn)行輸入），最后+1是為了保證 wk為一個(gè)奇數(shù)，窗口對(duì)稱(chēng)。但是，如果一個(gè)研究區(qū)域具有非常大的地物，這種增加窗口尺寸速度太慢則會(huì)耗費(fèi)較多時(shí)間。因此，可以使用第二種方式，通過(guò)指數(shù)增長(zhǎng)來(lái)改變窗口大小，計(jì)算如下式：

同樣的，式中k為迭代次數(shù)，b是初始窗口大?。ㄓ捎脩暨M(jìn)行輸入），這種方式的增長(zhǎng)速度較第一種方式快很多。

3.1.2 高程差閾值

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

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

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

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

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

3.2.1 參數(shù)輸入

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

3.2.1 參數(shù)輸出

• 地面點(diǎn)
• 非地面點(diǎn)

 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);
}

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