【Homework】【1--4】Learning resources for DQ Robotics in MATLAB

Learning resources for DQ Robotics in MATLAB

Lesson 1

在这里插入图片描述

代码


% Step 2: Define the real numbers a1 and a2
a1 = 123;
a2 = 321;% Step 3: Calculate and display a3 = a1 + a2
a3 = a1 + a2;
disp(['a3 (a1 + a2) = ', num2str(a3)])% Step 4: Calculate and display a3 = a1 * a2
a3 = a1 * a2;
disp(['a3 (a1 * a2) = ', num2str(a3)])% Step 5: Calculate and display a3 = a1 / a2
a3 = a1 / a2;
disp(['a3 (a1 / a2) = ', num2str(a3)])% Step 6: Define the real numbers theta1 and theta2
theta1 = pi / 12;
theta2 = 3 * pi;% Step 7: Calculate and display a3 = cos(theta1 + theta2)
a3 = cos(theta1 + theta2);
disp(['a3 (cos(theta1 + theta2)) = ', num2str(a3)])% Step 8: Calculate and display a3 = cos^2(theta1 + theta2) + sin^2(theta1 + theta2)
a3 = cos(theta1 + theta2)^2 + sin(theta1 + theta2)^2;
disp(['a3 (cos^2(theta1 + theta2) + sin^2(theta1 + theta2)) = ', num2str(a3)])% Step 9: Calculate and display a3 = exp(a1) + log(a2)
a3 = exp(a1) + log(a2);
disp(['a3 (exp(a1) + log(a2)) = ', num2str(a3)])% Step 10: Display the base of the natural logarithm, e, with 15 significant digits
format long
e_value = exp(1);
disp(['The base of the natural logarithm, e = ', num2str(e_value, 15)])

结果

a3 (a1 + a2) = 444
a3 (a1 * a2) = 39483
a3 (a1 / a2) = 0.38318
a3 (cos(theta1 + theta2)) = -0.96593
a3 (cos^2(theta1 + theta2) + sin^2(theta1 + theta2)) = 1
a3 (exp(a1) + log(a2)) = 2.619517318749063e+53
The base of the natural logarithm, e = 2.71828182845905

Lesson 2

在这里插入图片描述

代码

% quaternion_basics_homework.m
clear all;
include_namespace_dq% 1. Display message
disp('Starting quaternion basics homework');% 2. Store r1 as a rotation of π/8 about the x-axis
phi1 = pi / 8;
r1 = cos(phi1 / 2) + i_*sin(phi1 / 2)% 3. Store r2 as a rotation of π/16 about the y-axis
phi2 = pi / 16;
r2 = cos(phi2 / 2) + j_*sin(phi2 / 2)% 4. Store r3 as a rotation of π/32 about the z-axis
phi3 = pi / 32;
r3 = cos(phi3 / 2) + k_*sin(phi3 / 2)% 5. Calculate the sequential rotation and store in r4
r4 = r1 * r2 * r3% Plot r4
figure;
plot(r4);
title('Sequential Rotation Result (r4)');% 6. Find the reverse rotation of r4 and store in r5
r5 = conj(r4)% Plot r5
figure;
plot(r5);
title('Reverse Rotation (r5)');% 7. Rotate r5 by 360 degrees about the x-axis and store in r6
phi4 = 2 * pi;  % 360 degrees in radians
r6 = cos(phi4 / 2) + i_*sin(phi4 / 2);
disp('Rotation of r5 by 360 degrees about the x-axis (r6):');
r6 = r5 * r6% Plot r5 and r6 to confirm they represent the same rotation
figure;
subplot(1, 2, 1);
plot(r5);
title('Rotation r5');subplot(1, 2, 2);
plot(r6);
title('Rotation r6');

结果

Starting quaternion basics homework
r1 =
0.98079 + 0.19509i
r2 =
0.99518 + 0.098017j
r3 =
0.9988 + 0.049068k
r4 =
0.97395 + 0.19863i + 0.086491j + 0.066992k

r5 =
0.97395 - 0.19863i - 0.086491j - 0.066992k

Rotation of r5 by 360 degrees about the x-axis (r6):
r6 =
- 0.97395 + 0.19863i + 0.086491j + 0.066992k

手写作业

在这里插入图片描述

题目 1

