\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
		TABLE OF CONTENTS
//////////////////////////////

%% COMPILATION

       ~~ compilation for CPU
            >>>> hardware depending optimization for the CPU 

       ~~ compilation for GPU

%% INSTALLATION

%% PROGRAM CALL

%% OPTIONS IN THE INPUT FILE

        ~~~~~~~~~~~  RECONSTRUCTION FROM RADIO
        ~~~~~~~~~~~  RECONSTRUCTION FROM SINOGRAM
	~~~~~~~~~~~  OTHER OPTIONS 



%% MORE ABOUT FILTERS
     
        ~~~~~~~~~~~ Modification to the Ramp Filter

  	~~~~~~~~~~~ A ring artifact removal filter by Dominique Bernard





%%%%%%%%%%%%%%%% COMPILATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 the distributed sources contains the subdirectory SOURCES. In this directory, the script setup.py must be executed, providing an optional argument to build either the GPU binaries or the CPU ones.


python setup.py install --install-lib ${PWD}/dummy --install-scripts ${PWD}/dummyscripts

python setup.py install --install-lib "%cd%"\dummy --install-scripts  "%cd%"\dummyscripts


~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ compilation of the CPU version ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

To compile you need to have python installed, and the Numeric module.
To compile the CPU binaries : 

python setup.py  build_ext --define USECPU

This command, on my installation, creates the module : build/lib.linux-x86_64-2.5/PyHST_c.so.
In your case, It will be created in a possibly different location depending on you installation.

This module must be copied or linked to the directory  PyHST_c_CPU/

The PyHST script look for the binary module first in the  PyHST_c_GPU/ directory and, by default, in the PyHST_c_CPU/ one.


 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> hardware depending optimization for the CPU  >>>>>>>>>>>>>>>>>>>>>>>>

The file c_hst_recon_1over.c uses an optimized instruction for float-to-int conversion. In the case of a 64bit processor the concerned lines are the following :

#include<emmintrin.h>
/* #include<ieeefp.h>
 */	


#define FLOAT_TO_INT(in,out)  \
     out=_mm_cvtss_si32(_mm_load_ss(&(in)));
  /*      __asm__ __volatile__ ("fistpl %0" : "=m" (out) : "t" (in) : "st") ; */


On a 32bit architecture, instead, you must include the ieeefp instead of  emmintrin, and use the asm line in place of the mm routine call.




~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ compilation of the GPU version ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


To compile you need to have python , the Numeric module and the CUDA environment installed.
The first step is compiling libgputomo.so, in the GPU directory. You find the commandline, to do that, in the instruction fileThe compilation proceeds as in the CPU case passing the USEGPU parameter to setup.py, and copying the obtained module to the appropriate directory.



%%%%%%%%%%%%%%%%%%%%%%%%% INSTALLATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 The program is launched through the PyHST wrapper script that you find in the WRAPPERS directory.
The installation consists then  in simply changing some variable values with the proper directory names so that the SOURCES directory, where you did the compilation, can be found.



%%%%%%%%%%%%%%%%%%% PROGRAM CALL    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Call the program in this way :

installation_directory/WRAPPERS/PyHST  input.par


Where input.par is the name of a file containing the parameters.
This file will be disccused in the following section 


%%%%%%%%%%%%%%%%%%%%%%%%%%%% OPTIONS IN THE INPUT FILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


~~~~~~~~~~~  RECONSTRUCTION FROM RADIO ~~~~~~~~~~~~~~~~~~~~~~~~~
 


We give here an example input file, and discuss it. Then we list all the possible variables that can be set and their default value.

Here the file, it begins after the following === line, explanation start after the second == line.
The input.par file can also be found in the Phantom directory.

X==================================================================================
! HST_SLAVE PARAMETER FILE

! Parameters defining the projection file series
FILE_PREFIX = phantom/proj
NUM_FIRST_IMAGE = 0 ! No. of first projection file
NUM_LAST_IMAGE = 1499 ! No. of last projection file
NUMBER_LENGTH_VARIES = NO
LENGTH_OF_NUMERICAL_PART = 4 ! No. of characters
FILE_POSTFIX = .edf
FILE_INTERVAL = 1 ! Interval between input files

! Parameters defining the projection file format
NUM_IMAGE_1 = 2048 ! Number of pixels horizontally
NUM_IMAGE_2 = 1    ! Number of pixels vertically
IMAGE_PIXEL_SIZE_1 = 10000.000000 ! Pixel size horizontally (microns)
IMAGE_PIXEL_SIZE_2 = 10000.000000 ! Pixel size vertically

