java - Vista OpenCV Birdseye sin pérdida de datos

image-processing 3d

Question

Estoy usando OpenCV para obtener una vista panorámica de los cuadros capturados. Esto se realiza proporcionando un patrón de tablero de ajedrez en el avión que formará la vista panorámica.

Aunque parece que la cámara ya está bastante por encima de este plano, necesito que sea perfecta para determinar la relación entre píxeles y centímetros.

En la siguiente fase, los marcos de captura se deforman. Da el resultado esperado:

Sin embargo, al realizar esta transformación, se pierden datos fuera del patrón del tablero de ajedrez. Lo que necesito es rotar la imagen en lugar de deformar un cuadrángulo conocido.

Pregunta: ¿Cómo rotar una imagen por el ángulo de una cámara para que quede de arriba abajo?

Algún código para ilustrar lo que estoy haciendo actualmente:

Size chessboardSize = new Size(12, 8); // Size of the chessboard

Size captureSize = new Size(1920, 1080); // Size of the captured frames

Size viewSize = new Size((chessboardSize.width / chessboardSize.height) * captureSize.height, captureSize.height); // Size of the view

MatOfPoint2f imageCorners; // Contains the imageCorners obtained in a earlier stage

Mat H; // Homography

El código que encuentra las esquinas:

Mat grayImage = new Mat();
//Imgproc.resize(source, temp, new Size(source.width(), source.height()));
Imgproc.cvtColor(source, grayImage, Imgproc.COLOR_BGR2GRAY);
Imgproc.threshold(grayImage, grayImage, 0.0, 255.0, Imgproc.THRESH_OTSU);
imageCorners = new MatOfPoint2f();
Imgproc.GaussianBlur(grayImage, grayImage, new Size(5, 5), 5); 
boolean found = Calib3d.findChessboardCorners(grayImage, chessboardSize, imageCorners, Calib3d.CALIB_CB_NORMALIZE_IMAGE + Calib3d.CALIB_CB_ADAPTIVE_THRESH + Calib3d.CALIB_CB_FILTER_QUADS);

if (found) {

    determineHomography();
}

El código que determina la homografía:

Point[] data = imageCorners.toArray();

if (data.length < chessboardSize.area()) {
    return;
}

Point[] roi = new Point[] {

    data[0 * (int)chessboardSize.width - 0], // Top left
    data[1 * (int)chessboardSize.width - 1], // Top right
    data[((int)chessboardSize.height - 1) * (int)chessboardSize.width - 0], // Bottom left
    data[((int)chessboardSize.height - 0) * (int)chessboardSize.width - 1], // Bottom right
};

Point[] roo = new Point[] {
    new Point(0, 0),
    new Point(viewSize.width, 0),
    new Point(0, viewSize.height),
    new Point(viewSize.width, viewSize.height)
};

MatOfPoint2f objectPoints = new MatOfPoint2f(), imagePoints = new MatOfPoint2f();

objectPoints.fromArray(roo);
imagePoints.fromArray(roi);

Mat H = Imgproc.getPerspectiveTransform(imagePoints, objectPoints);

Finalmente, los fotogramas capturados se están deformando:

Imgproc.warpPerspective(capture, view, H, viewSize);

Answer 1

[Editar2] progreso actualizado

Puede haber más que solo rotación presente, así que probaría esto en su lugar:

1. imagen de preproceso

   Puede aplicar muchos filtros para eliminar el ruido de la imagen o normalizar las condiciones de iluminación (parece que su imagen publicada no lo necesita). Luego, simplemente binarise la imagen para simplificar más pasos. ver relacionado:

   • OpenCV para OCR: cómo calcular los niveles de umbral para OCR de imagen gris

2. detectar puntos de esquina cuadrados

   y almacenar sus coordenadas en alguna matriz con su topología

   double pnt[col][row][2];

   donde (col,row) es el índice del tablero de ajedrez y [2] almacena (x, y). Puede usar int pero double/float evitará conversiones innecesarias y redondeos durante el ajuste ...

   Las esquinas se pueden detectar (a menos que la inclinación / rotación esté cerca de 45 grados) al escanear los píxeles adyacentes diagonales de esta manera:

   Una diagonal debe estar en un color y la otra en diferente. Este patrón detectará un grupo de puntos alrededor del cruce, así que encuentre puntos cercanos y calcule su promedio.

   Si escanea toda la imagen, la parte superior for eje del ciclo también ordenará la lista de puntos, por lo que no es necesario ordenarla más. Después de promediar ordenar / ordenar los puntos a la topología de la cuadrícula (por ejemplo, por dirección entre 2 puntos más cercanos)

