使用 MATLAB 将 f(x,y,z)=0 拟合到一组 3D 点
Fit f(x,y,z)=0 to a set of 3D points using MATLAB
声明如下：
给定一组 3D 点 (x,y,z)，拟合由一定数量的参数定义的表面 IMPLICITLY
我绝对不是编程方面的专家，但我需要以某种方式完成这项工作。我考虑过其他程序，例如 OriginPro，它可以很容易地解决这个问题，但我想在 MATLAB 中完成它。
曲面定义为：
A*x^2+B*y^2+C*z^2+D*x+E*y+F*z+G*xy+H*xz+I*yz+J=0
考虑到曲线拟合工具箱只能拟合显式函数，你们有什么建议？
重要提示：我不是在寻求解决方案，只是就如何进行提供建议
这可以归结为求解线性方程组，其中每个点在您的系统中形成一个约束或方程。因此，您会在满足所有点的曲面方程中找到一组正确的系数。按原样使用问题中的方程式，可以找到满足曲面方程式的线性系统的 null space 。你需要做的是给定一组包含 x、y 和 z 坐标的 m 点，我们可以将上面的等式重新表述为 matrix-vector 与第一个参数相乘在技术上是一行矩阵，向量是适合您的平面的系数。在我们继续处理此问题的空 space 部分之前，这很重要。
特别是，你可以同意我的意见，我们可以用下面的matrix-vector乘法表示上面的内容：
[x^2 y^2 z^2 x y z xy xz yz 1][A] = [0]
[B]
[C]
[D]
[E]
[F]
[G]
[H]
[I]
[J]
我们的objective就是找到满足上述约束的系数A, B, ..., J。现在转到更一般的情况，因为你有 m 个点，我们可以构建我们的线性系统，从而在这个表达式的左侧构建一个系数矩阵：
[x_1^2 y_1^2 z_1^2 x_1 y_1 z_1 x_1*y_1 x_1*z_1 y_1*z_1 1][A] = [0]
[x_2^2 y_2^2 z_2^2 x_2 y_2 z_2 x_2*y_2 x_2*z_2 y_2*z_2 1][B] = [0]
[x_3^2 y_3^2 z_3^2 x_3 y_3 z_3 x_3*y_3 x_3*z_3 y_3*z_3 1][C] = [0]
... [D] = [0]
... [E] = [0]
... [F] = [0]
... [G] = [0]
... [H] = [0]
... [I] = [0]
[x_m^2 y_m^2 z_m^2 x_m y_m z_m x_m*y_m x_m*z_m y_m*z_m 1][J] = [0]
我们现在构建这个线性系统，并求解以找到我们的系数。诀窍是构建您在该线性系统左侧看到的矩阵，我将其称为 M。每一行都是这样的，你创建 [x_i^2 y_i^2 z_i^2 x_i y_i z_i x_i*y_i x_i*z_i y_i*z_i 1] 时 x_i、y_i 和 z_i 是第 ith (x,y,z) 坐标在你的数据集中。
一旦你构建了这个，你就会发现这个系统的空space。在 MATLAB 中有很多方法可以做到这一点。一种方法是根据需要简单地使用 null function on the matrix you build above and it will return to you a matrix where each column is a potential solution to the surface you are fitting above. That is, each column directly corresponds to the coefficients A, B, ..., J that would fit your data to the surface. You can also try using the singular value decomposition or QR decomposition，但是 null 函数是一个很好的起点，因为它已经使用了奇异值分解。
我想指出，只有当您提供的矩阵 不是 完整 rank 时，以上内容才有效。为了简化事情，如果您拥有的点数少于您拥有的参数数，就会发生这种情况。因此，此方法仅在您的数据中有多达 9 个点时才有效。如果你正好有 9，那么这个方法会很有效。如果少于 9，那么随着自由度数的增加，可能的解决方案会更多。具体来说，您将有 10 - m 个可能的解决方案，并且这些解决方案中的任何一个都是有效的。如果你有超过 10 个点并且它们都是 unique，这将被认为是一个满秩矩阵，因此对 null space 的唯一解决方案是具有系数全部设置为 0.
为了避免null space全为0的可能性，或者null space提供多个解决方案的可能性，你可能只需要一个 解决方案，并且您很可能有 10 个或更多可能的点要与您的数据相匹配。我可以提供的另一种方法只是上述方法的扩展，但我们不需要找到 null space。具体来说，您放宽其中一个系数，比如 J，您可以将其设置为您想要的任何值。例如，将其设置为 J = 1。因此，方程组现在发生变化，其中 J 从组合中消失，现在出现在方程组的右侧：
[x_1^2 y_1^2 z_1^2 x_1 y_1 z_1 x_1*y_1 x_1*z_1 y_1*z_1][A] = [-1]
[x_2^2 y_2^2 z_2^2 x_2 y_2 z_2 x_2*y_2 x_2*z_2 y_2*z_2][B] = [-1]
[x_3^2 y_3^2 z_3^2 x_3 y_3 z_3 x_3*y_3 x_3*z_3 y_3*z_3][C] = [-1]
... [D] = [-1]
... [E] = [-1]
... [F] = [-1]
... [G] = [-1]
... [H] = [-1]
[x_m^2 y_m^2 z_m^2 x_m y_m z_m x_m*y_m x_m*z_m y_m*z_m][I] = [-1]
因此您可以使用 linear least squares where the solution can be solved using the pseudoinverse. The benefit with this approach is that because the matrix is full rank, there is one and only one solution, thus being unique 找到参数 A、B、...、I。此外，这个公式很好，因为如果线性系统有精确解，则使用伪逆求解将提供精确解。如果系统没有精确的解决方案，这意味着并非所有约束都得到满足，则提供的解决方案是最小化数据与适合该数据的参数之间的最小平方误差的解决方案。
MATLAB 已经有了一个很棒的实用程序，可以通过线性最小二乘法求解系统 - 事实上，MATLAB 的核心功能是求解线性代数问题（如果您还不知道的话）。您可以使用矩阵左除法来解决问题。简单地说，假设你在引入 J 的松弛后构建了上面的系数矩阵，也被称为 M，那么问题的解决方案就是 coeff = M\(-ones(m,1)); 和 m是点数，coeff 是适合你的点的表面方程的系数。代码中的 ones 语句创建了一个包含 m 个元素的负列向量。
使用最小二乘法有一个更稳定和独特的解决方案，因为您专门将其中一个系数 J 限制为 1。使用 null space 方法只有在您比你做参数的点少，只要系数跨越空 space，就可能会给你不止一个解决方案。具体来说，您将获得 10 - m 个解决方案，它们都同样擅长拟合您的数据。
我希望这足以让您入门，祝您好运！
声明如下：
给定一组 3D 点 (x,y,z)，拟合由一定数量的参数定义的表面 IMPLICITLY
我绝对不是编程方面的专家，但我需要以某种方式完成这项工作。我考虑过其他程序，例如 OriginPro，它可以很容易地解决这个问题，但我想在 MATLAB 中完成它。
曲面定义为：
A*x^2+B*y^2+C*z^2+D*x+E*y+F*z+G*xy+H*xz+I*yz+J=0
考虑到曲线拟合工具箱只能拟合显式函数，你们有什么建议？
重要提示：我不是在寻求解决方案，只是就如何进行提供建议
这可以归结为求解线性方程组，其中每个点在您的系统中形成一个约束或方程。因此，您会在满足所有点的曲面方程中找到一组正确的系数。按原样使用问题中的方程式，可以找到满足曲面方程式的线性系统的 null space 。你需要做的是给定一组包含 x、y 和 z 坐标的 m 点，我们可以将上面的等式重新表述为 matrix-vector 与第一个参数相乘在技术上是一行矩阵，向量是适合您的平面的系数。在我们继续处理此问题的空 space 部分之前，这很重要。
特别是，你可以同意我的意见，我们可以用下面的matrix-vector乘法表示上面的内容：
[x^2 y^2 z^2 x y z xy xz yz 1][A] = [0]
[B]
[C]
[D]
[E]
[F]
[G]
[H]
[I]
[J]
我们的objective就是找到满足上述约束的系数A, B, ..., J。现在转到更一般的情况，因为你有 m 个点，我们可以构建我们的线性系统，从而在这个表达式的左侧构建一个系数矩阵：
[x_1^2 y_1^2 z_1^2 x_1 y_1 z_1 x_1*y_1 x_1*z_1 y_1*z_1 1][A] = [0]
[x_2^2 y_2^2 z_2^2 x_2 y_2 z_2 x_2*y_2 x_2*z_2 y_2*z_2 1][B] = [0]
[x_3^2 y_3^2 z_3^2 x_3 y_3 z_3 x_3*y_3 x_3*z_3 y_3*z_3 1][C] = [0]
... [D] = [0]
... [E] = [0]
... [F] = [0]
... [G] = [0]
... [H] = [0]
... [I] = [0]
[x_m^2 y_m^2 z_m^2 x_m y_m z_m x_m*y_m x_m*z_m y_m*z_m 1][J] = [0]
我们现在构建这个线性系统，并求解以找到我们的系数。诀窍是构建您在该线性系统左侧看到的矩阵，我将其称为 M。每一行都是这样的，你创建 [x_i^2 y_i^2 z_i^2 x_i y_i z_i x_i*y_i x_i*z_i y_i*z_i 1] 时 x_i、y_i 和 z_i 是第 ith (x,y,z) 坐标在你的数据集中。
一旦你构建了这个，你就会发现这个系统的空space。在 MATLAB 中有很多方法可以做到这一点。一种方法是根据需要简单地使用 null function on the matrix you build above and it will return to you a matrix where each column is a potential solution to the surface you are fitting above. That is, each column directly corresponds to the coefficients A, B, ..., J that would fit your data to the surface. You can also try using the singular value decomposition or QR decomposition，但是 null 函数是一个很好的起点，因为它已经使用了奇异值分解。
我想指出，只有当您提供的矩阵 不是 完整 rank 时，以上内容才有效。为了简化事情，如果您拥有的点数少于您拥有的参数数，就会发生这种情况。因此，此方法仅在您的数据中有多达 9 个点时才有效。如果你正好有 9，那么这个方法会很有效。如果少于 9，那么随着自由度数的增加，可能的解决方案会更多。具体来说，您将有 10 - m 个可能的解决方案，并且这些解决方案中的任何一个都是有效的。如果你有超过 10 个点并且它们都是 unique，这将被认为是一个满秩矩阵，因此对 null space 的唯一解决方案是具有系数全部设置为 0.
为了避免null space全为0的可能性，或者null space提供多个解决方案的可能性，你可能只需要一个 解决方案，并且您很可能有 10 个或更多可能的点要与您的数据相匹配。我可以提供的另一种方法只是上述方法的扩展，但我们不需要找到 null space。具体来说，您放宽其中一个系数，比如 J，您可以将其设置为您想要的任何值。例如，将其设置为 J = 1。因此，方程组现在发生变化，其中 J 从组合中消失，现在出现在方程组的右侧：
[x_1^2 y_1^2 z_1^2 x_1 y_1 z_1 x_1*y_1 x_1*z_1 y_1*z_1][A] = [-1]
[x_2^2 y_2^2 z_2^2 x_2 y_2 z_2 x_2*y_2 x_2*z_2 y_2*z_2][B] = [-1]
[x_3^2 y_3^2 z_3^2 x_3 y_3 z_3 x_3*y_3 x_3*z_3 y_3*z_3][C] = [-1]
... [D] = [-1]
... [E] = [-1]
... [F] = [-1]
... [G] = [-1]
... [H] = [-1]
[x_m^2 y_m^2 z_m^2 x_m y_m z_m x_m*y_m x_m*z_m y_m*z_m][I] = [-1]
因此您可以使用 linear least squares where the solution can be solved using the pseudoinverse. The benefit with this approach is that because the matrix is full rank, there is one and only one solution, thus being unique 找到参数 A、B、...、I。此外，这个公式很好，因为如果线性系统有精确解，则使用伪逆求解将提供精确解。如果系统没有精确的解决方案，这意味着并非所有约束都得到满足，则提供的解决方案是最小化数据与适合该数据的参数之间的最小平方误差的解决方案。
MATLAB 已经有了一个很棒的实用程序，可以通过线性最小二乘法求解系统 - 事实上，MATLAB 的核心功能是求解线性代数问题（如果您还不知道的话）。您可以使用矩阵左除法来解决问题。简单地说，假设你在引入 J 的松弛后构建了上面的系数矩阵，也被称为 M，那么问题的解决方案就是 coeff = M\(-ones(m,1)); 和 m是点数，coeff 是适合你的点的表面方程的系数。代码中的 ones 语句创建了一个包含 m 个元素的负列向量。
使用最小二乘法有一个更稳定和独特的解决方案，因为您专门将其中一个系数 J 限制为 1。使用 null space 方法只有在您比你做参数的点少，只要系数跨越空 space，就可能会给你不止一个解决方案。具体来说，您将获得 10 - m 个解决方案，它们都同样擅长拟合您的数据。
我希望这足以让您入门，祝您好运！