! Parameters defining background treatment
SUBTRACT_BACKGROUND = NO ! No background subtraction
BACKGROUND_FILE = N.A.

! Parameters defining flat-field treatment
CORRECT_FLATFIELD = NO ! No flat-field correction
FLATFIELD_CHANGING = N.A.
FLATFIELD_FILE = N.A.
FF_PREFIX = N.A.
FF_NUM_FIRST_IMAGE = N.A.
FF_NUM_LAST_IMAGE = N.A.
FF_NUMBER_LENGTH_VARIES = N.A.
FF_LENGTH_OF_NUMERICAL_PART = N.A.
FF_POSTFIX = N.A.
FF_FILE_INTERVAL = N.A.

TAKE_LOGARITHM = NO ! Take log of projection values

! Parameters defining experiment
ANGLE_BETWEEN_PROJECTIONS = 0.120000 ! Increment angle in degrees
ROTATION_VERTICAL = YES
ROTATION_AXIS_POSITION = 1023.500000 ! Position in pixels

! Parameters defining reconstruction
OUTPUT_SINOGRAMS = NO ! Output sinograms to files or not

START_VOXEL_1 =      1 ! X-start of reconstruction volume
START_VOXEL_2 =      1 ! Y-start of reconstruction volume
START_VOXEL_3 =      1 ! Z-start of reconstruction volume
END_VOXEL_1 =   2048 ! X-end of reconstruction volume
END_VOXEL_2 =   2048 ! Y-end of reconstruction volume
END_VOXEL_3 =      1 ! Z-end of reconstruction volume
OVERSAMPLING_FACTOR = 4 ! 0 = Linear, 1 = Nearest pixel
ANGLE_OFFSET = 0.000000 ! Reconstruction rotation offset angle in degrees
SINOGRAM_MEGABYTES = 400 ! Maximum size of sinogram storage (megabytes)

! Parameters extra features PyHST
DO_CCD_FILTER = NO ! CCD filter (spikes)
CCD_FILTER = "CCD_Filter"
CCD_FILTER_PARA = {"threshold": 100.000000 }
DO_SINO_FILTER = NO ! Sinogram filter (rings)
SINO_FILTER = "SINO_Filter"
ar = Numeric.ones(1748,'f')
ar[0]=0.0
ar[2:18]=0.0
SINO_FILTER_PARA = {"FILTER": ar }
DO_AXIS_CORRECTION = NO ! Axis correction
AXIS_CORRECTION_FILE = 'filename'
OPTIONS= { 'padding':'E' , 'axis_to_the_center':'Y' } ! Padding and position axis

! Parameters defining output file / format
OUTPUT_FILE = axistocenter.vol

# Iterative regression with total variation minimization
BETA_TV = 700 # TV parameter
ITERATIVE_CORRECTIONS=200 # number of iterations
DETECTOR_DUTY_OVERSAMPLING = 2

X==================================================================================


FILE_PREFIX = phantom/proj

     this tells where the projection can be found. The projection name is obtained appending XXXX.edf to the prefix, where XXXX is the projection number

NUM_FIRST_IMAGE = 0 ! No. of first projection file

	the first projection is phantom/proj0000.edf

NUM_LAST_IMAGE = 1499 ! No. of last projection file

        the last projection is phantom/proj1499.edf 


NUMBER_LENGTH_VARIES = NO

	if this is set to YES, the projection names will be  : phantom/proj0.edf , phantom/proj1.edf ...  phantom/proj1499.edf 

LENGTH_OF_NUMERICAL_PART = 4 ! No. of characters

      This matters if NUMBER_LENGTH_VARIES is set to no. Increase it if you have more the 10**4 projections


FILE_POSTFIX = .edf

	self-explaining 

FILE_INTERVAL = 1 ! Interval between input files

       increase it to n if you want to skip n-1 projection over n 

NUM_IMAGE_1 = 2048 ! Number of pixels horizontally

          in the projection, should coincide with the Dim_1 field in the Edf file.

NUM_IMAGE_2 = 1 ! Number of pixels vertically

	should coincide with the Dim_2 field in the Edf file.

SUBTRACT_BACKGROUND = NO 

	if you want to subtract a background from the projections set this to YES and set the following variable

BACKGROUND_FILE = N.A.
      
  	if the above option is active, set it to a file name , for example dark.edf


CORRECT_FLATFIELD = NO ! No flat-field correction

	For this phantom example flat-field correction are not used ( it is a phantom). For experimental data however this option is always set to YES


