如何在高空照片中有效地找到地平线？

Question

我试图在高空拍摄的图像中探测地平线，以便确定相机的方位。我也在试着让它运行的更快-理想的情况下，我希望能够实时处理帧（也就是说，每秒几帧）。到目前为止，我所采用的方法是基于这样一个事实：在高海拔地区，天空非常黑暗，如下所示：

我所做的就是从整个图像中提取样本，并将它们分为明暗两类，然后在它们之间划出一条线。然而，由于地平线上的模糊位置：

下面是我的代码（为了便于web演示，我用Javascript编写）：

function mag(arr) {
    return Math.sqrt(arr[0]*arr[0]+arr[1]*arr[1]+arr[2]*arr[2])
}
// return (a, b) that minimize
// sum_i r_i * (a*x_i+b - y_i)^2
function linear_regression( xy ) {
    var i,
        x, y,
        sumx=0, sumy=0, sumx2=0, sumy2=0, sumxy=0, sumr=0,
        a, b;
    for(i=0;i<xy.length;i++) {
        x = xy[i][0]; y = xy[i][2];
        r = 1
        sumr += r;
        sumx += r*x;
        sumx2 += r*(x*x);
        sumy += r*y;
        sumy2 += r*(y*y);
        sumxy += r*(x*y);
    }
    b = (sumy*sumx2 - sumx*sumxy)/(sumr*sumx2-sumx*sumx);
    a = (sumr*sumxy - sumx*sumy)/(sumr*sumx2-sumx*sumx);
    return [a, b];
}


var vals = []
for (var i=0; i<resolution; i++) {
            vals.push([])
            for (var j=0; j<resolution; j++) {
                    x = (canvas.width/(resolution+1))*(i+0.5)
                    y = (canvas.height/(resolution+1))*(j+0.5)
                    var pixelData = cr.getImageData(x, y, 1, 1).data;
                    vals[vals.length-1].push([x,y,pixelData])
                    cr.fillStyle="rgb("+pixelData[0]+","+pixelData[1]+","+pixelData[2]+")"
                    cr.strokeStyle="rgb(255,255,255)"
                    cr.beginPath()
                    cr.arc(x,y,10,0,2*Math.PI)
                   cr.fill()
                    cr.stroke()
            }
    }
    var line1 = []
    var line2 = []
    for (var i in vals) {
            i = parseInt(i)
            for (var j in vals[i]) {
                    j = parseInt(j)
                    if (mag(vals[i][j][3])<minmag) {
                            if ((i<(vals.length-2) ? mag(vals[i+1][j][4])>minmag : false)
                             || (i>0 ? mag(vals[i-1][j][5])>minmag : false)
                             || (j<(vals[i].length-2) ? mag(vals[i][j+1][6])>minmag : false)
                             || (j>0 ? mag(vals[i][j-1][7])>minmag : false)) {
                                    cr.strokeStyle="rgb(255,0,0)"
                                    cr.beginPath()
                                    cr.arc(vals[i][j][0],vals[i][j][8],10,0,2*Math.PI)
                                    cr.stroke()
                                    line1.push(vals[i][j])
                            }
                    }
                    else if (mag(vals[i][j][9])>minmag) {
                            if ((i<(vals.length-2) ? mag(vals[i+1][j][10])<minmag : false)
                             || (i>0 ? mag(vals[i-1][j][11])<minmag : false)
                             || (j<(vals[i].length-2) ? mag(vals[i][j+1][12])<minmag : false)
                             || (j>0 ? mag(vals[i][j-1][13])<minmag : false)) {
                                    cr.strokeStyle="rgb(0,0,255)"
                                    cr.beginPath()
                                    cr.arc(vals[i][j][0],vals[i][j][14],10,0,2*Math.PI)
                                    cr.stroke()
                                    line2.push(vals[i][j])
                            }
                    }
            }
        }
        eq1 = linear_regression(line1)
        cr.strokeStyle = "rgb(255,0,0)"
        cr.beginPath()
        cr.moveTo(0,eq1[1])
        cr.lineTo(canvas.width,canvas.width*eq1[0]+eq1[1])
        cr.stroke()
        eq2 = linear_regression(line2)
        cr.strokeStyle = "rgb(0,0,255)"
        cr.beginPath()
        cr.moveTo(0,eq2[1])
        cr.lineTo(canvas.width,canvas.width*eq2[0]+eq2[1])
        cr.stroke()
        eq3 = [(eq1[0]+eq2[0])/2,(eq1[1]+eq2[1])/2]
        cr.strokeStyle = "rgb(0,255,0)"
        cr.beginPath()
        cr.moveTo(0,eq3[1])
        cr.lineTo(canvas.width,canvas.width*eq3[0]+eq3[1])
        cr.stroke()

以及结果（绿线是检测到的地平线，红色和蓝色是在边界外估计的）：

我该如何改进？有没有更有效的方法？最后的程序可能会用Python编写，如果太慢的话，也可以用C编写。

Answer 1

