基于经度坐标校正鱼眼图像

最近开始鱼眼图像校正方面的研究，在这个过程中阅读博主元气少女缘结神的相关博客让我受益匪浅，在此对她表示感谢,另外所有代码在Github。

提取有效区域

在研究中仅仅考虑圆形的鱼眼图像，其他形状，如长方形，不在目前的研究范围。在校正鱼眼图像之前需要找到有效的图像区域，即圆形区域。借鉴张伟等人的《鱼眼图像校正算法研究》，在其3.5节改进的算法中提出了兼顾精度和效率的提取方法，大意是分别从图象的上下左右进行扫描，当扫描到大于灰度阈值并且向里一个像素仍然大于灰度阈值的时候认定扫描到有效鱼眼区域的边界。当然为了提取效率，论文中提到的类似将椭圆转换为圆形的部分不进行考虑。

参考博客 matlab/python+opencv提取圆形鱼眼图片的有效区域以及快速扫描算法提取鱼眼图像有效区域，略微改进其中提到的程序就可以形成下面的提取程序，该程序经过验证具有很快的提取效率，并且适用于大多数鱼眼图像。运行的时候只要通过imread读入图象到img，然后阈值T设置为40即可。

% 参考http://blog.csdn.net/dengxf01/article/details/53374014
% 参考http://blog.csdn.net/wd1603926823/article/details/45672741
% 张伟等 《鱼眼图像校正算法研究》
function [img_valid, R] = imageEffectiveAreaInterception(img, T)
% input
% img: rgb image
% T: gray threshold
% output
% img_valid: effective image area
% R: effective image area radius

% rgb to gray
% gray = 0.299r + 0.587g + 0.114b
img_gray = rgb2gray(img);

% image size
[m, n, k] = size(img_gray);

% scan from top to bottom
for i = 1:m
    flag = 0;
    for j = 1:n
        if(img_gray(i, j) >= T)
           if(img_gray(i+1, j) >= T)
               top = i;
               flag = 1;
               break;
           end
        end
    end
    if flag == 1; break; end
end

% scan from bottom to top
for i = m:-1:top
    flag = 0;
    for j = 1:n
        if(img_gray(i, j) >= T)
           if(img_gray(i-1, j) >= T)
               bottom = i;
               flag = 1;
               break;
           end
        end
    end
    if flag == 1; break; end
end

% scan from left to right
for j = 1:n
    flag=0;    
    for i = top:bottom
        if(img_gray(i, j) >= T)
           if(img_gray(i, j+1) >= T)
               left = j;
               flag = 1;               
               break;
           end
        end
    end
    if flag == 1; break; end
end

% scan from right to left
for j = n:-1:left
    flag = 0;    
    for i = top:bottom
        if(img_gray(i, j) >= T)
           if(img_gray(i, j-1) >= T)
               right = j;
               flag = 1;               
               break;
           end
        end
    end
    if flag == 1; break; end
end

% effictive area radius
R = max((right - left)/2, (bottom - top)/2);

% create effictive area 
img_valid = imcrop(img, [left, top, 2*R, 2*R]);
end

原图如下

效果图如下

基于经度坐标校正鱼眼图像

获取到鱼眼图像的有效区域后就可以进行下一步的矫正工作，张伟等人的***《鱼眼图像校正算法研究》***提到了基于经度坐标校正鱼眼图像，这是一种非常简单的校正方法。在论文中也提到了相关原理，但是阅读中因为作者没有提到太多细节，开始理解起来的时候比较费力，主要是没有提到计算过程中使用的坐标系以及相关证明。

通过阅读汪丹等人的***《一种不断重定位圆心的鱼眼图像校正方法》***才明白其中的原理。所谓经度坐标校正就是把鱼眼图像中像素的横坐标变换到原来的位置，而纵坐标不变，通过这样的变换会把圆形的鱼眼区域变换成正方形。下图（来源于《一种不断重定位圆心的鱼眼图像校正方法》）明确阐述了其中的原理，鱼眼图像中心坐标为O(x0,y0)，半径为R，A(x,y)是在圆形畸变区域边上的一点，B(u,v)是其经过校正后的坐标，其中B点肯定会位于校正后正方形的右边上。