3. Topología

   Para hacerlo robusto, uso una imagen girada y sesgada, por lo que la detección de topología es un poco complicada. Después de un rato de elaboración llego a esto:

   1. encontrar el punto p0 cerca del centro de la imagen

      Eso debería garantizar que haya vecinos para ese punto.

   2. encontrar el punto p más cercano

      Pero ignore los puntos diagonales ( |x/y| -> 1 +/- escala de cuadrados). Desde este punto, calcule el vector de primera base, llámelo u por ahora.

   3. encontrar el punto más cercano p

      De la misma manera que el n . ° 2, pero esta vez también ignora los puntos en la dirección +/- u ( |(uv)|/(|u|.|v|) -> 1 +/- sesgo / rotaciones). Desde este punto, calcule el vector de segunda base, llámelo v por ahora.

   4. normalizar u, v

      Elegí que el vector u apunta a +x y v a +y dirección. Entonces vector base con mayor |x| el valor debe ser u y con mayor |y| debería ser v . Por lo tanto, pruebe e intercambie si es necesario. Entonces solo niega si el signo equivocado. Ahora tenemos vectores de base para el centro de la pantalla (más lejos podrían cambiar).

   5. calcular topología

      Establezca el punto p0 como (u=0,v=0) como punto de inicio. Ahora recorra todos los puntos aún no coincidentes p . Para cada cálculo, predice la posición de los vecinos al sumar / restar vectores base de su posición. Luego encuentre el punto más cercano a esta ubicación y, si se encuentra, debe ser vecino, por lo tanto, configure su coordenada (u,v) en +/-1 del punto original p . Ahora actualice los vectores base para estos puntos y repita todo hasta que no se encuentre una nueva coincidencia. El resultado debería ser que la mayoría de los puntos deberían haber calculado sus coordenadas (u,v) que es lo que necesitamos.

   Después de esto, puede encontrar el min(u),min(v) y cambiarlo a (0,0) para que los índices no sean negativos si es necesario.

4. ajustar un polinomio para los puntos de esquina

   por ejemplo algo como:

   pnt[i][j][0]=fx(i,j)
   pnt[i][j][1]=fy(i,j)

   donde fx,fy son funciones polinómicas. Puedes probar cualquier proceso de adaptación. Intenté el ajuste polinómico cúbico con el uso de búsqueda de aproximación, pero el resultado no fue tan bueno como la interpolación bicúbica nativa (posiblemente debido a una distorsión no uniforme de la imagen de prueba), así que cambié a la interpolación bicúbica en lugar de la adaptación. Eso es más simple pero hace que la computación inversa sea muy difícil, pero se puede evitar a costa de la velocidad. Si necesita calcular el inverso de todos modos, vea

   • Tabla de búsqueda 2D compleja inversa

   Estoy usando una simple interpolación cúbica como esta:

   d1=0.5*(pp[2]-pp[0]);
   d2=0.5*(pp[3]-pp[1]);
   a0=pp[1];
   a1=d1;
   a2=(3.0*(pp[2]-pp[1]))-(2.0*d1)-d2;
   a3=d1+d2+(2.0*(-pp[2]+pp[1])); }
   coordinate = a0+(a1*t)+(a2*t*t)+(a3*t*t*t);

   donde pp[0..3] son 4 puntos de control conocidos consecuentes (nuestros cruces de cuadrícula), a0..a3 son coeficientes polinómicos calculados y la coordinate es punto en curva con el parámetro t . Esto se puede ampliar a cualquier cantidad de dimensiones.

   Las propiedades de esta curva son simples, es continua, comenzando en pp[1] y terminando en pp[2] mientras que t=<0.0,1.0> . La continuidad con los segmentos vecinos está asegurada con una secuencia común a todas las curvas cúbicas.

