如何在高海拔照片中有效地找到地平线?
我试图检测从高空拍摄的图像中的地平线,以确定相机的方向。 我也试图快速运行 – 理想情况下,我希望能够在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]
以下是频道混音的一些选项:
- 原版的
- 扁平单声道混合R + G + B.
- 红色通道
- 2 * R – B.
- R + G – B.
它看起来像#4是最清晰的地平线(感谢@Spektre让我更仔细地检查它),以比例混合颜色[红色2:绿色0:蓝色-1],你得到这个单色图像:
设置蓝色阴性意味着地平线上的蓝色阴霾用于消除那里的模糊性。 事实certificate这比使用红色和/或绿色更有效(尝试使用GIMP中的通道混合器)。
然后我们可以进一步澄清,如果你愿意,可以通过阈值处理(虽然你可以在采样后这样做),这里的灰度为25%:
使用Spektre垂直采样图像的方法,只需向下扫描,直到您看到值超过25%。 有3条线,你应该获得3 x,y对,从而重建曲线,知道它是抛物线。
为了更加稳健,请采集3个以上的样本并丢弃exception值。
我会这样做:
-
转换为BW
-
像这样从每一侧扫描水平和垂直线
垂直线从顶部扫描
黑线表示线的位置。 对于选定的一个, 绿色箭头显示扫描方向(向下)和颜色强度可视化方向(右)。 白色曲线是颜色强度图(所以你实际看到发生了什么)
- 选择一些网格步骤这是行之间的64点
- create temp array
int p[];
存储线 -
预处理每一行
-
p[0]
是线的第一像素的强度 -
p[1,...]
是由x
导出的H
和由y
导出的V
行(只是减去相邻的行像素)
模糊
p[1,...]
几次以避免噪音问题(从两侧避免位置偏移)。 -
-
扫描+将其整合回来
积分只是求和
c(i)=p[ 0 ] + p[ 1 ] + ... + p[ i ]
。 如果c
低于阈值,则您处于大气层外,因此开始扫描,如果从行开始是您从右侧扫描的区域。 记住达到阈值A-point
并继续扫描,直到达到峰值C-point
(第一个负值推导值或实际峰值…局部最大值)。 -
计算
B-point
为了方便你可以做到
B = 0.5*(A+C)
但如果你想要精确,那么大气强度呈指数增长,因此扫描从A
到C
推导并从中确定指数函数。 如果推导开始与它不同,则达到B-point
因此请记住所有B-points
(对于每一行)。
-
现在你有一套
B-points
因此,删除所有无效的
B-points
(你应该每行2个……从开始到结束)所以具有更大气氛的区域通常是正确的,除非你在视图中有一些黑暗的无缝近距离物体。 -
通过剩余的
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 }
,用于识别扫描方向
这是输出图像(即使对于旋转图像也能很好地工作)
除非您添加一些重型过滤,否则这对太阳光图像不起作用