从激光雷达点云转换为 2D 图像时 2D 图像发生变化

Question

我正在将激光雷达点云映射到二维图像。我已成功将点云映射到二维图像，但重叠时，二维图像与激光雷达点云有 0.5 米到 1 米的偏移。这种重叠对我来说很重要，因为我将使用由激光雷达数据创建的多边形，多边形和图像也不会正确重叠。二维图像的范围是激光雷达数据的范围。我无法弄清楚为什么存在小偏差，即使它们位于相同的坐标系中并且图像与激光雷达数据具有相同的范围。

这是我正在处理的数据。

https://drive.google.com/drive/folders/1cBNDJmoh8i3Sb4AoU5FLDqPXjy7OYfaW?usp=sharing

这是我已经厌倦的代码。

import laspy
import rasterio
import numpy as np
from PIL import Image
from scipy.interpolate import griddata
import matplotlib.pyplot as plt

import geopandas as gpd
class SingleTree:
    def __init__(self, las_data):
        self.points = np.vstack((las_data.x, las_data.y, las_data.z, las_data.red, las_data.green, 
                                 las_data.blue)).T

        
    def point_cloud_to_image(self, output_size):
        
        # Extract x, y, and color values
        x = self.points[:, 0]
        y = self.points[:, 1]
        red = self.points[:, 3]
        green = self.points[:, 4]
        blue = self.points[:, 5]
        
        # Normalize coordinates to range [0, 1]
        x_normalized = (x - x.min()) / (x.max() - x.min())
        y_normalized = (y - y.min()) / (y.max() - y.min())
        
        # Scale normalized coordinates to image size, ensuring they are within the bounds
        x_img = np.clip((x_normalized * (output_size[1] - 1)).astype(int), 0, output_size[1] - 1)
        y_img = np.clip((y_normalized * (output_size[0] - 1)).astype(int), 0, output_size[0] - 1)

        # Convert color values to 8-bit
        r = np.interp(red, (red.min(), red.max()), (0, 255)).astype(np.uint8)
        g = np.interp(green, (green.min(), green.max()), (0, 255)).astype(np.uint8)
        b = np.interp(blue, (blue.min(), blue.max()), (0, 255)).astype(np.uint8)

        # Initialize an empty image
        img = np.zeros(output_size + (3,), dtype=np.uint8)
        
        
        # Assign colors to image pixels safely
        for xi, yi, ri, gi, bi in zip(x_img, y_img, r, g, b):
            img[yi, xi] = [ri, gi, bi]
        
        return img

    def polygon_bound(self, poly_path):
        # Load the polygon data
        gdf = gpd.read_file(poly_path)

        # Assuming you're interested in the first polygon if there are multiple
        polygon = gdf.geometry[0]

        # Get the bounding box of the polygon
        minx, miny, maxx, maxy = polygon.bounds
        
        polygon_width = maxx - minx
        polygon_height = maxy - miny

        # The top-left corner coordinates
        top_left_x = minx
        top_left_y = miny
        
        return top_left_x, top_left_y, polygon_width, polygon_height

import os
from rasterio.transform import from_origin
# Assuming 'path' is the directory containing your '.las' files
path = r"lasdata path"
output_path = r"path"
number_of_files = len([name for name in os.listdir("directory of lasdata")])
discarded_lasdata = []
# Loop through your LAS files
for i in range(1, number_of_files):
    file = f"tree_i{i}.las"
    poly_file = "polygons\\poly{}.shp".format(i)
    image_path = os.path.join(path, file)
    
    # Read the LAS file
    las_data = laspy.read(image_path)
    
    # Process the point cloud to generate an image
    tree = SingleTree(las_data)
    
    top_left_x, top_left_y, polygon_width, polygon_height = tree.polygon_bound(poly_file)
    # Get metadata from LAS file for georeferencing
    scale = las_data.header.scale
    offset = las_data.header.offset
    
    # Assuming las_data is a LasData object from laspy
    min_x, min_y, min_z = las_data.header.min
    max_x, max_y, max_z = las_data.header.max
    
    lidar_width = round((max_x - min_x), 0)
    lidar_height = round((max_y - min_y), 0)
    
    if (lidar_width <= 0 or lidar_height <= 0):
        discarded_lasdata.append(i)
        continue
        
    ltop_left_x = min_x
    ltop_left_y = min_y
    
    resolution = 0.5
    
    image_width_pixels = int((lidar_width * 2 + 5) / resolution)
    image_height_pixels = int((lidar_height * 2 + 5) / resolution)
    
    img = tree.point_cloud_to_image((int(lidar_height*2), int(lidar_width*2)))
    
    # Create the geotransform
    transform = from_origin(ltop_left_x, ltop_left_y, resolution, -resolution)
    
    # Define the CRS (coordinate reference system)
    crs = "EPSG:25832"  # Replace with the correct EPSG code for your data
    
    # Define the rasterio metadata dictionary
    meta = {
        'driver': 'GTiff',
        'dtype': 'uint8',
        'nodata': None,
        'width': img.shape[1],
        'height': img.shape[0],
        'count': 3,
        'crs': crs,
        'transform': transform,
        'compress': 'lzw'
    }
    
    output_filename = os.path.join(output_path, f"lidar_test{i}.tif")
    
    # Write the image data and metadata to a GeoTIFF
    with rasterio.open(output_filename, 'w', **meta) as dst:
        for k in range(img.shape[2]):
            dst.write(img[:, :, k], k+1)

print("GeoTIFFs have been saved.")

您可以看到多边形没有与图像正确重叠，而所有激光雷达点都落在多边形内

编辑：第二次查看时，我发现激光雷达数据和图像的范围不相同，并且偏离了 0 到 1 米。我认为这是因为激光雷达数据的宽度和高度是浮点数，我们传递一个整数作为宽度和高度。但我仍然无法找到解决方案

Answer 1

我怀疑你的问题可能出在线路上

x_img = np.clip((x_normalized * (output_size[1] - 1)).astype(int), 0, output_size[1] - 1)
y_img = np.clip((y_normalized * (output_size[0] - 1)).astype(int), 0, output_size[0] - 1)

as

astype

只是删除小数，而不是四舍五入到更接近的整数。这与您的图像向左移动的事实是一致的。

您可以在打字之前尝试四舍五入，例如与

x_img = np.clip(np.round((x_normalized * (output_size[1]))), 0, output_size[0]-1).astype(int)
y_img = np.clip(np.round((y_normalized * (output_size[0]))), 0, output_size[1]-1).astype(int)