反转旋转插值 [英] Inversing an interpolation of rotation
问题描述
对于一个项目，我有一个矩阵< float> ，它旋转几度。我无法控制这个过程（假设它正在使用最近邻居），我想要反转这个旋转操作并获得初始矩阵（或非常接近它的矩阵）。
我最初的假设是，如果我用 -angle 旋转旋转的矩阵并裁剪中间部分，我会有原始矩阵但结果表明质量急剧下降。
考虑我的原始矩阵（图中的第一张图像）是10x10矩阵从1到100.我将它旋转+10度，然后旋转-10度。图中的第二张图像是我的结果矩阵。然后我从第二个矩阵的中间进行裁剪，并将其与初始矩阵相关联。
我用1000 1000随机矩阵测试了这个1000;当我用 bicubic 或 bilinear 插值旋转-10度时，平均相关结果约为0.37而最近邻居是0.25。
如果两个插值都是 bilinear 或 bicubic ，那么相关结果约为0.45-0.5。
我想知道是否有办法减少插值造成的损失。请注意，在实际实验中我没有原始图像，我只是估计旋转角度，因此旋转角度估计的精度会导致另一个性能下降。我在网上搜索但找不到任何关于它的信息。
这是我在matlab中的简单测试代码，
res = 0; 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）;
结束
res / 1000
我在 P1.csv 750x1000 矩阵中对 C ++ 进行了小测试。我通过 + 10deg 将其旋转然后返回 -10deg ，并在矩阵中心周围进行双线性插值。
产生的相关性（在 749x749 结果的中间平方）是 0.8275936 所以要么你没有关联相同的数据（可能是矩阵之间的某些偏移），或者你是以某种方式截断结果。例如，我从整数矩阵旋转代码中进行此操作，并且忘记删除整数截断，相关性大约是 0.3 ，这与您的声明类似。
由于我不在 Matlab 这里使用 C ++ 来源，您可以尝试移植或检查您的实施：
// ------------------------ -------------------------------------------------- -
const float deg = M_PI / 180.0;
const float rad = 180.0 / M_PI;
int x0，y0，r0;
矩阵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（x0，y0）CW如果ang> 0
{
int x，y，ix0，iy0，ix1，iy1;
float xx，yy，fx，fy，c，s，q;
B.resize（A.xs，A.ys）;
//圈内核
c = cos（-ang）; S = SIN（-Ang）;
//旋转
for（y = 0; y< A.ys; y ++）
for（x = 0; x< A.xs; x ++）
{
// offset so（0,0）是旋转中心
xx = x-x0;
yy = y-y0;
//旋转（fx，fy）ang
fx = float（（xx * c） - （yy * s））;
fy = float（（xx * s）+（yy * c））;
//抵消并转换为整数和权重
fx + = x0; IX0 =地板（FX）; FX-= IX0; IX1 = IX0 + 1; if（ix1> = A.xs）ix1 = ix0;
fy + = y0; iy0 =地板（年度）; fy- = iy0; iy1 = iy0 + 1; if（iy1> = A.ys）iy1 = iy0;
//双线性插值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;
}否则q = 0;
B [x] [y] = q;
}
}
// ---------------------------------- -----------------------------------------
float correl（matrix& ; A，矩阵& 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）返回0.0;
返回a / b;
}
// --------------------------------------- ------------------------------------
矩阵A 只是动态2D数组（我为此破坏了），如 float A [A.xs] [A.ys]; 其中 xs，ys 是大小。 A.resize（xs，ys）会将矩阵 A 调整为新大小。来源：
// --------------------- -------------------------------------------------- ----
类矩阵
{
public：
int xs，ys;
float ** a; // float a [xs] [ys]
matrix（）{a = NULL; XS = 0; YS = 0; }
矩阵（矩阵和q）{* this = q; }
~matrix（）{free（）; }
矩阵*运算符=（const矩阵* q）{* this = * q;归还这个; }
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];归还这个; }
float * operator []（int x）{return a [x]; };
void free（）{if（a）{if（a [0]）delete [] a [0];删除[] a; } a = NULL; XS = 0; YS = 0; }
void resize（int _xs，int _ys）
{
free（）;
如果（_xs< = 0）返回;
如果（_ys< = 0）返回;
a = new float * [_ xs]; if（a == NULL）返回;
float * aa = new float [_xs * _ys]; if（aa == NULL）返回;
xs = _xs; YS = _ys;
for（int x = 0; x< xs; x ++，aa + = ys）a [x] = aa;
}
};
// -------------------------------------------- -------------------------------
测试如下：
x0 = A.xs>> 1; //旋转中心
y0 = A.ys>> 1;
if（x0< y0）r0 = x0-1;否则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）;
由于采用双线性插值，旋转的单元格会流向相邻的单元格，因此如果需要旋转多次（例如，找出未知角度）然后你应该总是旋转原始矩阵，而不是在子结果矩阵上多次旋转。
这里预览 P1
如您所见，有些功能可用于检测展位的旋转角度和中心旋转...只需找到/交叉匹配 A 和 B 中的孔，然后计算差异角度。旋转计算偏移后你应该得到你所需要的一切...
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
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 ...
这篇关于反转旋转插值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
