Help for HORIZON

PURPOSE:
HORIZON detects strong and weak horizon borders for martian images.

The program currently supports only grayscale format. The input images
to the program must be a 1-band image. The program makes two assumptions
for input images. If input images don't fullfil the assumptions, the 
program will fail to detect horizon borders. The first assumption is 
that sky region is above ground region. The second assumption is that 
there is no object in sky region.

The output of this program will be a Comma-separated values (CSV) file. 
The length of the output CSV file should always be exactly the the same as 
input image's width. Each integer value in the CSV file represents the 
horizon border's LINE position in that specific SAMP. 
For example:
  324,325,-1,325,326,327, ... ...
This means that at column(SAMP) 0, the horizon border is at row(LINE) 324;
at column(SAMP) 1, the horizon border is at row(LINE) 325. At column(SAMP) 2, 
the horizon border value is -1, it means there is no horizon border detected 
or it is an outlier.

EXECUTION:

The program is quite sensitive to the parameters. Parameters' default values
are chosen experimentally. In the METHOD section, the algorithm used in this
program is introduced. Knowing exactly how the algorithm works will help
choose parameters' values. 

  1) HORIZON INP=input.img OUT=output.horizon
Use all default parameters' values, the program takes an input image, and then
generate a CSV file that contains the horizon borders.

  2) HORIZON INP=input.img OUT=output.horizon WH_ITER_MAX=0
By setting parameter WH_ITER_MAX to be 0, the program will turn off weak horizon
detection and refinement process. This is useful when salt-and-pepper noise was
not completely removed from preprocessing, and it is also useful when vignetting
effects present in the input image.

  3) HORIZON INP=input.img OUT=output.horizon FILTER_SIZE=9
By default the program uses a filter window size of 7 to blend salt-and-pepper
noise. This should work for most of cases. However, if a lot of noisy pixel 
appear in a given N x N space, then increase FILTER_SIZE should be helpful to
blend these noisy pixels to a degree which they don't affect the sky and ground
busyness levels calculation.

  4) HORIZON INP=input.img OUT=output.horizon SH_THRESH_MAX=700
The parameter SH_THRESH_MAX defines the upper searching bound for gradient 
magnitude image in preliminary horizon borders detection process. By decreasing
SH_THRESH_MAX value, the program will search more precisely for weaker horizon 
borders. This parameter is useful when strong horizon borders appear to be weaker 
than normal, or in another words, the sky pixels intensity is very close to ground
pixels intensity.

  5) HORIZON INP=input.img OUT=output.horizon BAR_THRESH=5 BAA_THRESH=10
The parameters BAR_THRESH and BAA_THRESH are used to identify and remove columns 
without sky region. If part of detected horizon borders are too close to the top
of image, then these horizon borders will be directly marked as false horizon,
and the value of these horizon borders will be set to -1. (Note: based on 
the algorithm in this program, the false horizon borders will always be very
close to the top of image, and will appear to be a zigzag pattern.) The parameter
BAR_THRESH and BAA_THRESH are useful when horizon borders occupy only part of the 
image.


METHOD:

The horizon detection algorithm consists of the following steps:

  1) Trim top, bottom, left, and right borders of the input image.
  2) Use a min, or max, or median filter with filter window size defined by 
     FILTER_SIZE parameter to smooth out salt-and-pepper noise in input images.
  3) Use Sobel operator to compute gradient magnitude image and gradient 
     direction image based on the smoothed image.
  4) Preliminary horizon borders detection by optimizing the difference of sky 
     and ground busyness levels. Horizon borders are a 1 dimensional array that 
     has the same width as the original input image. 
  5) Outliers and weak horizon borders detection and refinement.
  6) Remove the outliers that cannot be refined from pervious step. Remove the 
     columns without sky region. 

Depending on input image types, sometimes there are one or two lines appear on 
top, bottom, left, or right of input images that are significantly brighter or 
darker than their neighboring pixels. If these lines are not trimmed, it will 
result in very strong edges when calculating gradient magnitude images. These 
strong edges contribute to sky and ground busyness level calculation, and the 
accuracy of detecting horizon borders decrease significantly. 

Salt-and-pepper noise presents itself as sparsely occurring white and black pixels. 
Median filter should be effective to remove salt-and-pepper noise; max filter should 
be effective to remove pepper noise; min filter should be effective to remove salt 
noise. Filter window size is crucial when removing noise, and it is usually an odd 
number. In general, the larger the filter window size, the more blurry the resulted 
image is.

Sobel operator has one 3x3 horizontal kernel and one 3x3 vertial kernel. 

  horizontal kernel:  | -1  0  1 |       vertial kernel:  | -1  -2  -1 |
                      | -2  0  2 |                        |  0   0   0 |
                      | -1  0  1 |                        |  1   2   1 |

