matlab - suptitle - Inversión de una interpolación de rotación

suptitle matlab

Question

Para un proyecto, tengo una matrix<float> que se gira unos pocos grados. No tengo control sobre este proceso (supongo que está usando el vecino más cercano), quiero revertir esta operación de rotación y obtener la matriz inicial (o una matriz muy cercana a ella).

Mi suposición inicial fue que si rotar la matriz girada con un -angle y recortar la parte media, tendría la matriz original, pero los resultados indican que la calidad disminuye drásticamente.

Considere que mi matriz original (la primera imagen en la figura ) es una matriz de 10x10 de 1 a 100. La giro +10 grados, luego -10 grados. La segunda imagen en la figura es mi matriz resultante. Luego recojo desde el medio de la segunda matriz y la correlaciono con la matriz inicial .

Probé esto con 1000 matrices aleatorias de 1000 * 1000; cuando giro -10 grados con interpolación bicubic o bilinear , el resultado de la correlación promedio es alrededor de 0.37 mientras que el nearest neighbor es 0.25.

Si ambas interpolaciones son bilinear o bicubic , el resultado de la correlación es alrededor de 0.45-0.5.

Me pregunto si hay una manera de minimizar la pérdida causada por la interpolación. Tenga en cuenta que en el experimento real no tengo la imagen original, solo estoy estimando el ángulo de rotación, por lo que hay otra caída de rendimiento causada por la precisión de la estimación del ángulo de rotación. Busqué en línea pero no pude encontrar nada al respecto.

Aquí está mi código de prueba simple en matlab,

res = 0;
for i = 0:1000
a = uint8(rand(1000,1000)*255);
arr = imrotate(imrotate(a,10, ''bicubic''), -10, ''bicubic'');
[r1,c1] = size(a);
[r2,c2] = size(arr);
rd = ceil((c2-c1)/2);
cd = ceil((r2-r1)/2);
c_arr = arr(cd:end-cd, rd:end-rd);
res = res+corr2(a, c_arr);
end
res/1000

Answer 1

Hice una pequeña prueba en C ++ en su matriz P1.csv 750x1000 . Lo +10deg luego retrocedí en -10deg con la interpolación bilineal alrededor del centro de la matriz.

La correlación resultante (en el 749x749 medio de resultados de 0.8275936 ) es 0.8275936 Entonces, o bien no está correlacionando los mismos datos (quizás algún desplazamiento entre matrices) o está truncando el resultado de alguna manera. Por ejemplo, hago esto desde mi código de rotación de la matriz entera y, si bien olvido eliminar el entero que trunca, la correlación fue de alrededor de 0.3 que es similar a sus afirmaciones.

Como no utilizo Matlab aquí, mi fuente de C ++ puede intentar realizar un puerto o verificar con sus implementaciones:

//---------------------------------------------------------------------------
const float deg=M_PI/180.0;
const float rad=180.0/M_PI;
int x0,y0,r0;
matrix A,B,C;
float c=0.0,ang=10.0*deg;
//---------------------------------------------------------------------------
void rotcw(matrix &B,matrix &A,int x0,int y0,float ang) // rotate A -> B by angle ang around (x0,y0) CW if ang>0
{
int x,y,ix0,iy0,ix1,iy1;
float xx,yy,fx,fy,c,s,q;
B.resize(A.xs,A.ys);
// circle kernel
c=cos(-ang); s=sin(-ang);
// rotate
for (y=0;y<A.ys;y++)
for (x=0;x<A.xs;x++)
{
// offset so (0,0) is center of rotation
xx=x-x0;
yy=y-y0;
// rotate (fx,fy) by ang
fx=float((xx*c)-(yy*s));
fy=float((xx*s)+(yy*c));
// offset back and convert to ints and weights
fx+=x0; ix0=floor(fx); fx-=ix0; ix1=ix0+1; if (ix1>=A.xs) ix1=ix0;
fy+=y0; iy0=floor(fy); fy-=iy0; iy1=iy0+1; if (iy1>=A.ys) iy1=iy0;
// bilinear interpolation A[fx][fy] -> B[x][y]
if ((ix0>=0)&&(ix0<A.xs)&&(iy0>=0)&&(iy0<A.ys))
{
xx=float(A[ix0][iy0])+(float(A[ix1][iy0]-A[ix0][iy0])*fx);
yy=float(A[ix0][iy1])+(float(A[ix1][iy1]-A[ix0][iy1])*fx);
xx=xx+((yy-xx)*fy); q=xx;
} else q=0;
B[x][y]=q;
}
}
//---------------------------------------------------------------------------
float correl(matrix &A,matrix &B,int x0,int y0,int x1,int y1)
{
int x,y;
float sxy=0.0,sx=0.0,sy=0.0,sxx=0.0,syy=0.0,n=(x1-x0+1)*(y1-y0+1),a,b;
for (x=x0;x<=x1;x++)
for (y=y0;y<=y1;y++)
{
a=A[x][y];
b=B[x][y];
sx+=a; sxx+=a*a;
sy+=b; syy+=b*b;
sxy+=a*b;
}
a=(n*sxy)-(sx*sy);
b=sqrt((n*sxx)-(sx*sx))*sqrt((n*syy)-(sy*sy));
if (fabs(b)<1e-10) return 0.0;
return a/b;
}
//---------------------------------------------------------------------------