FLATFIELD_CHANGING = N.A.
	if it is set to NO just one file is used for the flat-field, and its name is given by the following parameter.
	If it is set to YES the flatfield changes across the projection
	
FLATFIELD_FILE = N.A.

	Set it to your FF file that will be used for all projection

FF_PREFIX  = N.A.

	This parameters is used in case FLATFIELD_CHANGING=YES to compose the FF file names of the files

FF_NUM_FIRST_IMAGE = N.A
	
	The number XXXX such that file FF_PREFIXXXXX.FF_POSTFIX is the  FF corresponding  to the first projection
	
.
FF_NUM_LAST_IMAGE = N.A.

	The number XXXX corresponding to the last available  FF

FF_NUMBER_LENGTH_VARIES = N.A.

	same meaning as for projections

FF_LENGTH_OF_NUMERICAL_PART = N.A.

	same meaning as above.

FF_POSTFIX = N.A.

	for example .edf

FF_FILE_INTERVAL = N.A.

	The interval length, expressed in  number of projections, between a FF and the following. Interpolation is done for FF at the interior of the intervals.


TAKE_LOGARITHM = NO 

	! Take log of projection values
	The value NO is used in the phantom example because the
projections    contain already the logharithm of the absorption. For real
data this must be YES unless you have already preprocessed the projections.


ANGLE_BETWEEN_PROJECTIONS = 0.120000 ! Increment angle in degrees

       	this should be clear

ROTATION_VERTICAL = YES
	
	set to NO if it is horizontal

ROTATION_AXIS_POSITION = 1023.500000 

	! Position in pixels
	the center of a pixel in its middle. 


OUTPUT_SINOGRAMS = NO 
	! Output sinograms to files or not

        
     

START_VOXEL_1 =      1 ! X-start of reconstruction volume

	This value and the others define the  reconstructed region	

START_VOXEL_2 =      1 ! Y-start of reconstruction volume
START_VOXEL_3 =      1 ! Z-start of reconstruction volume
END_VOXEL_1 =   2048 ! X-end of reconstruction volume
END_VOXEL_2 =   2048 ! Y-end of reconstruction volume
END_VOXEL_3 =      1 ! Z-end of reconstruction volume




OVERSAMPLING_FACTOR = 4 ! 0 = Linear, 1 = Nearest pixel

	This option is used in the CPU version. To speed up the CPU, instead of intepolation, oversampled values are precalculated, and the nearest is used.


ANGLE_OFFSET = 0.000000 ! Reconstruction rotation offset angle in degrees

	this can rotate your volume.
	

SINOGRAM_MEGABYTES = 400 ! Maximum size of sinogram storage (megabytes)

	When the reconstructed region consists of several slices, the  data for a  given number of slices are read at once, trying to stay within the above memory limit.


DO_CCD_FILTER = NO ! CCD filter (spikes)

	This option can activate a filter which acts on the projections.


CCD_FILTER = "CCD_Filter"

	This give the name of the module which is loaded for projection filtering. 
	The filters are compiled with the setup_filters.py script.
	An example is the  CCD_Filter_Template.c file , which generates  the CCD_Filter.so    module



CCD_FILTER_PARA = {"threshold": 100.000000 }

	This variable is a python dictionary, containing parameters for the CDD_Filter. Every filter has its own set of keyword that it expects to find in this dictionary


==========================================================================================
======================  OLD KEYWORD FOR COMPATIBILITY WITH OLD PAR FILES =================
DO_SINO_FILTER = NO ! Sinogram filter (rings)	

SINO_FILTER = "SINO_Filter"
	This give the name of the method which is used for sinogram filtering.
	( ring artefact removal ) 

ar = Numeric.ones(1748,'f')
ar[0]=0.0
ar[2:18]=0.0
SINO_FILTER_PARA = {"FILTER": ar }
	This variable is a python dictionary, containing parameters for the SINO_Filter. Every filter has its own set of keyword that it expects to find in this dictionary
================================ END OF THE OLD KEYWORD DESCRIPTION ========================
===========================================================================================




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  RING_FILTER alias SINO_FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DO_RING_FILTER = YES/NO

~~~~~~~~~~  If DO_RING_FILTER = YES  
            there are several options

----- OPTION A)

RING_FILTER = "RING_Filter"
ar = Numeric.ones(1300,'f')
ar[0]=0.0
RING_FILTER_PARA = {"FILTER": ar}