5. remapear píxeles

   simplemente use i,j como valores flotantes con un paso alrededor del 75% del tamaño de píxel para evitar huecos. Luego, simplemente recorra todas las posiciones (i,j) calcule (x,y) y copie píxeles de la imagen de origen en (x,y) a (i*sz,j*sz)+/-offset donde se desea la cuadrícula sz tamaño en píxeles

Aquí el C ++ :

//---------------------------------------------------------------------------
picture pic0,pic1;                          // pic0 - original input image,pic1 output
//---------------------------------------------------------------------------
struct _pnt
    {
    int x,y,n;
    int ux,uy,vx,vy;
    _pnt(){};
    _pnt(_pnt& a){ *this=a; };
    ~_pnt(){};
    _pnt* operator = (const _pnt *a) { x=a->x; y=a->y; return this; };
    //_pnt* operator = (const _pnt &a) { ...copy... return this; };
    };
//---------------------------------------------------------------------------
void vision()
    {
    pic1=pic0;                              // copy input image pic0 to pic1
    pic1.enhance_range();                   // maximize dynamic range of all channels
    pic1.treshold_AND(0,127,255,0);         // binarize (remove gray shades)
    pic1&=0x00FFFFFF;                       // clear alpha channel for exact color matching

    pic1.save("out_binarised.png");

    int i0,i,j,k,l,x,y,u,v,ux,uy,ul,vx,vy,vl;
    int qi[4],ql[4],e,us,vs,**uv;

    _pnt *p,*q,p0;
    List<_pnt> pnt;
    // detect square crossings point clouds into pnt[]
    pnt.allocate(512); pnt.num=0;
    p0.ux=0; p0.uy=0; p0.vx=0; p0.vy=0;
    for (p0.n=1,p0.y=2;p0.y<pic1.ys-2;p0.y++)   // sorted by y axis, each point has usage n=1
     for (      p0.x=2;p0.x<pic1.xs-2;p0.x++)
      if (pic1.p[p0.y-2][p0.x+2].dd==pic1.p[p0.y+2][p0.x-2].dd)
      if (pic1.p[p0.y-1][p0.x+1].dd==pic1.p[p0.y+1][p0.x-1].dd)
      if (pic1.p[p0.y-1][p0.x+1].dd!=pic1.p[p0.y+1][p0.x+1].dd)
      if (pic1.p[p0.y-1][p0.x-1].dd==pic1.p[p0.y+1][p0.x+1].dd)
      if (pic1.p[p0.y-2][p0.x-2].dd==pic1.p[p0.y+2][p0.x+2].dd)
       pnt.add(p0);
    // merge close points (deleted point has n=0)
    for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
     if (p->n)                              // skip deleted points
      for (p0=*p,j=i+1,q=p+1;j<pnt.num;j++,q++) // scan all remaining points
       if (q->n)                            // skip deleted points
        {
        if (q->y>p0.y+4) continue;          // scan only up do y distance <=4 (clods are not bigger then that)
        x=p0.x-q->x; x*=x;                  // compute distance^2
        y=p0.y-q->y; y*=y; x+=y;
        if (x>25) continue;                 // skip too distant points
        p->x+=q->x;                         // add coordinates (average)
        p->y+=q->y;
        p->n++;                             // increase ussage
        q->n=0;                             // mark current point as deleted
        }
    // divide the average coordinates and delete marked points
    for (p=pnt.dat,i=0,j=0;i<pnt.num;i++,p++)
     if (p->n)                              // skip deleted points
        {
        p->x/=p->n;
        p->y/=p->n;
        p->n=1;
        pnt.dat[j]=*p; j++;
        } pnt.num=j;
    // n is now encoded (u,v) so set it as unmatched (u,v) first
    #define uv2n(u,v) ((((v+32768)&65535)<<16)|((u+32768)&65535))
    #define n2uv(n) { u=n&65535; u-=32768; v=(n>>16)&65535; v-=32768; }
    for (p=pnt.dat,i=0;i<pnt.num;i++,p++) p->n=0;
    // p0,i0 find point near middle of image
    x=pic1.xs>>2;
    y=pic1.ys>>2;
    for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
     if ((p->x>=x)&&(p->x<=x+x+x)
       &&(p->y>=y)&&(p->y<=y+y+y)) break;
    p0=*p; i0=i;
    // q,j find closest point to p0
    vl=pic1.xs+pic1.ys; k=0;
    for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
     if (i!=i0)
        {
        x=p->x-p0.x;
        y=p->y-p0.y;
        l=sqrt((x*x)+(y*y));
        if (abs(abs(x)-abs(y))*5<l) continue;   // ignore diagonals
        if (l<=vl) { k=i; vl=l; }               // remember smallest distance
        }
    q=pnt.dat+k; j=k;
    ux=q->x-p0.x;
    uy=q->y-p0.y;
    ul=sqrt((ux*ux)+(uy*uy));
    // q,k find closest point to p0 not in u direction
    vl=pic1.xs+pic1.ys; k=0;
    for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
     if (i!=i0)
        {
        x=p->x-p0.x;
        y=p->y-p0.y;
        l=sqrt((x*x)+(y*y));
        if (abs(abs(x)-abs(y))*5<l) continue;   // ignore diagonals
        if (abs((100*ux*y)/((x*uy)+1))>75) continue;// ignore paralel to u directions
        if (l<=vl) { k=i; vl=l; }               // remember smallest distance
        }
    q=pnt.dat+k;
    vx=q->x-p0.x;
    vy=q->y-p0.y;
    vl=sqrt((vx*vx)+(vy*vy));
    // normalize directions u -> +x, v -> +y
    if (abs(ux)<abs(vx))
        {
        x=j ; j =k ; k =x;
        x=ux; ux=vx; vx=x;
        x=uy; uy=vy; vy=x;
        x=ul; ul=vl; vl=x;
        }
    if (abs(vy)<abs(uy))
        {
        x=ux; ux=vx; vx=x;
        x=uy; uy=vy; vy=x;
        x=ul; ul=vl; vl=x;
        }
    x=1; y=1;
    if (ux<0) { ux=-ux; uy=-uy; x=-x; }
    if (vy<0) { vx=-vx; vy=-vy; y=-y; }
    // set (u,v) encoded in n for already found points
    p0.n=uv2n(0,0);         // middle point
    p0.ux=ux; p0.uy=uy;
    p0.vx=vx; p0.vy=vy;
    pnt.dat[i0]=p0;
    p=pnt.dat+j;            // p0 +/- u basis vector
    p->n=uv2n(x,0);
    p->ux=ux; p->uy=uy;
    p->vx=vx; p->vy=vy;
    p=pnt.dat+k;            // p0 +/- v basis vector
    p->n=uv2n(0,y);
    p->ux=ux; p->uy=uy;
    p->vx=vx; p->vy=vy;

    // qi[k],ql[k] find closest point to p0
    #define find_neighbor                                                       /
    for (ql[k]=0x7FFFFFFF,qi[k]=-1,q=pnt.dat,j=0;j<pnt.num;j++,q++)             /
        {                                                                       /
        x=q->x-p0.x;                                                            /
        y=q->y-p0.y;                                                            /
        l=(x*x)+(y*y);                                                          /
        if (ql[k]>=l) { ql[k]=l; qi[k]=j; }                                     /
        }

    // process all matched points
    for (e=1;e;)
    for (e=0,p=pnt.dat,i=0;i<pnt.num;i++,p++)
     if (p->n)
        {
        // prepare variables
        ul=(p->ux*p->ux)+(p->uy*p->uy);
        vl=(p->vx*p->vx)+(p->vy*p->vy);
        // find neighbors near predicted position p0
        k=0; p0.x=p->x-p->ux; p0.y=p->y-p->uy; find_neighbor; if (ql[k]<<1>ul) qi[k]=-1;    // u-1,v
        k++; p0.x=p->x+p->ux; p0.y=p->y+p->uy; find_neighbor; if (ql[k]<<1>ul) qi[k]=-1;    // u+1,v
        k++; p0.x=p->x-p->vx; p0.y=p->y-p->vy; find_neighbor; if (ql[k]<<1>vl) qi[k]=-1;    // u,v-1
        k++; p0.x=p->x+p->vx; p0.y=p->y+p->vy; find_neighbor; if (ql[k]<<1>vl) qi[k]=-1;    // u,v+1
        // update local u,v basis vectors for found points (and remember them)
        n2uv(p->n); ux=p->ux; uy=p->uy; vx=p->vx; vy=p->vy;
        k=0; if (qi[k]>=0) { q=pnt.dat+qi[k]; if (!q->n) { e=1; q->n=uv2n(u-1,v); q->ux=-(q->x-p->x); q->uy=-(q->y-p->y); } ux=q->ux; uy=q->uy; }
        k++; if (qi[k]>=0) { q=pnt.dat+qi[k]; if (!q->n) { e=1; q->n=uv2n(u+1,v); q->ux=+(q->x-p->x); q->uy=+(q->y-p->y); } ux=q->ux; uy=q->uy; }
        k++; if (qi[k]>=0) { q=pnt.dat+qi[k]; if (!q->n) { e=1; q->n=uv2n(u,v-1); q->vx=-(q->x-p->x); q->vy=-(q->y-p->y); } vx=q->vx; vy=q->vy; }
        k++; if (qi[k]>=0) { q=pnt.dat+qi[k]; if (!q->n) { e=1; q->n=uv2n(u,v+1); q->vx=+(q->x-p->x); q->vy=+(q->y-p->y); } vx=q->vx; vy=q->vy; }
        // copy remembered local u,v basis vectors to points where are those missing
        k=0; if (qi[k]>=0) { q=pnt.dat+qi[k]; if (!q->vy) { q->vx=vx; q->vy=vy; }}
        k++; if (qi[k]>=0) { q=pnt.dat+qi[k]; if (!q->vy) { q->vx=vx; q->vy=vy; }}
        k++; if (qi[k]>=0) { q=pnt.dat+qi[k]; if (!q->ux) { q->ux=ux; q->uy=uy; }}
        k++; if (qi[k]>=0) { q=pnt.dat+qi[k]; if (!q->ux) { q->ux=ux; q->uy=uy; }}
        }
    // find min,max (u,v)
    ux=0; uy=0; vx=0; vy=0;
    for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
     if (p->n)
        {
        n2uv(p->n);
        if (ux>u) ux=u;
        if (vx>v) vx=v;
        if (uy<u) uy=u;
        if (vy<v) vy=v;
        }
    // normalize (u,v)+enlarge and create topology table
    us=uy-ux+1;
    vs=vy-vx+1;
    uv=new int*[us];
    for (u=0;u<us;u++) uv[u]=new int[vs];
    for (u=0;u<us;u++)
     for (v=0;v<vs;v++)
      uv[u][v]=-1;
    for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
     if (p->n)
        {
        n2uv(p->n);
        u-=ux; v-=vx;
        p->n=uv2n(u,v);
        uv[u][v]=i;
        }
    // bi-cubic interpolation
    double a0,a1,a2,a3,d1,d2,pp[4],qx[4],qy[4],t,fu,fv,fx,fy;
    // compute cubic curve coefficients a0..a3 from 1D points pp[0..3]
    #define cubic_init { d1=0.5*(pp[2]-pp[0]); d2=0.5*(pp[3]-pp[1]); a0=pp[1]; a1=d1; a2=(3.0*(pp[2]-pp[1]))-(2.0*d1)-d2; a3=d1+d2+(2.0*(-pp[2]+pp[1])); }
    // compute cubic curve cordinates =f(t)
    #define cubic_xy (a0+(a1*t)+(a2*t*t)+(a3*t*t*t));
    // safe access to grid (u,v) point copies it to p0
    // points utside grid are computed by mirroring
    #define point_uv(u,v)                                                       /
        {                                                                       /
        if ((u>=0)&&(u<us)&&(v>=0)&&(v<vs)) p0=pnt.dat[uv[u][v]];               /
        else{                                                                   /
            int uu=u,vv=v;                                                      /
            if (uu<0) uu=0;                                                     /
            if (uu>=us) uu=us-1;                                                /
            if (vv<0) vv=0;                                                     /
            if (vv>=vs) vv=vs-1;                                                /
            p0=pnt.dat[uv[uu][vv]];                                             /
            uu=u-uu; vv=v-vv;                                                   /
            p0.x+=(uu*p0.ux)+(vv*p0.vx);                                        /
            p0.y+=(uu*p0.uy)+(vv*p0.vy);                                        /
            }                                                                   /
        }

    //----------------------------------------
    //--- Debug draws: -----------------------
    //----------------------------------------

    // debug recolor white to gray to emphasize debug render
    pic1.recolor(0x00FFFFFF,0x00404040);

    // debug draw basis vectors
    for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
        {
        pic1.bmp->Canvas->Pen->Color=clRed;
        pic1.bmp->Canvas->Pen->Width=1;
        pic1.bmp->Canvas->MoveTo(p->x,p->y);
        pic1.bmp->Canvas->LineTo(p->x+p->ux,p->y+p->uy);
        pic1.bmp->Canvas->Pen->Color=clBlue;
        pic1.bmp->Canvas->MoveTo(p->x,p->y);
        pic1.bmp->Canvas->LineTo(p->x+p->vx,p->y+p->vy);
        pic1.bmp->Canvas->Pen->Width=1;
        }

    // debug draw crossings
    AnsiString s;
    pic1.bmp->Canvas->Font->Height=12;
    pic1.bmp->Canvas->Brush->Style=bsClear;
    for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
        {
        n2uv(p->n);
        if (p->n)
            {
            pic1.bmp->Canvas->Font->Color=clWhite;
            s=AnsiString().sprintf("%i,%i",u,v);
            }
        else{
            pic1.bmp->Canvas->Font->Color=clGray;
            s=i;
            }
        x=p->x-(pic1.bmp->Canvas->TextWidth(s)>>1);
        y=p->y-(pic1.bmp->Canvas->TextHeight(s)>>1);
        pic1.bmp->Canvas->TextOutA(x,y,s);
        }
    pic1.bmp->Canvas->Brush->Style=bsSolid;

    pic1.save("out_topology.png");

    // debug draw of bi-cubic interpolation fit/coveradge with half square step
    pic1=pic0;
    pic1.treshold_AND(0,200,0x40,0);            // binarize (remove gray shades)
    pic1.bmp->Canvas->Pen->Color=clAqua;
    pic1.bmp->Canvas->Brush->Color=clBlue;
    for (fu=-1;fu<double(us)+0.01;fu+=0.5)
     for (fv=-1;fv<double(vs)+0.01;fv+=0.5)
        {
        u=floor(fu);
        v=floor(fv);
        // 4x cubic curve in v direction
        t=fv-double(v);
        for (i=0;i<4;i++)
            {
            point_uv(u-1+i,v-1); pp[0]=p0.x;
            point_uv(u-1+i,v+0); pp[1]=p0.x;
            point_uv(u-1+i,v+1); pp[2]=p0.x;
            point_uv(u-1+i,v+2); pp[3]=p0.x;
            cubic_init; qx[i]=cubic_xy;
            point_uv(u-1+i,v-1); pp[0]=p0.y;
            point_uv(u-1+i,v+0); pp[1]=p0.y;
            point_uv(u-1+i,v+1); pp[2]=p0.y;
            point_uv(u-1+i,v+2); pp[3]=p0.y;
            cubic_init; qy[i]=cubic_xy;
            }
        // 1x cubic curve in u direction on the resulting 4 points
        t=fu-double(u);
        for (i=0;i<4;i++) pp[i]=qx[i]; cubic_init; fx=cubic_xy;
        for (i=0;i<4;i++) pp[i]=qy[i]; cubic_init; fy=cubic_xy;
        t=1.0;
        pic1.bmp->Canvas->Ellipse(fx-t,fy-t,fx+t,fy+t);
        }
    pic1.save("out_fit.png");

    // linearizing of original image
    DWORD col;
    double grid_size=32.0;  // linear grid square size in pixels
    double grid_step=0.01;  // u,v step <= 1 pixel

    pic1.resize((us+1)*grid_size,(vs+1)*grid_size); // resize target image
    pic1.clear(0);                                  // clear target image
    for (fu=-1;fu<double(us)+0.01;fu+=grid_step)    // copy/transform source image to target
     for (fv=-1;fv<double(vs)+0.01;fv+=grid_step)
        {
        u=floor(fu);
        v=floor(fv);
        // 4x cubic curve in v direction
        t=fv-double(v);
        for (i=0;i<4;i++)
            {
            point_uv(u-1+i,v-1); pp[0]=p0.x;
            point_uv(u-1+i,v+0); pp[1]=p0.x;
            point_uv(u-1+i,v+1); pp[2]=p0.x;
            point_uv(u-1+i,v+2); pp[3]=p0.x;
            cubic_init; qx[i]=cubic_xy;
            point_uv(u-1+i,v-1); pp[0]=p0.y;
            point_uv(u-1+i,v+0); pp[1]=p0.y;
            point_uv(u-1+i,v+1); pp[2]=p0.y;
            point_uv(u-1+i,v+2); pp[3]=p0.y;
            cubic_init; qy[i]=cubic_xy;
            }
        // 1x cubic curve in u direction on the resulting 4 points
        t=fu-double(u);
        for (i=0;i<4;i++) pp[i]=qx[i]; cubic_init; fx=cubic_xy; x=fx;
        for (i=0;i<4;i++) pp[i]=qy[i]; cubic_init; fy=cubic_xy; y=fy;
        // here (x,y) contains source image coordinates coresponding to grid (fu,fv) so copy it to col
        col=0; if ((x>=0)&&(x<pic0.xs)&&(y>=0)&&(y<pic0.ys)) col=pic0.p[y][x].dd;
        // compute liner image coordinates (x,y) by scaling (fu,fv)
        fx=(fu+1.0)*grid_size; x=fx;
        fy=(fv+1.0)*grid_size; y=fy;
        // copy col to it
        if ((x>=0)&&(x<pic1.xs)&&(y>=0)&&(y<pic1.ys)) pic1.p[y][x].dd=col;
        }
    pic1.save("out_linear.png");

    // release memory and cleanup macros
    for (u=0;u<us;u++) delete[] uv[u]; delete[] uv;
    #undef uv2n
    #undef n2uv
    #undef find_neighbor
    #undef cubic_init
    #undef cubic_xy
    #undef point_uv(u,v)
    }