问题：四元数乘法的通用形式是什么？将 $h_1 = a_1 + b_1\hat{i} + c_1\hat{j} + d_1\hat{k}$ 和 $h_2 = a_2 + b_2\hat{i} + c_2\hat{j} + d_2\hat{k}$ 相乘，找到结果 $h_3 = a_3 + b_3\hat{i} + c_3\hat{j} + d_3\hat{k}$ 的形式。

解答：
四元数乘法的通用公式是基于以下规则：

$\hat{i}^2 = \hat{j}^2 = \hat{k}^2 = -1$
$\hat{i}\hat{j} = \hat{k}$ , $\hat{j}\hat{k} = \hat{i}$ , $\hat{k}\hat{i} = \hat{j}$
$\hat{j}\hat{i} = -\hat{k}$ , $\hat{k}\hat{j} = -\hat{i}$ , $\hat{i}\hat{k} = -\hat{j}$

将 $h_1$ 和 $h_2$ 相乘得到 $h_3 = h_1 \cdot h_2$ ，结果可以通过展开每一项并应用以上乘法规则计算。

$h_3 = (a_1a_2 - b_1b_2 - c_1c_2 - d_1d_2) + (a_1b_2 + b_1a_2 + c_1d_2 - d_1c_2)\hat{i} + (a_1c_2 - b_1d_2 + c_1a_2 + d_1b_2)\hat{j} + (a_1d_2 + b_1c_2 - c_1b_2 + d_1a_2)\hat{k}$

题目 2

问题：四元数范数的通用形式是什么？简化 $\sqrt{h_1 h_1^*}$ 。

解答：
四元数 $h_1 = a_1 + b_1\hat{i} + c_1\hat{j} + d_1\hat{k}$ 的范数定义为：

$\| h_1 \| = \sqrt{h_1 h_1^*}$

其中 $h_1^* = a_1 - b_1\hat{i} - c_1\hat{j} - d_1\hat{k}$ 是 $h_1$ 的共轭。

计算 $h_1 h_1^*$ ：

$h_1 h_1^* = (a_1 + b_1\hat{i} + c_1\hat{j} + d_1\hat{k})(a_1 - b_1\hat{i} - c_1\hat{j} - d_1\hat{k}) = a_1^2 + b_1^2 + c_1^2 + d_1^2$

因此，四元数的范数为：
$\| h_1 \| = \sqrt{a_1^2 + b_1^2 + c_1^2 + d_1^2}$

题目 3

问题：证明每个单位四元数 $\cos\left(\frac{\varphi}{2}\right) + \mathbf{v}\sin\left(\frac{\varphi}{2}\right)$ 都具有单位范数。

解答：
单位四元数的定义是具有单位范数的四元数。设 $\cos\left(\frac{\varphi}{2}\right) + \mathbf{v}\sin\left(\frac{\varphi}{2}\right)$ ，其中 $\mathbf{v}$ 是一个单位向量（即 $\| \mathbf{v} \| = 1$ ）。

四元数 $r$ 的范数是其实部和虚部平方和的平方根。在这里，我们需要将 $r$ 写成实部和虚部的形式，然后计算它们的平方和。

对于四元数 $\cos\left(\frac{\varphi}{2}\right) + \mathbf{v} \sin\left(\frac{\varphi}{2}\right)$ ，其实部为 $\cos\left(\frac{\varphi}{2}\right)$ ，虚部为 $\mathbf{v} \sin\left(\frac{\varphi}{2}\right)$ 。

为了计算四元数的范数 $\|r\|$ ，我们取实部和虚部的平方和的平方根：

$\|r\| = \sqrt{\left(\cos\left(\frac{\varphi}{2}\right)\right)^2 + \left(\sin\left(\frac{\varphi}{2}\right) \|\mathbf{v}\|\right)^2}$

其中：

$\cos\left(\frac{\varphi}{2}\right)$ 是实部。
$\sin\left(\frac{\varphi}{2}\right) \|\mathbf{v}\|$ 是虚部的大小。

因为 $\mathbf{v}$ 是单位向量，意味着 $\|\mathbf{v}\| = 1$ ，代入上式可得：

$\|r\| = \sqrt{\cos^2\left(\frac{\varphi}{2}\right) + \sin^2\left(\frac{\varphi}{2}\right)}$

利用三角恒等式 $cos^2 x + \sin^2 x = 1$ ，进一步简化得到：

