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

我试图检测从高空拍摄的图像中的地平线，以确定相机的方向。 我也试图快速运行 – 理想情况下，我希望能够在Raspberry Pi上实时处理帧（即每秒几帧）。 到目前为止我采取的方法是基于这样一个事实：在高海拔地区，天空非常暗，如下：

来自太空的地球

我试过的是从整个图像中取样并将它们分成浅色和深色样本，并在它们之间画一条线。 然而，由于大气层的模糊性，这将地平线置于其实际位置之上：

这是我的代码（在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 ((iminmag : false) || (i>0 ? mag(vals[i-1][j][5])>minmag : false) || (jminmag : 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])0 ? mag(vals[i-1][j][11])<minmag : false) || (j<(vals[i].length-2) ? mag(vals[i][j+1][12])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编写。

考虑一些基本的通道混合和阈值处理，然后是@Spektre建议的垂直样本。 [根据@ Spektre的评论编辑改为2 * RB而不是R + GB]

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

多项选择

  1. 原版的
  2. 扁平单声道混合R + G + B.
  3. 红色通道
  4. 2 * R – B.
  5. R + G – B.

它看起来像＃4是最清晰的地平线（感谢@Spektre让我更仔细地检查它），以比例混合颜色[红色2：绿色0：蓝色-1]，你得到这个单色图像：

混合通道单声道

设置蓝色阴性意味着地平线上的蓝色阴霾用于消除那里的模糊性。 事实certificate这比使用红色和/或绿色更有效（尝试使用GIMP中的通道混合器）。

然后我们可以进一步澄清，如果你愿意，可以通过阈值处理（虽然你可以在采样后这样做），这里的灰度为25％：

25pc阈值

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

为了更加稳健，请采集3个以上的样本并丢弃exception值。

我会这样做：

  1. 转换为BW

  2. 像这样从每一侧扫描水平和垂直线

    垂直线从顶部扫描

    垂直线扫描

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

    1. 选择一些网格步骤这是行之间的64点
    2. create temp array int p[]; 存储线

    3. 预处理每一行

      • p[0]是线的第一像素的强度
      • p[1,...]是由x导出的H和由y导出的V行（只是减去相邻的行像素）

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

    4. 扫描+将其整合回来

      积分只是求和c(i)=p[ 0 ] + p[ 1 ] + ... + p[ i ] 。 如果c低于阈值，则您处于大气层外，因此开始扫描，如果从行开始是您从右侧扫描的区域。 记住达到阈值A-point并继续扫描，直到达到峰值C-point （第一个负值推导值或实际峰值…局部最大值）。

    5. 计算B-point

      为了方便你可以做到B = 0.5*(A+C)但如果你想要精确，那么大气强度呈指数增长，因此扫描从A到C推导并从中确定指数函数。 如果推导开始与它不同，则达到B-point因此请记住所有B-points （对于每一行）。

  3. 现在你有一套B-points

    因此，删除所有无效的B-points （你应该每行2个……从开始到结束）所以具有更大气氛的区域通常是正确的，除非你在视图中有一些黑暗的无缝近距离物体。

  4. 通过剩余的B-points近似一些曲线

[笔记]

您不能根据高度移动B-point位置，因为大气的视觉厚度也取决于观察者和光源（太阳）位置。 此外，您应该过滤剩余的B-points因为视野中的某些星星可能会弄得一团糟。 但我认为曲线近似应该足够了。

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

所以我在BDS2006 C ++ VCL中做到了 ……所以你必须改变对环境的图像访问 

 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>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=tr0) break; // non homogenous bump } pl=y; // - for end of exponential increasing intensity part for (i=c1-c0,y++;ytr1) { 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 } pl=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+tr0tr1) { pnt[pnts]=p; pnts++; n[p.id]++; } } // process H-lines for (y=da>>1;y>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=tr0) break; // non homogenous bump } pl=x; // - for end of eyponential increasing intensitx part for (i=c1-c0,x++;xtr1) { 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 } pl=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+tr0tr1) { pnt[pnts]=p; pnts++; n[p.id]++; } } pic1=pic0; // get the original image // chose id with max horizont points j=0; if (n[j]Canvas->Pen->Color=0x000000FF; // Red for (i=0;iCanvas->MoveTo(pnt[i].x,pnt[i].y); break; } for ( ;iCanvas->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 } ，用于识别扫描方向

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

正常

除非您添加一些重型过滤，否则这对太阳光图像不起作用

Interesting Posts

使用Bluez的C / C ++ BLE读/写示例

错误：没有这样的文件或目录 – C.

multithreadingudp服务器

为什么strcpy（）导致缓冲区溢出时没有分段错误？

将有符号整数除以2的幂

C中的无限递归

将预处理程序指令缩进为emacs中的C代码

SSE加载和添加

在int数组上使用memcmp严格符合吗？

头文件包含静态分析工具？