寫在前面

跨域描述符LCD可以實(shí)現(xiàn)二維圖片特征點(diǎn)到三維點(diǎn)云特征點(diǎn)的配準(zhǔn)，是個(gè)具有通用性的深度學(xué)習(xí)特征描述子。（圖片來源于論文LCD: Learned Cross-Domain Descriptors for 2D-3D Matching）
在Github開源的源碼里面給出了利用LCD進(jìn)行三維點(diǎn)云配準(zhǔn)的例程。align_point_cloud.py，這里對例程如何使用已經(jīng)訓(xùn)練好的模型來進(jìn)行三維點(diǎn)云配準(zhǔn)進(jìn)行解析。

運(yùn)行環(huán)境

python版本3.6.0以上
pytorch非CPU版本（可選）
Open3D
numpy及其它庫，自行下載
需要注意的是，官方的源碼中使用的Open3D版本較舊，在運(yùn)行程序時(shí)回出現(xiàn)新版本對應(yīng)函數(shù)不匹配的報(bào)錯(cuò)，詳情可以參考我之前發(fā)的解決辦法。

順帶一提，例程是不能在pycharm里面直接運(yùn)行的，需要輸入?yún)?shù)，具體是啥參照官方Github里面給的：

Aligning two point clouds with LCD This demo aligns two 3D colored
point clouds using our pre-trained LCD descriptor with RANSAC. How to
run:

$ python -m apps.align_point_cloud samples/000.ply samples/002.ply --logdir logs/LCD-D256/
For more information, use the --help option.

在代碼的根目錄下運(yùn)行python -m apps.align_point_cloud samples/000.ply samples/002.ply --logdir logs/LCD-D256/即可運(yùn)行調(diào)試。在pycharm里面需要帶參數(shù)運(yùn)行，具體方法自行搜索。

正文

1.庫

把該用的庫引一引：

import os
import json
import open3d
import torch
import argparse
import numpy as np

from lcd.models import *

但是pycharm提示我這個(gè)例程里面似乎并沒有用到pytorch，刪了也能正常運(yùn)行。

這部分里面的from lcd.models import * 對應(yīng)的lcd.models文件夾下有
patchnet和pointnet兩種分別對應(yīng)了二維和三維描述子生成網(wǎng)絡(luò)，例程中使用了pointnet下的PointNetAutoencoder類，之后會展開說說。

2.參數(shù)輸入和模型建立

parser = argparse.ArgumentParser()
parser.add_argument("source", help="path to the source point cloud")
parser.add_argument("target", help="path to the target point cloud")
parser.add_argument("--logdir", help="path to the log directory")
parser.add_argument("--voxel_size", default=0.1, type=float)
parser.add_argument("--radius", default=0.15, type=float)
parser.add_argument("--num_points", default=1024, type=int)
args = parser.parse_args()

(1) source → 輸入點(diǎn)云路徑
(2) target → 目標(biāo)點(diǎn)云路徑
(3) --logdir → 模型文件路徑
(4) --voxel_size → 體素大小，默認(rèn)為0.1
(5) --radius → 最近鄰搜索半徑，默認(rèn)為0.15
(6) --num_points → 采樣點(diǎn)對應(yīng)patch內(nèi)特征點(diǎn)個(gè)數(shù)，默認(rèn)為1024

得到輸入args為：

接著是讀取網(wǎng)絡(luò)配置文件config.json

logdir = args.logdir
config = os.path.join(logdir, "config.json")
config = json.load(open(config))

device = config["device"]

logdir 是 Log/LCD-D256目錄，保存了訓(xùn)練好的模型和配置等：

接下來

fname = os.path.join(logdir, "model.pth")
print("> Loading model from {}".format(fname))
model = PointNetAutoencoder(
    config["embedding_size"],
    config["input_channels"],
    config["output_channels"],
    config["normalize"],
)
model.load_state_dict(torch.load(fname)["pointnet"])
model.to(device)
model.eval()

fname為模型文件model.pth的路徑
我們來看一下 PointNetAutoencoder類的輸入
程序這邊model直接和config保持了一致。
在獲得配置文件后進(jìn)行模型加載：

model.load_state_dict(torch.load(fname)["pointnet"])

fname路徑對應(yīng)的model文件里面有 pointnet 和 patchnet 兩批參數(shù)，這邊因?yàn)橹挥悬c(diǎn)云的網(wǎng)絡(luò)所以只需要讀取[“pointnet”]。

model.to(device)

將模型加載在配置讀取出來的設(shè)備上（GPU）

model.eval()

eval()挺重要的，在模型預(yù)測階段我們需要Dropout層和batch normalization層設(shè)置到預(yù)測模式。

3.extract_uniform_patches

定義函數(shù)extract_uniform_patches，獲得點(diǎn)云Patch。

def extract_uniform_patches(pcd, voxel_size, radius, num_points):
    kdtree = open3d.geometry.KDTreeFlann(pcd)
    downsampled = pcd.voxel_down_sample(voxel_size)
    points = np.asarray(downsampled.points)
    patches = []
    for i in range(points.shape[0]):
        k, index, _ = kdtree.search_hybrid_vector_3d(points[i], radius, num_points)
        if k < num_points:
            index = np.random.choice(index, num_points, replace=True)
        xyz = np.asarray(pcd.points)[index]
        rgb = np.asarray(pcd.colors)[index]
        xyz = (xyz - points[i]) / radius  # normalize to local coordinates
        patch = np.concatenate([xyz, rgb], axis=1)
        patches.append(patch)
    patches = np.stack(patches, axis=0)
    return downsampled, patches

輸入pcd, voxel_size, radius, num_points

kdtree = open3d.geometry.KDTreeFlann(pcd)
downsampled = pcd.voxel_down_sample(voxel_size)
points = np.asarray(downsampled.points)

kdtree一會最近鄰搜索用的；
對點(diǎn)云進(jìn)行體素降采樣得到downsampled，是一個(gè)Open3D點(diǎn)云對象;
得到采樣點(diǎn)points。

 patches = []
    for i in range(points.shape[0]):
        k, index, _ = kdtree.search_hybrid_vector_3d(points[i], radius, num_points)
        if k < num_points:
            index = np.random.choice(index, num_points, replace=True)
        xyz = np.asarray(pcd.points)[index]
        rgb = np.asarray(pcd.colors)[index]
        xyz = (xyz - points[i]) / radius  # normalize to local coordinates
        patch = np.concatenate([xyz, rgb], axis=1)
        patches.append(patch)
    patches = np.stack(patches, axis=0)
    return downsampled, patche

最近鄰搜索kdtree.search_hybrid_vector_3d對每個(gè)降采樣之后的點(diǎn)進(jìn)行最近鄰搜索。
返回一個(gè)Tuple[int， open3d.utility.IntVector， open3d.utility.DoubleVector]，程序中的K為搜索點(diǎn)的個(gè)數(shù)，index為搜索點(diǎn)對應(yīng)的索引，是個(gè)一維數(shù)組。

后面的

 if k < num_points:
            index = np.random.choice(index, num_points, replace=True)

是當(dāng)k的值小于給定目標(biāo)點(diǎn)的個(gè)數(shù)(1024)時(shí)將進(jìn)行隨機(jī)填充，將index填充到目標(biāo)長度。

xyz = np.asarray(pcd.points)[index]
rgb = np.asarray(pcd.colors)[index]

取xyz為采樣點(diǎn)的位置信息，rgb為采樣點(diǎn)的顏色信息。

xyz = (xyz - points[i]) / radius

點(diǎn)云位置歸一化處理。

		patch = np.concatenate([xyz, rgb], axis=1)
        patches.append(patch)
    patches = np.stack(patches, axis=0)

將位置和顏色信息合并，對所有采樣點(diǎn)處理完畢后將列表patches堆疊為一個(gè)[points.shape[0], num_points, 6]（例程中為[1469, 1024, 6]）維度的矩陣。

