身处相机内外参之间(EG3D/NeRF/3D Gaussian Splatting)

文摘   2024-11-19 07:00   上海  

作者 | 八氨合氯化钙  编辑 | 自动驾驶之心

点击下方卡片,关注“3D视觉之心”公众号

第一时间获取3D视觉干货

>>点击进入→3D视觉之心技术交流群

此文记录个人学习过程,如有错误欢迎私信交流,十分感谢!

“相机系统是一个非常令我头疼的事情,我真的很不擅长这个。”

这篇blog是为了“解析”我遇到的一些仓库中的关于相机系统的代码,这么说可能会很奇怪,因为对于比较专业的人来说,这些根本不能算是问题;但由于我是一个只知道PyTorch的文盲,所以这对我来说,是问题,而且很大。

哪怕在这里语境下的“相机”只是最最简单的理想模型,但即使这样,想修改与其相关的代码对没有受过相关训练的人来说还是有些困难的。尤其当需要对相机系统做一些更“定制化”的操作时候,浅显且不够直接的理解就不太能支持继续推进下去了。

我曾将我的窘况诉说给一位学长,他听完非常震惊“不是吧这都不会”(此处自动脑补“虾头座椅电脑”表情包)。他建议我去回炉重造,去看GAMES101或者洋人的公开课。但我没看几分钟就睡着了,睡醒以后我不禁陷入思考:在那个田园的时代,人们或许是通过看一些关于3DV的课或者书从而掌握了其理论基础,然后可能在OpenCV,OpenGL,Unity等库里进行实践,由于那些库写的都很严谨,所以很快就可以将理论与实际匹配上,然后熟悉这一套东西。

但在2024年,在新时代图形学的背景下,我如果再这么做可能有点缘木求鱼了。所以有没有一种顺应新时代的掌握这些的方法呢?大部分DL+3D项目都是缝来缝去,作者可能自己都不知道我的coordinates convention是什么。所以我觉得更适合新时代的方法是头铁的去硬看那些代码,把这些都看明白了,自然就会了。

这篇blog以这样的组织进行:

  • 会先从EG3D的相机系统出发,这里提供了一种构造相机位姿矩阵(下文简称为c2w)的方法,同时可以让我们对NeRF-based的方法有个可视化的认识。在这个部分我们可以对c2w有一个直观的认识。
  • 然后我们会返回去解决初学NeRF-PyTorch时不太关心的那些相机代码。相当于对上一部分的练习,同时温习一下透视投影。
  • 之后我们会阅读PanoHead的预处理部分,其基于3DDFA_V2。这个预处理本质上和EG3D的预处理做的事情是一样的,只是这个基于3DDFA的预处理被PanoHead的作者整理的更自洽。在这套代码里,我们会正好遇到相机外参(下文简称为w2c)和c2w的关系。这一套预处理比较复杂,涉及正交投影和一些琐碎的数字图像处理。相当于某种程度的过关考核。
  • 再然后我们会解读3D Gaussian Splatting的代码。这个如今大热的显式表示可以让我们复习一下在NeRF里常常被忽略的透视投影。
  • 最后,我们会从更深入的角度下讨论在代码中经常出现的旋转操作,用李群和李代数的视角进行一些浅显的讨论,来从某种意义上“升华”一下我们的理解。

所以这篇blog至少需要读者大概了解上述工作,可能咿呀学语般的玩过其中一些的代码,大概能稀里糊涂的说出什么相机内外参,刚体变换之类的名词,对情况有大致的了解。(天哪这简直是我.jpg)

EG3D

在EG3D那一套代码里,我们经常可以看到c2w是这么被召唤出来的,先是有一个地方定义了一个:

camera_lookat_point = torch.tensor([0, 0, 0.2], device=self.device)

然后在需要c2w的地方,会有:

cam2world_pose = LookAtPoseSampler.sample(3.14/2 + 2 * 3.14 * frame_idx / num_frames, 3.14/2,
camera_lookat_point, radius=2.75, device=self.device)

类似这样的操作,然后就神奇的得到一个有效的4×4矩阵了。例如:

>>> cam2world = LookAtPoseSampler.sample(math.pi/4, math.pi/4, torch.tensor([0, 0, 0]), radius=2.5)
>>> tensor([[[ 0.7071, -0.3536,  0.6124, -1.5309],
             [ 0.0000, -0.8660, -0.5000,  1.2500],
             [ 0.7071,  0.3536, -0.6124,  1.5309],
             [ 0.0000,  0.0000,  0.0000,  1.0000]]])

现在我们关心LookAtPoseSampler.sample的实现,这个实现展示了一个很经典的用look, at, up构造相机外参的方法。EG3D中,整套代码是在这样的一个settings下:

这里用的是左手系,不过影响不大。整套代码的实现是:

class LookAtPoseSampler:
    """
    Same as GaussianCameraPoseSampler, except the
    camera is specified as looking at 'lookat_position', a 3-vector.

    Example:
    For a camera pose looking at the origin with the camera at position [0, 0, 1]:
    cam2world = LookAtPoseSampler.sample(math.pi/2, math.pi/2, torch.tensor([0, 0, 0]), radius=1)
    "
""

    @staticmethod
    def sample(horizontal_mean, vertical_mean, lookat_position, horizontal_stddev=0, vertical_stddev=0, radius=1, batch_size=1, device='cpu'):
        h = torch.randn((batch_size, 1), device=device) * horizontal_stddev + horizontal_mean
        v = torch.randn((batch_size, 1), device=device) * vertical_stddev + vertical_mean
        v = torch.clamp(v, 1e-5, math.pi - 1e-5)

        theta = h
        v = v / math.pi
        phi = torch.arccos(1 - 2*v)

        camera_origins = torch.zeros((batch_size, 3), device=device)

        camera_origins[:, 0:1] = radius*torch.sin(phi) * torch.cos(math.pi-theta)
        camera_origins[:, 2:3] = radius*torch.sin(phi) * torch.sin(math.pi-theta)
        camera_origins[:, 1:2] = radius*torch.cos(phi)

        # forward_vectors = math_utils.normalize_vecs(-camera_origins)
        forward_vectors = math_utils.normalize_vecs(lookat_position - camera_origins)
        return create_cam2world_matrix(forward_vectors, camera_origins)

最开始,是先构造出球坐标系下的。然后,我们会根据半径r给出相机的位置camera_origins

由于与标准的右手系下的球坐标系不太一样,所以这里变成了-。之后,我们会计算一个forward_vector,这个其实就是LookAt。通过将输入的lookat_positioncamera_origins的差归一化,我们会得到一个从相机位置指向look at位置的一个方向。然后就进入了create_cam2world_matrix(forward_vectors, camera_origins)

def create_cam2world_matrix(forward_vector, origin):
    """
    Takes in the direction the camera is pointing and the camera origin and returns a cam2world matrix.
    Works on batches of forward_vectors, origins. Assumes y-axis is up and that there is no camera roll.
    "
""

    forward_vector = math_utils.normalize_vecs(forward_vector)
    up_vector = torch.tensor([0, 1, 0], dtype=torch.float, device=origin.device).expand_as(forward_vector)

    right_vector = -math_utils.normalize_vecs(torch.cross(up_vector, forward_vector, dim=-1))
    up_vector = math_utils.normalize_vecs(torch.cross(forward_vector, right_vector, dim=-1))

    rotation_matrix = torch.eye(4, device=origin.device).unsqueeze(0).repeat(forward_vector.shape[0], 1, 1)
    rotation_matrix[:, :3, :3] = torch.stack((right_vector, up_vector, forward_vector), axis=-1)

    translation_matrix = torch.eye(4, device=origin.device).unsqueeze(0).repeat(forward_vector.shape[0], 1, 1)
    translation_matrix[:, :3, 3] = origin
    cam2world = (translation_matrix @ rotation_matrix)[:, :, :]
    assert(cam2world.shape[1:] == (4, 4))
    return cam2world

这个矩阵很直接的描述了如何从相机坐标系下的某个点,变换回世界坐标系。我们可以带入[0.0.0.1],[1.0.0.1]这两个特殊的点:

相机外参实际上是这个矩阵的逆,即w2c。从坐标系变换来理解c2w和w2c,总是有些困难的。更倾向于用这种方法理解了c2w,然后将相机外参当作这个矩阵的逆就好了。可以做一些简单的可视化来加深印象,下图画一个对于人头前半部分和整个部分的相机:

上述脚本中设置了lookat_position为[0, 0, 0.2],所以可以看到这些相机的延长线汇聚到了z轴上方的一个点。

然后是关于相机内参的事情,在EG3D里相机内参是在预处理阶段直接给定的。其本身没有什么特殊的,只是选取了一个视场角比较小的内参,在EG3D的代码中,可以看到:

intrinsics = torch.tensor([[4.2647, 0, 0.5], [0, 4.2647, 0.5], [0, 0, 1]], device=device)

上面的4.2647是归一化后的内参:

所以:

内参在这套代码里的目的,是用来决定发射光线的方向的。具体来说是起到“视口变换”的逆变换,即从归一化平面坐标系转移到相机坐标系:

def raysampler(cam2world_matrix, intrinsics, resolution=16):
    """
    Create batches of rays and return origins and directions.

    cam2world_matrix: (N, 4, 4)
    intrinsics: (N, 3, 3)
    resolution: int

    ray_origins: (N, M, 3)
    ray_dirs: (N, M, 2)
    "
""
    N, M = cam2world_matrix.shape[0], resolution**2
    cam_locs_world = cam2world_matrix[:, :3, 3]
    fx = intrinsics[:, 0, 0]
    fy = intrinsics[:, 1, 1]
    cx = intrinsics[:, 0, 2]
    cy = intrinsics[:, 1, 2]
    sk = intrinsics[:, 0, 1]

    uv = torch.stack(torch.meshgrid(torch.arange(resolution, dtype=torch.float32, device=cam2world_matrix.device), torch.arange(resolution, dtype=torch.float32, device=cam2world_matrix.device), indexing='ij')) * (1./resolution) + (0.5/resolution)
    uv = uv.flip(0).reshape(2, -1).transpose(1, 0)
    uv = uv.unsqueeze(0).repeat(cam2world_matrix.shape[0], 1, 1)
    
    x_cam = uv[:, :, 0].view(N, -1)
    y_cam = uv[:, :, 1].view(N, -1)
    z_cam = torch.ones((N, M), device=cam2world_matrix.device)

    x_lift = (x_cam - cx.unsqueeze(-1) + cy.unsqueeze(-1)*sk.unsqueeze(-1)/fy.unsqueeze(-1) - sk.unsqueeze(-1)*y_cam/fy.unsqueeze(-1)) / fx.unsqueeze(-1) * z_cam
    y_lift = (y_cam - cy.unsqueeze(-1)) / fy.unsqueeze(-1) * z_cam

    cam_rel_points = torch.stack((x_lift, y_lift, z_cam, torch.ones_like(z_cam)), dim=-1)

    world_rel_points = torch.bmm(cam2world_matrix, cam_rel_points.permute(0, 2, 1)).permute(0, 2, 1)[:, :, :3]

    ray_dirs = world_rel_points - cam_locs_world[:, None, :]
    ray_dirs = torch.nn.functional.normalize(ray_dirs, dim=2)

    ray_origins = cam_locs_world.unsqueeze(1).repeat(1, ray_dirs.shape[1], 1)

    return ray_origins, ray_dirs

sk是一个用于校正x和y轴不垂直的系数,我们可以认为其是0。我们可以看出整个代码的逻辑是创建一个[0, 1]的uv,实际上就是一个1×1的grid上的坐标。他们现在需要被换算到相机坐标系:

大多数时候,平面z都是1。然后这些相机坐标系下的点左乘刚才的c2w,就可以变换到世界坐标系。继而求出ray_dirs

求出每个像素上光线的方向后,下面就是要考虑沿光线方向采样。在EG3D中,其对于FFHQ数据集,经验性质的选择了ray_start=2.25ray_end=3.3。不同的数据集这个采样的上下限不一样,我推测是试+大概齐估计的。所以对于一个相机,它发出的光线大概就是下面这样:

我们可以看到采样出的那个光锥,以及视场角确实很小。给定一个预训练模型后,这里能改动的地方其实不多。首先更改ray_startray_end大概率会有一些无效的值,感觉能做的操作只剩下修改采样点数,和换成那个“倒数线性插值”。

以及ray_start=2.25ray_end=3.3也和相机半径r在EG3D训练FFHQ被指定为2.75有关,可以看到图中的坐标基本是在-0.5~0.5,因为后续需要在triplane里去插值。

最后画一个多个相机发出光线的样子:

这个动图还是挺有意思的。所有EG3D生成出来的图,都来自于中心那蓝色区域,蓝色区域会由不同的z所导引出的不同的特征图f按照triplane的方法来产生不同的响应,仿佛一个在涌现着什么的水晶球,从某个角度看过去就会有一张又一张的人脸。

“视锥之内,即是NeRF。”

所以在这个settings下,“透视变换”的概念就比较弱,相机内参给定焦距,焦距给定光线方向,除此以外就没有“透视投影”的事情了。

NeRF-PyTorch

上述EG3D的管线中,并没有使用“透视变换”。但实际上在原版NeRF实现中,我们不应忘记,是有NDC空间的使用的。而从正常的欧式空间变换到NDC空间用的就是透视变换。

对于拍摄的facing-forward的数据集,我们应该在训练时变换到NDC空间,这样深度就被缩放到 [-1.1]了,有利于NN去拟合。同时在此基础上对z进行均匀采样,在NDC空间里自动变成了对近处多采一些对远处少采一些的倒数线性采样。

在这篇博客中,详细阐述了NeRF之中的NDC和lindisp。我们在这里不作赘述,重点放在解析其他的一些相机工具代码。

我们主要讨论load_llff.py中的一些函数,这个版本的NeRF代码中传递的poses一般为[N, 3, 5],N为数量,5那儿多出来的那一维是cat了图像的宽高,估计的焦距。中间是3而不是4,是因为存储时舍去了[0.0.0.1]这一行。这里我们使用的poses是由LLFF中的imgs2poses.py导出的,其中储存的就是c2w。

一个常用的函数是求取这些位姿的平均姿态pose_avg:

def viewmatrix(z, up, pos):
    vec2 = normalize(z)
    vec1_avg = up
    vec0 = normalize(np.cross(vec1_avg, vec2))
    vec1 = normalize(np.cross(vec2, vec0))
    m = np.stack([vec0, vec1, vec2, pos], 1)
    return m

def poses_avg(poses):

    hwf = poses[0, :3, -1:]

    center = poses[:, :3, 3].mean(0)
    vec2 = normalize(poses[:, :3, 2].sum(0))
    up = poses[:, :3, 1].sum(0)
    c2w = np.concatenate([viewmatrix(vec2, up, center), hwf], 1)
    
    return c2w

这个是很直接的,我们求取了每个轴方向向量的平均,相机位置的平均,最后得到一个表示相机平均位置的c2w,如下图所示:

一个在处理facing-forward的数据时的一个很有用的操作是recenter_poses,通过对全体姿态左乘一个平均姿态的逆,可以让任意settings下的相机都转变到一个“跟世界坐标系一致”的样子下:

def recenter_poses(poses):

    poses_ = poses+0
    bottom = np.reshape([0,0,0,1.], [1,4])
    c2w = poses_avg(poses)
    c2w = np.concatenate([c2w[:3,:4], bottom], -2)
    bottom = np.tile(np.reshape(bottom, [1,1,4]), [poses.shape[0],1,1])
    poses = np.concatenate([poses[:,:3,:4], bottom], -2)

    poses = np.linalg.inv(c2w) @ poses
    poses_[:,:3,:4] = poses[:,:3,:4]
    poses = poses_
    return poses

这么说可能有点抽象,下面是图示:

我们可以很轻松的看出哪些是recenter后的相机,可以看到,recenter消除了后面做NDC时的歧义,同时方便了后面生成螺旋轨迹的相机轨迹。

另一个比较复杂的就是spherify_poses,它分为两个部分,先对相机轨迹进行“球化”,然后再其基础上构造新的相机轨迹(用于可视化):

def spherify_poses(poses, bds):
    
    p34_to_44 = lambda p : np.concatenate([p, np.tile(np.reshape(np.eye(4)[-1,:], [1,1,4]), [p.shape[0], 1,1])], 1)
    
    rays_d = poses[:,:3,2:3]
    rays_o = poses[:,:3,3:4]

    def min_line_dist(rays_o, rays_d):
        A_i = np.eye(3) - rays_d * np.transpose(rays_d, [0,2,1])
        b_i = -A_i @ rays_o
        pt_mindist = np.squeeze(-np.linalg.inv((np.transpose(A_i, [0,2,1]) @ A_i).mean(0)) @ (b_i).mean(0))
        return pt_mindist

    pt_mindist = min_line_dist(rays_o, rays_d)
    
    center = pt_mindist
    up = (poses[:,:3,3] - center).mean(0)

    vec0 = normalize(up)
    vec1 = normalize(np.cross([.1,.2,.3], vec0))
    vec2 = normalize(np.cross(vec0, vec1))
    pos = center
    c2w = np.stack([vec1, vec2, vec0, pos], 1)

    poses_reset = np.linalg.inv(p34_to_44(c2w[None])) @ p34_to_44(poses[:,:3,:4])

    rad = np.sqrt(np.mean(np.sum(np.square(poses_reset[:,:3,3]), -1)))
    
    sc = 1./rad
    poses_reset[:,:3,3] *= sc
    bds *= sc
    rad *= sc
    
    centroid = np.mean(poses_reset[:,:3,3], 0)
    zh = centroid[2]
    radcircle = np.sqrt(rad**2-zh**2)
    new_poses = []
    
    for th in np.linspace(0.,2.*np.pi, 120):

        camorigin = np.array([radcircle * np.cos(th), radcircle * np.sin(th), zh])
        up = np.array([0,0,-1.])

        vec2 = normalize(camorigin)
        vec0 = normalize(np.cross(vec2, up))
        vec1 = normalize(np.cross(vec2, vec0))
        pos = camorigin
        p = np.stack([vec0, vec1, vec2, pos], 1)

        new_poses.append(p)

    new_poses = np.stack(new_poses, 0)
    
    new_poses = np.concatenate([new_poses, np.broadcast_to(poses[0,:3,-1:], new_poses[:,:3,-1:].shape)], -1)
    poses_reset = np.concatenate([poses_reset[:,:3,:4], np.broadcast_to(poses[0,:3,-1:], poses_reset[:,:3,-1:].shape)], -1)
    
    return poses_reset, new_poses, bds

代码的前半部分是一个有趣的线性代数问题(min_line_dist),给定N条直线,求出一个点到这些直线的距离和最小。这个问题可以直接直接得到解析解。考虑光线原点o,光线方向d,直线上的点可以表述为o+td,考虑直线外一点p,其点到直线距离为(均为列向量):

最后一个等号是因为方向向量d的内积d是1,所以I-d长成了一个幂等矩阵。

考虑有N条光线时,有:

这其实就是一个线性方程组:

于是就有了上面代码里的子函数。接下来,这个函数用不同的叉乘顺序计算了一个类似c2w的矩阵,其位置是刚才估计的直线中心。然后对所有位姿左乘其逆阵,这样就让位姿的lookat对准原点了。同时将相机位置的半径缩放为1,完成“球化”:

(但这个“球化”并不意味着相机姿态真的处理成了一个球面上均匀分布的样子,不要混淆。)

最后一个需要解读的函数是render_path_spiral,这个函数可以生成一个螺旋型的轨迹,来可视化出NeRF在facing-forward时的novel pose:

def render_path_spiral(c2w, up, rads, focal, zdelta, zrate, rots, N):
    render_poses = []
    rads = np.array(list(rads) + [1.])
    hwf = c2w[:,4:5]
    
    for theta in np.linspace(0., 2. * np.pi * rots, N+1)[:-1]:
        c = np.dot(c2w[:3,:4], np.array([np.cos(theta), -np.sin(theta), -np.sin(theta*zrate), 1.]) * rads) 
        z = normalize(c - np.dot(c2w[:3,:4], np.array([0,0,-focal, 1.])))
        render_poses.append(np.concatenate([viewmatrix(z, up, c), hwf], 1))
    return render_poses

可以看出是先创建这样的螺旋线:

与解析几何中熟悉的螺旋线不同,这里的z是用sin(.)来有确保有界的。这个函数的其他输入来自于load_llff_data函数的一个分支:

else:

    c2w = poses_avg(poses)
    print('recentered', c2w.shape)
    print(c2w[:3,:4])

    ## Get spiral
    # Get average pose
    up = normalize(poses[:, :3, 1].sum(0))

    # Find a reasonable "focus depth" for this dataset
    close_depth, inf_depth = bds.min()*.9, bds.max()*5.
    dt = .75
    mean_dz = 1./(((1.-dt)/close_depth + dt/inf_depth))
    focal = mean_dz

    # Get radii for spiral path
    shrink_factor = .8
    zdelta = close_depth * .2
    tt = poses[:,:3,3] # ptstocam(poses[:3,3,:].T, c2w).T
    rads = np.percentile(np.abs(tt), 90, 0)
    c2w_path = c2w
    N_views = 120
    N_rots = 2
    if path_zflat:
        #             zloc = np.percentile(tt, 10, 0)[2]
        zloc = -close_depth * .1
        c2w_path[:3,3] = c2w_path[:3,3] + zloc * c2w_path[:3,2]
        rads[2] = 0.
        N_rots = 1
        N_views/=2

    # Generate poses for spiral path
    render_poses = render_path_spiral(c2w_path, up, rads, focal, zdelta, zrate=.5, rots=N_rots, N=N_views)

通过场景的界限bds,我们可以知道一个类似近远平面的度量,将两者进行倒数线性插值,作为一个估计,来给出focal,focal在def render_path_spiral用于计算lookat的那个固定点。当path_zflat为真时,计算出的rads的z值会归0,于是刚才的螺旋线将退化为圆。

这就是NeRF中的一些关于相机的工具函数了。

3DDFA

在这个部分,我们其实关心的是EG3D中对人脸图片的预处理。这个预处理在处理2D人脸数据时常用的“FFHQ alignment”上更进了一步。但由于EG3D给的预处理代码写的比较分散,不太适合作为“教具”。而基于EG3D的工作PanoHead里提供了一个更自洽和规整的脚本,来让我们可以体会对齐这一步。

PanoHead的预处理用到了3DDFA,这篇工作的全称是Face Alignment in Full Pose Range: A 3D Total Solution,它的本意是为了为了提升大角度偏转情况下人脸对齐的效果,在这个预处理中用到它的目的其实和它的本意有点出入,但这个问题不大,我们也可以把这个当作一个熟悉陌生工作中相机settings的例子。我这里用一张3DDFA GitHub仓库里的teaser来图示一下这个工作具体做了什么事情。

括号里的那一项是一个3DMM重建出的人脸,然后这个人脸上的每一个点,都被R旋转,然后被Pr做正交投影,这个概念非常的简单,正交投影矩阵只是:

你可以看出它直接忽略掉了z,等价于一个焦距非常大的相机,没有近大远小的性质,在这个任务的情景下是非常合适的简化。而它又乘了一个缩放系数f,那么实际上f* Pr在做的操作有一个专门的名字:弱透视投影

在一个比较滑稽的情境下,我被指出分不清强投影和弱投影。实际上这个概念非常简单,只是我过于土鳖,没听过。弱投影应该就是上面的“弱透视投影”,我后来又去搜了一下,发现并没有“强投影”这个专有名词,大概率这个强投影只是指普通的透视投影。

这里的f,你可以把它按照弱透视投影里的那个scaling来理解,或者更简单地,单纯的认为它只是将标准3DMM里那些顶点的量级(大约几百)调整为更符合图像中的量级(可能-1~1)。

所以这里就有了一个很有意思的比较,再看刚才的式子:

这个过程被总结为MVP变换,但其实GAMES101上讲的MVP变换描述了一个更正规情况下的事情。比如在GAMES101中的视图变换,其变换矩阵就是用上面说的lookat的事情构造的,并不是单纯的只有旋转R,透视变换时其作的也是正常的透视变换而非简化后的弱透视变换。

根源是我们这里在单纯的描述“从一张图片中重建出一个人脸mesh”时,我们不需要那么详细的描述。我们可以假定焦距无限远,然后不考虑近大远小,因为人头本身就不是一个很长的东西。

只不过在这里,我们需要“意识到”投影变换的存在。正如前文所说,一个滑稽的事情就是:假如你入门3DV只是跟我一样通过NeRF入门的,那么你大概率“意识不到”投影变换的存在,你会把它弱化,结构为初二物理里平面镜成像规律中就能学到的“近大远小”。你虽然知道有这么个东西,但你大概率很难对这个有什么“intuition”。因为在NeRF里每个像素都是ray casting出来的,想象成是某种“投影”并不直接。因为从体渲染的角度来看,NeRF是backward mapping, 是发出光线去查点,你并不需要像forward mapping般的主动的去投影。

意识到这些以后,我们可以看一下代码了。

在recrop_images.py之前,其实每张输入的人脸图片是通过dlib库检测了一波人脸关键点的。如下图所示:

从FFHQ数据集的对齐开始,就一直沿用一种方法:从这些关键点中,提取左眼,右眼,嘴部的关键点,然后整合出一个quad,其包含4个顶点,在图上可以连成一个类似平行四边形的图案,如上图中的红色圆点(灰色区域表示他们在原图外围)。FFHQ对齐的精髓在于以这4个点对图片进行仿射变换,从而将图片拉正。这样可以保证GAN学习到一个更规整的人脸分布,减少一些掉san的伪影。

但在这里,我们不仅要对齐这个人头,我们还想用3DDFA估计出这个人头此时的位姿,这要求我们对quad进行进一步的处理。同时我们希望位姿符合EG3D的相机系统,所以接下来会有更复杂的一些处理。首先,利用已有的,通过关键点计算出的quad对图片进行仿射变换:

bound = np.array([[0, 0], [0, size-1], [size-1, size-1], [size-1, 0]], dtype=np.float32)
mat = cv2.getAffineTransform(quad[:3], bound[:3])
img = crop_image(img_orig, mat, size, size)

然后用3DDFA扫一下仿射变换后的图片img:

param_lst, roi_box_lst = tddfa(img, boxes)
box_idx = find_center_bbox(roi_box_lst, w, h)

这里find_center_bbox就是单纯的在所有roi_box_lst里,最接近图片中心的那个。但大部分时候图片里只有一张人脸,我们可以认为box_idx就是0。接下来我们会利用第box_idx个param_lst和roi_box_lst来进行处理:

param = param_lst[box_idx]
P = param[:12].reshape(3, -1)  # camera matrix
s_relative, R, t3d = P2sRt(P)

我们可以打印一个P出来,会发现:

(Pdb) print(P)
[[ 4.7336455e-04  3.7309933e-06  1.8318256e-05  5.8912811e+01]
 [ 1.1534430e-06  4.8227943e-04  1.3704226e-05  6.9054771e+01]
 [-1.9893454e-05 -1.7727274e-05  4.7678972e-04 -6.6671005e+01]]
def P2sRt(P):
    """ decompositing camera matrix P.
    Args:
        P: (3, 4). Affine Camera Matrix.
    Returns:
        s: scale factor.
        R: (3, 3). rotation matrix.
        t2d: (2,). 2d translation.
    "
""
    t3d = P[:, 3]
    R1 = P[0:1, :3]
    R2 = P[1:2, :3]
    s = (np.linalg.norm(R1) + np.linalg.norm(R2)) / 2.0
    r1 = R1 / np.linalg.norm(R1)
    r2 = R2 / np.linalg.norm(R2)
    r3 = np.cross(r1, r2)

    R = np.concatenate((r1, r2, r3), 0)
    return s, R, t3d

注意t3d中,实际上有效的只有前两个数,在tddfa中t3d的最后一位永远是-66.67,因为这里没有任何深度之类的事情(焦距无限远)。然后我们提取第一行和第二行,将其归一化,叉乘出第三行,得到新的正交阵作为旋转R。同时计算第一行和第二行模场的平均值,我们就可以估计出1/f,代码中记成了s(scale)。根据3DDFA用的BFM model的量级,其实f可以估计出大致是2000。

(Pdb) print(s_relative)
0.0004781045136041939
(Pdb) print(R)
[[ 0.9992211   0.00787572  0.03866785]
 [ 0.00239068  0.9995937   0.02840398]
 [-0.03842843 -0.02828942  0.9987962 ]]
(Pdb) print(t3d)
[ 58.91281   69.05477  -66.671005]

然后,需要调整一下t3d,在代码中称为"Adjust z-translation in object space",我推测这么做的原因是3DDFA估计出的结果,是隐含在BFM标准人头的那个坐标系下的:

如上图所示,蓝色的人脸是标准的BFM,橙色的是将z轴减去0.5mean(z)以后的BFM。因为我们想要那个人头(是某种程度上)在原点附近的,0.5mean(z)作为一个经验值就被作者采用了。可以看到右侧的橙色人脸,如果我们脑补一个完整的脑壳出来,那么差不多就是一个在原点处的人头了。

这样做虽然是调整"z-translation",但我们这里其实是没有真正的深度的,所以这样做的目的其实是修正t3d中的t2d,即:

R_ = param[:12].reshape(3, -1)[:, :3]
u = tddfa.bfm.u.reshape(3, -1, order='F')
trans_z = np.array([ 0, 0, 0.5*u[2].mean() ]) # Adjust the object center
trans = np.matmul(R_, trans_z.reshape(3,1))
t3d += trans.reshape(3)

在正脸的时候,trans大概是:

[[ 0.6798876 ]
 [ 0.50863648]
 [17.69619383]]

由于是正脸,所以修正项很小。

在一些侧着头的情况下,大概是:

[[-11.60984414]
 [ -0.08104343]
 [ 14.03402536]]

这时候修正项就不能忽略了。

接下来,我们需要对t3d进行归一化,这一步基本算是数字图像处理课程里的习题,t3d本身是在3DDFA的尺寸下的,其中tddfa.size一般是120,我们需要通过之前Face_Boxes得到的roi_box_lst把t3d先缩放回原始图像上:

sx, sy, ex, ey = roi_box_lst[0]
scale_x = (ex - sx) / tddfa.size
scale_y = (ey - sy) / tddfa.size

由于t3d的坐标系和图像坐标系的y轴是反的(如下图):

所以进一步换算为:

t3d[0] = (t3d[0]-1) * scale_x + sx
t3d[1] = (tddfa.size-t3d[1]) * scale_y + sy

然后再根据图片的w和h进行归一化:

t3d[0] = (t3d[0] - 0.5*(w-1)) / (0.5*(w-1)) # Normalize to [-1,1]
t3d[1] = (t3d[1] - 0.5*(h-1)) / (0.5*(h-1)) # Normalize to [-1,1], y is flipped for image space

同时要对缩放系数s_relative也进行“重缩放”:

s_relative = s_relative * 2000
scale_x = (ex - sx) / (w-1)
scale_y = (ey - sy) / (h-1)
s = (scale_x + scale_y) / 2 * s_relative

完成这些以后,我们终于可以对quad进行加工了:

quad_c = quad_c + quad_x * t3d[0]
quad_c = quad_c - quad_y * t3d[1]
quad_x = quad_x * s
quad_y = quad_y * s
c, x, y = quad_c, quad_x, quad_y
quad = np.stack([c - x - y, c - x + y, c + x + y, c + x - y]).astype(np.float32)

我们可以看出,对quad进行修改的目的,是为了让quad圈定的图片中的人头的中心更贴近(TDDFA认为的)图像中心。即争取让每一张图片在对齐后,输入给TDDFA后都能输出近似为[60, 60, -66.7]的t3d(tddfa.size为120)。

最后一步就是要得到c2w矩阵来做训练,由于现在的图像已经被裁剪成没有“translation”了,我们只需要拿来旋转R,就可以产生一个符合PanoHead的c2w了。

这里有一个比较不自然的事情,TDDFA输出的R其实是将BFM从canonical变换成图像中的姿态。在坐标系规定合适的情况下,它是等价于w2c中的R的,但许多时候不一定合适,这样可能x或y方向就会差个负号。我们下面的讨论里就当这里的R是等价于w2c里的R的。

但与之前不同,我们并不是很清楚相机的位置,不能直接用角度定义。我们可以先给出相机外参w2c,然后求逆。当谈及相机外参时,我们往往会提到"translation vector",即 。一个关于t的表述是:世界坐标系下的原点在相机坐标系下的表示。因为当你将[0.0.0.1]左乘w2c后,你会得到[t.1]。

考虑旋转R和平移t构成的w2c,我们可以构造出其逆阵:

所以考虑c2w中的旋转和相机位置C,我们有:

那么直接代入t=-RC就会得到(注意向量r和向量u对于向量c都是正交的):

