%%?鏡頭桶形失真校正（短焦鏡頭）

img_origin1?=?imread（‘E:\360data\重要數(shù)據(jù)\桌面\2.2.jpg‘）;
img_origin?=?rgb2gray（img_origin1）;

k1?=?-1.4?*?10^（-6）;???%?形變參數(shù)，根據(jù)實際情況調(diào)整
k2?=?-1.4?*?10^（-6）;

cx?=?116;
cy?=?26;

img_size?=?size（?img_origin?）;
img_undist?=?zeros（?img_size（1）+cy*2img_size（2）+cx*2?）;
new_size?=?size（img_undist）;
img_undist?=?uint8（?img_undist?）;

%%?校正
tic
counter?=?zeros（img_size）;
counter?=?uint8（counter）;

error?=?zeros（new_size）;
error?=?double（error）;

yy?=?zeros（img_size）;
xx?=?zeros（img_size）;

for?l1?=?1:（img_size（1）+cy*2）??%?垂直方向
????
????y?=?l1?-?img_size（1）/2-cy;
????
????for?l2?=?1:（img_size（2）+cx*2）??%?水平方向
????????
????????x?=?l2?-?img_size（2）/2-cx;
????????
????????x1?=?round（?x?*?（?1?+?k1?*?x?*?x?+?k1?*?y?*?y?）?）;
????????x1_err?=??x?*?（?1?+?k1?*?x?*?x?+?k1?*?y?*?y?）?-?x1;
????????y1?=?round（?y?*?（?1?+?k2?*?x?*?x?+?k2?*?y?*?y?）?）;
????????y1_err?=??y?*?（?1?+?k2?*?x?*?x?+?k2?*?y?*?y?）?-?y1;
????????err?=?x1_err?*?x1_err?+?y1_err?*?y1_err;
????????y1?=?y1?+?img_size（1）/2;
????????x1?=?x1?+?img_size（2）/2;
????????
????????
????????if?y1>0?&&?y1<（img_size（1）+1）?&&?x1>0?&&?x1<（img_size（2）+1）
????????????yy（y1x1）?=?l1;
????????????xx（y1x1）?=?l2;
????????????error（l1l1）?=?err;????
????????????counter（y1x1）=counter（y1x1）+1;
????????????img_undist（l1l2）?=?img_origin（y1?x1）;
????????end
????end
end
toc
%%?原圖用的地方的像素
img_lost?=?img_origin;
for?i=1:img_size（1）
????for?j=1:img_size（2）
????????if?xx（ij）>0?&&?yy（ij）>0?
????????????img_lost（ij）?=?img_origin（ij）;
????????else
????????????img_lost（ij）?=?255;
????????end
????end
end
figure（1）imshow（img_lost）;title（‘原圖用的地方‘）;

%%?校正以后的圖上重復(fù)的地方（多個點對應(yīng)原圖上的一個點）
img_lost2?=?img_undist;
for?i=1:img_size（1）
????for?j=1:img_size（2）
????????if?xx（ij）>0?&&?yy（ij）>0
????????????img_lost2（yy（ij）xx（ij））?=?255;
????????end
????end
end
figure（2）imshow（img_lost2）;title（‘校正后的圖上重復(fù)的地方‘）;
%%?找到校正以后圖上最接近正確的點（有很多重復(fù)）
%?因為有很多點重復(fù)，所以找到最能一一對應(yīng)的那個店

error_last?=?error;

for?l1?=?1:（img_size（1）+cy*2）??%?垂直方向
????
????y?=?l1?-?img_size（1）/2-cy;
????
????for?l2?=?1:（img_size（2）+cx*2）??%?水平方向
????????
????????x?=?l2?-?img_size（2）/2-cx;
????????
????????x1?=?round（?x?*?（?1?+?k1?*?x?*?x?+?k1?*?y?*?y?）?）;
???????