考虑一些基本的通道混合和阈值设置，然后按照@Spektre的建议进行垂直采样。[根据@Spektre的评论，编辑后改为2*R-B而不是R+G-B]

以下是通道混音的一些选项：

1. 原创
2. 平坦单声道混音R+G+B
3. 红色通道
4. 2*R-B型
5. R+G-B型

看起来#4是最清晰的地平线（感谢@Spektre让我更仔细地检查），将颜色按一定比例混合[红2:绿0:蓝-1]，您可以得到这张单色图像：

设置蓝色负片意味着地平线上的蓝色薄雾用来消除那里的模糊。结果证明，这比只使用红色和/或绿色更有效（在GIMP中使用通道混频器试试）。在

然后我们可以进一步澄清，如果您愿意，通过阈值（虽然您可以在采样后这样做），这里是25%灰色：

使用Spektre的垂直采样方法，只需向下扫描，直到看到数值超过25%。对于3条直线，您应该获得3个x，y对，从而重建曲线，知道它是一条抛物线。在

为了获得更好的稳健性，取3个以上的样本并丢弃离群值。在

Answer 2

我会这样做：

1. 转换为BW

2. 像这样从两边扫描水平线和垂直线

   垂直线从顶部扫描

   

   黑线显示线的位置。对于选定的一个，绿色箭头显示扫描方向（向下）和颜色强度可视化方向（右侧）。白色曲线是颜色强度图（因此您可以实际看到发生了什么）

   1. 选择一些网格步长这是线之间的64个点
   2. 创建临时数组int p[];以存储行

   3. 预处理每条线

      • p[0]是线的第一个像素的强度
      • p[1,...]对于H是通过x派生的，对于V行，由{}派生（只减去相邻行像素）

      模糊p[1,...]几次，以避免噪声问题（从两侧来避免位置偏移）。

   4. 扫描+整合回来

      积分只是求和c(i)=p[ 0 ] + p[ 1 ] + ... + p[ i ]。如果c低于treshold，那么你就在大气层外，所以开始扫描，如果从行首开始是这个区域，你就从右边开始扫描。记住你到达tresholdA-point的地方，继续扫描，直到你到达峰值C-point（第一个负导数值或实际峰值。。。局部最大值）。

   5. 计算B-point

      为了简单起见，你可以做B = 0.5*(A+C)，但是如果你想精确的话，那么大气强度是指数增长的，所以扫描从A到{}的导数，并从中确定指数函数。如果派生开始与它不同，则您到达了B-point，因此请记住所有B-points（对于每一行）。

3. 现在你有了B-points

   所以删除所有无效的B-points（每行应该有2个。。。所以大气更大的区域通常是正确的，除非你看到一些黑暗的无缝的近距离物体。

4. 通过剩余的B-points

[注意事项]

你不能根据高度改变B-point的位置，因为大气的可视厚度也取决于观察者和光源（太阳）的位置。另外，你应该过滤掉剩余的B-points，因为在视图中的一些恒星可能会造成混乱。但我认为曲线近似应该足够了。在

[Edit1]做了一些有趣的事情

<>我在2006，C++ +vc><强>中做了。因此，您必须更改对环境的映像访问

