云计算百科
云计算领域专业知识百科平台

点云配准之-ColorICP中公式详细推导

相关内容: 1、点云配准之-Colored Point Cloud Registration(ColorICP)

一、计算梯度

计算梯度

d

p

d_p

dp 的最小二乘目标函数为:

L

(

d

p

)

=

p

N

p

(

C

p

(

f

(

p

)

p

)

C

(

p

)

)

2

p

N

p

(

C

(

p

)

+

d

p

(

f

(

p

)

p

)

C

(

p

)

)

2

,

(10)

L(d_p) = \\sum_{p' \\in N_p} \\big(C_p(f(p') – p) – C(p')\\big)^2 \\approx \\sum_{p' \\in N_p} \\big(C(p) + d_p^\\top(f(p') – p) – C(p')\\big)^2, \\tag{10}

L(dp)=pNp(Cp(f(p)p)C(p))2pNp(C(p)+dp(f(p)p)C(p))2,(10)

同时约束

d

p

n

p

=

0

d_p^\\top n_p = 0

dpnp=0

实现代码片段如下:

for (int i = 0; i < frame::size(*frame); i++) {
// Estimate color gradient
const auto& point = frame::point(*frame, i);
const auto& normal = frame::normal(*frame, i);
const double intensity = frame::intensity(*frame, i);

Eigen::Matrix<double, 1, 4> A = Eigen::Matrix<double, 1, 4>::Zero(k_photo_neighbors, 4);
Eigen::VectorXd b = Eigen::VectorXd::Zero(k_photo_neighbors);

// dp^T np = 0
A.row(0) = normal;
b[0] = 0.0;

// Intensity gradient in the tangent space
for (int j = 1; j < k_photo_neighbors; j++) {
const int index = neighbors[k * i + j];
const auto& point_ = frame::point(*frame, index);
const double intensity_ = frame::intensity(*frame, index);
const Eigen::Vector4d projected = point_ (point_ point).dot(normal) * normal;
A.row(j) = projected point;
b(j) = (intensity_ intensity);
}

Eigen::Matrix3d H = (A.transpose() * A).block<3, 3>(0, 0);
Eigen::Vector3d e = (A.transpose() * b).head<3>();
gradients->intensity_gradients[i] << H.inverse() * e, 0.0;
}

一、符号与背景先说明

上面公式10是 ColorICP中强度梯度估计或光度一致性建模的标准形式。ColorICP

1 点与邻域

  • p

    R

    3

    p \\in \\mathbb{R}^3

    pR3:当前点

  • N

    p

    N_p

    Np:点

    p

    p

    p 的邻域

  • p

    N

    p

    p' \\in N_p

    pNp:邻居点


2 强度函数

  • C

    (

    )

    C(\\cdot)

    C():强度(Intensity / Color)函数

    • LiDAR:反射强度
    • RGB-D:像素亮度

3 映射函数

f

(

)

f(\\cdot)

f()

  • f

    (

    p

    )

    f(p')

    f(p):将邻居点映射到 以

    p

    p

    p 为参考的局部坐标系

  • 在实现的代码中,本质是:

f

(

p

)

邻居点在切平面上的投影

f(p') \\approx \\text{邻居点在切平面上的投影}

f(p)邻居点在切平面上的投影


4 梯度变量

  • d

    p

    R

    3

    d_p \\in \\mathbb{R}^3

    dpR3:点

    p

    p

    p 处的强度梯度

d

p

=

C

(

p

)

d_p = \\nabla C(p)

dp=C(p)


二、原始代价函数(非线性)

给出的第一行是真实但不可直接解的形式:

L

(

d

p

)

=

p

N

p

(

C

p

(

f

(

p

)

p

)

C

(

p

)

)

2

L(d_p)= \\sum_{p' \\in N_p} \\Big( C_p(f(p') – p) – C(p') \\Big)^2

L(dp)=pNp(Cp(f(p)p)C(p))2

解释:

  • C

    p

    (

    )

    C_p(\\cdot)

    Cp():在点

    p

    p

    p 处定义的局部强度函数

  • 希望:

    p

    p

    p 处的局部模型,去拟合邻居点

    p

    p'

    p 的强度


三、一阶泰勒展开(关键步骤)

目标

近似:

C

p

(

f

(

p

)

p

)

C_p(f(p') – p)

Cp(f(p)p)


p

p

p 处做一阶泰勒展开

C

(

p

+

Δ

p

)

C

(

p

)

+

C

(

p

)

Δ

p

C(p + \\Delta p) \\approx C(p) + \\nabla C(p)^\\top \\Delta p

C(p+Δp)C(p)+C(p)Δp

其中:

Δ

p

=

f

(

p

)

p

\\Delta p = f(p') – p

Δp=f(p)p


代入得到

C

p

(

f

(

p

)

p

)

C

(

p

)

+

d

p

(

f

(

p

)

p

)

C_p(f(p') – p) \\approx C(p) + d_p^\\top (f(p') – p)

Cp(f(p)p)C(p)+dp(f(p)p)

其中:

d

p

=

C

(

p

)

d_p = \\nabla C(p)

dp=C(p)


四、代入原始代价函数

将近似结果代入原式:

L

(

d

p

)

=

p

N

p

(

C

p

(

f

(

p

)

p

)

C

(

p

)

)

2

 

p

N

p

(

C

(

p

)

+

d

p

(

f

(

p

)

p

)

C

(

p

)

)

2

\\begin{aligned} L(d_p) &= \\sum_{p' \\in N_p} \\Big( C_p(f(p') – p) – C(p') \\Big)^2 \\ &\\approx \\sum_{p' \\in N_p} \\Big( C(p) + d_p^\\top (f(p') – p) – C(p') \\Big)^2 \\end{aligned}

L(dp)=pNp(Cp(f(p)p)C(p))2 pNp(C(p)+dp(f(p)p)C(p))2

这正是原文中给出的 公式 (10)。

推导完成


五、把公式 (10) 改写成标准最小二乘形式

这是理解代码和实现的关键。


定义残差

对每个邻点

p

p'

p

r

p

(

d

p

)

=

C

(

p

)

+

d

p

(

f

(

p

)

p

)

C

(

p

)

r_{p'}(d_p)= C(p) + d_p^\\top (f(p') – p) – C(p')

rp(dp)=C(p)+dp(f(p)p)C(p)


堆叠成线性系统

记:

A

p

=

(

f

(

p

)

p

)

R

1

×

3

A_{p'} = (f(p') – p)^\\top \\in \\mathbb{R}^{1 \\times 3}

Ap=(f(p)p)R1×3

b

p

=

C

(

p

)

C

(

p

)

b_{p'} = C(p') – C(p)

bp=C(p)C(p)

则:

r

p

=

A

p

d

p

b

p

r_{p'} = A_{p'} d_p – b_{p'}

rp=Apdpbp


最小二乘问题

min

d

p

p

N

p

A

p

d

p

b

p

2

\\min_{d_p} \\sum_{p' \\in N_p} \\left| A_{p'} d_p – b_{p'} \\right|^2

dpminpNpApdpbp2


写成矩阵形式

min

d

p

A

d

p

b

2

\\min_{d_p} \\left| A d_p – b \\right|^2

dpminAdpb2


正规方程

A

A

d

p

=

A

b

A^\\top A d_p = A^\\top b

AAdp=Ab


六、切平面约束

在实现的代码和论文中,一定还有这一条:

n

p

d

p

=

0

n_p^\\top d_p = 0

npdp=0

原因

  • 强度变化 只在表面切平面内有意义
  • 法向方向是深度方向,光度变化不可靠

合并约束后的系统

[

n

p

 

(

f

(

p

1

)

p

)

]

d

p

=

[

0

C

(

p

1

)

C

(

p

)

]

\\begin{bmatrix} n_p^\\top \\ (f(p_1) – p)^\\top \\\\ \\vdots \\end{bmatrix} d_p= \\begin{bmatrix} 0 \\\\ C(p_1) – C(p) \\\\ \\vdots \\end{bmatrix}

[np (f(p1)p)]dp=

0C(p1)C(p)

对应代码中:

A.row(0) = normal;
b(0) = 0;


七、代码中的一一对应关系

数学代码

d

p

d_p

dp

intensity_gradient[i].head<3>()

f

(

p

)

p

f(p') – p

f(p)p

projected – point

C

(

p

)

C

(

p

)

C(p') – C(p)

C(p)C(p)

intensity_ – intensity
切平面约束 A.row(0) = normal

八、模型的物理意义

用邻域点的强度差,在线性近似下拟合出“表面上的强度变化方向”,这就是光度梯度。


二、光度雅克比计算

原文对公式如下:

r

C

(

p

,

q

)

(

T

)

=

(

C

p

f

s

)

ξ

i

=

C

p

(

f

)

,

J

f

(

s

)

,

J

s

(

ξ

)

(28–29)

\\nabla r_C^{(p,q)}(T) = \\frac{\\partial (C_p \\circ f \\circ s)}{\\partial \\xi_i} = \\nabla C_p(f), J_f(s), J_s(\\xi) \\tag{28–29}

rC(p,q)(T)=ξi(Cpfs)=Cp(f),Jf(s),Js(ξ)(28–29)


一、问题回顾与符号统一

光度残差定义为(公式 18):

r

C

(

p

,

q

)

(

T

)

=

C

p

(

f

(

s

(

q

,

T

)

)

)

C

(

q

)

r_C^{(p,q)}(T) = C_p\\big(f(s(q,T))\\big) – C(q)

rC(p,q)(T)=Cp(f(s(q,T)))C(q)

其中:

  • q

    Q

    q \\in Q

    qQ:源点云中的点

  • p

    P

    p \\in P

    pP:目标点云中的对应点

  • s

    (

    q

    ,

    T

    )

    s(q,T)

    s(q,T):刚性变换,将

    q

    q

    q 变换到

    P

    P

    P 坐标系

  • f

    (

    )

    f(\\cdot)

    f():将 3D 点投影到

    p

    p

    p 的切平面

  • C

    p

    (

    )

    C_p(\\cdot)

    Cp():定义在切平面上的连续颜色函数

  • C

    (

    q

    )

    C(q)

    C(q):常数,对

    T

    T

    T 无关

优化变量只有

T

T

T(或其最小参数

ξ

\\xi

ξ

因此求导时:

r

C

(

p

,

q

)

(

T

)

=

ξ

C

p

(

f

(

s

(

q

,

T

)

)

)

\\nabla r_C^{(p,q)}(T) = \\frac{\\partial}{\\partial \\xi}C_p\\big(f(s(q,T))\\big)

rC(p,q)(T)=ξCp(f(s(q,T)))


二、函数复合结构(关键)

我们先把残差写成复合函数形式:

r

C

(

p

,

q

)

(

T

)

=

(

C

p

f

s

)

(

q

,

T

)

r_C^{(p,q)}(T) = (C_p \\circ f \\circ s)(q,T)

rC(p,q)(T)=(Cpfs)(q,T)

对应的函数链是:

ξ

s

s

(

q

,

ξ

)

f

f

(

s

(

q

,

ξ

)

)

C

p

C

p

(

f

(

s

(

q

,

ξ

)

)

)

\\xi \\xrightarrow{ s } s(q,\\xi) \\xrightarrow{ f } f(s(q,\\xi)) \\xrightarrow{ C_p } C_p(f(s(q,\\xi)))

ξs

s(q,ξ)f

f(s(q,ξ))Cp

Cp(f(s(q,ξ)))

这是一个典型的三层复合函数,因此直接使用链式法则。


三、第一步:最外层求导(颜色函数)

1 局部颜色函数的一阶形式

在 上 节中,

C

p

(

u

)

C_p(u)

Cp(u) 被一阶近似为:

C

p

(

u

)

C

(

p

)

+

d

p

u

C_p(u) \\approx C(p) + d_p^\\top u

Cp(u)C(p)+dpu

因此其梯度为常量向量:

C

p

(

u

)

=

d

p

R

1

×

2

\\nabla C_p(u) = d_p \\in \\mathbb{R}^{1 \\times 2}

Cp(u)=dpR1×2

注意:

u

u

u 是切平面坐标(2D),但嵌入在 3D 空间中


2 对中间变量求导

设:

y

=

f

(

s

(

q

,

ξ

)

)

y = f(s(q,\\xi))

y=f(s(q,ξ))

则:

y

C

p

(

y

)

=

C

p

(

y

)

=

d

p

\\frac{\\partial}{\\partial y} C_p(y) = \\nabla C_p(y) = d_p

yCp(y)=Cp(y)=dp


四、第二步:切平面投影函数

f

f

f 的 Jacobian

投影函数回顾(公式 9)

f

(

s

)

=

s

n

p

(

s

p

)

n

p

f(s) = s – n_p (s – p)^\\top n_p

f(s)=snp(sp)np

整理为:

f

(

s

)

=

(

I

n

p

n

p

)

s

+

n

p

n

p

p

f(s) = \\big(I – n_p n_p^\\top\\big)s + n_p n_p^\\top p

f(s)=(Inpnp)s+npnpp

其中:

  • n

    p

    n_p

    np 是常量(点

    p

    p

    p 的法向)

  • p

    p

    p 是常量


Jacobian

J

f

(

s

)

J_f(s)

Jf(s)

s

s

s 求导:

J

f

(

s

)

=

f

s

=

I

n

p

n

p

J_f(s) = \\frac{\\partial f}{\\partial s} = I – n_p n_p^\\top

Jf(s)=sf=Inpnp

这是一个

3

×

3

3\\times3

3×3 的投影矩阵,将向量投影到切平面。


五、第三步:刚性变换

s

(

q

,

ξ

)

s(q,\\xi)

s(q,ξ) 的 Jacobian

1 刚性变换形式

刚性变换写作:

s

(

q

,

ξ

)

=

R

(

ξ

)

q

+

t

(

ξ

)

s(q,\\xi) = R(\\xi)q + t(\\xi)

s(q,ξ)=R(ξ)q+t(ξ)

在 Gauss–Newton 中,

R

(

ξ

)

R(\\xi)

R(ξ)

T

k

T_k

Tk 附近线性化:

R

(

ξ

)

I

+

[

ω

]

×

R(\\xi) \\approx I + [\\omega]_\\times

R(ξ)I+[ω]×

因此:

s

(

q

,

ξ

)

q

+

ω

×

q

+

t

s(q,\\xi) \\approx q + \\omega \\times q + t

s(q,ξ)q+ω×q+t


2 对

ξ

=

(

α

,

β

,

γ

,

a

,

b

,

c

)

\\xi = (\\alpha,\\beta,\\gamma,a,b,c)

ξ=(α,β,γ,a,b,c) 求导

Jacobian 为:

J

s

(

ξ

)

=

s

(

q

,

ξ

)

ξ

[

0

q

z

q

y

1

0

0

q

z

0

q

x

0

1

0

q

y

q

x

0

0

0

1

]

J_s(\\xi) = \\frac{\\partial s(q,\\xi)}{\\partial \\xi} \\begin{bmatrix} 0 & q_z & -q_y & 1 & 0 & 0 \\\\ -q_z & 0 & q_x & 0 & 1 & 0 \\\\ q_y & -q_x & 0 & 0 & 0 & 1 \\end{bmatrix}

Js(ξ)=ξs(q,ξ)

0qzqyqz0qxqyqx0100010001

这是标准的 SE(3) 点 Jacobian。


六、完整链式法则推导

将三部分串起来:

r

C

(

p

,

q

)

(

T

)

=

C

p

f

f

s

s

ξ

=

C

p

(

f

)

J

f

(

s

)

J

s

(

ξ

)

\\begin{aligned} \\nabla r_C^{(p,q)}(T) &= \\frac{\\partial C_p}{\\partial f} \\frac{\\partial f}{\\partial s} \\frac{\\partial s}{\\partial \\xi} &= \\nabla C_p(f) J_f(s) J_s(\\xi) \\end{aligned}

rC(p,q)(T)=fCpsfξs=Cp(f)Jf(s)Js(ξ)

即得到论文中的公式:

r

C

(

p

,

q

)

(

T

)

=

C

p

(

f

)

J

f

(

s

)

J

s

(

ξ

)

\\boxed{ \\nabla r_C^{(p,q)}(T) = \\nabla C_p(f) J_f(s)J_s(\\xi) }

rC(p,q)(T)=Cp(f)Jf(s)Js(ξ)


七、维度检查

项维度

C

p

(

f

)

=

d

p

\\nabla C_p(f)=d_p

Cp(f)=dp

1

×

3

1 \\times 3

1×3

J

f

(

s

)

J_f(s)

Jf(s)

3

×

3

3 \\times 3

3×3

J

s

(

ξ

)

J_s(\\xi)

Js(ξ)

3

×

6

3 \\times 6

3×6

最终 Jacobian

1

×

6

1 \\times 6

1×6

完全符合 Gauss-Newton 所需的残差 Jacobian。


八、总结

光度残差对位姿的导数 = 颜色变化方向 × 切平面投影 × 位姿扰动对点位置的影响


九、工程实现提示

在代码中通常这样实现:

// dp: 1×3
// Jf: 3×3 = I – np*np^T
// Js: 3×6
JrC = dp * Jf * Js;

其中:

  • dp:预处理阶段计算
  • Jf:每个 correspondence 固定
  • Js:依赖当前点 q

三、几何雅克比计算

原文对应公式如下:

r

G

(

p

,

q

)

(

T

)

=

n

p

J

s

(

ξ

)

(30)

\\nabla r_G^{(p,q)}(T) = n_p^\\top J_s(\\xi) \\tag{30}

rG(p,q)(T)=npJs(ξ)(30)


一、从几何残差的定义出发

几何残差定义(公式 19)为:

r

G

(

p

,

q

)

(

T

)

=

(

s

(

q

,

T

)

p

)

n

p

r_G^{(p,q)}(T) = \\big(s(q,T) – p\\big)^\\top n_p

rG(p,q)(T)=(s(q,T)p)np

其中:

  • s

    (

    q

    ,

    T

    )

    s(q,T)

    s(q,T):将点

    q

    q

    q 通过刚性变换

    T

    T

    T 变换到

    P

    P

    P 坐标系

  • p

    p

    p:目标点(常量)

  • n

    p

    n_p

    np:点

    p

    p

    p 的法向量(常量)

优化变量只有

T

T

T(或其最小参数

ξ

\\xi

ξ


二、去掉与优化无关的常量项

将残差展开:

r

G

(

p

,

q

)

(

T

)

=

n

p

s

(

q

,

T

)

n

p

p

\\begin{aligned} r_G^{(p,q)}(T) &= n_p^\\top s(q,T) – n_p^\\top p \\end{aligned}

rG(p,q)(T)=nps(q,T)npp

由于

p

p

p

n

p

n_p

np 都是常量:

ξ

(

n

p

p

)

=

0

\\frac{\\partial}{\\partial \\xi}\\big(n_p^\\top p\\big) = 0

ξ(npp)=0

因此:

r

G

(

p

,

q

)

(

T

)

=

ξ

(

n

p

s

(

q

,

T

)

)

\\nabla r_G^{(p,q)}(T) = \\frac{\\partial}{\\partial \\xi}\\big(n_p^\\top s(q,T)\\big)

rG(p,q)(T)=ξ(nps(q,T))


三、明确函数结构

这是一个线性形式 + 刚性变换的组合:

r

G

(

p

,

q

)

(

T

)

=

n

p

s

(

q

,

ξ

)

r_G^{(p,q)}(T) = n_p^\\top s(q,\\xi)

rG(p,q)(T)=nps(q,ξ)

可视为复合函数:

ξ

s

s

(

q

,

ξ

)

n

p

(

)

n

p

s

(

q

,

ξ

)

\\xi \\xrightarrow{ s } s(q,\\xi) \\xrightarrow{ n_p^\\top(\\cdot) } n_p^\\top s(q,\\xi)

ξs

s(q,ξ)np()

nps(q,ξ)


四、对刚性变换

s

(

q

,

ξ

)

s(q,\\xi)

s(q,ξ) 求导

1 刚性变换回顾

局部线性化后(公式 20):

s

(

q

,

ξ

)

R

(

ξ

)

q

+

t

(

ξ

)

s(q,\\xi) \\approx R(\\xi)q + t(\\xi)

s(q,ξ)R(ξ)q+t(ξ)

T

k

T_k

Tk 附近:

s

(

q

,

ξ

)

q

+

ω

×

q

+

t

s(q,\\xi) \\approx q + \\omega \\times q + t

s(q,ξ)q+ω×q+t

其中:

ξ

=

(

α

,

β

,

γ

,

a

,

b

,

c

)

\\xi = (\\alpha,\\beta,\\gamma,a,b,c)

ξ=(α,β,γ,a,b,c)


2 Jacobian

J

s

(

ξ

)

J_s(\\xi)

Js(ξ)

ξ

\\xi

ξ 求导得到:

J

s

(

ξ

)

=

s

(

q

,

ξ

)

ξ

[

0

q

z

q

y

1

0

0

q

z

0

q

x

0

1

0

q

y

q

x

0

0

0

1

]

J_s(\\xi) = \\frac{\\partial s(q,\\xi)}{\\partial \\xi} \\begin{bmatrix} 0 & q_z & -q_y & 1 & 0 & 0 \\\\ -q_z & 0 & q_x & 0 & 1 & 0 \\\\ q_y & -q_x & 0 & 0 & 0 & 1 \\end{bmatrix}

Js(ξ)=ξs(q,ξ)

0qzqyqz0qxqyqx0100010001

维度:

3

×

6

3 \\times 6

3×6


五、应用链式法则

从:

r

G

(

p

,

q

)

(

T

)

=

n

p

s

(

q

,

ξ

)

r_G^{(p,q)}(T) = n_p^\\top s(q,\\xi)

rG(p,q)(T)=nps(q,ξ)

直接对

ξ

\\xi

ξ 求导:

r

G

(

p

,

q

)

(

T

)

=

ξ

(

n

p

s

(

q

,

ξ

)

)

=

n

p

s

(

q

,

ξ

)

ξ

=

n

p

J

s

(

ξ

)

\\begin{aligned} \\nabla r_G^{(p,q)}(T) &= \\frac{\\partial}{\\partial \\xi} \\big(n_p^\\top s(q,\\xi)\\big) &= n_p^\\top \\frac{\\partial s(q,\\xi)}{\\partial \\xi} &= n_p^\\top J_s(\\xi) \\end{aligned}

rG(p,q)(T)=ξ(nps(q,ξ))=npξs(q,ξ)=npJs(ξ)


六、最终结果

r

G

(

p

,

q

)

(

T

)

=

n

p

J

s

(

ξ

)

\\boxed{ \\nabla r_G^{(p,q)}(T) = n_p^\\top J_s(\\xi) }

rG(p,q)(T)=npJs(ξ)

这正是论文中的公式 (30)。


七、维度检查

项维度

n

p

n_p^\\top

np

1

×

3

1 \\times 3

1×3

J

s

(

ξ

)

J_s(\\xi)

Js(ξ)

3

×

6

3 \\times 6

3×6

结果

1

×

6

1 \\times 6

1×6

符合 Gauss–Newton 对单残差的 Jacobian 要求。


八、几何意义

几何残差对位姿的敏感度 = 沿法向方向的位移变化 × 位姿扰动对点位置的影响

这正是 point-to-plane ICP 的核心思想。


九、与 point-to-plane ICP 的关系

当只保留几何项(

σ

=

1

\\sigma = 1

σ=1)时:

E

(

T

)

=

(

p

,

q

)

(

(

s

(

q

,

T

)

p

)

n

p

)

2

E(T) = \\sum_{(p,q)} \\big((s(q,T)-p)^\\top n_p\\big)^2

E(T)=(p,q)((s(q,T)p)np)2

其 Jacobian 正是:

J

=

n

p

J

s

(

ξ

)

J = n_p^\\top J_s(\\xi)

J=npJs(ξ)

** Colored ICP 在几何部分与经典 point-to-plane ICP 完全一致**。


十、工程实现示例

// np: Vector3d
// Js: Matrix<double,3,6>
RowVectorXd JrG = np.transpose() * Js;


赞(0)
未经允许不得转载:网硕互联帮助中心 » 点云配准之-ColorICP中公式详细推导
分享到: 更多 (0)

评论 抢沙发

评论前必须登录!