最后得到的结果即是相机位置的模长,由于相机位置是在球坐标系下采样得到的,所以就是半径。

于是我们直接取t=[0,0,-radius],即可得到w2c,继而求逆得到c2w:

def eg3dcamparams(R_in):
    camera_dist = 2.7
    intrinsics = np.array([[4.2647, 0, 0.5], [0, 4.2647, 0.5], [0, 0, 1]])
    # assume inputs are rotation matrices for world2cam projection
    R = np.array(R_in).astype(np.float32).reshape(4,4)
    # add camera translation
    t = np.eye(4, dtype=np.float32)
    t[2, 3] = - camera_dist

    # convert to OpenCV camera
    convert = np.array([
        [1, 0, 0, 0],
        [0, -1, 0, 0],
        [0, 0, -1, 0],
        [0, 0, 0, 1],
    ]).astype(np.float32)

    # world2cam -> cam2world
    P = convert @ t @ R
    cam2world = np.linalg.inv(P)

    # add intrinsics
    label_new = np.concatenate([cam2world.reshape(16), intrinsics.reshape(9)], -1)
    return label_new

到这里,才完成了全部。我们终于从一张wild-image出发,得到了对齐后的图像和对应的符合EG3D格式的相机位姿。

3D Gaussian Splatting

在3DGS中,有着更加“practical”的相机系统。由于3DGS是一种显式的表达方法,其跟NeRF在相机实现上的最大不同就是因为其是forward mapping,所以必须手动实现投影的过程。除了这一点以外,由于3DGS的官方代码开发周期比较长,所以里面很多令人困惑的地方,下面我们一并过一下。

3DGS中将相机的一些属性组织进一个Camera类中:

class Camera(nn.Module):
    def __init__(self, colmap_id, R, T, FoVx, FoVy, image, gt_alpha_mask,
                 image_name, uid,
                 trans=np.array([0.0, 0.0, 0.0]), scale=1.0, data_device = "cuda"
                 ):
        super(Camera, self).__init__()

        self.uid = uid
        self.colmap_id = colmap_id
        self.R = R
        self.T = T
        self.FoVx = FoVx
        self.FoVy = FoVy
        self.image_name = image_name

        try:
            self.data_device = torch.device(data_device)
        except Exception as e:
            print(e)
            print(f"[Warning] Custom device {data_device} failed, fallback to default cuda device" )
            self.data_device = torch.device("cuda")

        self.original_image = image.clamp(0.0, 1.0).to(self.data_device)
        self.image_width = self.original_image.shape[2]
        self.image_height = self.original_image.shape[1]

        if gt_alpha_mask is not None:
            self.original_image *= gt_alpha_mask.to(self.data_device)
        else:
            self.original_image *= torch.ones((1, self.image_height, self.image_width), device=self.data_device)

        self.zfar = 100.0
        self.znear = 0.01

        self.trans = trans
        self.scale = scale

        self.world_view_transform = torch.tensor(getWorld2View2(R, T, trans, scale)).transpose(0, 1).cuda()
        self.projection_matrix = getProjectionMatrix(znear=self.znear, zfar=self.zfar, fovX=self.FoVx, fovY=self.FoVy).transpose(0,1).cuda()
        self.full_proj_transform = (self.world_view_transform.unsqueeze(0).bmm(self.projection_matrix.unsqueeze(0))).squeeze(0)
        self.camera_center = self.world_view_transform.inverse()[3, :3]

具体的逻辑是,当在Scene中加载数据集中的每一张图片时,旋即读取相应的COLMAP或synthesis数据集的transform,对应实例化一个Camera出来。然后会来到scene/dataset_readers.py里,可以看到其为这两种形式的数据集layout提供了两种回调函数:

sceneLoadTypeCallbacks = {
    "Colmap": readColmapSceneInfo,
    "Blender" : readNerfSyntheticInfo
}

但无论是哪种,我们都会发现,程序赋值的R,都是相机外参的旋转的转置:

def readColmapCameras(cam_extrinsics, cam_intrinsics, images_folder):

 ...
 
        R = np.transpose(qvec2rotmat(extr.qvec))
    ...
    
def readCamerasFromTransforms(path, transformsfile, white_background, extension=".png"):

 ...

            # get the world-to-camera transform and set R, T
            w2c = np.linalg.inv(c2w)
            R = np.transpose(w2c[:3,:3])  # R is stored transposed due to 'glm' in CUDA code
            
            ...

这个其实是历史遗留问题,我们返回来看Camera实现里关于变换的计算:

self.world_view_transform = torch.tensor(getWorld2View2(R, T, trans, scale)).transpose(0, 1).cuda()
self.projection_matrix = getProjectionMatrix(znear=self.znear, zfar=self.zfar, fovX=self.FoVx, fovY=self.FoVy).transpose(0,1).cuda()
self.full_proj_transform = (self.world_view_transform.unsqueeze(0).bmm(self.projection_matrix.unsqueeze(0))).squeeze(0)
self.camera_center = self.world_view_transform.inverse()[3, :3]

查看getWorld2View2:

def getWorld2View2(R, t, translate=np.array([.0, .0, .0]), scale=1.0):
    Rt = np.zeros((4, 4))
    Rt[:3, :3] = R.transpose()
    Rt[:3, 3] = t
    Rt[3, 3] = 1.0

    C2W = np.linalg.inv(Rt)
    cam_center = C2W[:3, 3]
    cam_center = (cam_center + translate) * scale
    C2W[:3, 3] = cam_center
    Rt = np.linalg.inv(C2W)
    return np.float32(Rt)

也就是说输入进Camera的R本身就是w2c转置过的,也就是c2w里的旋转,然后这里又取个转置得到Rt,那么说明Rt应该是w2c,然后这里留了一个修正姿态的废案,给Rt取逆得到c2w。然后再给c2w求逆得到w2c,所以这个函数返回的是相机外参/w2c,这没什么问题。

但在给self.world_view_transform赋值的时候,刚传出来的w2c会再转置一下。所以我们得到的其实是:

同样的,在创建透视投影矩阵时,最后也是将结果转置一下,得到。这里如果我们检视getProjectionMatrix(),会发现其和我们熟悉的透视投影矩阵不太一样,代码中实现的是:

而我们更熟悉的投影矩阵是:

于是投影矩阵P即:

你或许发现第一行的第三个元素和第二行的第三个元素,跟代码实现不符。这应该是代码实现的BUG,但由于函数内定义的视锥是对称的(r=-l,t=-b),所以这一项为0,不会影响结果。

为了对这个过程有更加透彻的理解,这里展示了一个透视矩阵视场角,近平面,远平面变化时,视锥体(黑)和canonical(红)的样子:

注意这里z轴是被缩放到[0.1],所以其看起来是个长方体而非正方体。

我们继续回到代码,之后的self.full_proj_transform是,刚才的self.world_view_transform是,在render时我们传入CUDA侧的变换都是转置过的。这其实是因为CUDA代码中是按照glm的规范,即列主序实现的矩阵乘。例如在cuda_rasterizer/auxiliary.h中:

__forceinline__ __device__ float3 transformPoint4x3(const float3& p, const float* matrix)
{
 float3 transformed = {
  matrix[0] * p.x + matrix[4] * p.y + matrix[8] * p.z + matrix[12],
  matrix[1] * p.x + matrix[5] * p.y + matrix[9] * p.z + matrix[13],
  matrix[2] * p.x + matrix[6] * p.y + matrix[10] * p.z + matrix[14],
 };
 return transformed;
}

__forceinline__ __device__ float4 transformPoint4x4(const float3& p, const float* matrix)
{
 float4 transformed = {
  matrix[0] * p.x + matrix[4] * p.y + matrix[8] * p.z + matrix[12],
  matrix[1] * p.x + matrix[5] * p.y + matrix[9] * p.z + matrix[13],
  matrix[2] * p.x + matrix[6] * p.y + matrix[10] * p.z + matrix[14],
  matrix[3] * p.x + matrix[7] * p.y + matrix[11] * p.z + matrix[15]
 };
 return transformed;
}

可以看到这些索引并不是按照我们以为的:

...
matrix[0] * p.x + matrix[1] * p.y + matrix[2] * p.z + matrix[3],
...

索引排列。

接下来我们简要讨论下为什么要将这两个变换传入CUDA侧。传入CUDA的两个变换矩阵viewmatrix和projmatrix有不同的作用,首先通过这两个变换,可以判断每个线程处理的高斯核是否在关心的视锥内,如果不在可以直接结束该线程:

__forceinline__ __device__ bool in_frustum(int idx,
 const float* orig_points,
 const float* viewmatrix,
 const float* projmatrix,
 bool prefiltered,
 float3& p_view)
{
 float3 p_orig = { orig_points[3 * idx], orig_points[3 * idx + 1], orig_points[3 * idx + 2] };

 // Bring points to screen space
 float4 p_hom = transformPoint4x4(p_orig, projmatrix);
 float p_w = 1.0f / (p_hom.w + 0.0000001f);
 float3 p_proj = { p_hom.x * p_w, p_hom.y * p_w, p_hom.z * p_w };
 p_view = transformPoint4x3(p_orig, viewmatrix);

 if (p_view.z <= 0.2f)// || ((p_proj.x < -1.3 || p_proj.x > 1.3 || p_proj.y < -1.3 || p_proj.y > 1.3)))
 {
  if (prefiltered)
  {
   printf("Point is filtered although prefiltered is set. This shouldn't happen!");
   __trap();
  }
  return false;
 }
 return true;
}

projmatrix计算出的结果p_proj,会联合cov2d,用于计算radii,从而得到该高斯核影响那些像素:

// Compute extent in screen space (by finding eigenvalues of
// 2D covariance matrix). Use extent to compute a bounding rectangle
// of screen-space tiles that this Gaussian overlaps with. Quit if
// rectangle covers 0 tiles. 
float mid = 0.5f * (cov.x + cov.z);
float lambda1 = mid + sqrt(max(0.1f, mid * mid - det));
float lambda2 = mid - sqrt(max(0.1f, mid * mid - det));
float my_radius = ceil(3.f * sqrt(max(lambda1, lambda2)));
float2 point_image = { ndc2Pix(p_proj.x, W), ndc2Pix(p_proj.y, H) };
uint2 rect_min, rect_max;
getRect(point_image, my_radius, rect_min, rect_max, grid);
if ((rect_max.x - rect_min.x) * (rect_max.y - rect_min.y) == 0)
return;

计算cov2d时,需要用viewmatrix得到在相机下高斯点的位置,然后用EWA sampling中仿射变换的一阶近似来得到J,才能得到cov2d:

// The following models the steps outlined by equations 29
// and 31 in "EWA Splatting" (Zwicker et al., 2002). 
// Additionally considers aspect / scaling of viewport.
// Transposes used to account for row-/column-major conventions.
float3 t = transformPoint4x3(mean, viewmatrix);

const float limx = 1.3f * tan_fovx;
const float limy = 1.3f * tan_fovy;
const float txtz = t.x / t.z;
const float tytz = t.y / t.z;
t.x = min(limx, max(-limx, txtz)) * t.z;
t.y = min(limy, max(-limy, tytz)) * t.z;

glm::mat3 J = glm::mat3(
focal_x / t.z, 0.0f, -(focal_x * t.x) / (t.z * t.z),
0.0f, focal_y / t.z, -(focal_y * t.y) / (t.z * t.z),
0, 0, 0);

glm::mat3 W = glm::mat3(
viewmatrix[0], viewmatrix[4], viewmatrix[8],
viewmatrix[1], viewmatrix[5], viewmatrix[9],
viewmatrix[2], viewmatrix[6], viewmatrix[10]);

glm::mat3 T = W * J;

glm::mat3 Vrk = glm::mat3(
cov3D[0], cov3D[1], cov3D[2],
cov3D[1], cov3D[3], cov3D[4],
cov3D[2], cov3D[4], cov3D[5]);

glm::mat3 cov = glm::transpose(T) * glm::transpose(Vrk) * T;

这就是传入这两个变换的逻辑,如果你对这其中的一些推导感兴趣,可以查阅这篇博客。

在讨论中我们忽略了一个传入CUDA侧的参数:相机位置。它的作用是计算点和相机的方向,从而计算球谐系数。在Python侧,相机位置通过计算:

self.camera_center = self.world_view_transform.inverse()[3, :3]

来实现,我们在先前已经讨论过为什么对w2c求逆可以得到相机位置了。

Lie theory of rotations

我们现在已经结合许多实际的项目代码分析并可视化了一些内容,在最后我们需要将目光放在一直以来都被认为习以为常的旋转R上。

了解到这一层对于搬运代码来说完全足够,但这会萌生更多新的问题,比如罗德里格斯公式除了可以几何意义上推导,还可以怎么来?为什么有时候这个公式左侧不叫R而叫exp(.)?为什么实际的旋转是四元数表示的旋转的两倍?以及四元数这一堆是什么?

我们需要用一个更抽象但合适的理论来概括一下,即需要借助李群(Lie Group)和李代数(Lie algebra)。由于我不是数学或物理专业出身,所以下面的叙述以建立直觉和认识为主。对于像我这样没有学过数学的人来说,阅读《视觉SLAM十四讲》的第四章是最合适的切入方式,但仅仅阅读那一段还是太过于管中窥豹了,接下来我们先以《视觉SLAM》十四讲的部分引入,然后再此基础上追加一些内容。

