1.原理簡述：
根據兩幅圖像中的一些已知對應點（控制點對），建立函數關系式，通過坐標變換，實現失真圖像的幾何校正。
設兩幅圖像坐標系統之間畸變關系能用解析式來描述：
根據上述的函數關系，可以依次計算畸變圖像每個像素的矯正坐標值，保持各像素值不變，這樣生成一幅矯正圖像。
2.實現過程：
（1）因此首先要得到多項式，matlab提供了擬合多項式的函數 Isqcurvefit,
格式:lsqcurvefit(f,a,x,y)
f:符號函數句柄
a:最開始預估的值(預擬合的未知參數的估計值)。如上面的問題如果我們預估A為1,B為2,則a=[1 2]
x:我們已經獲知的x的值
y:我們已經獲知的x對應的y的值
函數的返回值為對應多項式系數組成的一維數組。
示例（二次多項式）：建立由校正圖像到畸變圖像的函數
function [F] = fun(a,b) x = b(:,1); y = b(:,2); F = a(1)+a(2)*x+a(3)*y+a(4)*x.^2+a(5)*x.*y+a(6)*y.^2; end x0 = fixedPoints(:,1); y0 = fixedPoints(:,2); x1 = movingPoints(:,1); y1 = movingPoints(:,2); data = [x1,y1]; a = [1 1 1 1 1 1]; a1 = lsqcurvefit('fun',a,data,x0); a2 = lsqcurvefit('fun',a,data,y0);
（2）根據得到的二項式，由校正圖像每個像素坐標（x，y）出發，算出在已知畸變圖像上的對應坐標（x'，y'），使像元一一對應，賦予校正圖像對應畸變圖像的像元
的像素值，最終得到校正圖像。
示例：
J2 = uint8(zeros(size(J))); for rgb = 1:3 for i = 1:m for j = 1:n 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 J2(i,j,rgb) = J1(round(fun(a1,[i,j])),round(fun(a2,[i,j])),rgb); end end end end
這樣生成的圖像像素分布不規則，會出現像素擠壓、疏密不均的現象，因此最后還需對不規則圖像進行灰度內插生成規則的柵格圖像。
源碼：
I = imread('sp.tif'); J = imread('tm.tif'); [m,n] = size(J); [o,p] = size(I); %cpselect(J,I); %xlswrite('data1.xls',fixedPoints); %xlswrite('data2.xls',movingPoints); %%重采樣 J1 = J(1:m/o:end,1:n/p:end,1:3); [c,d,q]= size(J1); fixedPoints = xlsread('data1.xls'); movingPoints = xlsread('data2.xls'); x0 = fixedPoints(:,1); y0 = fixedPoints(:,2); x1 = movingPoints(:,1); y1 = movingPoints(:,2); data = [x1,y1]; a = [1 1 1 1 1 1]; a1 = lsqcurvefit('fun',a,data,x0); a2 = lsqcurvefit('fun',a,data,y0); %多項式幾何校正 J2 = uint8(zeros(size(J))); for rgb = 1:3 for i = 1:m for j = 1:n 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 J2(i,j,rgb) = J1(round(fun(a1,[i,j])),round(fun(a2,[i,j])),rgb); end % J2(round(fun(a1,[i,j])),round(fun(a2,[i,j]))) = J(i,j); % end end end end [x,y] = size(J2); %根據數據游標取截取區域的左上方點 J3 = imcrop(I,[98 180 60*o/x 60*p/y]); J4 = imcrop(J2,[41 80 60 60]); [k,y,z] = size(J3); [h,t,e] = size(J4); %%重采樣 %雙線性內插法 u = h/k; v = t/y; itp = uint8(zeros(k,y,3)); for rgb = 1:3 for i = ceil(1/u):k-1 iu = floor(i*u); for j = ceil(1/v):y-1 jv = floor(j*v); itp(i,j,rgb) = (1-(i*u-iu))*(1-(j*v-jv))*J4(iu,jv,rgb)... +(1-(i*u-iu))*(j*v-jv)*J4(iu,jv+1,rgb)... +(i*u-iu)*(1-(j*v-jv))*J4(iu+1,jv,rgb)... +(i*u-iu)*(j*v-jv)*J4(iu+1,jv+1,rgb); end end end %去黑邊 for rgb = 1:3 for i = 1:3 for j = 1:175 itp(i,j,rgb) = J4(ceil(i*u),ceil(j*v),rgb); itp(145,j,rgb) = J4(43,ceil(j*v),rgb); end end for j = 1:2 for i = 1:145 itp(i,j,rgb) = J4(ceil(i*u),ceil(j*v),rgb); itp(i,175,rgb) = J4(ceil(i*u),61,rgb); end end end imwrite(J3,'Core1.bmp','bmp'); imwrite(itp,'Core2.bmp','bmp'); subplot(231),imshow(J),title('待配准圖像'); subplot(232),imshow(I),title('基准圖像'); subplot(233),imshow(J2),title('多項式幾何校正后'); subplot(234),imshow(J3),title('基准影像裁剪后'); subplot(235),imshow(itp),title('校正影像裁剪重采樣后'); % %基准圖重采樣 % kh = zuixiaogongbeishu(k,h); % yt = zuixiaogongbeishu(y,t); % u = h/kh;v = t/yt; % J5 = J3(1:k/kh:end,1:y/yt:end,1:3); % %配准圖 雙線性內插法重采樣 % itp = uint8(zeros(kh,yt,3)); % for rgb = 1:3 % for i = floor(1/u):kh-1 % iu = floor(i*u); % for j = floor(1/v):yt-1 % jv = floor(j*v); % itp(i,j,rgb) = (1-(i*u-iu))*(1-(j*v-jv))*J4(iu,jv,rgb)... % +(1-(i*u-iu))*(j*v-jv)*J4(iu,jv+1,rgb)... % +(i*u-iu)*(1-(j*v-jv))*J4(iu+1,jv,rgb)... % +(i*u-iu)*(j*v-jv)*J4(iu+1,jv+1,rgb); % end % end % end % %去黑邊 % for rgb = 1:3 % for i = 1:144 % for j = 1:10675 % itp(i,j,rgb) = J4(ceil(i*u),ceil(j*v),rgb); % end % end % for j = 1:175 % for i = 1:6235 % itp(i,j,rgb) = J4(ceil(i*u),ceil(j*v),rgb); % end % end % end % % itp1 = uint8(zeros(6193,10615,3)); % itp1(1:6193,1:10615,1:3) = itp(1:6193,1:10615,1:3); % imwrite(J5,'Crop1.bmp','bmp'); % J5 = imread('Crop1.bmp'); % imwrite(itp1,'Crop2.bmp','bmp'); % J6 = imread('Crop2.bmp'); % subplot(231),imshow(J),title('待配准圖像'); % subplot(232),imshow(I),title('基准圖像'); % subplot(233),imshow(J2),title('多項式幾何校正后'); % subplot(234),imshow(J5),title('基准影像裁剪重采樣后'); % subplot(235),imshow(J6),title('校正影像裁剪重采樣后'); % a = imread('基准.bmp'); % b = imread('重采樣后圖像.bmp'); % c = imcrop(a,[1,100,100,100]); % d = imcrop(b,[1,100,100,100]); % imwrite(c,'Core3.bmp','bmp'); % imwrite(d,'Core4.bmp','bmp');