How to find horizon line efficiently in a high-altitude photo?

Question

I am trying to detect the horizon in images taken from high altitude, so as to determine the orientation of the camera. I am also trying to have this run fast - ideally, I'd like to be able to process frames in real time (that is, a few frames per second) on a Raspberry Pi. The approach I've been taking so far is based on the fact that at high altitudes the sky is very dark, like so:

What I've tried is taking samples from all over the image and separating them into light and dark samples, and drawing a line between them. However, this places the horizon above its actual location, due to the fuzzyness of the amosphere:

And here's my code (in Javascript for ease of web demo):

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()

And the result (green line is the detected horizon, red and blue are estimated outside bounds):

How can I improve this? And is there a more efficient way to do it? The final program will probably be written in Python, or C if that's too slow.

Answer 1

Consider some basic channel mixing and thresholding, followed by vertical samples as @Spektre suggests. [Edited to change to 2*R-B instead of R+G-B following @Spektre's comment]

Here are some options on the channel mixing:

1. Original
2. Flat mono mix R+G+B
3. Red channel
4. 2*R - B
5. R + G - B

It looks like #4 is the clearest horizon (thanks @Spektre for making me check this more carefully), mixing the colours in a ratio [Red 2: Green 0: Blue -1], you get this monochrome image:

Setting blue negative means that the blue haze over the horizon is used to kill off the fuzziness there. This turns out to be more effective than just using red and/or green (try it with the Channel Mixer in the GIMP).

Then we can clarify further, if you like, by thresholding (although you could do this after sampling), here at 25% grey:

Using Spektre's approach of vertically sampling the image, just scan down until you see the value go over 25%. With 3 lines, you should gain 3 x,y pairs and thus reconstruct the curve knowing that it is a parabola.

For more robustness, take more than 3 samples and discard outliers.

Answer 2

I would do it like this:

1. convert to BW

2. scan Horizontal and Vertical lines from each side like this

   vertical lines scan from top

   

   The black line shows the line position. For selected one the green arrows shows direction of scan (down) and direction of color intensity visualization (right). The white curve is color intensity graph (so you actually see what is happening)

   1. select some grid step this is 64 points between lines
   2. create temp array int p[]; to store line

   3. pre process each line

      • p[0] is intensity of first pixel of line
      • p[1,...] is derivation by x for H and by y for V lines (just substract neighboring line pixels)

      blur p[1,...] few times to avoid noise problems (from both sides to avoid position shift).

   4. scan + integrate it back

      integration is just summing c(i)=p[ 0 ] + p[ 1 ] + ... + p[ i ]. If c is below treshold you are outside atmosphere so start scanning and if from start of line is this area you are scanning from the right side. Remember where you reach treshold A-point and continue scanning until you reach peak C-point (first negative derivation value or real peak ... local max value).

   5. compute B-point

      for ease you can do B = 0.5*(A+C) but if you want be precise then atmosphere intensity grows exponentially so scan derivations from A to C and determine exponential function from it. If derivation start differs from it you reached B-point so remember all B-points (for each line).

3. now you have set of B-points

   So delete all invalid B-points (you should have 2 per each line ... from start and from end) so area with bigger atmosphere is often the right one unless you have some dark seamless close object in view.

4. approximate some curve through remaining B-points

[Notes]

You cannot shift B-point position based on altitude because the visual thickness of atmosphere also depends on observer and light source (Sun) positions. Also you should filter remaining B-points because some stars in view could do a mess. But I think the curve approximation should be enough.

[Edit1] did some stuff around for fun

so I did it in BDS2006 C++ VCL ... so you have to change image access to your environment

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);
}

input image is pic0, output image is pic1 they are my classes so some members are:

• xs,ys size of image in pixels
• p[y][x].dd is pixel at (x,y) position as 32 bit integer type
• bmp is GDI bitmap
• rgb2i() converts all RGB pixels to intensity integer values <0-765> (r+g+b)

As you can see all horizon points are in pnt[pnts] array where:

• x,y is position of horizon point
• l is atmosphere thickness (exponential part)
• id is { 0,1,2,3 } which identify scan direction

Here is output image (works well even for rotated images)

This will not work for sun glow images unless you add some heavy duty filtering