Fit f(x,y,z)=0 to a set of 3D points using MATLAB

Question

The statement is the following:

Given a 3D set of points (x,y,z), fit a surface defined by a certain number of parameters IMPLICITLY

I'm definitely no expert in programming, but I need to get this done in one way or another. I've considered other programs, such as OriginPro, which can solve this problem pretty easily, but I want to have it done in MATLAB.

The surface is defined by:

A*x^2+B*y^2+C*z^2+D*x+E*y+F*z+G*xy+H*xz+I*yz+J=0

Considering that the Curve Fitting Toolbox can only fit explicit functions, what would you guys suggest?

IMPORTANT REMARK: I'm not asking for a solution, just advice on how to proceed

Answer 1

This can be chalked up to solving a linear system of equations where each point forms a constraint or equation in your system. You would thus find the right set of coefficients in your surface equation that satisfies all points. Using the equation in your question as is, one would find the null space of the linear system that satisfies the surface equation. What you'll have to do is given a set of m points that contain x, y and z coordinates, we can reformulate the above equation as a matrix-vector multiplication with the first argument technically being a matrix of one row and the vector being the coefficients that fit your plane. This is important before we proceed to the null space part of this problem.

In particular, you can agree with me that we can represent the above in the following matrix-vector multiplication:

[x^2 y^2 z^2 x y z xy xz yz 1][A]    =    [0]
                              [B]
                              [C]
                              [D]
                              [E]
                              [F]
                              [G]
                              [H]
                              [I]
                              [J]

Our objective is to find the coefficients A, B, ..., J that would satisfy the constraint above. Now moving onto the more general case, since you have m points, we can build our linear system and thus a matrix of coefficients on the left side of this expression:

[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]

We now build this linear system, and solve to find our coefficients. The trick is to build the matrix that you see on the left hand side of this linear system, which I will call M. Each row is such that you create [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] with x_i, y_i and z_i being the i^th (x,y,z) coordinate in your dataset.

Once you build this, you would thus find the null space of this system. There are many methods to do this in MATLAB. One way is to simply use the 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 if you like, but the null function is a good place to start as it uses singular value decomposition already.

---

I would like to point out that the above will only work if you provide a matrix that is not full rank. To simplify things, this can happen if the number of points you have is less than the number of parameters you have. Therefore, this method would only work if you had up to 9 points in your data. If you have exactly 9, then this method will work very well. If you have less than 9, then there will be more potential solutions as the number of degrees of freedom increases. Specifically, you will have 10 - m possible solutions and any of those solutions is valid. If you have more than 10 points and they are all unique, this would be considered a full rank matrix and thus the only solution to the null space is the trivial one with the coefficients being all set to 0.

In order to escape the possibility of the null space being all 0, or the possibility of the null space providing more than one solution, you probably just want one solution, and you most likely have 10 or more possible points that you want to fit your data with. An alternative method that I can provide is simply an extension of the above but we don't need to find the null space. Specifically, you relax one of the coefficients, say J and you can set that to any value you wish. For example, set it to J = 1. Therefore, the system of equations now changes where J disappears from the mix and it now appears on the right side of the system:

[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]

You can thus find the parameters A, B, ..., I using 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. Additionally, this formulation is nice because if there is an exact solution to the linear system, solving with the pseudoinverse will provide the exact solution. If there is no exact solution to the system, meaning that not all constraints are satisfied, the solution provided is one that minimizes the least squared error between the data and the parameters that were fit with that data.

MATLAB already has an awesome utility to solve a system through linear least squares - in fact, the core functionality of MATLAB is to solve linear algebra problems (if you didn't know that already). You can use matrix left division to solve the problem. Simply put, given that you built your matrix of coefficients above after you introduce the relaxation of J also being called M, the solution to the problem is simply coeff = M\(-ones(m,1)); with m being the number of points and coeff being the coefficients for the surface equation that fit your points. The ones statement in the code creates a column vector of ones that are negative that has m elements.

Using the least squares approach has a more stable and unique solution as you are specifically constraining one of the coefficients, J, to be 1. Using the null space approach only works if you have less points than you do parameters and will possibly give you more than one solution so long as the coefficients span the null space. Specifically, you will get 10 - m solutions and they are all equally good at fitting your data.

I hope this is enough to get you started and good luck!