Backprojection algorithm assignment 8 ICG
=========================================

The original 2D image is placed in an NxN matrix ORIG. 
It's practical to set N to a power of 2, e.g. a square of 16x16. 

The following matrix displays a non-filled square 
with 0.0 indicating background and 1.0 indicating foreground:

0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 
0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 
0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0 
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 

This original is scanned from different angles PHI (in radians), for each
angle PHI delivering a 1D projection PROJ (array of length M, with M=2*N).

Scanning with angle PHI:
At the start of each projection, PROJ is cleared.
Each element (x,y)>0.0 from ORIG is added to an element of PROJ by:
  index in PROJ = center of PROJ + x*cos(PHI) + y*sin(PHI) 
  PROJ[index] += ORIG[x][y] 
So PROJ is filled for each PHI.
(step 3 of algorithm below)

A filter can be applied now, e.g. convolution- or fourier method. 
If no filter is applied, skip this part.
(step 4 of algorithm below)

Reconstruction with backprojection by adding PROJ to RECO for each angle PHI. 
(step 4 of algorithm below)

Different slices are generated by adding results of more projections with
angle PHI.
The first, slice 1, uses angles PHI=0.0 and PHI=pi/2.0 radians. 
The second, slice 2, adds the projections PHI=pi/4.0 and PHI=3*pi/4.0 to
the previous result in RECO.
Etcetera...

The filters (convolution- and fourier method) are described in section 
12.2.4 of the Handleiding ICG.
The algorithm for the fourier method is given at the end of this section:
  step 1 is projection with function fourn.c,
  step 2 is filtering, 
  step 3 is backprojection with fourn.c. 
See Appendix D of Handleiding ICG for a description of the FFT function.


Backprojection algorithm (from Handleiding ICG section 12.2.3.1)
================================================================

1. fill matrix ORIG(x,y) with original figure
   initialize matrix RECO(x,y) to 0.0
2. loop over angles PHI:
   -  initialize array PROJ(phi,t) to 0.0
   3. project ORIG(x,y) with angle PHI onto PROJ(phi,t)
      this may be done by nested for-loops over x and y:
        index in PROJ = center of PROJ + x*cos(PHI) + y*sin(PHI) 
        PROJ[index] += ORIG[x][y] 
   4. apply filter to PROJ(phi,t)
   5. backprojection, add PROJ to RECO by nested loop over x and y:
        x*cos(PHI) + y*sin(PHI) gives index in array PROJ
        RECO[x][y] += PROJ[index]
6. normalize RECO


Practice
========

suppose N = power of 2
ORIG  =  original (matrix NxN)
PROJ  =  projection (array M, e.g. M=2*N)
RECO  =  reconstruction (matrix MxM)
NORM  =  normalized reconstruction (matrix MxM)
PHI   =  scanning angle in radians

initialize
- fill ORIG
- clear RECO

first slice (for PHI=0.0  and  PHI=PI/2.0)
- for each PHI
  - make projection: scan ORIG with angle PHI -> PROJ
  - apply filter to PROJ (optional step)
  - backprojection: add PROJ with angle PHI to RECO
- after processing all projections
  - normalize RECO, result in NORM
  - display NORM

second slice (add two more angles  PHI=PI/4.0  and  PHI=PI*3.0/4.0)

next slices ...


March 1, 2011,
November 3, 2011,
Peter Klok.