>>>>>>>>  this is the same as SINO_FILTER = "SINO_Filter".
          The Filter array is used to filter the average projection
           ( which is then subtracted )
           

----- OPTION B)

RING_FILTER = "RING_Filter_THRESHOLDED"
ar = Numeric.ones(1300,'f')
ar[0]=0.0
RING_FILTER_PARA = {"FILTER": ar, threshold=0.15}


>>>>>>> as above but the filtering is done on an averaged projection
        from which one has beforehand removed the projection 
        of a thresholded reconstruction



DO_AXIS_CORRECTION = NO ! Axis correction
AXIS_CORRECTION_FILE = 'filename'

	this file is used if 	DO_AXIS_CORRECTION = YES.
   It can have one or two columns. The first columns is the horizontal correction. The second , if present, the vertical correction.



OPTIONS= { 'padding':'E' , 'axis_to_the_center':'Y' } ! Padding and position axis

	padding =='E' pads with extrema, '0' ( zero ) with zeros , defaults to 'E' 
	    in the filtered backprojection algorithm, the fourier transform is taken of the slice. To avoid border effect the data are padded.
	
        axis_to_the_center :  'Y' or 'N' 
 		If 'Y' is chosen the reconstructed region is centered over the rotation axis.



OUTPUT_FILE = axistocenter.vol
	
	the name of the output volume file.



~~~~~~~~~~~  RECONSTRUCTION FROM SINOGRAM ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


The par file for reconstructon from sinogram is sligthly different .
Here is an example


X===========================================================================================

RECONSTRUCT_FROM_SINOGRAMS= YES
NUM_FIRST_IMAGE = 0     ! No. of first projection file
NUM_LAST_IMAGE =  900   ! No. of last projection file


NUM_FIRST_SINOGRAM=0
NUM_LAST_SINOGRAM=0

SINOGRAM_PREFIX=avsino

SINOGRAM_POSTFIX= .edf
LENGTH_OF_NUMERICAL_PART =4
NUMBER_LENGTH_VARIES=NO

OUTPUT_FILE = abs.vol

OVERSAMPLING_FACTOR = 4 ! 0 = Linear, 1 = Nearest pixel


START_VOXEL_1 =      1 ! X-start of reconstruction volume
START_VOXEL_2 =      1 ! Y-start of reconstruction volume
END_VOXEL_1 =   1336 ! X-end of reconstruction volume
END_VOXEL_2 =   1336 ! Y-end of reconstruction volume

ROTATION_AXIS_POSITION = 526.0 ! Position in pixels

ANGLE_OFFSET=0.0

ANGLE_BETWEEN_PROJECTIONS= 0.3995560488346282


NO_SINOGRAM_FILTERING=0



IMAGE_PIXEL_SIZE_1 = 30.000000 ! Pixel size horizontally (microns)
IMAGE_PIXEL_SIZE_2 = 30.000000 ! Pixel size vertically

OPTIONS= { 'padding':'E' , 'axis_to_the_center':'Y' } ! Padding and position axis
X=============================================================================================








~~~~~~~~~~~  OTHER OPTIONS  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



    RENORMALIZE_FOR_TALBOT=0
    TALBOT_P2=0
    TALBOT_TD=0

Here with their default values. If you set RENORMALIZE_FOR_TALBOT=1, and the others according 
to your experiment, there will be an additional factor
     TALBOT_P2/ TALBOT_TD/2/PI *  pixelsize

	
    CUMSUM_INTEGRAL=0   # 0 no, 1 yes, 2 subract average, 3 subtract slope 
	
	



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  MORE ABOUT FILTERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Modification to the Ramp Filter ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

the default filter is a linear ramp, equal to 0 for the zero frequency and equal to 0.5 for the two highest frequencies.
Fourier filtering can be disactivated 
		
NO_SINOGRAM_FILTERING = 1

It cas also be modified by a multiplicative function :

def  FOURIER_FILTER(self, x):
	return 0.5+(1-0.5)*math.cos( math.pi*  x * 2 )


in this function x takes the value 0 for the zero frequency and 0.5 ( -0.5) for the highest ones.
The returned value is used to multiply the ramp filter.
The parameters given in the input file are accessible as members of the self object, if you need them.

This gives you some freedom. However in the general case you want to take care also of the ramp, for example if you dont want the filter go to zero at zero.

The most general case is obtained with the following option :



FOURIER_FILTER_HAS_RAMP=YES


and then defining the FOURIER_FILTER with arguments (self,  i, dim_fft).
The index i being the number of the fourier component, and dim_fft the numbers of data to fft.
The index i is 0 for the lowest frequency, 1 for the highest. For the other  values  the quotient (i/2) is proportional to the absolute value of the frequency.