return downsampled, patches

返回值為Open3D點(diǎn)云對象downsampled和patches?？梢钥闯龊瘮?shù)功能就是將輸入的點(diǎn)云降采樣并得到所有采樣點(diǎn)對應(yīng)的patch。

4.compute_lcd_descriptors

計(jì)算LCD特征描述子

def compute_lcd_descriptors(patches, model, batch_size, device):
    batches = torch.tensor(patches, dtype=torch.float32)
    batches = torch.split(batches, batch_size)
    descriptors = []
    with torch.no_grad():
        for i, x in enumerate(batches):
            x = x.to(device)
            z = model.encode(x)
            z = z.cpu().numpy()
            descriptors.append(z)
    return np.concatenate(descriptors, axis=0)

輸入為生成的patch和model以及選擇參數(shù)batch_size（例程中為128）、device（例程中為CUDA）。

batches = torch.tensor(patches, dtype=torch.float32)
batches = torch.split(batches, batch_size)

將輸入的patch轉(zhuǎn)換為張量并將生成的張量按照batch_size進(jìn)行分割。例程中的batch為若干(128, 1024, 6)的張量塊。

	with torch.no_grad():
        for i, x in enumerate(batches):
            x = x.to(device)
            z = model.encode(x)
            z = z.cpu().numpy()
            descriptors.append(z)

只是想看一下訓(xùn)練的效果，并不是想通過驗(yàn)證集來更新網(wǎng)絡(luò)時(shí)，可以使用with torch.no_grad()。
對于z = model.encode(x)，在類PointNetAutoencoder中有：

    def encode(self, x):
        z = self.encoder(x)
        if self.normalize:
            z = F.normalize(z)
        return z

然后pointnet.py前面有寫

self.normalize = normalize
self.input_channels = input_channels
self.embedding_size = embedding_size
self.encoder = PointNetEncoder(embedding_size, input_channels)

其中類PointNetEncoder如下：

class PointNetEncoder(nn.Module):
    def __init__(self, embedding_size, input_channels=3):
        super(PointNetEncoder, self).__init__()
        self.input_channels = input_channels
        self.stn1 = STN3D(input_channels)
        self.stn2 = STN3D(64)
        self.mlp1 = nn.Sequential(
            nn.Conv1d(input_channels, 64, 1),
            nn.BatchNorm1d(64),
            nn.ReLU(),
            nn.Conv1d(64, 64, 1),
            nn.BatchNorm1d(64),
            nn.ReLU(),
        )
        self.mlp2 = nn.Sequential(
            nn.Conv1d(64, 64, 1),
            nn.BatchNorm1d(64),
            nn.ReLU(),
            nn.Conv1d(64, 128, 1),
            nn.BatchNorm1d(128),
            nn.ReLU(),
            nn.Conv1d(128, 1024, 1),
            nn.BatchNorm1d(1024),
            nn.ReLU(),
        )
        self.fc = nn.Linear(1024, embedding_size)

正是訓(xùn)練得到的三維點(diǎn)云編碼器模型,得到一個(gè)CUDA tensor格式的數(shù)據(jù)、

z = z.cpu().numpy()

先將其轉(zhuǎn)換成cpu float-tensor隨后再轉(zhuǎn)到numpy格式。 numpy不能讀取CUDA tensor 需要將它轉(zhuǎn)化為 CPU tensor，得到一個(gè)(128,256)的描述子矩陣。

return np.concatenate(descriptors, axis=0)

將得到的所有描述子堆疊后輸出，例程中的輸出格式為(1469，256)。

5.點(diǎn)云讀取與處理

source = open3d.io.read_point_cloud(args.source)
target = open3d.io.read_point_cloud(args.target)

讀取源點(diǎn)云和目標(biāo)點(diǎn)云。

source_points, source_patches = extract_uniform_patches(
    source, args.voxel_size, args.radius, args.num_points)

