libigl第三章:矩阵和线性代数
时间:2023-02-26 14:00:01
Libigl 严重依赖密集和稀疏的线性代数 Eigen 库。 除几何处理例程外,libigl 还具有引导 Eigen 线性代数例程让它感觉更像 Matlab 等高级代数库。
切片
Matlab 数组切片是一个非常熟悉和强大的例程。 这允许读取或写入可能不连续的子矩阵。 让我们考虑一下 Matlab 代码:
B = A(R,C);
假如A是一个 m × n m\times n m×n矩阵,R是长度为j的行索引的向量,C是长度为k的列索引的向量,所以B是 j × k j\times k j×k根据R和C从A中提取的的索引。libigl通过相同的操作slice函数给出:
VectorXi R,C; MatrixXd A,B; ... igl::slice(A,R,C,B);
请注意,A 和 B 也可以是稀疏矩阵。
同样,考虑 Matlab 代码:
A(R,C) = B;
这个代码的意思是一个 j × k j\times k j×k通过R和C提供的索引写入A。libigl通过相同的操作slice_into函数给出:
igl::slice_into(B,R,C,A);
排序
Matlab 与其他高级语言相比,提取排序和比较范围的索引变得非常容易。 例如在 Matlab 可以这样写:
[Y,I] = sort(X,1,
'ascend'
)
;
所以,如果X是一个 m × n m\times n m×n的矩阵,那么Y也是一个 m × n m\times n m×n的矩阵,其中Y是由X每一列元素按照降序排列后的结果。第二个输出记录了X是如何排序生成Y的,它们存在着这样的对应关系 Y ( i , j ) = X ( I ( i , j ) , j ) Y(i,j)=X(I(i,j),j) Y(i,j)=X(I(i,j),j).I记录的是一个索引矩阵
同样的功能在libigl中也被支持
igl::sort(X,1,true,Y,I);
类似地,可以在 Matlab 中使用以下命令对整行进行排序:
[Y,I] = sortrows(X,'ascend');
其中,I为长度为m的行索引的向量,其中 Y = X ( I , : ) Y=X(I,:) Y=X(I,:),
在 libigl 中,由以下的代码实现
igl::sortrows(X,true,Y,I);
其中,I同样表示的是排序的索引,可以通过igl::slice(X,I,1,Y),获得相同的Y;
libigl 中提供了类似的函数:max、min 和 unique。
其他的matlab形式的函数
函数名 | 描述 |
---|---|
igl::all | 是否矩阵中的所有元素均为非零元素 |
igl::any | 是否矩阵中的任意元素都是非零元素 |
igl::cat | 连接两个矩阵(在处理稀疏矩阵时很有用) |
igl::ceil | 舍入元素为最接近的整数 |
igl::cumsum | 矩阵元素的累积和 |
igl::colon | 像 Matlab 的 :, 类似于 Eigen 的 LinSpaced |
igl::components | 连接图的部件(和matlab的graphconncomp比较) |
igl::count | 对矩阵行或是列中的非零元素进行计数 |
igl::cross | 每一行的叉乘 |
igl::cumsum | 累计和 |
igl::dot | 每一行的点积 |
igl::eigs | 解决稀疏特征值的问题 |
igl::find | 寻找非零元素的下标 |
igl::floor | 对矩阵的元素进行舍入以取整 |
igl::histc | 对构建直方图时的出现次数进行计数 |
igl::hsv_to_rgb | 将hsv表示的颜色转换为rgb表示(类似于matlab的hsv2rgb) |
igl::intersect | 设置矩阵元素的交集 |
igl::isdiag | 判断矩阵是否为对角矩阵 |
igl::ismember | 判断A中的元素是否出现在B中 |
igl::jet | 沿着彩虹的量化颜色 |
igl::max | 计算每一行或是列的最大的元素 |
igl::median | 计算每一列的中位数 |
igl::min | 计算每一行或是列的最小值 |
igl::mod | 计算每一列的模 |
igl::null | 计算矩阵的零空间基 |
igl::nchoosek | 计算长度为n的向量的所有长度为k的组合 |
igl::orth | 基的正交化 |
igl::parula | 生成从蓝色到黄色的量化颜色图 |
igl::pinv | 计算 Moore-Penrose 伪逆 |
igl::rgb_to_hsv | 将 RGB 颜色转换为 HSV(参见 Matlab 的 rgb2hsv) |
igl::repmat | 沿列和行重复矩阵 |
igl::round | 每元素四舍五入到整数 |
igl::setdiff | 设置矩阵元素的差异 |
igl::setunion | 设置矩阵元素的并集 |
igl::setxor | 设置矩阵元素的异“或” |
igl::slice | 使用索引列表对矩阵的部分进行切片:(对比Matlab的 B = A(I,J)) |
igl::slice_mask | 使用布尔掩码对矩阵的部分进行切片:(对比matlab的B = A(M,N)) |
igl::slice_into | 使用索引列表对矩阵赋值的左侧进行切片(参见 Matlab 的 B(I,J) = A) |
igl::sort | 对矩阵的元素或行进行排序 |
igl::speye | 单位初始化稀疏矩阵 |
igl::sum | 沿列或行求和(稀疏矩阵) |
igl::unique | 提取矩阵的唯一元素或行 |
拉普拉斯方程(这块内容需要后续深入研究)
几何处理中常见的线性系统是拉普拉斯方程:
Δ z = 0 \Delta z=0 Δz=0
受限于一些边界条件,例如 Dirichlet (狄利克雷)边界条件(固定值)
z ∣ ∂ S = z b c z|_{\partial S}=z_{bc} z∣∂S=zbc
在离散设置中,线性系统可以写成:
L z = 0 \mathbf{Lz=0} Lz=0
其中 L \mathbf L L是一个 n × n n\times n n×n的拉普拉斯算子, z \mathbf z z则是由每个顶点值所组成的向量,大多数 z \mathbf z z对应于内部顶点并且是未知的,但有些 z \mathbf z z表示边界顶点处的值。它们的值是已知的,因此我们可以将它们的对应项移到右侧。
从概念上讲,如果我们已经将 z \mathbf z z排序,很容易实现在 z \mathbf z z中内部顶点首先出现,然后是边界顶点。
( L i n , i n L i n , b L b , i n L b , b ) ( z i n z b ) = ( 0 i n z b c ) \begin{pmatrix}\mathbf{L}_{in,in} & \mathbf{L}_{in,b} \\\mathbf{L}_{b,in} & \mathbf{L}_{b,b}\end{pmatrix}\begin{pmatrix}\mathbf{z}_{in} \\\mathbf{z}_{b} \end{pmatrix}=\begin{pmatrix}\mathbf{0}_{in} \\\mathbf{z}_{bc} \end{pmatrix} (Lin,inLb,inLin,bLb,b)(zinzb)=(0inzbc)
方程的底部块不再有意义,所以我们只考虑顶部块:
( L i n , i n L i n , b ) ( z i n z b ) = 0 i n \begin{pmatrix}\mathbf{L}_{in,in}& \mathbf{L}_{in,b} \end{pmatrix}\begin{pmatrix}\mathbf{z}_{in} \\\mathbf{z}_{b} \end{pmatrix}=\mathbf{0}_{in} (Lin,inLin,b)(zinzb)=0in
我们可以将已知值移动到右侧:
L i n , i n z i n = − L i n , b z b \mathbf{L}_{in,in}\mathbf{z}_{in}=\mathbf{-L}_{in,b}\mathbf{z}_b Lin,inzin=−Lin,bzb
最后,我们可以求解内部顶点处未知值的方程。
但是,我们的顶点通常不会以这种方式排序。一种选择是对 V 进行排序,然后按上述方法进行操作,然后对解 Z 进行排序以匹配 V。但是,该解不是很通用。
使用数组切片不需要显式排序。相反,我们可以切出子矩阵块(例如 L i n , i n , L i n , b \mathbf{L}_{in,in},\mathbf{L}_{in,b} Lin,in,Lin,b)并直接遵循上面的线性代数。然后我们可以将解分割成 Z 对应于内部顶点的行
二次能量最小化
相同的拉普拉斯方程可以通过在相同边界条件下最小化狄利克雷能量来等效地导出:
m i n i n i z e z 1 2 ∫ S ∣ ∣ ∇ z ∣ ∣ 2 d A mininize_z\frac{1}{2}\int_S||\nabla z||^2dA mininizez21∫S∣∣∇z∣∣2dA
在我们的离散网格上,回想一下,这变成
m i n i m i z e z 1 2 z T G T G z → m i n i m i z e z z T L z minimize_z\frac{1}{2}\mathbf{z}^T\mathbf{G}^T\mathbf{Gz}\to minimize_z\mathbf{z}^T\mathbf{Lz} minimizez21zTGTGz→minimizezzTLz
在受固定值边界条件约束的网格上最小化一些能量的一般问题是如此广泛,以至于 libigl 有一个专用的 api 来解决此类系统。
让我们考虑一个受到不同公共约束的一般二次最小化问题:
m i n i m i z e z 1 2 z T Q z + z T B + c o n s t a n t minimize_z\frac{1}{2}\mathbf{z}^T\mathbf{Qz}+\mathbf{z}^T\mathbf{B}+constant minimizez21zTQz+zTB+constant
约束于:
z b = z b c a n d A e q z = B e q \mathbf{z}_b=\mathbf{z}_{bc}and\mathbf{A}_{eq}\mathbf{z}=\mathbf{B}_{eq} zb=zbcandAeqz=Beq
其中:
- Q \mathbf{Q} Q通常是一个 n × n n\times n n×n的稀疏矩阵,为二次系数的半正定矩阵(Hessian)
- B \mathbf B B是一个 n × 1 n\times1 n×1的线性系数向量
- z b \mathbf{z}_b zb是 z \mathbf z z的 ∣ b ∣ × 1 |b|\times 1 ∣b∣×1的一部分,对应于边界或是固定的顶点