How To Find Horizon Line Efficiently In A High-altitude Photo?

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 t

Solution 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.

Solution 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-linesfor (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 leftfor (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 areafor (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 partfor (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;  // peakif (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 areafor (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 partfor (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;  // peakif (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-linesfor (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 leftfor (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 areafor (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 partfor (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;  // peakif (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 areafor (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 partfor (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;  // peakif (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;    // Redfor (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