Inversing an interpolation of rotation

Question

For a project, I have a matrix<float> which is rotated few degrees. I have no control over this process (assume it is using nearest neighbour), I want to reverse this rotation operation and obtain the initial matrix (or a matrix very close to it).

My initial assumption was if I rotate the rotated matrix with -angle and crop the middle part, I'd have the original matrix but the results indicate the quality drops dramatically.

Consider my original matrix (the first image in the figure) is 10x10 matrix from 1 to 100. I rotate it +10 degrees, then -10 degrees. The second image in the figure is my resulting matrix. Then I crop from the middle of the second matrix and correlate it with the initial matrix.

I tested this with 1000 random matrix of 1000*1000; when I rotate -10 degrees with bicubic or bilinear interpolation, the average correlation result is around 0.37 whereas nearest neighbor is 0.25.

If both interpolations are bilinear or bicubic, then the correlation result is around 0.45-0.5.

I'm wondering if there is a way to minimize the loss caused by interpolation. Note that in the real experiment I don't have the original image, I'm just estimating rotation angle, so there is another performance drop caused by the precision of the rotation angle estimation. I searched online but couldn't find anything about it.

Here is my simple test code in 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

I did small test in C++ on your P1.csv 750x1000 matrix. I rotated it by +10deg then back by -10deg with bilinear interpolation around matrix center.

Resulting correlation (on the 749x749 mid square of result) is 0.8275936 So either you are not correlating the same data (perhaps some offset between matrices) or you are truncating the result somehow. For example I make this from my integer matrix rotation code and while forget to remove integer truncating the correlation was around 0.3 which is similar to your claims.

As I do not use Matlab here my C++ source you can try to port or check with your implementations:

//---------------------------------------------------------------------------
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 is just dynamic 2D array (I busted for this) like float A[A.xs][A.ys]; where xs,ys is the size. A.resize(xs,ys) will resize matrix A to new size. Here source:

//---------------------------------------------------------------------------
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;
        }
    };
//---------------------------------------------------------------------------

The test looks like this:

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);

Due to bilinear interpolation the rotated cells are bleeding to neighboring cells so if you need to rotate many times (for example to find out the unknown angle) then you should always rotate the original matrix instead of applying rotation multiple times on sub-result matrix.

Here preview for your P1

on the left original matrix A in the middle rotated matrix B by +10deg CW and on right matrix C rotated back by -10deg CW. Blue pixels are positive and red pixels are negative values. The green rectangle is correlated area (sqrt of square overlapped area)

[Edit1] I play with the coloring a bit

let a0=-13.487; a1=9.3039; be the min and max values from your A matrix. Then to compute RGB color from any value from A,B or C I used this:

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;
    }

And here the recolored result:

As you can see there are Features that could be used to detect booth the rotation angle and center of rotation ... Just locate/cross match the holes in A and B and then compute the difference angle. After rotation compute offset and you should get all you need ...