【ATL03学习】

介绍

关于四叉树应用于点云数据的分割应用相对较少，网上有少量的四叉树应用于图像分割的文章，我自己写了一个基础的四叉树分割点云，欢迎大家讨论指正

建立四叉树

{ D k − 1 ′ = { X ∣ X ⊆ D k − 1 , s u m ( X ) > max ⁡ l i m i t } D k = s p l i t ( D k − 1 ′ ) \left\{\begin{matrix} D'_{k-1}=\{X|X\subseteq D_{k-1},sum(X)>\max_{limit}\} \\ D_k=split(D'_{k-1}) \end{matrix}\right. {Dk−1′​={X∣X⊆Dk−1​,sum(X)>maxlimit​}Dk​=split(Dk−1′​)​
对整片光子区域建立划分为矩形片区，对于矩形片区中当矩形内光子的数量大于${\max }_{limit}$时继续划分为小矩阵，否则停止。

def dgs(k, x, h, maxl, l):
    # 构建四叉树
    ks = {}
    k2 = np.array([[k[0], k[0] + (k[1] - k[0]) / 2, k[2], k[2] + (k[3] - k[2]) / 2],
                   [k[0] + (k[1] - k[0]) / 2, k[1], k[2], k[2] + (k[3] - k[2]) / 2],
                   [k[0], k[0] + (k[1] - k[0]) / 2, k[2] + (k[3] - k[2]) / 2, k[3]],
                   [k[0] + (k[1] - k[0]) / 2, k[1], k[2] + (k[3] - k[2]) / 2, k[3]]])
    for i in range(4):
        # print(i)
        ns = tj(k2[i], x, h)
        # print(k2[i])
        if ns > maxl:
            ks[str(i)] = dgs(k2[i], x, h, maxl, l)
        else:
            k2s = k2[i]
            cens = np.around(np.log2(l / (k2s[1] - k2s[0] + 0.1)))
            # print(cens)
            k2s = np.append(k2s, [ns, int(cens)], axis=0)
            ks[str(i)] = k2s
            # print(ns)
    return ks

大津法求阈值

利用每一级的光子数量构建频率直方图，求类间方差，获取阈值。

def otsu(jz):
    #otsu算法求阈值
    jz = np.array(jz)
    t = jz[:, 5]
    s = jz[:, 4]
    ns = np.array(hfh(t, s))
    t = list(ns[:, 0])
    s = ns[:, 1]
    p = s / np.sum(s)
    w1 = []
    w2 = []
    miu1 = []
    miu2 = []
    miu = []
    I = list(range(len(t)))
    sig2 = []
    for i in range(len(t)):
        w1.append(sum(p[:i]))
        w2.append(sum(p[i:]))
        miu1.append(sum(I[:i]) * sum(p[:i]) / (w1[i] + 0.01))
        miu2.append(sum(I[i:]) * sum(p[i:]) / (w2[i] + 0.01))
        miu.append(w1[i] * miu1[i] + w2[i] * miu2[i])
        sig2.append(w1[i] * (miu2[i] - miu[i]) + w2[i] * (miu2[i] - miu[i]))

    return t[sig2.index(max(sig2))]

结果展示

完整代码

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def tj(k, x, h):
    #统计方块内的光子数量
    return np.sum((x >= k[0]) & (x <= k[1]) & (h >= k[2]) & (h <= k[3]))


def area(k):
    #求面积
    return (k[1] - k[0]) * (k[3] - k[2])


def dgs(k, x, h, maxl, l):
    # 构建四叉树
    ks = {}
    k2 = np.array([[k[0], k[0] + (k[1] - k[0]) / 2, k[2], k[2] + (k[3] - k[2]) / 2],
                   [k[0] + (k[1] - k[0]) / 2, k[1], k[2], k[2] + (k[3] - k[2]) / 2],
                   [k[0], k[0] + (k[1] - k[0]) / 2, k[2] + (k[3] - k[2]) / 2, k[3]],
                   [k[0] + (k[1] - k[0]) / 2, k[1], k[2] + (k[3] - k[2]) / 2, k[3]]])
    for i in range(4):
        # print(i)
        ns = tj(k2[i], x, h)
        # print(k2[i])
        if ns > maxl:
            ks[str(i)] = dgs(k2[i], x, h, maxl, l)
        else:
            k2s = k2[i]
            cens = np.around(np.log2(l / (k2s[1] - k2s[0] + 0.1)))
            # print(cens)
            k2s = np.append(k2s, [ns, int(cens)], axis=0)
            ks[str(i)] = k2s
            # print(ns)
    return ks



def zhuanghua(x, h, k):
    #判断方块内是否含有光子
    return (x >= k[0]) & (x <= k[1]) & (h >= k[2]) & (h <= k[3])


def qz(x, h, jz, thr):
    #去噪函数，jz为转化矩阵，thr为阈值
    label = np.zeros(len(x))
    for i in jz:
        if i[5] >= thr:
            label[zhuanghua(x, h, i)] = 1
    return label


def zh(ks, kf):
    #字典转化为,矩阵
    for i in ks:
        if isinstance(ks[i], dict):
            zh(ks[i], kf)
        else:
            kf.append(list(ks[i]))
    return kf


def hfh(t, s):
    #统计各级的光子数量
    ns = dict()
    for i in range(len(t)):
        if t[i] in ns.keys():
            ns[t[i]] += s[i]
        else:
            ns[t[i]] = s[i]
    ns = sorted(ns.items(), key=lambda x: x[0])
    return ns


def otsu(jz):
    #otsu算法求阈值
    jz = np.array(jz)
    t = jz[:, 5]
    s = jz[:, 4]
    ns = np.array(hfh(t, s))
    t = list(ns[:, 0])
    s = ns[:, 1]
    p = s / np.sum(s)
    w1 = []
    w2 = []
    miu1 = []
    miu2 = []
    miu = []
    I = list(range(len(t)))
    sig2 = []
    for i in range(len(t)):
        w1.append(sum(p[:i]))
        w2.append(sum(p[i:]))
        miu1.append(sum(I[:i]) * sum(p[:i]) / (w1[i] + 0.01))
        miu2.append(sum(I[i:]) * sum(p[i:]) / (w2[i] + 0.01))
        miu.append(w1[i] * miu1[i] + w2[i] * miu2[i])
        sig2.append(w1[i] * (miu2[i] - miu[i]) + w2[i] * (miu2[i] - miu[i]))

    return t[sig2.index(max(sig2))]

f = 'ATL03_20200327094144_00130706_005_01_gt3r.csv'
data = pd.read_csv(f)
print(data.keys())
#预处理
x, h, lat = np.array(data['x']), np.array(data['h']), np.array(data['lat'])
r1 = np.argsort(x)
x, h, lat = x[r1], h[r1], lat[r1]
r_id = np.where((lat < 44.9) & (lat > 44.8))
print(r_id[-1])
x = x[r_id]
h = h[r_id]
lat = lat[r_id]
#去噪前点云
plt.figure()
plt.scatter(x,h,marker='.',c='black')
plt.xlabel('Along-track distance/m')
plt.ylabel('Heights/m')
plt.show()

xiankuang = np.array([np.min(x), np.max(x), np.min(h), np.max(h)])
s = np.sum((x >= xiankuang[0]) & (x <= xiankuang[1]) & (h >= xiankuang[2]) & (h <= xiankuang[3]))
# k=xiankuang
l = xiankuang[1] - xiankuang[0]
ks = dgs(xiankuang, x, h, 2, l)
# print(ks.values())
jz = zh(ks, [])

#jz2 = np.array(jz)
thr = otsu(jz)
lab = qz(x, h, jz, thr)
plt.figure()
plt.scatter(x[lab == 0], h[lab == 0], marker='.', c='black', label='noise')
plt.scatter(x[lab == 1], h[lab == 1], marker='.', c='red', label='signal')
plt.xlabel('Along-track distance/m')
plt.ylabel('Heights/m')
plt.legend()
plt.show()
# s1 = (x < 2) & (h < 1000)
# lab = np.zeros(len(x))
# lab = qz(x, h, ks, 0.1, lab)

还是对字典结构不熟悉，后续又将字典结构转化为矩阵做后续处理，欢迎大家批评指正。

参考文献

[1]刘翔,张立华,戴泽源,陈秋,周寅飞.一种无输入参数的强噪声背景下ICESat-2点云去噪方法[J].光子学报,2022,51(11):354-364.
[2]Xie H, Sun Y, Xu Q, Li B, Guo Y, Liu X, Huang P, Tong X. Converting along-track photons into a point-region quadtree to assist with ICESat-2-based canopy cover and ground photon detection[J]. International Journal of Applied Earth Observation and Geoinformation, 2022, 112: 102872