通過exract_uniform_patches得到降采樣后的源點(diǎn)云和其對應(yīng)的patch。

source_descriptors = compute_lcd_descriptors(
    source_patches, model, batch_size=128, device=device)

通過compute_lcd_descriptors得到源點(diǎn)云的特征描述子source_descriptors。

source_features = open3d.pipelines.registration.Feature()
source_features.data = np.transpose(source_descriptors)  # T
print("Extracted {} features from source".format(len(source_descriptors)))

open3d.pipelines.registration.Feature()在官方文檔中的解釋如下：

大概可以理解為在Open3D的piplines里面用于儲存用于配準(zhǔn)的特征數(shù)據(jù)的類，需要將生成的特征描述子轉(zhuǎn)置后賦給source_features中的data——source_features.data = np.transpose(source_descriptors)

目標(biāo)點(diǎn)云的處理方式一致：

target_points, target_patches = extract_uniform_patches(
    target, args.voxel_size, args.radius, args.num_points
)
target_descriptors = compute_lcd_descriptors(
    target_patches, model, batch_size=128, device=device
)
target_features = open3d.pipelines.registration.Feature()
target_features.data = np.transpose(target_descriptors)
print("Extracted {} features from target".format(len(target_descriptors)))

6. 特征點(diǎn)配準(zhǔn)

threshold = 0.075
result = open3d.pipelines.registration.registration_ransac_based_on_feature_matching(
    source_points,
    target_points,
    source_features,
    target_features,
    True,
    threshold,
    open3d.pipelines.registration.TransformationEstimationPointToPoint(False),
    4,
    [open3d.pipelines.registration.CorrespondenceCheckerBasedOnDistance(threshold)],
    open3d.pipelines.registration.RANSACConvergenceCriteria(4000000, 500),
)

注意：這部分的代碼是在新版本Open3D的基礎(chǔ)上修改的，并非原代碼。
修改原因見我之前發(fā)的Open3D 15.1 報(bào)錯(cuò) module ‘open3d‘ has no attribute ‘registration‘（跑LCD代碼時(shí)報(bào)錯(cuò)）
open3d.pipelines.registration.registration_ransac_based_on_feature_matching函數(shù)新版本官方給出的解釋如下：

對兩點(diǎn)云進(jìn)行特征配準(zhǔn)得到結(jié)果result,

if result.transformation.trace() == 4.0:
    success = False

當(dāng)變換矩陣的跡為4時(shí)（沒有產(chǎn)生旋轉(zhuǎn)）認(rèn)定未完成配準(zhǔn)。

information = open3d.pipelines.registration.get_information_matrix_from_point_clouds(
    source_points, target_points, threshold, result.transformation
)
n = min(len(source_points.points), len(target_points.points))
if (information[5, 5] / n) < 0.3:  # overlap threshold
    success = False

通過獲得信息矩陣來判定是否配準(zhǔn)成功。

if not success:
    print("Cannot align two point clouds")
    exit(0)

如果配準(zhǔn)未成功，直接結(jié)束程序。

print("Success!")
print("Visualizing alignment result...")
source.estimate_normals(open3d.geometry.KDTreeSearchParamHybrid(radius=0.2, max_nn=30))
target.estimate_normals(open3d.geometry.KDTreeSearchParamHybrid(radius=0.2, max_nn=30))
source.paint_uniform_color([1, 0.706, 0])
target.paint_uniform_color([0, 0.651, 0.929])
source.transform(result.transformation)
open3d.visualization.draw_geometries([source, target])
print(result.transformation)

如果配準(zhǔn)成功完成，將配準(zhǔn)后的點(diǎn)云變換后著色進(jìn)行可視化輸出得到如下結(jié)果：
文章來源地址http://www.zghlxwxcb.cn/news/detail-481623.html

到了這里，關(guān)于[點(diǎn)云配準(zhǔn)]LCD(2D-3D特征配準(zhǔn)算法)例程align_point_cloud.py解析的文章就介紹完了。如果您還想了解更多內(nèi)容，請?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！