The following is an example of implementation of a band limited filter

( http://www.mathematica-journal.com/issue/v6i2/article/murrell/murrell.pdf)


htransfdictio={}
def FOURIER_FILTER(self,  i, dim_fft):
   return   self.Htransform( dim_fft)[i]
def Htransform(self, N):
    if self. htransfdictio.has_key(N):
        return self.htransfdictio[N]
    h=Numeric.zeros(N,"d")
    for i in range(N/2):
        if (i%2):
            h[i]=-1.0/i/i/Numeric.pi/Numeric.pi
            h[N-i]=-1.0/i/i/Numeric.pi/Numeric.pi
    h[0]=1.0/4
    import FFT
    h=fft(h).real
    self.htransfdictio[N]=h
    return h




 This is an example of a step filter which integrates the phase derivative from a talbot radiography

FOURIER_FILTER_HAS_RAMP=YES
def FOURIER_FILTER(self,  i, dim_fft):
    if i==0:
            return 0.0
    elif i==dim_fft/2:
            return 0.0
    elif i >0  :
            return 1.0/2/3.1415
    else :
            return -1.0/2/3.1415



~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ A ring artifact removal filter by Dominique Bernard ~~~~~~~~~~~


Description in progress.
Look
 Filtre_DBRINGS.doc                  

# -------------- Total variation regression ----------------------

This section deals with an iterative reconstruction algorithm (hence,
different from the filtered back projection), that solves the
minimization problem

x = argmin 1/2||H x - y||^2  + beta_tv TV(x)

where:

* x is the image to be reconstructed
* H is the projection matrix
* y is the sinogram
* TV(x) is the total variation semi-norm of x, that is, the l1 norm of
  its gradient
* beta_tv is a parameter controlling the relative importance of the two
  terms in the minimization

This algorithm tends to reconstruct piecewise-constant images. It is
therefore suitable for images with a limited number of phases, such as
many images in materials science. It is not suited, however, for images
with a very irregular texture, or for images with smooth large-scale
gradients.

BETA_TV = 700

The parameter controlling the relative importance of the total variation
penalization. The greater BETA_TV, the more piecewise-constant will be
the image.

ITERATIVE_CORRECTIONS=200 # number of iterations

The number of iterations of the algorithm.


The optional parameter ITERATIVE_CORRECTIONS_NOPREC is by default set to None.
In the default caseall iteration are done without preconditionning the fidelity term.
If  ITERATIVE_CORRECTIONS_NOPREC is assigned an integer number N, the first iterations are done
without preconditionning, and after the Nth the preconditionning is done.

DO_LOCAL_TOMO=0
by default,
if set to one the reconstructed slice beyond the always-in-view limits 
follows the fidelity-term gradient without regularisation.

LOCAL_TOMO_ROI_RADIUS=1000000
by default huge.
This defines a roi which includes  by default the always-in-view limits.
If one changes the default setting it to smaller value , this parameter becomes  active :
   the local-tomo regularisation is applied in the region 
given by the intersection of the ROI circle and the always-in-view region. 

DETECTOR_DUTY_OVERSAMPLING = 2

If this parameter is used, one projection acquired over a small angular
cone (for continuous acquisition) is considered to be the sum of several
projections corresponding to directions regularly spaced inside the small
angular cone. The number of different directions is given by
DETECTOR_DUTY_OVERSAMPLING. Use this parameter if the number of
projections is small compared to the detector size, in the case of
continuous acquisition.


#----------------- Non Local Means Regularisation
The default of DENOISING_TYPE is 1 ( tv regression with fast gradient and fista).

A different type of denoising is used if DENOISING_TYPE=2
The minimal set of keys is

DENOISING_TYPE=2
CALM_ZONE_LEN = 31
NLM_NOISE_GEOMETRIC_RATIO = 2.0
NLM_NOISE_INITIAL_FACTOR  =1.4


The number of ITERATIVE_CORRECTIONS is the number of subsequent nlm_denoising steps.
CALM_ZONE_LEN is the diameter of the area where a typical variance (noise) is calculated.
This area is calculated  in the initial filtered back projected reconstruction 
at a point such that  the noise is the least.
At the first iteration the nlm averaging is applied with a noise*NLM_NOISE_INITIAL_FACTOR
as parameter for the nlm method.
Subsequently the typical noise is divided by NLM_NOISE_GEOMETRIC_RATIO at each iteration.