Fitting

Least Square

Many implementations shown in NURBS book didn’t work, like knot_insert, knot_remove, degree_increase and degree_decrease. Another approach was made to contour this problem, which is by using the Least Square method.

Note

Mathematically speaking, we use the Galerkin’s method, which reduces the integral of the residual square.

Formulation

Let \(\mathbf{C}\) and \(\mathbf{D}\) be two curves defined by

\[\mathbf{C}(u) = \sum_{i=0}^{n-1} f_{i}(u) \cdot \mathbf{P}_i\]
\[\mathbf{D}(u) = \sum_{i=0}^{m-1} g_{i}(u) \cdot \mathbf{Q}_i\]

With respectivelly knot vectors \(\mathbf{U}\) and \(\mathbf{V}\)

The objective is to find the values of \(\mathbf{Q}\) such \(\mathbf{D}\) keeps as near as possible to \(\mathbf{C}\). It’s made by reducing the residual square \(J\)

\[J\left(\mathbf{Q}\right) = \dfrac{1}{2}\int_{a}^{b} \|\mathbf{C}(u)-\mathbf{D}(u)\|^2 \ du\]
\[J\left(\mathbf{Q}\right) = \dfrac{1}{2}\mathbf{Q}^{T} \cdot \left[GG\right] \cdot \mathbf{Q} - \mathbf{Q}^{T} \cdot \left[GF\right] \cdot \mathbf{P} + \dfrac{1}{2}\mathbf{P}^{T} \cdot \left[FF\right] \cdot \mathbf{P}\]

With

\[\left[GG\right]_{ij} = \int_{a}^{b} g_{i}(u) \cdot g_{j}(u) \ du\]
\[\left[GF\right]_{ij} = \int_{a}^{b} g_{i}(u) \cdot f_{j}(u) \ du\]
\[\left[FF\right]_{ij} = \int_{a}^{b} f_{i}(u) \cdot f_{j}(u) \ du\]

Minimizing

Since \(J\) is convex, minimizing \(J\) is obtained by talking its gradient and setting it to zero

\[\nabla J = \mathbf{0} \Rightarrow \left[GG\right] \cdot \mathbf{Q}^{\star} = \left[GF\right] \cdot \mathbf{P}\]

Now it depends only on the computation of the terms \([GG]\) and \([GF]\) and solving the linear system to find \(\mathbf{Q}^{\star}\) from \(\mathbf{P}\). It’s a linear transformation by using matrix \(T\):

\[\mathbf{Q}^{\star} = \left[T\right] \cdot \mathbf{P}\]
\[\left[T\right] = \left[GG\right]^{-1}\left[GF\right]\]

Error

An important mesure is about the error, which is in fact the value of \(J_{min} = J\left(\mathbf{Q}^{\star}\right)\):

\[J\left(\mathbf{Q}^{\star}\right) = \dfrac{1}{2}\left(\left[T\right] \cdot \mathbf{P}\right)^{T} \cdot \left[GG\right] \cdot \left(\left[T\right] \cdot \mathbf{P}\right) -\left(\left[T\right] \cdot \mathbf{P}\right)^{T} \cdot \left[GF\right] \cdot \mathbf{P} + \dfrac{1}{2}\mathbf{P} \cdot \left[FF\right] \cdot \mathbf{P}\]

Expanding and simplifying we get

\[J\left(\mathbf{Q}^{\star}\right) = \dfrac{1}{2} \mathbf{P}^{T} \cdot [E] \cdot \mathbf{P}\]
\[\left[E\right] = \left[FF\right] - \left[GF\right]^{T} \cdot \left[GG\right]^{-1} \cdot \left[GF\right]\]

Essentially the error comes from the transcription from the set of basis functions \(\ \Upsilon_{U}\) to the other set of basis functions \(\Upsilon_{V}\). If \(\Upsilon_{U} \subset \Upsilon_{V}\), that means, the basis functions of \(\mathbf{D}\) describes completelly the basis functions of \(\mathbf{C}\), then the error is zero.

For example, if \(\mathbf{V}\) (knotvector from \(\mathbf{D}\)) is obtained from a knot_insert or degree_increase from \(\mathbf{U}\) (knotvector from \(\mathbf{C}\)), then \(\Upsilon_{U} \subset \Upsilon_{V}\) hence \(J\left(\mathbf{Q}^{\star}\right) = 0\).

But if \(\mathbf{V}\) is obtained from a knot_remove or degree_decrease from \(\mathbf{U}\), then \(\mathbf{D}\) may not be equal to \(\mathbf{C}\), giving a non-zero error.

Least Square with restrictions

Sometimes, it’s needed that the new curve (after transformation) interpolate exactly at some nodes. For example, the extremities points of a curve don’t move when applying knot_remove (left image bellow) or degree_decrease, while the standard least square given above would give a result like in the right figure.

pic1 pic2

Formulation

Let \(\mathbf{C}\) and \(\mathbf{D}\) be two curves defined by