$\|r\| = \sqrt{1} = 1$

因此，每个单位四元数确实具有单位范数。

Lesson 3

在这里插入图片描述

代码

% dual_quaternion_basics_homework_part1.m
clear all;
close all;
include_namespace_dq% 1. Store x_1 representing a translation of 1 meter along the x-axis,
%    followed by a rotation of π/8 about the z-axis.
t1 = i_; % Translation of 1 meter along the x-axis
r1 = cos(pi/16) + k_*sin(pi/16); % Rotation of π/8 about the z-axis
x_1 = r1 + 0.5 * E_ * t1 * r1% 2. Store x_2 representing a translation of -1 meter along the z-axis,
%    followed by a rotation of -π/2 about the x-axis.
t2 = -k_; % Translation of -1 meter along the z-axis
r2 = cos(-pi/4) + i_*sin(-pi/4); % Rotation of -π/2 about the x-axis
x_2 = r2 + 0.5 * E_ * t2 * r2% 3. Calculate the sequential rotation of the neutral reference frame,
%    first applying x_0 = 1, then x_1, followed by x_2.
x_0 = DQ(1); % Neutral reference frame
disp('Sequential rotation result (x_3):');
x_3 = x_0 * x_1 * x_2 % Sequential rotation% Plot x_3
figure;
plot(x_3);
title('Sequential rotation result x_3');% 4. Obtain the rotation part of x_3 without using rotation()
disp('Rotation part of x_3 (without using rotation()):');
rotation_part_x3 = x_3 - E_ * D(x_3)
rotation_part_x3_2 = P(x_3)% 5. Obtain the translation part of x_3 without using translation()
disp('Translation part of x_3 (without using translation()):');
translation_part_x3 = 2 * D(x_3) * conj(rotation_part_x3)% 6. Obtain the pose transformation of a rotation of π/8 about the z-axis,
%    followed by a translation of 1 meter along the x-axis, and store it in x_4
t41 = 0;
r41 = cos(pi/16) + k_*sin(pi/16); % Rotation of π/8 about the z-axis
x41 = r41 + E_ * 0.5 * t41 * r41;t42 = i_; % Translation of 1 meter along the x-axis
r42 = 1;
x42 = r42 + E_ * 0.5 * t42 * r42;
x_4 = x41 * x42% Check if x_4 is equal to x_1
% If different, explain the reason
if ~(x_4 == x_1)disp('x_4 and x_1 are different because the order of rotation and translation affects the final pose.');
end

结果

x_1 = (0.98079 + 0.19509k) + E*(0.49039i - 0.097545j)
x_2 = (0.70711 - 0.70711i) + E*(0.35355j - 0.35355k)
Sequential rotation result (x_3):
x_3 = (0.69352 - 0.69352i - 0.13795j + 0.13795k) + E*(0.41573 + 0.27779i + 0.27779j - 0.41573k)

在这里插入图片描述

Rotation part of x_3 (without using rotation()):
rotation_part_x3 = 0.69352 - 0.69352i - 0.13795j + 0.13795k
rotation_part_x3_2 = 0.69352 - 0.69352i - 0.13795j + 0.13795k
Translation part of x_3 (without using translation()):
translation_part_x3 = 1i - 1k
x_4 = (0.98079 + 0.19509k) + E*(0.49039i + 0.097545j)
x_4 and x_1 are different because the order of rotation and translation affects the final pose.

$\underline{x}_4$ 与 $\underline{x}_1$ 不同的原因在于旋转和平移的顺序影响了最终的姿态变换结果。

在三维空间中，旋转和平移是非交换的操作，也就是说，先进行旋转然后平移，与先平移然后旋转，通常会得到不同的结果。因为在三维空间中，旋转会改变坐标系的方向，进而影响随后的平移方向。

只有在特定情况下会得到相同的结果，包括：

旋转角度为零：

如果旋转角度为零，则没有实际的旋转效果，平移的方向不会受到影响。因此，不论平移和旋转的顺序，最终的结果都是相同的。

旋转轴与平移方向平行：

如果旋转是绕某个轴（例如 $x$ 轴）进行，而平移也是沿该轴的方向（即 $x$ 轴方向），那么旋转不会改变平移的方向。因此，无论先平移再旋转还是先旋转再平移，结果都会相同。

平移向量为零：

如果平移量为零，则没有实际的平移效果。这种情况下，无论先旋转还是先平移，最终位置都相同。

单位对偶四元数的公式是先平移后旋转
在这里插入图片描述

Lesson 4

在这里插入图片描述

代码

%% dual_quaternion_basics_part2_homework.m
clear all;
close all;
include_namespace_dq;% Define the direction vector for the y-axis
l_y = j_; % 1.Define four different points on the y-axis
p1 = DQ([0, 1, 0]);    % (0, 1, 0)
p2 = DQ([0, -1, 0]);   % (0, -1, 0)
p3 = DQ([0, 2, 0]);    % (0, 2, 0)
p4 = DQ([0, -2, 0]);   % (0, -2, 0)% Calculate the moment vectors
m1 = cross(p1, l_y);
m2 = cross(p2, l_y);
m3 = cross(p3, l_y);
m4 = cross(p4, l_y);% Define the Plücker lines for the y-axis
l1_dq = l_y + E_ * m1
l2_dq = l_y + E_ * m2
l3_dq = l_y + E_ * m3
l4_dq = l_y + E_ * m4% Check if all Plücker lines are the same
disp('Are all Plücker lines the same?');
are_lines_same = isequal(l1_dq, l2_dq) && isequal(l2_dq, l3_dq) && isequal(l3_dq, l4_dq);
disp(are_lines_same);%% 2. Define the point p1 and calculate its distance to l1_dq
p1 = 2 * i_ + j_ + k_; % Point (2, 1, 1)
disp('Distance from p1 to l1_dq:');
d_p1_l1 = norm(cross(p1, l_y) - m1)
sqrt(DQ_Geometry.point_to_line_squared_distance(p1,l1_dq))%% 3. Find the y-z plane using two possible plane normals, n_pi1 and n_pi2
n_pi1 = i_  % Normal in x direction
n_pi2 = -i_  % Opposite normal in -x direction% Define a point on the y-z plane
p_pi = DQ(0); % Origin point on the plane% Calculate distances for each normal
d_pi1 = dot(p_pi, n_pi1);
d_pi2 = dot(p_pi, n_pi2);% Define planes
pi1 = n_pi1 + E_ * d_pi1
pi2 = n_pi2 + E_ * d_pi2%% 4. Calculate the distance between p1 and pi1, and between p1 and pi2
disp('Distance from p1 to pi1:');
d_p1_pi1 = dot(p1, n_pi1) - d_pi1
DQ_Geometry.point_to_plane_distance(p1,pi1)
disp('Distance from p1 to pi2:');
d_p1_pi2 = dot(p1, n_pi2) - d_pi2
DQ_Geometry.point_to_plane_distance(p1,pi2)%% 5. Build a 5mx5mx5m cubic region centered at the origin with 6 planes
% Define the half side length of the cube
half_side = 2.5;% Define the planes using normals pointing towards the origin
% Each plane is represented by its normal vector and a distance from the origin
% Plane x = 2.5 (normal pointing towards origin)
pi3 = -i_ + E_ * (-half_side);     
% Plane x = -2.5
pi4 = i_ + E_ * (-half_side);    
% Plane y = 2.5
pi5 = -j_ + E_ * (-half_side);     
% Plane y = -2.5
pi6 = j_ + E_ * (-half_side);    
% Plane z = 2.5
pi7 = -k_ + E_ * (-half_side);     
% Plane z = -2.5
pi8 = k_ + E_ * (-half_side);    % Plot the cubic region
figure;
hold on;% Plot each plane with the specified size and center location
plot(pi3, 'plane', 10, 'location', DQ([2.5, 0, 0]));   % Plane at x = 2.5
plot(pi4, 'plane', 10, 'location', DQ([-2.5, 0, 0]));  % Plane at x = -2.5
plot(pi5, 'plane', 10, 'location', DQ([0, 2.5, 0]));   % Plane at y = 2.5
plot(pi6, 'plane', 10, 'location', DQ([0, -2.5, 0]));  % Plane at y = -2.5
plot(pi7, 'plane', 10, 'location', DQ([0, 0, 2.5]));   % Plane at z = 2.5
plot(pi8, 'plane', 10, 'location', DQ([0, 0, -2.5]));  % Plane at z = -2.5% Set axis limits to show only the range from -2.5 to +2.5 (size of the cube)
xlim([-2.5 2.5]);
ylim([-2.5 2.5]);
zlim([-2.5 2.5]);% Draw the cube's edge lines
% Define the 8 vertices of a cube with side length 5 (from -2.5 to 2.5)
vertices = [-2.5, -2.5, -2.5;2.5, -2.5, -2.5;2.5, 2.5, -2.5;-2.5, 2.5, -2.5;-2.5, -2.5, 2.5;2.5, -2.5, 2.5;2.5, 2.5, 2.5;-2.5, 2.5, 2.5
];% Define the edges connecting the vertices to form the cube's edges
edges = [1, 2; 2, 3; 3, 4; 4, 1;   % Bottom square5, 6; 6, 7; 7, 8; 8, 5;   % Top square1, 5; 2, 6; 3, 7; 4, 8    % Vertical edges connecting top and bottom
];% Plot the edges
for i = 1:size(edges, 1)plot3([vertices(edges(i, 1), 1), vertices(edges(i, 2), 1)], ...[vertices(edges(i, 1), 2), vertices(edges(i, 2), 2)], ...[vertices(edges(i, 1), 3), vertices(edges(i, 2), 3)], 'k-', 'LineWidth', 2);
end% Add grid and labels
grid on
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');% Set the view to a 3D perspective
view(3);
hold off;%% 6. Find the closest plane to p1 among pi3, pi4, pi5, pi6, pi7, pi8
distances = [DQ_Geometry.point_to_plane_distance(p1,pi3),  % Distance to pi3DQ_Geometry.point_to_plane_distance(p1,pi4), % Distance to pi4DQ_Geometry.point_to_plane_distance(p1,pi5),  % Distance to pi5DQ_Geometry.point_to_plane_distance(p1,pi6), % Distance to pi6DQ_Geometry.point_to_plane_distance(p1,pi7),  % Distance to pi7DQ_Geometry.point_to_plane_distance(p1,pi8) % Distance to pi8
] % 数据类型为doubledistances_dq= [dot(p1, -i_) - (-2.5), ...  Distance to pi3dot(p1, i_) - (-2.5), ...   Distance to pi4dot(p1, -j_) - (-2.5), ...  Distance to pi5dot(p1, j_) - (-2.5), ...   Distance to pi6dot(p1, -k_) - (-2.5), ...  Distance to pi7dot(p1, k_) - (-2.5) ...    Distance to pi8
] % 数据类型为DQ% Find the minimum distance and corresponding plane
[min_distance, closest_plane_index] = min(distances);
disp('The closest plane to p1 is:');
disp(['pi', num2str(closest_plane_index + 2)]);
disp('The distance to the closest plane is:');
disp(min_distance);

