WLOP算法python实现源码

由于查了github上没有python版本的wlop算法，于是自己看论文按照公式实现了下。

各个参数，公式设置严格按照论文

输入为numpy （N*3）形状的点云

输出为降采样的点云

没有考虑太多计算复杂度

class WLOP:
    def __init__(self, points, mu=0.45):
        self.original_points = np.array(points)
        min_coords = np.min(self.original_points, axis=0)
        max_coords = np.max(self.original_points, axis=0)
        d_bb = np.linalg.norm(max_coords - min_coords)  
        self.h = 4 * np.sqrt(d_bb)  
        self.radius = self.h / 4 
        self.mu = mu  
        self.sample_points = None  

    def generate_initial_samples(self, n=1024):

        min_coords = np.min(self.original_points, axis=0)
        max_coords = np.max(self.original_points, axis=0)
        return np.random.rand(n, 3) * (max_coords - min_coords) + min_coords

    def compute_density_weights(self, points, tree):
       
        weights = np.zeros(len(points))
        radius2 = self.radius ** 2
        iradius16 = -1.0 / radius2  # exp(-r²/(h/4)²) 

        for i in range(len(points)):
            neighbors = tree.query_ball_point(points[i], self.radius)
            distances = np.linalg.norm(points[neighbors] - points[i], axis=1)
            total = 1.0  # 避免除以0
            for d in distances:
                if d > 1e-8:
                    total += np.exp(d * iradius16)
            weights[i] = 1.0 / total
        return weights

    def run(self, iterations=25, samples_num=1024):
        # 初始化
        self.sample_points = self.generate_initial_samples(samples_num)  # 固定1024点
        original_tree = cKDTree(self.original_points)
        sample_tree = cKDTree(self.sample_points)

        # 预计算原始点密度v（迭代前计算）
        v = self.compute_density_weights(self.original_points, original_tree)

        for _ in range(iterations):
            # 计算样本点密度w
            w = self.compute_density_weights(self.sample_points, sample_tree)

            # 更新每个样本点
            new_samples = np.zeros_like(self.sample_points)
            for i in range(len(self.sample_points)):
                # 计算吸引项（原始点贡献）
                neighbors = original_tree.query_ball_point(self.sample_points[i], self.radius)
                dists = np.linalg.norm(self.original_points[neighbors] - self.sample_points[i], axis=1)
                alpha_weights = np.exp(-dists ** 2 / (self.radius ** 2)) * v[neighbors] / (dists + 1e-8)
                sum_alpha = np.sum(alpha_weights)
                if sum_alpha > 1e-8:
                    new_samples[i] = np.dot(alpha_weights, self.original_points[neighbors]) / sum_alpha
                else:
                    new_samples[i] = self.sample_points[i]

                # 计算排斥项（样本点间贡献）
                sample_neighbors = sample_tree.query_ball_point(self.sample_points[i], self.radius)
                sample_dists = np.linalg.norm(self.sample_points[sample_neighbors] - self.sample_points[i], axis=1)
                valid = (sample_neighbors != i) & (sample_dists > 1e-8)
                beta_weights = -np.exp(-sample_dists ** 2 / (self.radius ** 2)) * w[sample_neighbors] / (
                            sample_dists + 1e-8)
                sum_beta = np.sum(beta_weights[valid])
                if sum_beta < -1e-8:
                    diff = self.sample_points[i] - self.sample_points[sample_neighbors][valid]
                    new_samples[i] += self.mu * np.dot(beta_weights[valid], diff) / sum_beta
                    print("new_samples[i] = ", new_samples[i])
            self.sample_points = new_samples
            print('iteration done!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
            sample_tree = cKDTree(self.sample_points)  # 更新KD树

    def get_samples(self):
        return self.sample_points