Aproximación de curvas con BSpline

Comienzo este Blog con uno de los temas fundamentales de la Computación Gráfica, las curvas y superficies están en el corazón de la disciplina y entre las representaciones más poderosas para curvas y superficies tenemos a la BSpline.

Como programador de gráficos creo que no hay mejor forma de conocer este hermoso tema que poner las manos en la masa. En este post voy a mostrar cómo aproximar e interpolar datos n-dimensionales con BSpline usando el método de los mínimos cuadrados. Mostraré tanto las ecuaciones como el código fuente y trataré de hacer algunos paralelos entre ellos, las optimizaciones de código las dejaré como ejercicio para los lectores :).

BSpline Fitting
Aproximación de BSpline

¿Para qué aproximar datos con BSpline? en la práctica la aproximación BSpline se usa para:

  • Poder suavizar datos ruidosos (e.g. series de tiempo)
  • Generar curvas C^n continuas (cuyas n diferenciales existan) lo cual es importante para optimización, si se quiere encontrar los mínimos y máximos por ejemplo
  • Comprimir datos reduciéndolos a unos pocos puntos de control
  • Poder modificar los datos manteniendo continuidad C^n a través de puntos de control
  • Animar objetos a través de trayectorias suaves y continuas

Es necesario distinguir la aproximación de la interpolación. Una curva aproximada no necesariamente pasa por los todos los puntos (pero trata de acercarse), mientras que una curva interpolada si pasa por todos los puntos. La interpolación BSpline se utiliza para computar todos los “datos intermedios” conociendo sólo unas pocas muestras, por ejemplo permite computar colores intermedios en una fotografía, animaciones intermedias en una secuencia de key frames, deformaciones intermedias en una secuencia de deformaciones, transformaciones intermedias, etc.

Pero antes de entrar de lleno a la aproximación repasemos los conceptos básicos de BSplines y su formulación.

La BSpline

La curva BSpline es una función paramétrica Q(t) \in \mathbb{R}^{n}. El parámetro t es un escalar que está en el intervalo [0; 1] (el intervalo puede ser cualquiera, pero elegimos este por convención).

Grado y Continuidad

El grado de la BSpline es \emph{D}, lo que quiere decir que la máxima potencia de t es D. Por ejemplo, si D = 2 la BSpline es cuadrática, si D = 3 la BSpline es cúbica y así sucesivamente.

La continuidad paramétrica de la BSpline está dada por el número de veces que podemos diferenciar la curva con respecto al parámetro t. Es el mismo concepto de continuidad que usamos en análisis matemático. Para la BSpline, la continuidad depende del grado de la curva, por ejemplo la BSpline cúbica tiene continuidad C^{2} lo que quiere decir que posee continuidad posicional, tangencial y de curvatura.

Formulación

La BSpline está definida por N+1 puntos de control P_{i} \in \mathbb{R}^{n}, los cuales corresponden a los coeficientes de la función Q(t).

La BSpline es una combinación lineal de los puntos de control P_{i} con las funciones base B_{i}(t) \in \mathbb{R} las cuales son polinomios formados a partir del parámetro t:

Q(t) = \sum_{i=0}^{N}{P_{i} B_{i}(t)} \quad ,\forall t \in [0 ; 1]

Si codificamos la función Q(t) directamente en C#, quedaría como sigue (para el caso Q(t) \in \mathbb{R}^{2}):

            List curve = new List();
            for (double t = 0; t < 1.01; t += 0.01)
            {
                double[] Bt = BSplineBasis(t, knots, D, N);

                Point Qt = new Point();
                Qt.X = 0; Qt.Y = 0;
                for (int i = 0; i <= N; ++i)
                {
                    Qt.X += P[i].X * Bt[i];
                    Qt.Y += P[i].Y * Bt[i];
                }
                curve.Add(Qt);
            }

B-Spline
Ejemplo de BSpline generada con el código de arriba

Además de los puntos de control P_{i}, es necesario especificar una secuencia no decreciente de escalares llamados “knots”: \{ t_0, t_1, ..., t_N,..., t_{N+D}, t_{N+D+2} \}. El vector de knots influye en la forma de la curva. Existen muchas formas de especificar el vector de knots, básicamente hay dos grupos: knots uniformes y knots no uniformes.

Por ultimo, las funciones base B_{i}(t) \in \mathbb{R} se definen recursivamente según el algoritmo de Cox-DeBoor:

N_{i,j}(t) = \frac{t - t_{i}}{t_{i+j} - t_{i}} N_{i,j-1}(t) + \frac{t_{i+j+1} - t}{t_{i+j+1} - t_{i+1}} N_{i+1,j-1}(t)

para 0 \leq i \leq N + D - j y 1 \leq j \leq D. La recursión termina con la definición de N_{i,0}:

