GEO4060 Project Assignment 2015

A simple library for image processing

• What is a PGM image?
• Convert to a binary image
• Edge detection
• Recontruct initial image from its edges
• Unit testing framework
• Performance analysis: timing your application
• The project assignement
• The final time of delivery
• Setup of your git repository

What is a PGM image?

1. A "magic number" for identifying the file type. A pgm image's magic number is the two characters "P5".
2.  Whitespace (blanks, TABs, CRs, LFs).
3. A Width, formatted as ASCII characters in decimal.
4. Whitespace.
5. A Height, again in ASCII decimal.
6. Whitespace.
7. The maximum gray value (Maxval), again in ASCII decimal. Must be less than 65536, and more than zero.
8. A single whitespace character (usually a newline).
9. A raster of Height rows, in order from top to bottom. Each row consists of Width gray values, in order from left to right. Each gray value is a number from 0 through Maxval, with 0 being black and Maxval being white. Each gray value is represented in pure binary by either 1 or 2 bytes. If the Maxval is less than 256, it is 1 byte. Otherwise, it is 2 bytes. The most significant byte is first.

Convert to a binary image

Edge detection

• If a pixel has the same value as its four surrounding neighbours (i.e. no edge) then the value of \(edge(i, j)\) will be zero.
• If the pixel is very different from its four neighbours (i.e. a possible edge) then \(edge(i, j)\) will be large in magnitude. If you are familiar with the discretization of partial differential equations, you will recognise that edge is the second derivative of \(pic\).


• We will always consider \(i\) and \(j\) to lie in the range \(1, 2, . . . Width\) and \(1, 2, . . . Height\) respectively.



• Pixels that lie outside this range (e.g. \(pic(i, 0)\) or \(pic(Width + 1, j)\)) are set to zero.



1. To implement this algorithm, you will first need to be able to read in a picture.
2. The result of your edge detection will be stored in a PGM image too. Therefore you also need to implement a routine to write a PGM image.
3. Note that the arrays have to be extended in each dimension to accommodate the boundary conditions. We can implement the boundary conditions by setting these halos to zero. You should take care that computation only take place on the interior of the pictures, e.g. loops should start at 1 and not 0.
4. Implement the edge detection method above to compute the edges in newpic based on the input stored in oldpic, and look at the output picture.
5. How can you validate your edge detection method?
6. Run the code on multiple images and check that it works correctly.




Recontruct initial image from its edges

It is actually possible to do the inverse operation, i.e. to reconstruct the initial picture from its edges. This requires multiple iterations and you should be able to reconstruct the original picture as follows:
• define a new data structure (the simplest one would be an array; but it may not be the best one...) called edge
• read the initial edges data file into bigpic
• zero the arrays oldpic, newpic and edge
• scatter bigpic to edge and set oldpic = edge
• repeat for many iterations:
• loop over \(i = 1, 2 . . . Width; j = 1, 2, . . . Height\)


\[ newpic(i,j) = {1 \over 4} \{ oldpic(i-1,j) + oldpic(i+1,j) + \\ \qquad \qquad \quad oldpic(i,j-1) + oldpic(i,j+1) - edge(i,j) \} \]


• end loop
• set oldpic = newpic
• end loop over iterations
• write out the final picture as before
For information, this is the Jacobi algorithm for solving the 2D Poisson equation \( \nabla ^2 pic = edge\).

Implement this reverse algorithm and explain how you can verify your code. For instance, can you easily create a "fake" image? You should start to see the original picture emerging after around 1000 iterations; after 100000 it should be practically indistinguishable from the original. Rather than running for a fixed number of iterations, it is better to stop when the output picture has converged to some tolerance, i.e. when it changes very little between each iteration. We can quantify this by computing \(\Delta\) defined by: \[ \Delta ^2 = {1 \over Width \times Height} \sum _{i=1;j=1} ^{Width;Height} (newpic(i,j) -oldpic(i,j))^2 \] and stopping when \(\Delta\) is less than some tolerance, e.g. 0.1. Use a larger picture (as supplied) and time your program to measure the time (CPU and elapsed).



Unit testing framework


Each tool (main program) and subroutine/function must be carefully tested and these tests may need to be reproduced if a developer changes the implementation of a routine/function. This is what we called unit testing.
For our project, you need to make sure that every function and main program can be tested we "fake images". A "Fake" image is an image that has been created "manually" or with a known algorithm and for which it is straightforward to get the result of the tested routine. When testing a routine, you need to compare the generated result with a reference that we have previously produced.



Performance analysis: timing your application


It is important to be able to time your application and have a detailed overview of the performance of your code. For this, you will implement a module containing routines to time each subroutine/functions of your tools. We wish to have a standardized outputs such as:

      Copyright (C) 2015, UIO

     Process CPU Time (s) | Process Elapsed Time (s) 
     =====================|==========================
               0.030      |           0.030               
     =====================|==========================

     Started on 04/03/2015 at 21:10:18 MET +01:00 from GMT
     Stopped on 04/03/2015 at 21:10:18 MET +01:00 from GMT
 


CPU time is the time for which the CPU was busy executing the task. It does not take into account the time spent in waiting for I/O (disk IO or network IO). Since I/O operations, such as reading files from disk, are performed by the OS, these operations may involve noticeable amount of time in waiting for I/O subsystems to complete their operations. This waiting time will be included in the elapsed time, but not CPU time. Hence CPU time is usually less than the elapsed time.
To measure the CPU time of your code, you can use CPU_TIME Fortran intrinsic subroutine:

   REAL(kind=8)  :: start_time, end_time
   ...
   call CPU_time(start_time)
   ... 
   call CPU_time(end_time)
   
   print*, 'CPU time = ', end_time - start_time
 

To measure the elasped time:

   INTEGER      :: ir
   real(kind=8) :: elapsed_time_start, elapsesd_time_end 
   ...
  CALL SYSTEM_CLOCK(COUNT=elapsed_time_start, count_rate=ir)
  ...
  
  CALL SYSTEM_CLOCK(COUNT=elapsed_time_end)
  print*, 'Elapsed time = ', (elapsed_time_end - elapsed_time_start ) / real(ir)
 

When running a code in "production" i.e. to get scientific results (not during the development phase), you would like to disable these timings as they may slow down your code.

For this you can for instance define a logical called USE_TIMING which call timing routines when true; for instance:

 LOGICAL :: USE_TIMING
 ...
 if (USE_TIMING) call timeAppInit()
 ...
 if (USE_TIMING) call timeAppEnd()
 


The value of USE_TIMING will be an input parameter of your tool (value chosen by users).



The project assignement


The purpose of the project assignment is to create a set of tools to manipulate and process PGM images. You will have to:
1. define an appropriate data structure (object) for PGM images, taking into account that your library could be extended in the future to handle various format. You must justify your choice in a separate document.
2. a routine to read PGM image
3. a routine to write PGM image
4. a routine to convert a PGM image into a binary (0 or 1) image with a threshold chosen by the user.
5. a routine for edge detection
6. routines to time your application: we should be able to measure and print information concerning the CPU and elapsed time
7. program(s) for users to run the implemented algorithms on a PGM image of their choice
8. a makefile or a simple Unix shell script for compiling and linking the application must be provided. You must provide a simple way to compile your application in debug mode i.e. with debug compiler options and ``standard'' mode.
9. a log file should be created by default with information on the run (which subroutines were called, time spent in each subroutine, etc.) and a silent mode (no log file) must be provided too.
10. your application and each tool much be tested (unit testing framework) and it must be easy to redo these tests if required. A test can be for instance reading a PGM file and write it back into another file and then compare (diff) the two files.
11. a documentation must be provided, explaining how to use the set of tools you have implemented and how to extend them (a cookbook for implementing new image processing algorithms, new data formats, new tools)
12. your source code must be stored in a repository (github, bitbucket, etc.) and I need to be able to access it.




The final time of delivery


The final time of delivery is Friday 22.05.2015 at 23:59



Setup of your git repository



On sverdrup in the directory of your choice:

mkdir ProjectFortran
cd ProjectFortran
git config --global user.name "Your Name"
git config --global user.email "you@some.domain"
git config --global color.ui "auto"
git config --global core.editor "your_editor"
git init

Display the status of the repository:

git status

Add specific files to the staging area:

git add filename_1 filename_2

Commit changes in the staging area to the repository's history: (Without -m and a message, this command runs a text editor.)

git commit -m "Some message explaining your changes"



If you need to create a Makefile to compile your code, start from one example provided in github/Fortran.
Once you have finished your project, you can either use github or bitbucket to store it and make it accessible (at least to me) or you can send me a tarball containing your project:

cd ..
tar cvf ProjectFortran.tar ProjectFortran
gzip ProjectFortran.tar

This command is to be done in the same directory you have created ProjectFortran directory. The resulting file is then called ProjectFortran.tar.gz. Please note that you do not need to send the data directory (containg PGM images).