在大一的时候,我们学过如何计算向量叉乘(外积):

右边的结果可以进一步拆成矩阵与向量的乘法:

接下来我们这里按SLAM十四讲里的方法来引入“李代数”,考虑旋转矩阵R,我们假设其代表某个相机的旋转,即随时间变化。那么根据旋转矩阵的正交性:

非正式地,我们对两边关于时间求导,得到:

不失一般性的,我们取=0且R(0)=I,于是右边可以写成:

我们看到,向量(t0)某种程度上反应了R的变化。除此以外,由于:

我们假设在t0附近(t0)=,那么上式就是一个普通的常微分方程,初始值为R(0)=I:

这里t只是我们为了进行上述推理所虚设的一个值,他其实没必要是0,只要在R(0)=I,那么在的邻域内就可以有这样的关系。

我们其实还没有定义这个情景下的指数运算。但我们能感觉到,冥冥之中,好像我如果有一个向量,然后把这个向量写成反对称阵,然后做一下这个“exp(.)”就可以得到R了。除此以外,我们其实也并不了解“为什么上面要这么引入”,为什么要构造一个似是而非的R(t)”和凑个微分方程,感觉很生硬。

但其实,李群、李代数的初衷好像就是为了求解微分方程。不过好像不是上面这样。

下面我们先引入群:我们刚才是习惯上想用R(t)来描述一组R,我们最好引入另一个结构来描述一组R,即群。群的定义是:

则称(G,.)为一个群。

所以我们一直在处理的旋转R,它和矩阵乘法就组成了一个群。显然旋转变换的复合还是旋转变换,且矩阵乘法满足结合律,旋转变换存在幺元即单位旋转,每个旋转变换都存在逆阵。我们将行列式为1的这些旋转矩阵,与矩阵乘法组成的群,称为特殊正交群(Special Orthogonal Group),记作SO(n):

我们关心的往往是二维和三维上的旋转,即SO(2)SO(3)。更重要的是,旋转是可以连续变化的。这种具有连续(光滑)性质的群,称为李群。我们可以简单的认为,由于李群的连续性,使得我们比较熟悉的微积分可以作为分析李群的工具。

这个指数映射是具有一般性的,我们考虑幺元g(0)=1附近的无穷小的群:

那么根据重要极限,我们就有:

所以刚才的例子就是生成元取J=i时的情况。所以从生成元和幺元出发,我们可以得到整个李群。这就是为什么,在刚才微分方程背景下的引入中,我们有“只要在R(t0)=I,那么在的邻域内就可以有这样的关系。”。

有了对李群和李代数的更直观的了解,现在我们引入李代数的定义:

其中二元运算[,]称为李括号(Lie bracket),这个运算只要求自反,不要求结合律,用于表达两个元素的差异。

形式上即泰勒展开。对于指数映射出的矩阵exp(A),这里有一个比较重要的性质:

这是由于任何复方阵A都可以作相似变换,得到一个三角阵:

那么两边作指数映射:

于是我们就得到了罗德里格斯公式,这个公式在SMPL和FLAME都用到了。这事实上说明, 里的so(3)其实就是几何意义上的旋转向量。相应地,我们也可以定义对数映射,来将SO(3)中的元素对应到so(3)上:

这个式子比较难处理,更明智的做法是通过指数映射的结论两边取迹:

现在我们从指数映射上解答了罗德里格斯公式,下面我们开始讨论一下四元数,我们知道四元数是所谓:

我们会发现,这个结构可以用矩阵和矩阵乘法来等价描述,下面我们定义这样的复矩阵:

这样,通过用a,b,c,d线性组合这四个矩阵,我们也可以得到一个复矩阵U:

这也是个特殊的群,我们定义:

称之为二阶特殊幺正(unitary)群,实际上,单位四元数以及其乘法构成了一个群,而2X2的幺正矩阵以及矩阵乘法也构成了一个群,并且每个单位四元数都可以对应到某个幺正矩阵中。这种关系正是群的同构,我们可以通过研究SU(2)来研究四元数。

现在,我们将SO(3)和SU(2)放在一起,考察他们的生成元,在SO(3)的定义中,我们有:

如果我们用矩阵来表达生成元,那么满足上面两个条件的一组基矩阵可以是:

以及,我们发现,他们之间的李括号存在这样的关系:

我们再考虑SU(2)上的生成元:

由于此时是复矩阵,所以:

可以验证,如下三个基矩阵可以用来表达生成元:

所以我们其实也得到了su(2),即:

此时,这几个基矩阵之间的李括号为:

我们将变换作用上去:

用二项式定理展开,并整理,可以得到:

这里:

结合先前的式子,这实际上给出了一种对旋转后的球谐函数进行计算的方法:

上式给出的复球谐,而我们实际用的其实是实球谐,这里的推导只是为了让我们对球谐函数是怎么来的有一个除开调和分析以外的概念,实际上指导意义不大。

写这一部分的时候还是很艰难的,因为我一直以来代数就不好。以及这一部分讨论的这些理论内容,大部分出自数学或物理教材,他们用的符号标记和说法工科出身很难适应,我在图书馆中找了很久,一个很“elementary”的教材是《李群》(邵丹 邵亮 郭紫 著),这本书帮助很大。

至于为什么要坚持写这个部分,这些推导肯定都是不严谨而且以后也用不到的。但通过学习没学过的东西,然后能有选择性的学和跳过哪些内容,从而构建一个逻辑自洽的roadmap,是很重要的一个能力。这个部分只是为了这碟醋包的饺子。我有时就在想,如果当年高等代数课不逃课,事情会不会有所不同。所以这也算是克服心魔了吧。

End


【3D视觉之心】技术交流群
3D视觉之心是面向3D视觉感知方向相关的交流社区,由业内顶尖的3D视觉团队创办!聚焦维重建、Nerf、点云处理、视觉SLAM、激光SLAM、多传感器标定、多传感器融合、深度估计、摄影几何、求职交流等方向。扫码添加小助理微信邀请入群,备注:学校/公司+方向+昵称(快速入群方式)


扫码添加小助理进群

【3D视觉之心】知识星球

3D视觉之心知识星球主打3D感知全技术栈学习,星球内部形成了视觉/激光/多传感器融合SLAM、传感器标定、点云处理与重建、视觉三维重建、NeRF与Gaussian Splatting、结构光、工业视觉、高精地图等近15个全栈学习路线,每天分享干货、代码与论文,星球内嘉宾日常答疑解惑,交流工作与职场问题。



3D视觉之心
3D视觉与SLAM、点云相关内容分享
 最新文章