MATLAB桶形变换

%% 镜头桶形失真校正（短焦镜头）

img_origin1 = imread（‘E:\360data\重要数据\桌面\2.2.jpg‘）;
img_origin = rgb2gray（img_origin1）;

k1 = -1.4 * 10^（-6）;   % 形变参数，根据实际情况调整
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（‘原图用的地方‘）;

%% 校正以后的图上重复的地方（多个点对应原图上的一个点）
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（‘校正后的图上重复的地方‘）;
%% 找到校正以后图上最接近正确的点（有很多重复）
% 因为有很多点重复，所以找到最能一一对应的那个店

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