通过作者列出下面的公式，可以方便地通过圆内的像素转换为正方形中的像素。

但是我有一个疑问就是作者是通过计算圆边上的点得来的，如果是圆内的弧会映射成正方形内的竖线吗？带着这样的疑问我阅读了黄有度等人的《一种鱼眼图象到透视投影图象的变换模型》，里面的定理1明确证明了空间中的直线经过球面投影后会变成投影面上的椭圆弧。这里可以自己证明，把点A当成椭圆弧上的点，椭圆长轴为2*R，短轴长为椭圆弧与X轴两个交点间的距离，列出椭圆公式带入坐标点即可证明。

得到通过圆内像素坐标计算校正后像素的公式后就可以写代码了，首先参考了博客鱼眼图像经度坐标校正（2D），里面的程序可以顺利运行，但是有一个很大的问题就是当图片很大的时候运行速度会很慢，因为程序中用了多层for循环，导致最终的复杂度为O(n^2)。将其中的代码略加修改之后如下，去掉了后面的插值程序，并且改成了我之前提到的提取有效区域函数imageEffectiveAreaInterception，这跟博客中使用的kuaisusaomiao是同样的作用，

function C = fisheye_longitude_correction_other(A, T)
tic;
[A,R]=imageEffectiveAreaInterception(A,T);
[m,n,k]=size(A);
C=[];
x=n/2;
y=m/2;
for u=1:m
    for v=1:n
        i=u;
        j=round(sqrt(R^2-(y-u)^2)*(v-x)/R+x);
        if(R^2-(y-u)^2<0)
            continue;
        end
        C(u,v,1)=A(i,j,1);
        C(u,v,2)=A(i,j,2);
        C(u,v,3)=A(i,j,3);
    end
end
C=uint8(C);
toc;
end

下面是该程序的运行结果，一共88秒多。其中图片A就是本篇博客最开始用到的测试图片，这张图片大小为3744*5616，是一张比较大的图片，另外将阈值T设置为40，

Elapsed time is 88.339909 seconds.

为了改进这个程序，我将多层for循环改成了矩阵化的运算，其中参考了博客for循环矩阵化和论坛两个向量分别保存行列坐标，提取相应元素，程序如下,

% 参考http://blog.csdn.net/wd1603926823/article/details/45841841
% 张伟等 《鱼眼图像校正算法研究》
% 黄有度等 《一种鱼眼图象到透视投影图象的变换模型》
% 汪丹等 《一种不断重定位圆心的鱼眼图像校正方法》
% for循环矩阵化 http://blog.csdn.net/myj0513/article/details/6871482
% 两个向量分别保存行列坐标，提取相应元素 http://www.ilovematlab.cn/thread-63242-1-1.html
function result = fisheye_longitude_correction(img, T)
% input
% image: rgb image
% T: gray threshold
% output
% result: corrected image

tic; % record run time
[img_valid, R] = imageEffectiveAreaInterception(img, T); % effective image area 
[m, n, k] = size(img_valid);
result = zeros(m, n, 3); % preallocation
x0 = n/2; y0 = m/2; % image center

% core code
[v, u] = meshgrid(1:m, 1:n);
tmp = round(sqrt(abs(R^2-(u-y0).^2)).*(v-x0)/R+x0);

% reshape
for i = 1:k
   img_gray = img_valid(:, :, i);
   coordinate = u(:) + (tmp(:) - 1)*m;
   channel = img_gray(coordinate);
   channel_reshape = reshape(channel, m, n);
   result(:, :, i) = channel_reshape;
end

% transform
result=uint8(result);
toc;
end

其中最核心的部分是将两层for循环中转换坐标的部分进行了矩阵化，复杂度变成了O(n)，通过meshgrid生成坐标矩阵，然后对坐标矩阵进行变换，这种方法设置跟之前完全相同，下面是程序的运行时间，共约1.5秒，大大缩短了运行时间。

Elapsed time is 1.496467 seconds.

当然最后了两种方法的结果都是一样的，如下图,