\[\mathbf{C}(u) = \sum_{i=0}^{n-1} f_{i}(u) \cdot \mathbf{P}_i\]
\[\mathbf{D}(u) = \sum_{i=0}^{m-1} g_{i}(u) \cdot \mathbf{Q}_i\]

With respectivelly knot vectors \(\mathbf{U}\) and \(\mathbf{V}\)

The objective is to find the values of \(\mathbf{Q}\) such \(\mathbf{D}\) keeps as near as possible to \(\mathbf{C}\) and satisfies the interpolation restrictions:

\[\mathbf{D}(z_{k}) = \mathbf{C}(z_{k})\]

For \(1 \le k \le m\) nodes \(z_{i} \in \left[a, \ b\right]\)

The same way, we reduce the residual square \(J\), but we add lagrange multipliers \(\lambda\) related to constraint functions \(h(\mathbf{Q})\)

\[J\left(\mathbf{Q}\right) = \dfrac{1}{2}\int_{a}^{b} \|\mathbf{C}(u)-\mathbf{D}(u)\|^2 \ du\]
\[\bar{J}\left(\mathbf{Q}, \lambda\right) = J\left(\mathbf{Q}\right) + \sum_{i=0}^{k-1} \lambda_{i} \cdot h_{i}\left(\mathbf{Q}\right)\]
\[\bar{J}\left(\mathbf{Q}, \lambda\right) = J\left(\mathbf{Q}\right) + \lambda^{T} \cdot \left(\left[G\right]\cdot \mathbf{Q} - \left[F\right] \cdot \mathbf{P}\right)\]

With

\[\left[G\right]_{ij} = g_{j}(z_i)\]
\[\left[F\right]_{ij} = f_{j}(z_i)\]

Minimizing

The same way as before, but with two variables

\[\dfrac{\partial \bar{J}}{\partial \mathbf{Q}} = \mathbf{0} \Rightarrow \left[GG\right] \cdot \mathbf{Q}^{\star} + \left[G\right]^{T} \cdot \lambda = \left[GF\right] \cdot \mathbf{P}\]
\[\dfrac{\partial \bar{J}}{\partial \lambda} = \mathbf{0} \Rightarrow \left[G\right] \cdot \mathbf{Q} = \left[F\right] \cdot \mathbf{P}\]

We mount the expanded matrix with these two equations

\[\begin{split}\begin{bmatrix}\left[GG\right] & \left[G\right]^{T} \\ \left[G\right] & \left[0\right] \end{bmatrix}\begin{bmatrix}\mathbf{Q} \\ \lambda \end{bmatrix} = \begin{bmatrix}\left[GF\right] \\ \left[F\right] \end{bmatrix} \cdot \mathbf{P}\end{split}\]

If it’s solvable (the matrix is not singular), it has the solution for \(\mathbf{Q}\):

\[\left[LL\right] = \left[G\right] \left[GG\right]^{-1}\left[G\right]^{T}\]
\[\left[LG\right] = \left[LL\right]^{-1} \left[G\right] \left[GG\right]^{-1}\]
\[\left[QG\right] = \left[GG\right]^{-1} \left( \left[I\right] - \left[G\right]^{T}\cdot \left[LG\right]\right)\]
\[\left[QF\right] = \left[GG\right]^{-1} \left[G\right]^{T} \left[LL\right]^{-1}\]
\[\left[T\right] = \left[QG\right] \left[GF\right] + \left[QF\right] \left[F\right]\]
\[\mathbf{Q} = \left[T\right] \cdot \mathbf{P}\]

For the error, the expression of the matrix in terms of basis matrix is too complex. We use the computed matrix \(\left[T\right]\) to this expression

Warning

Repeted nodes makes the expanded matrix become singular (det = 0)

Dimension of matrices

Number rows

\(k\)

\(m\)

\(n\)

\(k\)

\(LL\)

\(G, LG\)

\(F\)

\(m\)

\(QF\)

\(GG, QG\)

\(T, GF\)

\(n\)

\(E, FF\)

Discrete Least Square

This type is used when you want to find \(\mathbf{D}(u)\) near \(\mathbf{C}(u)\), but you cannot express \(\mathbf{C}(u)\) as a linear combination of points and basis function.

That means, you want to find \(\mathbf{Q}\) from

  • \(m\) basis functions \(g\)

  • \(k\) nodes \(z_{i} \in \left[a, \ b\right]\)

  • \(k\) points \(\mathbf{Z}_{i} = \mathbf{C}(z_i)\)

\[\mathbf{D}(u) = \sum_{i=0}^{m-1} g_{i}(u) \cdot \mathbf{Q}_i\]

Note

There are three cases for \(k\). The first one gives error

  • \(k < m\) is a under-determined problem which has no unique solution

  • \(k = m\) is a interpolation problem

  • \(k > m\) is a over-determined problem which we use the least squares method

Then, define the residual function \(J\) and get its minimum

\[J\left(\mathbf{Q}\right) = \dfrac{1}{2} \sum_{i=0}^{k-1} \|\mathbf{D}(z_i) - \mathbf{Z}_i\|^2\]

We derivate and set it to zero

TO DO