Homework 1

Histogram equalization

byte Figure 1. Digital elevation map of Mount St. Helens area.

Figure 1 shows a digital elevation map of the Mount St. Helens area in Washington. Start by reproducing this figure on your screen.

1. Change directory to hw1/dem
2. Run

   scons byte.view
   

3. Examine the file byte.rsf which refers to the byte (unsigned character) numbers which get displayed on the screen.
   1. Open byte.rsf with a text editor to check its contents.
   2. Run

      sfin byte.rsf
      

      to check the data size and format.
   3. Run

      sfattr < byte.rsf
      

      to check data attributes. What is the maximum and minimum value? What is the mean value? For an explanation of different attributes, run sfattr without input.

Each image has a certain distribution of values (a histogram). The histogram for the elevation map values is shown in Figure 2. When different values in a histogram are not uniformly distributed, the image can have a low contrast. One way of improving the contrast is histogram equalization.

hist Figure 2. Normalized histogram (solid line) and cumulative histogram (dashed line) of the digital elevation data.

Let $f(x,y)$ be the original image. The equalized image will be $F(x,y)$. Let $h(f)$ be the histogram (probability distribution) of the original image values. Let $H(F)$ be the histogram of the modified image. The mapping of probabilities suggests

$$H(F) dF = h(f) df$$ (1)

or, if we want the modified histogram to be uniform,

$$\frac{d F}{d f} = C h(f)$$ (2)

where C is a constant. Solving equation 2, we obtain the mapping

$$F(f) = f_0 + C \int\limits_{f_0}^f h(\phi) d\phi\;,$$ (3)

where $f_0$ is the minimum value of $f$.

The algorithm for histogram equalization consists of the following three steps:

1. Taking an input image $f(x,y)$, compute its histogram $h(f)$.
2. Compute the cumulative histogram $F(f)$ according to equation (3). Choose an appropriate $C$ so that the range of $F$ is the same as the range of $f$.
3. Map every pixel $f(x,y)$ to the corresponding $F(x,y)$.

Your task:

1. Among the Madagascar programs, find a program that implements histogram equalization. Hint: you may find the sfdoc utility useful.
2. Edit the SConstruct file to add histogram equalization. Create a new figure and compare it with Figure 1.
3. Check the effect of equalization by recomputing the histogram in Figure 2 with equalized data. Run

   scons hist.view
   

   to display the figure on your screen.

from rsf.proj import *

# Download data
txt = 'st-helens_after.txt' 
Fetch(txt,'data',
      server='https://raw.githubusercontent.com',
      top='agile-geoscience/notebooks/master')
Flow('data.asc',txt,'/usr/bin/tail -n +6')

# Convert to RSF format
Flow('data','data.asc',
     '''
     echo in=$SOURCE data_format=ascii_float
     label=Elevation unit=m
     n1=979  o1=557.805  d1=0.010030675 label1=X
     n2=1400 o2=5107.965 d2=0.010035740 label2=Y |
     dd form=native | 
     clip2 lower=0 | lapfill grad=y niter=200 
     ''')

# Convert to byte form
Flow('byte','data',
     '''
     dd form=native |
     byte pclip=100 allpos=y bias=2231
     ''')

# Display
Result('byte',
       '''
       grey yreverse=n screenratio=1 
       title="Digital Elevation Map" 
       ''')

# Histogram
Flow('hist','byte',
     '''
     dd type=float |
     histogram n1=256 o1=0 d1=1 |
     dd type=float
     ''')
Plot('hist',
     'graph label1=Value label2=Occurence title=Histogram')

# Cumulative histogram
Flow('cumu','hist','causint')

Result('hist','hist cumu',
       '''
       cat axis=2 ${SOURCES[1]} | scale axis=1 |
       graph label1=Value label2="Normalized Occurence"
       title=Histogram dash=0,1
       ''')

# ADD HISTOGRAM EQUALIZATION

End()

Homework 1