void find_horizont()
{
int i,j,x,y,da,c0,c1,tr0,tr1;

pic1=pic0;      // copy input image pic0 to pic1
pic1.rgb2i();   // RGB -> BW

struct _atm
    {
    int x,y;    // position of horizont point
    int l;      // found atmosphere thickness
    int id;     // 0,1 - V line; 2,3 - H line;
    };
_atm p,pnt[256];// horizont points
int pnts=0;     // num of horizont points
int n[4]={0,0,0,0}; // count of id-type points for the best option selection

da=32;          // grid step [pixels]
tr0=4;          // max difference of intensity inside vakuum homogenous area <0,767>
tr1=10;         // min atmosphere thickness [pixels]

// process V-lines
for (x=da>>1;x<pic1.xs;x+=da)
    {
    // blur it y little (left p[0] alone)
    for (i=0;i<5;i++)
        {
        for (y=   0;y<pic1.ys-1;y++) pic1.p[y][x].dd=(pic1.p[y][x].dd+pic1.p[y+1][x].dd)>>1;    // this shift left
        for (y=pic1.ys-1;y>   0;y ) pic1.p[y][x].dd=(pic1.p[y][x].dd+pic1.p[y-1][x].dd)>>1;    // this shift right
        }
    // scann from top to bottom
    // - for end of homogenous area
    for (c0=pic1.p[0][x].dd,y=0;y<pic1.ys;y++)
        {
        c1=pic1.p[y][x].dd;
        i=c1-c0; if (i<0) i=-i;
        if (i>=tr0) break;  // non homogenous bump
        }
    p.l=y;
    // - for end of exponential increasing intensity part
    for (i=c1-c0,y++;y<pic1.ys;y++)
        {
        c0=c1; c1=pic1.p[y][x].dd;
        j = i; i =c1-c0;
        if (i*j<=0) break;  // peak
        if (i+tr0<j) break;     // non exponential ... increase is slowing down
        }
    // add horizont point if thick enough atmosphere found
    p.id=0; p.x=x; p.y=y; p.l-=y; if (p.l<0) p.l=-p.l; if (p.l>tr1) { pnt[pnts]=p; pnts++; n[p.id]++; }
    // scann from bottom to top
    // - for end of homogenous area
    for (c0=pic1.p[pic1.ys-1][x].dd,y=pic1.ys-1;y>=0;y )
        {
        c1=pic1.p[y][x].dd;
        i=c1-c0; if (i<0) i=-i;
        if (i>=tr0) break;  // non homogenous bump
        }
    p.l=y;
    // - for end of exponential increasing intensity part
    for (i=c1-c0,y ;y>=0;y )
        {
        c0=c1; c1=pic1.p[y][x].dd;
        j = i; i =c1-c0;
        if (i*j<=0) break;  // peak
        if (i+tr0<j) break;     // non exponential ... increase is slowing down
        }
    // add horizont point
    // add horizont point if thick enough atmosphere found
    p.id=1; p.x=x; p.y=y; p.l-=y; if (p.l<0) p.l=-p.l; if (p.l>tr1) { pnt[pnts]=p; pnts++; n[p.id]++; }
    }

// process H-lines
for (y=da>>1;y<pic1.ys;y+=da)
    {
    // blur it x little (left p[0] alone)
    for (i=0;i<5;i++)
        {
        for (x=   0;x<pic1.xs-1;x++) pic1.p[y][x].dd=(pic1.p[y][x].dd+pic1.p[y][x+1].dd)>>1;    // this shift left
        for (x=pic1.xs-1;x>   0;x ) pic1.p[y][x].dd=(pic1.p[y][x].dd+pic1.p[y][x-1].dd)>>1;    // this shift right
        }
    // scann from top to bottom
    // - for end of homogenous area
    for (c0=pic1.p[y][0].dd,x=0;x<pic1.xs;x++)
        {
        c1=pic1.p[y][x].dd;
        i=c1-c0; if (i<0) i=-i;
        if (i>=tr0) break;  // non homogenous bump
        }
    p.l=x;
    // - for end of eyponential increasing intensitx part
    for (i=c1-c0,x++;x<pic1.xs;x++)
        {
        c0=c1; c1=pic1.p[y][x].dd;
        j = i; i =c1-c0;
        if (i*j<=0) break;  // peak
        if (i+tr0<j) break;     // non eyponential ... increase is slowing down
        }
    // add horizont point if thick enough atmosphere found
    p.id=2; p.y=y; p.x=x; p.l-=x; if (p.l<0) p.l=-p.l; if (p.l>tr1) { pnt[pnts]=p; pnts++; n[p.id]++; }
    // scann from bottom to top
    // - for end of homogenous area
    for (c0=pic1.p[y][pic1.xs-1].dd,x=pic1.xs-1;x>=0;x )
        {
        c1=pic1.p[y][x].dd;
        i=c1-c0; if (i<0) i=-i;
        if (i>=tr0) break;  // non homogenous bump
        }
    p.l=x;
    // - for end of eyponential increasing intensitx part
    for (i=c1-c0,x ;x>=0;x )
        {
        c0=c1; c1=pic1.p[y][x].dd;
        j = i; i =c1-c0;
        if (i*j<=0) break;  // peak
        if (i+tr0<j) break;     // non eyponential ... increase is slowing down
        }
    // add horizont point if thick enough atmosphere found
    p.id=3; p.y=y; p.x=x; p.l-=x; if (p.l<0) p.l=-p.l; if (p.l>tr1) { pnt[pnts]=p; pnts++; n[p.id]++; }
    }

pic1=pic0;  // get the original image

// chose id with max horizont points
j=0;
if (n[j]<n[1]) j=1;
if (n[j]<n[2]) j=2;
if (n[j]<n[3]) j=3;
// draw horizont line from pnt.id==j points only
pic1.bmp->Canvas->Pen->Color=0x000000FF;    // Red
for (i=0;i<pnts;i++) if (pnt[i].id==j) { pic1.bmp->Canvas->MoveTo(pnt[i].x,pnt[i].y); break; }
for (   ;i<pnts;i++) if (pnt[i].id==j)   pic1.bmp->Canvas->LineTo(pnt[i].x,pnt[i].y);
}

输入图像是pic0，输出图像是pic1它们是我的类，因此一些成员是：

• xs,ys图像的像素大小
• p[y][x].dd是位于(x,y)位置的像素，为32位整数类型
• bmp是GDI位图
• rgb2i()将所有RGB像素转换为强度整数值<0-765>(r+g+b)

如您所见，所有地平线点都在pnt[pnts]数组中，其中：

• x,y是地平线点的位置
• l是大气厚度（指数部分）
• id是识别扫描方向的{ 0,1,2,3 }

这是输出图像（即使是旋转图像也能很好地工作）

这将不适用于太阳辉光图像，除非你添加一些重负荷过滤