matrix A es solo una matriz 2D dinámica (busqué esto) como float A[A.xs][A.ys]; donde xs,ys es el tamaño. A.resize(xs,ys) cambiará el tamaño de la matriz A a un nuevo tamaño. Aquí fuente:

//---------------------------------------------------------------------------
class matrix
{
public:
int xs,ys;
float **a; // float a[xs][ys]
matrix() { a=NULL; xs=0; ys=0; }
matrix(matrix& q) { *this=q; }
~matrix() { free(); }
matrix* operator = (const matrix *q) { *this=*q; return this; }
matrix* operator = (const matrix &q) { resize(q.xs,q.ys); for (int x=0;x<xs;x++) for (int y=0;y<ys;y++) a[x][y]=q.a[x][y]; return this; }
float* operator[] (int x) { return a[x]; };
void free() { if (a) { if (a[0]) delete[] a[0]; delete[] a; } a=NULL; xs=0; ys=0; }
void resize(int _xs,int _ys)
{
free();
if (_xs<=0) return;
if (_ys<=0) return;
a=new float*[_xs]; if (a==NULL) return;
float *aa=new float[_xs*_ys]; if (aa==NULL) return;
xs=_xs; ys=_ys;
for (int x=0;x<xs;x++,aa+=ys) a[x]=aa;
}
};
//---------------------------------------------------------------------------

La prueba se ve así:

x0=A.xs>>1; // center for rotation
y0=A.ys>>1;
if (x0<y0) r0=x0-1; else r0=y0-1; // mid square size for correltaion
rotcw(B,A,x0,y0,+ang);
rotcw(C,B,x0,y0,-ang);
c=correl(A,C,x0-r0,y0-r0,x0+r0,y0+r0);

Debido a la interpolación bilineal, las células giradas sangran a las células vecinas, por lo que si necesita girar muchas veces (por ejemplo, para conocer el ángulo desconocido), siempre debe girar la matriz original en lugar de aplicar la rotación varias veces en la matriz de resultados secundarios.

Aquí una vista previa para su P1

en la matriz original izquierda A en la matriz media girada B por +10deg CW y en la matriz derecha C girada hacia atrás por -10deg CW . Los píxeles azules son positivos y los rojos son negativos. El rectángulo verde es el área correlacionada (sqrt del área superpuesta cuadrada)

[Edit1] Juego un poco con el color

deje a0=-13.487; a1=9.3039; a0=-13.487; a1=9.3039; ser los valores mínimo y máximo de su matriz A Luego, para calcular el color RGB a partir de cualquier valor de A,B o C , utilicé esto:

DWORD col(float x)
{
DWORD c; int sh;
if (x>=0) { sh= 0; x/=a1; } // positive values in Blue
else { sh=16; x/=a0; } // negative values in Red
x*=255.0*50.0; // 50.0x saturated to emphasize used values
c=x; if (c>255) c=255; // clamp to 8bit per channel
return c<<sh;
}

Y aquí el resultado recolored:

Como puede ver, hay características que podrían usarse para detectar el ángulo de rotación y el centro de rotación de la cabina ... Simplemente ubique / cruce los agujeros en A y B y luego calcule el ángulo de diferencia. Después de la rotación, calcula el desplazamiento y deberás obtener todo lo que necesites ...