N_{i,0}(t) = \left \{ \begin{array}{ll} 1, & t_{i} \leq t < t_{i+1} \\ 0, & \mbox{otros casos} \end{array} \right.

para 0 \leq i \leq N + D. Cuya implementación en código C# sería la siguiente:

       double[] BSplineBasis(double t, double[] knots, int D, int N)
        {
            double[] basis = new double[N+1];

            for (int i = 0; i <= N; ++i)
            {
                basis[i] = Nij(i, D, knots, t);
            }

            return basis;
        }

        double Nij(int i, int j, double[] knots, double t)
        {
            if (j == 0)
            {
                if (t < knots[i + 1] && t >= knots[i])
                    return 1.0;
                else
                    return 0.0;
            }

            double coefA, coefB, denom;

            denom = knots[i + j] - knots[i];
            if (denom != 0.0)
                coefA = ((t - knots[i]) / denom) * Nij(i, j - 1, knots, t);
            else
                coefA = 0;

            denom = knots[i + j + 1] - knots[i + 1];
            if (denom != 0.0)
                coefB = ((knots[i + j + 1] - t) / denom) * Nij(i + 1, j - 1, knots, t);
            else
                coefB = 0;

            return coefA + coefB;
        }

Esta implementación, aunque fiel a la formula, puede sufrir de algunos errores al finalizar la recursión debido a la precisión de punto flotante, pero es simple subsanarlo de la siguiente forma:

            if (j == 0)
            {
                if (Math.Abs(t - 1.0) < 1e-6)
                {
                    if ((t < knots[i + 1] || Math.Abs(t - knots[i + 1]) < 1e-6) && t >= knots[i])
                        return 1.0;
                    else
                        return 0.0;
                }

                if (t < knots[i + 1] && t >= knots[i])
                    return 1.0;
                else
                    return 0.0;
            }

La BSpline uniforme (con knots uniformemente espaciados) nunca pasa por los puntos de control, mientras que la BSpline no uniforme puede interpolar algunos puntos. Una elección particular de knots no uniformes que interpola los extremos es la siguiente (usado para generar el ejemplo de arriba):

knots = \{ \underbrace{0, 0, ..., 0}_{D+1}, \frac{1}{N+1-D}, \frac{2}{N+1-D}, ..., \frac{N-D}{N+1-D}, \underbrace{1, 1, ..., 1}_{D+1} \}

Este es el vector que usaremos nosotros. En código C# lo podemos expresar de la siguiente forma:

            double[] knots = new double[N + D + 2];

            int k = 0;
            for (k = 0; k <= D; ++k)
                knots[k] = 0.0;

            for (; k <= N; ++k)
                knots[k] = (k - D) / (N + 1.0 - D);

            for (; k <= N+D+1; ++k)
                knots[k] = 1.0;

Aproximación de datos

Dado un conjunto ordenado de datos de entrada \{ Q_0, Q_1, ..., Q_N \} \quad, Q_i \in \mathbb{R}^{n}, tenemos que determinar un conjunto de puntos de control \{ P_0, P_1, ..., P_N \} \quad ,P_i \in \mathbb{R}^{n} tales que:

Q_k(t_k) = \sum_{i=0}^{N}{P_i B_i(t_k)}

Nota que usamos t_k en lugar de t y por consiguiente Q_k en lugar de Q. Ya que Q_k son los datos de entrada y tenemos que asumir que fueron muestreados en secuencia creciente \{ \{t_0, Q_0\}, \{t_1, Q_1\}, ... \{t_M, Q_M\} \} . Esto nos va a ayudar a plantear el sistema lineal de forma matricial de la siguiente manera:

Sea el vector columna \mathrm{\bar Q} = \{ Q_0, Q_1, ..., Q_M \}, el vector columna \mathrm{\bar P} = \{ P_0, P_1, ..., P_N \} y la matriz \mathrm{A}_{M+1 \times N+1}:

\mathrm{A}_{M+1 \times N+1} = \left[ \begin{array}{cccc} B_0(t_0) & B_1(t_0) & \cdots & B_N(t_0) \\ B_0(t_1) & B_1(t_1) & \cdots & B_N(t_1) \\ \vdots & \vdots & \ddots & \vdots \\ B_0(t_M) & B_1(t_M) & \cdots & B_N(t_M) \end{array} \right]

El sistema lineal entonces se puede reescribir de la siguiente forma:

\mathrm{\bar Q} = \mathrm{A} \ \mathrm{\bar P}

Podríamos resolver este sistema directamente encontrando la inversa de la matriz \mathrm{A}:

\mathrm{\bar P} = \mathrm{A}^{-1} \ \mathrm{\bar Q}

Para computar \mathrm{A}^{-1} es necesario que \mathrm{A} sea cuadrada, es decir, que el número de puntos de control P_i sea igual al número de datos Q_i (M = N). Esto puede funcionar siempre y cuando el número M sea bajo (una docena de puntos como mucho), ya que normalmente el numero de puntos de control debe ser menor que el número de datos M > N, de otro modo la curva resultante (de existir) se volvería totalmente salvaje.

La mejor solución, en el sentido de los mínimos cuadrados, sería encontrar el vector \mathrm{\bar P} que minimiza la función:

arg\min_{\mathrm{\bar P}} | \mathrm{A} \ \mathrm{\bar P} - \mathrm{\bar Q} |^2

Esta función es cuadrática en los componentes de \mathrm{\bar P} y por lo tanto su gráfica es un paraboloide cuyo mínimo global se encuentra en el vértice donde sus derivadas parciales con respecto de \mathrm{\bar P} se vuelven cero. Esa minimización se puede resolver eficientemente utilizando las ecuaciones normales:

\mathrm{A}^T \ \mathrm{\bar Q} = \mathrm{A}^T \ \mathrm{A} \ \mathrm{\bar P}
\mathrm{\bar P} = (\mathrm{A}^T \ \mathrm{A})^{-1} \ \mathrm{A}^T \ \mathrm{\bar Q}

Lamentablemente este sistema puede quedar mal condicionado (ill conditioned) debido a que algunos autovalores de la matriz (\mathrm{A}^T \ \mathrm{A}) son cercanos a cero. Sin embargo, un caso especial en que la matriz no queda mal condicionada es cuando el vector de knots es uniforme (excepto en los extremos), es decir:

knots = \{ \underbrace{0, 0, ..., 0}_{D+1}, \frac{1}{N+1-D}, \frac{2}{N+1-D}, ..., \frac{N-D}{N+1-D}, \underbrace{1, 1, ..., 1}_{D+1} \}

Si se elige un vector de knots inadecuado, la matriz puede quedar mal condicionada. No obstante, el sistema aún puede ser resuelto usando Singular Value Decomposition (es comúnmente usado para invertir matrices mal condicionadas), en nuestro caso resulta más adecuado resolverlo utilizando la factorización de Cholesky ya que la matriz (\mathrm{A}^T \ \mathrm{A}) es simétrica, esparcida y en banda (tiene entradas simétricas formando bandas a lo largo de la diagonal, el resto son ceros).

Para este post yo he utilizado el método de factorización LU, en particular la implementación de John Burkardt. La cual empaqueté dentro en la clase MatrixNxM.

El siguiente fragmento en C# muestra el proceso de solución:

        List LeastSquaresFit(List samplePoints)
        {
            if (samplePoints.Count < 4)
                throw new ArgumentException();

            int sampleQuantity = samplePoints.Count;

            //En una aproximación, el numero de puntos de control
            //tiene que ser menor que el numero de muestras
            int controlQuantity = (int)(sampleQuantity * 0.9); //entre [0.1 ; 0.9]

            if( controlQuantity < 4 )
                controlQuantity = 4;

            int N = controlQuantity - 1; //0 - N
            int M = sampleQuantity - 1; //0 - M
            int D = 3; //esta BSpline será cúbica.
            int NplusDplusOne = N + D + 1;

            //
            // Creo el vector de knots.
            //
            double[] knots = new double[NplusDplusOne + 1];

            int k = 0;
            for (k = 0; k <= D; ++k)
                knots[k] = 0.0;

            for (; k <= N; ++k)
                knots[k] = (k - D) / (N + 1.0 - D);

            for (; k <= NplusDplusOne; ++k)
                knots[k] = 1.0;

            //
            // Tenemos que resolver el sistema lineal
            // por minimos cuadrados:
            //
            // (A^T) Q = (A^T) A P
            //
            // para P
            //

            //Creo la matriz A
            MatrixNxM AMat = new MatrixNxM(M + 1, N + 1);
            for (int i = 0; i <= M; i++)
            {
                for (int j = 0; j <= N; j++)                 {                     double tk = i / (double)M;                     AMat[i, j] = Nij(j, D, knots, tk);                 }             }             //La transpuesta A^T             MatrixNxM ATMat = AMat.Transpose();             //La matriz (A^T) A ==> esta matriz es simétrica, esparcida y en banda!
            MatrixNxM ATtimesA = ATMat * AMat;

            //Hallamos la matriz inversa resolviendo el sistema (A^T)A = A^T
            MatrixNxM AInv = ATtimesA.SolveLU(ATMat);

            //P = (A^-1) Q
            MatrixNxM QX = new MatrixNxM(M + 1, 1);
            MatrixNxM QY = new MatrixNxM(M + 1, 1);
            for (int i = 0; i <= M; i++)
            {
                QX[i, 0] = samplePoints[i].X;
                QY[i, 0] = samplePoints[i].Y;
            }

            MatrixNxM PX = AInv * QX;
            MatrixNxM PY = AInv * QY;

            List controlPoints = new List();

            for (int i = 0; i <= N; i++)
                controlPoints.Add(new Point() { X = PX[i, 0], Y = PY[i, 0] });

            return controlPoints;
        }

Código Fuente

Espero que hayan disfrutado el post. Les dejo el código fuente para que se diviertan.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s