//---------------------------------------------------------------------------

Lo siento, sé que es mucho código, pero al menos lo comenté tanto como pude. El código no está optimizado en aras de la simplicidad y la capacidad de comprensión, la linealización de la imagen final se puede escribir mucho más rápido. También elegí el grid_size y grid_step en esa parte del código manualmente. En su lugar, debe calcularse a partir de la imagen y las propiedades físicas conocidas.

Utilizo mi propia clase de picture para imágenes, por lo que algunos miembros son:

• xs,ys tamaño de la imagen en píxeles
• p[y][x].dd es un píxel en la posición (x,y) como tipo entero de 32 bits
• clear(color) : borra toda la imagen
• resize(xs,ys) : redimensiona la imagen a una nueva resolución
• bmp - Mapa de bits GDI encapsulado VCL con acceso a Canvas

También uso la plantilla de lista dinámica mía, así que:

• List<double> xxx; es lo mismo que double xxx[];
• xxx.add(5); agrega 5 al final de la lista
• Elemento de matriz de acceso xxx[7] (seguro)
• xxx.dat[7] matriz de acceso xxx.dat[7] (acceso directo inseguro pero rápido)
• xxx.num es el tamaño real utilizado de la matriz
• xxx.reset() borra la matriz y establece xxx.num = 0
• xxx.allocate(100) preasigna espacio para 100 elementos

Aquí están las imágenes de salida del resultado secundario. Para hacer las cosas más robustas, cambié la imagen de entrada a una más distorsionada:

Para hacerlo visualmente más agradable, volví a colorear el blanco con el gris. Las líneas rojas son bases u locales y las azules son los vectores base v locales. Los números de vectores 2D blancos son coordenadas de topología (u,v) y los números escalares grises son índices de cruces en pnt[] para topología pero puntos no coincidentes.

[Notas]

Este enfoque no funcionará para rotaciones cercanas a 45 grados. Para tales casos, debe cambiar la detección de cruce de patrón cruzado a más y también las condiciones de topología y las ecuaciones cambian un poco. Sin mencionar u, v selección de dirección.