matlab - Ajuste de curva con puntos y en posiciones x repetidas(brazos Galaxy Spiral)

2d curve-fitting

Question

Actualmente tengo un programa MATLAB que toma imágenes RGB de los brazos en espiral trazados de las galaxias y selecciona el componente más grande del brazo y traza solo eso.

He intentado usar la herramienta de ajuste de curvas incorporada de matlab con una estría suavizante para ajustarla y obtengo el siguiente resultado:

He intentado utilizar interp1 con adaptación paramétrica para obtener solo malos resultados.

¿Hay alguna manera de adaptarse a este tipo de curva?

Answer 1

Su error se debe a que maneja su curva 2D como función, que no es el caso (tiene más valores y para la misma x ) y es por eso que el ajuste falla en el lado derecho (cuando toca el área sin función).

Para remediarlo, debe separar el ajuste de la curva en cada dimensión. De modo que puede ajustar cada eje como una función separada. Para eso necesita usar un parámetro de función diferente (no x). Si ordena sus puntos de alguna manera (por ejemplo, por la distancia de la curva desde el punto de inicio, o por el ángulo polar o por lo que sea), entonces puede usar el índice de punto como dicho parámetro de función.

Entonces has hecho algo como esto:

y(x) = fit((x0,y0),(x1,y1),(x2,y2)...)

que devuelve un polinomio para y(x) . En cambio, deberías hacer algo como esto:

x(t) = fit(( 0,x0),( 1,x1),( 2,x2)...)
y(t) = fit(( 0,y0),( 1,y1),( 2,y2)...)

donde t es su nuevo parámetro ajustado al orden del punto en la lista ordenada. La mayoría de las curvas usan parámetros en el rango t=<0.0,1.0> para facilitar el cálculo y el uso. Entonces, si obtuviste N puntos, puedes convertir el índice de puntos i=<0,N-1> al parámetro de curva t como este:

t=i/(N-1);

Cuando ploteas necesitas cambiar tu

plot(x,y(x))

A:

plot(x(t),y(t))

Hice un ejemplo simple de cúbica simple de interpolación simple en C ++ / VCL para su tarea, así que es mejor que vea lo que quiero decir:

    picture pic0,pic1;
        // pic0 - source
        // pic1 - output
    int x,y,i,j,e,n;
    double x0,x1,x2,x3,xx;
    double y0,y1,y2,y3,yy;
    double d1,d2,t,tt,ttt;
    double ax[4],ay[4];
    approx a0,a3; double ee,m,dm; int di;
    List<_point> pnt;
    _point p;

    // [extract points from image]
    pic0.load("spiral_in.png");
    pic1=pic0;
    // scan image
    xx=0.0; x0=pic1.xs;
    yy=0.0; y0=pic1.ys;
    for (y=0;y<pic1.ys;y++)
     for (x=0;x<pic1.xs;x++)
      // select red pixels
      if (DWORD(pic1.p[y][x].dd&0x00008080)==0)     // low blue,green
       if (DWORD(pic1.p[y][x].dd&0x00800000)!=0)    // high red
        {
        // recolor to green (just for visual check)
        pic1.p[y][x].dd=0x0000FF00;
        // add found point to a list
        p.x=x;
        p.y=y;
        p.a=0.0;
        pnt.add(p);
        // update bounding box
        if (x0>p.x) x0=p.x;
        if (xx<p.x) xx=p.x;
        if (y0>p.y) y0=p.y;
        if (yy<p.y) yy=p.y;
        }
    // center of bounding box for polar sort origin
    x0=0.5*(x0+xx);
    y0=0.5*(y0+yy);
    // draw cross (for visual check)
    x=x0; y=y0; i=16;
    pic1.bmp->Canvas->Pen->Color=clBlack;
    pic1.bmp->Canvas->MoveTo(x-i,y);
    pic1.bmp->Canvas->LineTo(x+i,y);
    pic1.bmp->Canvas->MoveTo(x,y-i);
    pic1.bmp->Canvas->LineTo(x,y+i);
    pic1.save("spiral_fit_0.png");
    // cpmpute polar angle for sorting
    for (i=0;i<pnt.num;i++)
        {
        xx=atan2(pnt[i].y-y0,pnt[i].x-x0);
        if (xx>0.75*M_PI) xx-=2.0*M_PI; // start is > -90 deg
        pnt[i].a=xx;
        }
    // bubble sort by angle (so index in point list can be used as curve parameter)
    for (e=1;e;)
     for (e=0,i=1;i<pnt.num;i++)
      if (pnt[i].a>pnt[i-1].a)
        {
        p=pnt[i];
        pnt[i]=pnt[i-1];
        pnt[i-1]=p;
        e=1;
        }
    // recolor to grayscale gradient (for visual check)
    for (i=0;i<pnt.num;i++)
        {
        x=pnt[i].x;
        y=pnt[i].y;
        pic1.p[y][x].dd=0x00010101*((250*i)/pnt.num);
        }
    pic1.save("spiral_fit_1.png");

    // [fit spiral points with cubic polynomials]
    n =6;                               // recursions for accuracy boost
    m =fabs(pic1.xs+pic1.ys)*1000.0;    // radius for control points fiting
    dm=m/50.0;                          // starting step for approx search
    di=pnt.num/25; if (di<1) di=1;      // skip most points for speed up
    // fit x axis polynomial
    x1=pnt[0          ].x;  // start point of curve
    x2=pnt[  pnt.num-1].x;  // endpoint of curve
    for (a0.init(x1-m,x1+m,dm,n,&ee);!a0.done;a0.step())
    for (a3.init(x2-m,x2+m,dm,n,&ee);!a3.done;a3.step())
        {
        // compute actual polynomial
        x0=a0.a;
        x3=a3.a;
        d1=0.5*(x2-x0);
        d2=0.5*(x3-x1);
        ax[0]=x1;
        ax[1]=d1;
        ax[2]=(3.0*(x2-x1))-(2.0*d1)-d2;
        ax[3]=d1+d2+(2.0*(-x2+x1));
        // compute its distance to points as the fit error e
        for (ee=0.0,i=0;i<pnt.num;i+=di)
            {
            t=double(i)/double(pnt.num-1);
            tt=t*t;
            ttt=tt*t;
            x=ax[0]+(ax[1]*t)+(ax[2]*tt)+(ax[3]*ttt);
            ee+=fabs(pnt[i].x-x);                   // avg error
//          x=fabs(pnt[i].x-x); if (ee<x) ee=x;     // max error
            }
        }
    // compute final x axis polynomial
    x0=a0.aa;
    x3=a3.aa;
    d1=0.5*(x2-x0);
    d2=0.5*(x3-x1);
    ax[0]=x1;
    ax[1]=d1;
    ax[2]=(3.0*(x2-x1))-(2.0*d1)-d2;
    ax[3]=d1+d2+(2.0*(-x2+x1));
    // fit y axis polynomial
    y1=pnt[0          ].y;  // start point of curve
    y2=pnt[  pnt.num-1].y;  // endpoint of curve
    m =fabs(y2-y1)*1000.0;
    di=pnt.num/50; if (di<1) di=1;
    for (a0.init(y1-m,y1+m,dm,n,&ee);!a0.done;a0.step())
    for (a3.init(y2-m,y2+m,dm,n,&ee);!a3.done;a3.step())
        {
        // compute actual polynomial
        y0=a0.a;
        y3=a3.a;
        d1=0.5*(y2-y0);
        d2=0.5*(y3-y1);
        ay[0]=y1;
        ay[1]=d1;
        ay[2]=(3.0*(y2-y1))-(2.0*d1)-d2;
        ay[3]=d1+d2+(2.0*(-y2+y1));
        // compute its distance to points as the fit error e
        for (ee=0.0,i=0;i<pnt.num;i+=di)
            {
            t=double(i)/double(pnt.num-1);
            tt=t*t;
            ttt=tt*t;
            y=ay[0]+(ay[1]*t)+(ay[2]*tt)+(ay[3]*ttt);
            ee+=fabs(pnt[i].y-y);                   // avg error
//          y=fabs(pnt[i].y-y); if (ee<y) ee=y;     // max error
            }
        }
    // compute final y axis polynomial
    y0=a0.aa;
    y3=a3.aa;
    d1=0.5*(y2-y0);
    d2=0.5*(y3-y1);
    ay[0]=y1;
    ay[1]=d1;
    ay[2]=(3.0*(y2-y1))-(2.0*d1)-d2;
    ay[3]=d1+d2+(2.0*(-y2+y1));
    // draw fited curve in Red
    pic1.bmp->Canvas->Pen->Color=clRed;
    pic1.bmp->Canvas->MoveTo(ax[0],ay[0]);
    for (t=0.0;t<=1.0;t+=0.01)
        {
        tt=t*t;
        ttt=tt*t;
        x=ax[0]+(ax[1]*t)+(ax[2]*tt)+(ax[3]*ttt);
        y=ay[0]+(ay[1]*t)+(ay[2]*tt)+(ay[3]*ttt);
        pic1.bmp->Canvas->LineTo(x,y);
        }
    pic1.save("spiral_fit_2.png");

Usé tu imagen de entrada que proporcionaste en OP. Aquí están las etapas de salida

selección de punto espiral:

orden de los puntos por el ángulo polar:

resultado de ajuste final:

Como puede ver, el ajuste no es muy bueno porque:

• Uso solo cúbico para todo el brazo (puede necesitar un mayor grado de polinomio o subdivisión en parches)
• Utilizo la extracción de puntos de la imagen y la curva es gruesa, así que tengo múltiples puntos por cada ángulo polar (tienes el punto original, así que no debería ser un problema) Debería usar thinning algo, pero perezoso para agregarlo ...

En el ejemplo de C ++ utilizo mi propia clase de imagen, así que aquí algunos miembros:

• xs,ys tamaño de la imagen en píxeles
• p[y][x].dd es píxel en la posición (x, y) como tipo entero de 32 bits
• p[y][x].db[4] es el acceso de píxeles por bandas de color (r, g, b, a)
• p.load(filename),p.save(filename) adivina qué ... carga / guarda la imagen
• p.bmp->Canvas es un acceso de mapa de bits GDI para que pueda usar GDI también

El ajuste lo realiza mi clase de búsqueda de aproximación desde:

• Cómo funciona la búsqueda de aproximación

Así que solo copie la class approx desde allí.

La plantilla List<T> es solo una matriz dinámica (lista):

• List<int> q; es lo mismo que int q[];
• q.num tiene la cantidad de elementos dentro
• q.add() agrega un nuevo elemento vacío al final de la lista
• q.add(10) agrega 10 como elemento nuevo al final de la lista

[Notas]

Como ya tiene la lista de puntos, no necesita escanear la imagen de entrada para obtener puntos ... por lo que puede ignorar esa parte del código ...

Si necesita BEZIER en lugar de polinomio de interpolación, puede convertir los puntos de control directamente en:

• Interpolación cúbica vs. BEZIER cúbico

Si la forma de la curva objetivo no es fija, entonces también puedes intentar ajustar directamente la ecuación espiral mediante un círculo paramétrico como la ecuación con centro cambiante y radio variable. Eso debería ser mucho más preciso y la mayoría de los parámetros se pueden calcular sin ajuste.

[Edit1] mejor descripción de la adaptación polinómica de la mina

Estoy usando interpolación cúbica desde el enlace de arriba con estas propiedades:

• para 4 puntos de entrada p0,p1,p2,p3 la curva comienza en p1 ( t=0.0 ) y termina en p2 ( t=1.0 ). los puntos p0,p3 son accesibles ( t=-1.0 t=2.0 ) y aseguran la condición de continuidad entre parches. Entonces, la derivación en p0 y p1 es la misma para todos los parches vecinos. Es lo mismo que fusionar parches BEZIER juntos.

El ajuste polinomial es fácil:

1. Configuré p1,p2 en los puntos finales en espiral

   entonces la curva comienza y termina donde debería

2. Busco p0,p3 cerca de p1,p2 hasta cierta distancia m

   recordando la coincidencia más cercana de la curva polinómica con los puntos originales. Puede usar una distancia promedio o máxima para esto. La clase approx hace todo el trabajo que necesita solo para calcular la distancia ee en cada iteración.

   para m uso múltiples de tamaño de imagen. Si es demasiado grande, perderá precisión (o necesitará más recursiones y ralentizará las cosas), si es demasiado baja puede delimitar el área donde deberían estar los puntos de control y el ajuste se deformará.

   La iteración que inicia el paso dm es parte de m y si el cálculo es demasiado pequeño será lento. Si está demasiado alto, puede perder el mínimo / máximo local, donde la solución está resultando en un ajuste incorrecto.

   Para acelerar el cálculo, estoy usando solo 25 puntos seleccionados de manera uniforme de los puntos (no es necesario usarlos todos), el paso está en di

La separación de las dimensiones x,y es la misma que acaba de cambiar todas las x por y contrario, el código es el mismo.