Convolve horizontal kernel with the original image will give us the edges horizontally 
denoted as Gx, and convolve vertial kernel with the original image will give us the 
edges vertically denoted as Gy. Based on the properties of the two kernels, Gx is 
defined as increasing in the "right" direction, and Gy is defined as increasing in the 
"down" direction. Combine Gx and Gy use the equation listed below will give us gradient 
magnitude image and gradient direction image.
  
  1. G = sqrt(Gx^2 + Gy^2)    gradient magnitude image.
  2. theta = arctan(Gy / Gx)  gradient direction image, depending on whether Gx and Gy 
                              are positive or negative numbers, theta belongs to different
                              quadrants. 

The preliminary horizon borders are calculated by maximizing the difference of sky and 
ground busyness levels. The optimization function Jn is defined below:
  
  Jn = Ground_busyness_level - Sky_busyness_level                               (1)

which
                           sum(ground pixels from gradient magnitude image)
  Ground_busyness_level = --------------------------------------------------    (2)
                          count(ground pixels from gradient magnitude image)

                        sum(sky pixels from gradient magnitude image)
  Sky_busyness_level = -----------------------------------------------          (3)
                       count(sky pixels from gradient magnitude image) 

The complete algorithm of finding maximum Jn value Jn_max is defined below:

       thresh_max - thresh_min
  n = ------------------------- + 1
       search_increment_factor

  for k = 1 to n
                       thresh_max - thresh_min
      t = thresh_min + ----------------------- x (k - 1)
                               n - 1

      for samp = 0 to image_width
          border_tmp[samp] = image_height

          for line = 0 to image_height  
              if gradient_magnitude_image(line, samp) >= t
                  border_tmp[samp] = line
                  break
              end if
          end for
      end for

      compute Ground_busyness_level using equation (2) listed above
      compute Sky_busyness_level using equation (3) listed above

      Jn = Ground_busyness_level - Sky_busyness_level

      if Jn > Jn_max
          Jn_max = Jn
          border = border_tmp
      end if
  end for

Outliers and weak horizon borders detection and refinement is an iterative process. 
The maximum iteration of detection and refinement is determined by either WH_ITER_MAX 
parameter or if there is no outlier and weak horizon detected. A horizon column is 
defined as outlier if its value is bigger or smaller than the average of its N neighbors 
by a threshold value. Determine whether a horizon column has weak horizon above it 
depends on two inputs from gradient magnitude image and gradient direction image. First, 
in gradient magnitude image, if there is a pixel value that is bigger or smaller than 
the average of all pixle values above the current detected position by a threshold value, 
then the program will further check gradient direction image. In gradient direction image, 
if a NxN space with the suspicious weak horizon pixel at the center contains more than 1/3 
pixels that are almost equal to the suspicious weak horizon pixel, then the suspicious weak 
horizon pixel will be defined as weak horizon pixel. The reason to use both gradient 
magnitude image and gradient direction image is to increase the accuracy of weak horizon 
detection. If only gradient magnitude image is used, the weak horizon detected can highly 
be a salt-and-pepper noise that was not successfully removed from preprocessing step. 
Once the horizon columns are marked as outlier or weak horizon existence, the optimizing
difference of sky and ground busyness levels algorithm introduced above will be performed 
on these columns individually.

In order to remove columns without sky region, the program must be able to first detect
them. The program analyzes the shape of the detected horizon borders to determine which 
part is true horizon and which part is columns without sky region. The false horizon borders
detected for the columns without sky region will satisfy the following two conditions:

  1) very close to the top of the image.
  2) the shape of the horizon borders appears to be in a "zigzag" pattern. 

For testing condition 1, the average of horizon border is calculated, and parameters BAA_THRESH
and BAR_THRESH threshold values are defined. For testing condition 2, the average of the absolute 
difference of horizon border is calculated, and parameters BAADA_THRESH and BAADR_THRESH threshold
valeus are defined. BAA_THRESH stands for border average accept threshold; BAR_THRESH stands for 
border average reject threshold; BAADA_THRESH stands for border average of absolute difference
accept threshold; BAADR_THRESH stands for border average of absolute difference reject threshold. 
The algorithm of determining true and false horizon is shown below:
  
  compute border_avg
  compute border_avg_absdiff

  if border_avg <= BAR_THRESH ||
     (border_avg <= BAA_THRESH && border_avg_absdiff >= BAADA_THRESH) ||
     border_avg_absdiff >= BAADR_THRESH
         marked as outliers or columns without sky region
  else if border_avg > BAA_THRESH && border_avg_absdiff < BAADA_THRESH
         marked as true horizon
  else 
         break the horizon borders in half, and recursively call this 
         algorithm until all the horizon borders are categorized.      
  end if

HISTORY:
2016-04-05  Steven Lu  Initial horizon by Steven Lu.

COGNIZANT PROGRAMMER: Steven Lu

PARAMETERS: