利用多项式实现图像几何校正(Matlab实现)

放弃爱

1.原理简述：

1.   根据两幅图像中的一些已知对应点（控制点对），建立函数关系式，通过坐标变换，实现失真图像的几何校正。
   
2. 
   
3.                                    ![1704593-20190607214039234-614794812.jpg][]
   
4. 
   
5. 
   
6. 
   
7.     设两幅图像坐标系统之间畸变关系能用解析式来描述：
   
8. 
   
9.                           ![1704593-20190607214412623-284507974.png][]
   
10. 
    
11.      根据上述的函数关系，可以依次计算畸变图像每个像素的矫正坐标值，保持各像素值不变，这样生成一幅矫正图像。

复制代码

2.实现过程：

1. （1）因此首先要得到多项式，matlab提供了拟合多项式的函数 Isqcurvefit,

复制代码

** 格式**:lsqcurvefit(f,a,x,y)

1.    f:符号函数句柄
   
2. 
   
3.             a:最开始预估的值(预拟合的未知参数的估计值)。如上面的问题如果我们预估A为1,B为2,则a=\[1 2\]
   
4. 
   
5.             x:我们已经获知的x的值
   
6. 
   
7.             y:我们已经获知的x对应的y的值
   
8. 
   
9.            函数的返回值为对应多项式系数组成的一维数组。

复制代码

示例（二次多项式）：建立由校正图像到畸变图像的函数

1. function [F] = fun(a,b)
   
2. x = b(:,1);
   
3. y = b(:,2);
   
4. F = a(1)+a(2)*x+a(3)*y+a(4)*x.^2+a(5)*x.*y+a(6)*y.^2;
   
5. end
   
6. 
   
7. 
   
8. x0 = fixedPoints(:,1);
   
9. y0 = fixedPoints(:,2);
   
10. x1 = movingPoints(:,1);
    
11. y1 = movingPoints(:,2);
    
12. data = [x1,y1];
    
13. a = [1 1 1 1 1 1];
    
14. a1 = lsqcurvefit('fun',a,data,x0);
    
15. a2 = lsqcurvefit('fun',a,data,y0);
    
16. 
    
17.      （2）根据得到的二项式，由校正图像每个像素坐标（x，y）出发，算出在已知畸变图像上的对应坐标（x'，y'），使像元一一对应，赋予校正图像对应畸变图像的像元

复制代码

的像素值，最终得到校正图像。
示例：

1. J2 = uint8(zeros(size(J)));
   
2. for rgb = 1:3
   
3.     for i = 1:m
   
4.         for j = 1:n
   
5.             if round(fun(a1,[i,j]))>=1 && round(fun(a1,[i,j]))<=c && round(fun(a2,[i,j]))>=1 && round(fun(a2,[i,j]))<=d
   
6.                  J2(i,j,rgb) = J1(round(fun(a1,[i,j])),round(fun(a2,[i,j])),rgb);
   
7.             end
   
8.         end
   
9.     end
   
10. end

复制代码

** ** 这样生成的图像像素分布不规则，会出现像素挤压、疏密不均的现象，因此最后还需对不规则图像进行灰度内插生成规则的栅格图像。
源码：

1. I = imread('sp.tif');
   
2. J = imread('tm.tif');
   
3. [m,n] = size(J);
   
4. [o,p] = size(I);
   
5. %cpselect(J,I);
   
6. %xlswrite('data1.xls',fixedPoints);
   
7. %xlswrite('data2.xls',movingPoints);
   
8. 
   
9. 
   
10. %%重采样
    
11. J1 = J(1:m/o:end,1:n/p:end,1:3);
    
12. [c,d,q]= size(J1);
    
13. 
    
14. fixedPoints = xlsread('data1.xls');
    
15. movingPoints = xlsread('data2.xls');
    
16. x0 = fixedPoints(:,1);
    
17. y0 = fixedPoints(:,2);
    
18. x1 = movingPoints(:,1);
    
19. y1 = movingPoints(:,2);
    
20. data = [x1,y1];
    
21. a = [1 1 1 1 1 1];
    
22. a1 = lsqcurvefit('fun',a,data,x0);
    
23. a2 = lsqcurvefit('fun',a,data,y0);
    
24. %多项式几何校正
    
25. J2 = uint8(zeros(size(J)));
    
26. for rgb = 1:3
    
27.     for i = 1:m
    
28.         for j = 1:n
    
29.             if round(fun(a1,[i,j]))>=1 && round(fun(a1,[i,j]))<=c && round(fun(a2,[i,j]))>=1 && round(fun(a2,[i,j]))<=d
    
30.                  J2(i,j,rgb) = J1(round(fun(a1,[i,j])),round(fun(a2,[i,j])),rgb);
    
31.             end
    
32.     %           J2(round(fun(a1,[i,j])),round(fun(a2,[i,j]))) = J(i,j);
    
33.     %           end
    
34.         end
    
35.     end
    
36. end
    
37. [x,y] = size(J2);
    
38. 
    
39. %根据数据游标取截取区域的左上方点
    
40. J3 = imcrop(I,[98 180 60*o/x 60*p/y]);
    
41. J4 = imcrop(J2,[41 80 60 60]);
    
42. [k,y,z] = size(J3);
    
43. [h,t,e] = size(J4);
    
44. 
    
45. %%重采样
    
46. %双线性内插法
    
47. u = h/k;
    
48. v = t/y;
    
49. itp = uint8(zeros(k,y,3));
    
50. for rgb = 1:3
    
51.     for i = ceil(1/u):k-1
    
52.         iu = floor(i*u);
    
53.         for j = ceil(1/v):y-1
    
54.             jv = floor(j*v);
    
55.             itp(i,j,rgb) = (1-(i*u-iu))*(1-(j*v-jv))*J4(iu,jv,rgb)...
    
56.                        +(1-(i*u-iu))*(j*v-jv)*J4(iu,jv+1,rgb)...
    
57.                        +(i*u-iu)*(1-(j*v-jv))*J4(iu+1,jv,rgb)...
    
58.                        +(i*u-iu)*(j*v-jv)*J4(iu+1,jv+1,rgb);
    
59.         end
    
60.     end
    
61. end
    
62. %去黑边
    
63. for rgb = 1:3
    
64.     for i = 1:3
    
65.         for j = 1:175
    
66.           itp(i,j,rgb) = J4(ceil(i*u),ceil(j*v),rgb);
    
67.           itp(145,j,rgb) = J4(43,ceil(j*v),rgb);
    
68.         end
    
69.     end
    
70.     for j = 1:2
    
71.         for i = 1:145
    
72.            itp(i,j,rgb) = J4(ceil(i*u),ceil(j*v),rgb);
    
73.            itp(i,175,rgb) = J4(ceil(i*u),61,rgb);
    
74.         end
    
75.     end
    
76. end
    
77. imwrite(J3,'Core1.bmp','bmp');
    
78. imwrite(itp,'Core2.bmp','bmp');
    
79. 
    
80. subplot(231),imshow(J),title('待配准图像');
    
81. subplot(232),imshow(I),title('基准图像');
    
82. subplot(233),imshow(J2),title('多项式几何校正后');        
    
83. subplot(234),imshow(J3),title('基准影像裁剪后');
    
84. subplot(235),imshow(itp),title('校正影像裁剪重采样后');
    
85. 
    
86. 
    
87. 
    
88. % %基准图重采样
    
89. % kh = zuixiaogongbeishu(k,h);
    
90. % yt = zuixiaogongbeishu(y,t);
    
91. % u = h/kh;v = t/yt;
    
92. % J5 = J3(1:k/kh:end,1:y/yt:end,1:3);
    
93. % %配准图 双线性内插法重采样
    
94. % itp = uint8(zeros(kh,yt,3));
    
95. % for rgb = 1:3
    
96. %     for i = floor(1/u):kh-1
    
97. %         iu = floor(i*u);
    
98. %         for j = floor(1/v):yt-1
    
99. %             jv = floor(j*v);
    
100. %             itp(i,j,rgb) = (1-(i*u-iu))*(1-(j*v-jv))*J4(iu,jv,rgb)...
     
101. %                        +(1-(i*u-iu))*(j*v-jv)*J4(iu,jv+1,rgb)...
     
102. %                        +(i*u-iu)*(1-(j*v-jv))*J4(iu+1,jv,rgb)...
     
103. %                        +(i*u-iu)*(j*v-jv)*J4(iu+1,jv+1,rgb);
     
104. %         end
     
105. %     end
     
106. % end
     
107. % %去黑边
     
108. % for rgb = 1:3
     
109. %     for i = 1:144
     
110. %         for j = 1:10675
     
111. %           itp(i,j,rgb) = J4(ceil(i*u),ceil(j*v),rgb);
     
112. %         end
     
113. %     end
     
114. %     for j = 1:175
     
115. %         for i = 1:6235
     
116. %            itp(i,j,rgb) = J4(ceil(i*u),ceil(j*v),rgb);
     
117. %         end
     
118. %     end
     
119. % end
     
120. %
     
121. % itp1 = uint8(zeros(6193,10615,3));
     
122. % itp1(1:6193,1:10615,1:3) = itp(1:6193,1:10615,1:3);
     
123. % imwrite(J5,'Crop1.bmp','bmp');
     
124. % J5 = imread('Crop1.bmp');
     
125. % imwrite(itp1,'Crop2.bmp','bmp');
     
126. % J6 = imread('Crop2.bmp');
     
127. 
     
128. % subplot(231),imshow(J),title('待配准图像');
     
129. % subplot(232),imshow(I),title('基准图像');
     
130. % subplot(233),imshow(J2),title('多项式几何校正后');        
     
131. % subplot(234),imshow(J5),title('基准影像裁剪重采样后');
     
132. % subplot(235),imshow(J6),title('校正影像裁剪重采样后');
     
133. 
     
134. % a = imread('基准.bmp');
     
135. % b = imread('重采样后图像.bmp');
     
136. % c = imcrop(a,[1,100,100,100]);
     
137. % d = imcrop(b,[1,100,100,100]);
     
138. % imwrite(c,'Core3.bmp','bmp');
     
139. % imwrite(d,'Core4.bmp','bmp');

复制代码