结果

l1_dq = 1j
l2_dq = 1j
l3_dq = 1j
l4_dq = 1j
Are all Plücker lines the same?1
Distance from p1 to l1_dq:
d_p1_l1 = 2.2361
ans = 2.2361
n_pi1 = 1i
n_pi2 = - 1i
pi1 = 1i
pi2 = - 1i
Distance from p1 to pi1:
d_p1_pi1 = 2
ans = 2
Distance from p1 to pi2:
d_p1_pi2 = - 2
ans = -2

在这里插入图片描述

distances = 6×1    0.50004.50001.50003.50001.50003.5000distances_dq = 0.5         4.5         1.5         3.5         1.5         3.5
The closest plane to p1 is:
pi3
The distance to the closest plane is:0.5000

第1问

Plücker 坐标的定义是基于方向向量和矩量向量的组合。只要这两个向量的组合满足几何上描述的同一条线，那么即使使用了不同的点，表示的 Plücker 线也是等价的。
在 Plücker 坐标中，不同点只要与同一方向组合成的矩量相同，所描述的几何线就是相同的。

第4问

符号的含义：如果 $d_{p_1, \pi1}$ 和 $d_{p_1, \pi2}$ 的符号不同，说明点 $p_1$ 位于两个平面的不同侧。对于一个平面来说，距离的符号代表点在平面法向量的正侧还是负侧。
不同符号的原因：因为这两个平面法向量方向相反，一个是 $+ i$ ，另一个是 $- i$ ，并且它们在不同的 $x$ 位置。因此，点 $p_1$ 与这两个平面的相对位置不同。结果上，点到这两个平面的距离会具有相反的符号，表示点在一个平面的一侧，同时在另一个平面的另一侧。

两个 y-z 平面的法向量相反，因此同一个点到它们的距离符号可能会不同。这表示该点位于两个平面的相对两侧。

如何绘制立方体的边框。

步骤 1: 定义立方体的顶点

立方体有8个顶点，在三维空间中，每个顶点都由 (x, y, z) 坐标表示。根据题目的要求，立方体的大小从 -2.5 到 2.5，所以每个顶点的坐标将在这个范围内。

vertices = [-2.5, -2.5, -2.5;2.5, -2.5, -2.5;2.5, 2.5, -2.5;-2.5, 2.5, -2.5;-2.5, -2.5, 2.5;2.5, -2.5, 2.5;2.5, 2.5, 2.5;-2.5, 2.5, 2.5
];

这段代码定义了8个顶点的坐标。每一行表示一个顶点的 (x, y, z) 坐标。比如第一个顶点 (-2.5, -2.5, -2.5) 就是立方体的一个角点，坐标值分别为 x = -2.5，y = -2.5，z = -2.5。

步骤 2: 定义边框的连接关系

要绘制立方体的边框，我们需要知道哪些顶点之间是相连的。我们通过定义一个 edges 数组来存储这些连接信息。每一行包含两个顶点的索引，这些顶点通过边相连。

edges = [1, 2;  % 1 <-> 2: 连接顶点 1 和 22, 3;  % 2 <-> 3: 连接顶点 2 和 33, 4;  % 3 <-> 4: 连接顶点 3 和 44, 1;  % 4 <-> 1: 连接顶点 4 和 15, 6;  % 5 <-> 6: 连接顶点 5 和 66, 7;  % 6 <-> 7: 连接顶点 6 和 77, 8;  % 7 <-> 8: 连接顶点 7 和 88, 5;  % 8 <-> 5: 连接顶点 8 和 51, 5;  % 1 <-> 5: 连接顶点 1 和 52, 6;  % 2 <-> 6: 连接顶点 2 和 63, 7;  % 3 <-> 7: 连接顶点 3 和 74, 8   % 4 <-> 8: 连接顶点 4 和 8
];

edges 数组中的每一行表示一条边的连接。共定义了立方体12条边（4条底边，4条顶边，以及4条竖直边）。

步骤 3: 绘制立方体的边框

接下来，通过 plot3 函数根据 edges 数组中的连接关系逐条绘制立方体的边。plot3 函数绘制的是三维空间中的线段，因此它需要三个坐标点数组：X、Y 和 Z。

for i = 1:size(edges, 1)plot3([vertices(edges(i, 1), 1), vertices(edges(i, 2), 1)], ...%x坐标[vertices(edges(i, 1), 2), vertices(edges(i, 2), 2)], ...%y坐标[vertices(edges(i, 1), 3), vertices(edges(i, 2), 3)], ...%z坐标'k-', 'LineWidth', 2);
end

这段代码逐条绘制边：

vertices(edges(i, 1), :) 获取当前边的第一个顶点的坐标。
vertices(edges(i, 2), :) 获取当前边的第二个顶点的坐标。
plot3 函数将这两个顶点连接成一条线段，并使用 'k-' 指定线的颜色为黑色，'LineWidth', 2 设置线的宽度为2。

在 plot3 函数内部，我们通过 edges(i, 1) 和 edges(i, 2) 来获取当前边所连接的两个顶点的索引。

edges(i, 1)：指当前第 i 条边的第一个顶点索引。
edges(i, 2)：指当前第 i 条边的第二个顶点索引。
vertices(edges(i, 1), 1)：获取当前边第一个顶点的 x 坐标。
vertices(edges(i, 2), 1)：获取当前边第二个顶点的 x 坐标。

将这两个 x 坐标值用数组形式 [x1, x2] 表示，传递给 plot3，用于定义边的起点和终点的 x 坐标。

类似地，第二个参数和第三个参数分别表示这两个顶点的 y 和 z 坐标.

整个 plot3 函数的调用，构成了一个从顶点 vertices(edges(i, 1), :) 到顶点 vertices(edges(i, 2), :) 的线段，表示立方体的其中一条边。通过循环遍历所有的 edges，这段代码依次绘制立方体的所有 12 条边，从而绘制出整个立方体的边框。