/tomo/pyhst

##  The PyHST program   is Copyright (C) 2002-2008 of the 
##  European Synchrotron Radiation Facility (ESRF).

##  You may use, distribute and copy the PyMCA XRF Toolkit under the terms of
##  GNU General Public License version 3 or (at your option) any later version.

import resource
import time
import string
import math
import re
import sys
import os

from copy import copy
from copy import deepcopy
from socket import htonl
from mmap import mmap, MAP_PRIVATE
from os.path import basename
from logger import logger

import fast_EdfFile
from FastEdf import extended_fread, extended_fread_2d
from PIL import Image
import numpy as Numeric

import build_config

if hasattr(build_config, 'DEFINES') and ("-DPYHST_IO_BENCHMARK" in build_config.DEFINES):
    reconstruct = False
else:
    reconstruct = True

if hasattr(build_config, 'DEFINES') and ("-DPYHST_RECON_BENCHMARK" in build_config.DEFINES):
    simulate_data = True
else:
    simulate_data = False
    
if hasattr(build_config, 'DEFINES') and ("-DHW_USE_PARALLEL_IO" in build_config.DEFINES):
    parallel_io = True
else:
    parallel_io = False

# This code should be after loading of PyHST, i.e. g_thread_init
# should be called before loading vipsCC module, otherwise
# everything will crash! (On operations with timers)
if hasattr(build_config, 'DISABLE_VIPS') and build_config.DISABLE_VIPS:
    use_vips = False
else:
    try:
	from vipsCC import *
	use_vips = True
    except:
	use_vips = False
    
edf_timer = 0			# Measures time needed for data loading
init_timer = 0			# Initalization and cleanup stages of C-code
c_timer = 0 			# Reconstruction and write-out stages of C-code
mem_timer = 0			# Transposition stage of C-code
pre_timer = 0			# Python-preprocessing

all_timer = time.time()		# Measures overall execution time

start_time = time.time()
import PyHST_c_CPU as PyHST_c
init_timer += time.time() - start_time


def ModulesForFF(  ip ,  INTERVALS):
    if type(INTERVALS)==type(1):
        return ip%INTERVALS, INTERVALS
    else:
        sum=0
        INTERVALS=Numeric.array(INTERVALS)
        if(INTERVALS[0] == 0):
            INTERVALS = INTERVALS[1:]-INTERVALS[:-1]
        for i in range(len(INTERVALS)):
            if sum+INTERVALS[i] > ip:
                return ip-sum, INTERVALS[i]
            elif sum+INTERVALS[i] == ip:
                return 0, INTERVALS[i]
                
            sum=sum+INTERVALS[i]
        raise Exception, " sum of INTERVALS for FF is smaller than projection index "

EdfFile=fast_EdfFile

DEG2RAD=math.pi/180.0

def treat_par_file(s):
    """
      pythonify the old format (fortran)
      input file given as a string s
    """
    logger.info("PYTHONIFYING INPUT FILE")
    s=string.replace(s,"!","#")
    lines=string.split(s,"\n")
    new_lines=""
    for line in lines:
        if( string.find(line,"FILE ")>=0 or string.find(line,"FILE=")>=0):
            line=add_quotes(line)
        elif( string.find(line,"PREFIX")>=0):
             line=add_quotes(line)
        elif( string.find(line,"POSTFIX")>=0):
             line=add_quotes(line)
        new_lines=new_lines+"\n"+line
    # try to execute the input file
    logger.info("CHECKING INPUTFILE FOR EXECUTION")
    try:
       NO =0
       YES=1
       new_lines= string.replace(new_lines, "\"N.A.\"", "\"N_A_\"")
       new_lines= string.replace(new_lines, "N.A.", "\"N_A_\"")
       exec(new_lines)
    except:
       print new_lines
       import traceback
       traceback.print_exception(sys.exc_type, sys.exc_value, sys.exc_traceback)

       raise Exception, " something went wrong interpreting input file "
    logger.info("INPUT FILE IS GOOD FOR ME")
    return new_lines


def add_quotes(s):
    pos=string.find(s,"=")
    if(pos==-1): return s
    new_s=s[:pos+1]
    pos=pos+1
    toggle=0
    while(pos<len(s) and s[pos]==" "):
        pos=pos+1
    if(pos==len(s)):
        return s
    new_s=new_s+"\""
    while(pos<len(s) and s[pos]!=" "):
        new_s=new_s+s[pos]
        pos=pos+1
    new_s=new_s+"\""
    return new_s


class Singleton(type):
    instance = None
    def __init__(self, name, bases, dict):
	super(Singleton, self).__init__(name, bases, dict)

    def __call__(self, *args, **Kw):
	if self.instance is None:
	    self.instance = super(Singleton, self).__call__(*args, **Kw)
	return self.instance

class EDFBuffer:
    __metaclass__ = Singleton
    def __init__(self):
	self.item = None
    def get(self, S):
	if ((self.item is None) or (S[0] * S[1] > self.item.size)):
	    self.item = Numeric.zeros(S, 'f')
	elif (S[0] * S[1] < self.item.size):
	    self.item.resize(S)

	return self.item    

class EDFBuffers:
    def __init__(self):
	self.item = None
    def get(self, S):
	if ((self.item is None) or (S[0] * S[1] > self.item.size)):
	    self.item = Numeric.zeros(S, 'f')
	elif (S[0] * S[1] < self.item.size):
	    self.item.resize(S)

	return self.item    


class SIMULATEReader:
    def __init__(self, Parameters, multipass):
	self.preload = True
    def open_file(self, name):
	return None
    def read_data(self, ipro, ima, P, S):
	if not(hasattr(self,"data")) or (self.S != S) or (self.P != P):
	    self.P = P
	    self.S = S
	    self.data = Numeric.zeros((S[0], S[1]), dtype=Numeric.float32)
	    
	return self.data

class EDFReader:
    def __init__(self, Parameters, multipass):
        self.edf_buffer = EDFBuffer()
	self.preload = True
    def open_file(self, name):
	return EdfFile.EdfFile(name, 1)
    def read_data(self, ipro, ima, P, S):
      sizex, sizey = ima.Images[0].Dim1, ima.Images[0].Dim2
      image = self.edf_buffer.get(S)
      ima.File.seek((((P[0]*sizex)+P[1])*4)+ima.Images[0].DataPosition,0)
      extended_fread_2d(image, 4 * S[1], S[0], 4 * sizex, ima.File,  string.upper(ima.SysByteOrder)!=string.upper(ima.Images[0].ByteOrder))
      return image


class EDFMultiReader:
    def __init__(self, Parameters, multipass):
        self.edf_buffers = []
	self.preload = True
    def open_file(self, name):
	return EdfFile.EdfFile(name, 1)
    def read_data(self, ipro, ima, P, S):
      sizex, sizey = ima.Images[0].Dim1, ima.Images[0].Dim2

      if len(self.edf_buffers) <= ipro:
        self.edf_buffers.append(EDFBuffers())
      image = self.edf_buffers[ipro].get(S)

      ima.File.seek((((P[0]*sizex)+P[1])*4)+ima.Images[0].DataPosition,0)
      extended_fread_2d(image, 4 * S[1], S[0], 4 * sizex, ima.File,  string.upper(ima.SysByteOrder)!=string.upper(ima.Images[0].ByteOrder))
      return image


class IMAGEReader:    
    def __init__(self, Parameters, multipass):
	self.lendian  = not (htonl(1) == 1)
	self.preload = False
    def open_file(self, name):
	return Image.open(name)
    def read_data(self, ipro, ima, P, S):
#	print ima.tile
#	ima.tile = [ima.tile[0][0], [P[1], P[0], P[1] + S[1], P[0] + S[0]], ima.tile[0][2], ima.tile[0][3]];
	image = Numeric.asarray(ima.crop([P[1], P[0], P[1] + S[1], P[0] + S[0]]));
	image.flags.writeable = True
	if self.lendian:
	    return image.byteswap(False)
        else:
    	    return image

class VIPSReader:    
    def __init__(self, Parameters, multipass):
	self.lendian  = not (htonl(1) == 1)
	self.preload = False
	if multipass:
	    if re.match('\\.(v|ppm|pbm|pgm)$', Parameters.FILE_POSTFIX):
		self.mmap = False
	    else:
		logger.warning("Using MMAPed I/O")
		self.mmap = True
	else:
	    self.mmap = False
    def open_file(self, name):
	if self.mmap:
	    info = VImage.VImage(name)
	    self.f = open(name, "rb")
	    self.data = mmap(self.f.fileno(), os.path.getsize(name), MAP_PRIVATE)
	    vim = VImage.VImage_frombuffer(self.data, info.Xsize(), info.Ysize(), info.Bands(), info.BandFmt())
	else:
	    vim = VImage.VImage (name)
	return vim
    def read_data(self, ipro, ima, P, S):
	area = ima.extract_area(P[1], P[0], S[1], S[0]);
	#pim = Image.frombuffer (self.mode, [S[1],S[0]], area.tobuffer(), 'raw', self.mode, 0, 1)
	#image = Numeric.asarray(pim)
	#image.flags.writeable = True
	image = Numeric.frombuffer(area.tobuffer(), dtype='f')
	image = image.reshape([S[0], S[1]])
	if self.mmap and self.lendian:
	    return image.byteswap(False)
	else:
	    return image.copy()


def open_image_file(reader, n, Parameters, rescaleN=1):
    name=Parameters.FILE_PREFIX

    if(rescaleN):
        n *= Parameters.FILE_INTERVAL

    if( Parameters.NUMBER_LENGTH_VARIES):
        name = name + ("%d"%n) + Parameters.FILE_POSTFIX
    else:
        number = "%d"%n
	if (len(number) < Parameters.LENGTH_OF_NUMERICAL_PART ):
	    number = ("0" *  (Parameters.LENGTH_OF_NUMERICAL_PART - len(number))) + number
        else:
    	    number = number[:Parameters.LENGTH_OF_NUMERICAL_PART]

        name= name + number + Parameters.FILE_POSTFIX

#    return EdfFile.EdfFile(name, 1)
    return reader.open_file(name)

def read_image_data(reader, ipro, ima, vertical, pos_edf, size_edf, BINNING=None):
   if BINNING is None:
       P = copy(pos_edf)
       S = copy(size_edf)
   else:
       P = [ pos_edf[0]*BINNING , pos_edf[1]*BINNING ] 
       S = [ size_edf[0]*BINNING, size_edf[1]*BINNING ] 
       
#   P.reverse()
#   S.reverse()

#   image=ima.GetData(0,Pos=P, Size=S, DataType="FloatValue")
   image = reader.read_data(ipro, ima, P, S)

   if BINNING is not None:
       image = Numeric.reshape(image, [size_edf[0], BINNING , size_edf[1], BINNING])
       image = Numeric.swapaxes(image, 1, 2)
       image = Numeric.sum(Numeric.sum(image, axis=-1), axis=-1)

   if(vertical==0):
       image = Numeric.transpose(image)
  
   return image


def extract_edf_N(prefix , n ,  lvaries, nlength, postfix, vertical , pos_edf, size_edf , Parameters, rescaleN=1, BINNING=None):
    name=prefix
    if(rescaleN):
        n=n*Parameters.FILE_INTERVAL
    if(lvaries):
        name=name+("%d"%n)+postfix
    else:
        number = "%d"%n
	if(len(number)< nlength ):
	        number = ("0" *  (nlength-len(number)))+ number
        else:
           number=number[:nlength]
        name=name+number+postfix
    return extract_edf(name , vertical    , pos_edf, size_edf  , BINNING  )
                  
def extract_edf(name , vertical  , pos_edf, size_edf   , BINNING=None   ):
   ima=EdfFile.EdfFile(name, 1)
   
   if BINNING is None:
       P=copy(pos_edf)
       S=copy(size_edf)
   else:
       P = [ pos_edf[0]*BINNING , pos_edf[1]*BINNING    ] 
       S = [ size_edf[0]*BINNING, size_edf[1]*BINNING   ] 
       
   P.reverse()
   S.reverse()
   image=ima.GetData(0,Pos=P, Size=S, DataType="FloatValue")

   if BINNING is not None:
       image=Numeric.reshape(image,[    size_edf[0], BINNING , size_edf[1], BINNING ])
       image = Numeric.swapaxes(image,1,2)
       image=Numeric.sum(Numeric.sum(image, axis=-1), axis=-1)

   if(vertical==0):
       image=Numeric.transpose(image)
       
   return image
       

def read_one_slice (name, sizex, sizey, size_pixel, sliceno, totalslices):
     File=open(name, "r+b")
     File.seek(0,0)
     Pos=[0,0]
     Size=sizex,sizey
     Size=list(Size)
     if Size[0]==0:Size[0]=sizex-Pos[0]
     if Size[1]==0:Size[1]=sizey-Pos[1]
     Data=Numeric.array([], "f")
     #for y in range(0,sizey):
     offset=size_pixel*sliceno*sizex*sizey
     #File.seek((((y*sizex)+Pos[0])*size_pixel)+offset,0)
     File.seek(offset,0)
     Data=Numeric.fromstring(File.read(sizey*sizex*size_pixel), "f")
	 #Data[y]=line
         #Data=Numeric.concatenate((Data,line))
     Data=Numeric.reshape(Data, (sizey,sizex))
     Data=Numeric.transpose(Data)
     File.close()

     return Data


## def  OverlappingLogic (pos_edf,  size_edf, axis, last_slice  , marge  ):
##         pos_edf_= copy(pos_edf)
##         if(pos_edf_[axis]>=marge):
## 	  pos_edf_[axis]=pos_edf_[axis]-marge
##         else:
##           pos_edf_[axis]=0
##         dsize=pos_edf[axis]-pos_edf_[axis]
        
##         if(size_edf[axis]+pos_edf[axis]+marge<=last_slice):
##            dsize=dsize+marge
##         else:
##            dsize=dsize+  last_slice - (size_edf[axis]+pos_edf[axis])

##         size_edf_=copy(size_edf)
##         size_edf_[axis]= size_edf[axis]+dsize

## 	return pos_edf_, size_edf_


def  OverlappingLogic (pos_edf,  size_edf, axis, last_slice  , marge  ):
        pos_edf_= copy(pos_edf)
        if(pos_edf_[axis]>=marge):
	  pos_edf_[axis]=pos_edf_[axis]-marge
        else:
          pos_edf_[axis]=0
        dsize=pos_edf[axis]-pos_edf_[axis]
        
        if(size_edf[axis]-1+pos_edf[axis]+marge<=(last_slice)):
           dsize=dsize+marge
        else:
           dsize=dsize+  last_slice - ((size_edf[axis]-1)+pos_edf[axis])

        size_edf_=copy(size_edf)
        size_edf_[axis]= size_edf[axis]+dsize
	return pos_edf_, size_edf_

def Filter_and_Trim(item, pos_edf_, size_edf_, pos_edf, size_edf, CCD_FILTER, CCD_AXIS_LONGITUDINAL_CORRECTION,  axis_correctionsL,  Xcorr, Ycorr, CCD_FILTER_PARA):
    if CCD_FILTER:
        CCD_FILTER(item, pos_edf_, size_edf_, pos_edf, size_edf,   CCD_FILTER_PARA     )
        
    if CCD_AXIS_LONGITUDINAL_CORRECTION:
        CCD_AXIS_LONGITUDINAL_CORRECTION(item, pos_edf_, size_edf_, pos_edf, size_edf,  axis_correctionsL,  Xcorr, Ycorr)
    
    begin0 = (pos_edf[0]- pos_edf_[0])
    begin1 = (pos_edf[1]- pos_edf_[1])  
    end0   = begin0 + size_edf[0]
    end1   = begin1 + size_edf[1]
    return item[  begin0  : end0 , begin1 : end1      ]
   

def simplereadxml(filename,what):
     e = open(filename,'r')
     value = None
     for line in e.readlines():
         posbeg = line.find('<'+what+'>')
         if posbeg > -1:
             posend = line.find('</'+what+'>')
             value = line[posbeg+len(what)+2:posend].strip()
             break
     e.close()
     return value

 
def  WriteInfo(name ,
               dim1,
               dim2,
               dim3
              ):
    f=open(name+".info","w")
    f.write("! PyHST_SLAVE VOLUME INFO FILE\n")
    f.write("NUM_X =  %d\n" % dim1)
    f.write("NUM_Y =  %d\n" % dim2)
    f.write("NUM_Z =  %d\n" % dim3)
    import sys
    if sys.byteorder=="big":
        f.write("BYTEORDER = HIGHBYTEFIRST\n")
    else:
        f.write("BYTEORDER = LOWBYTEFIRST\n")

    f.close()
    try:
      import os
      os.system("chmod og+r "+  name+".info"      )
      os.system("chmod og+r "+  name   )
    except:
      pass



def  WriteInfoXML(name ,Parameters):
    import sys
    if sys.byteorder=="big":
        BYTEORDER = "HIGHBYTEFIRST"
    else:
        BYTEORDER = "LOWBYTEFIRST"


    f=open(name+".xml","w")
    f.write("<!-- PyHST VOLUME XML FILE -->\n")
    f.write("<tomodb2>\n")
    f.write("<reconstruction>\n")
    f.write("<idAc>%s</idAc>\n" % Parameters.idAc)
    f.write("<listSubVolume>\n")
    f.write("<subVolume>\n")
    f.write("<SUBVOLUME_NAME>%s</SUBVOLUME_NAME>\n" % basename(name))

    sx =  Parameters.END_VOXEL_1 - Parameters.START_VOXEL_1 +1
    sy =  Parameters.END_VOXEL_2 - Parameters.START_VOXEL_2 +1
    sz =  Parameters.END_VOXEL_3 - Parameters.START_VOXEL_3 +1
    x =  Parameters.START_VOXEL_1 
    y =   Parameters.START_VOXEL_2 
    z =   Parameters.START_VOXEL_3 
    
    f.write("<SIZEX>%d</SIZEX>\n" % sx )
    f.write("<SIZEY>%d</SIZEY>\n" % sy )
    f.write("<SIZEZ>%d</SIZEZ>\n" % sz )
    f.write("<ORIGINX>%d</ORIGINX>\n" % x )
    f.write("<ORIGINY>%d</ORIGINY>\n" % y )
    f.write("<ORIGINZ>%d</ORIGINZ>\n" % z )
    f.write("<DIM_REC>%d</DIM_REC>\n" %  (sx*sy*sz) )


    f.write("<BYTE_ORDER>%s</BYTE_ORDER> \n" % BYTEORDER)


    f.write("</subVolume>\n")
    f.write("</listSubVolume>\n")


#    f.write("<complementaryInfoList>\n")
#    f.write("</complement>\n")
#
#    f.write("</complementaryInfoList>\n")
#
    
    f.write("</reconstruction>\n")
    f.write("</tomodb2>\n")


    try:
      import os
      os.system("chmod og+r "+  name+".xml"      )
      os.system("chmod og+r "+  name   )
    except:
      pass


module_mode = False

if len(sys.argv)!=2 :
  try:
    sys._getframe(1)
    module_mode = True
  except:
    raise Exception, " WRONG NUMBER OF ARGUMENTS, give par_file_name "



if module_mode:
    s = sys._getframe(1).f_globals['parameters']
else:
    filename=sys.argv[1]
    try:    
	f=open(filename,"r")
    except:
	logger.error(" problems reading file %s", filename)
	raise Exception, " EXITING "

    s=f.read()
    f.close()

s=treat_par_file(s)


NO=0
YES=1
class Parameters:
    METHOD = "FBP"
    FFT_OVERSAMPLING_FACTOR = 2
    DFI_KERNEL_SIZE = 7
    DFI_KERNEL_POINTS = 1023
    
    RECONSTRUCT_FROM_SINOGRAMS=0
    NO_SINOGRAM_FILTERING=0
    DO_CCD_FILTER=0
    DO_SINO_FILTER=0

    CCD_FILTER=""
    SINO_FILTER=""
    CCD_FILTER_PATH="/"
    SINO_FILTER_PATH="/"

    CCD_FILTER_PARA ={}
    SINO_FILTER_PARA={}

    DO_AXIS_CORRECTION=0
    AXIS_CORRECTION_FILE=""
    OPTIONS= { 'padding':'E' }    
    FOURIER_FILTER=None
    FOURIER_FILTER_HAS_RAMP = NO

    OUTPUT_SINOGRAMS=0

    FILE_INTERVAL =      1

    BICUBIC=0
    SUMRULE=0

    SAVE_JPEG_SLICES=0
    JPEG_QUALITY=100
    DO_HISTOGRAM=0
    DO_PROJECTION_MEDIAN=0
    DO_PROJECTION_MEAN=0

    DO_AXIS_LONGITUDINAL_CORRECTION=0

    PROJECTION_MEDIAN_FILENAME = "projectionmedian.edf"

    DOUBLEFFCORRECTION=0

    ZEROCLIPVALUE=1.0e-9
    ONECLIPVALUE =None

    OFFSETRADIOETFF=0.0
    BINNING = None
    XYCORRECTIONS=0
    ANGLES_FILE=None
    PENTEZONE=10.0

    ZEROOFFMASK=0


    
    exec(s)



if( Parameters.ANGLES_FILE is not None  ):
    angles_ = string.split( open( Parameters.ANGLES_FILE).read())
    angles=[]
    for a in angles_:
        try:
            angles.append(   string.atof(a)     )
        except:
            break
    angles=(Numeric.array(angles)*DEG2RAD).astype("f")
else:
    angles=0




PARALLEL_MACHINE = 0

logger.debug(Parameters.ANGLE_BETWEEN_PROJECTIONS*DEG2RAD)


if(len(sys.argv)==6):

  PARALLEL_MACHINE = 1

  nmachines = int(sys.argv[5])
  nmach     = int(sys.argv[4])
  tot_slices = Parameters.END_VOXEL_3 - Parameters.START_VOXEL_3 + 1
  partial_slices = int ( 0.999999 +  (  0.0 + tot_slices )/ nmachines   ) 
  START_VOXEL_3 =  partial_slices  * nmach
  END_VOXEL_3   =  min ( Parameters.END_VOXEL_3 , partial_slices  * ( nmach+1) -1 )  
  Parameters.START_VOXEL_3   = START_VOXEL_3 +1
  Parameters.END_VOXEL_3     = END_VOXEL_3 +1
  Parameters.OUTPUT_FILE = Parameters.OUTPUT_FILE  + "_" + str (nmach) 


  logger.debug("%i %i", START_VOXEL_3,  END_VOXEL_3)
 
if(( Parameters.RECONSTRUCT_FROM_SINOGRAMS) or (module_mode)):

################################################
 if module_mode:
    Parameters.NUM_FIRST_SINOGRAM = 0
    Parameters.NUM_LAST_SINOGRAM = len(sys._getframe(1).f_globals['sinograms']) - 1
    numpj = Parameters.NUM_LAST_SINOGRAM
 elif( "NUM_LAST_IMAGE" in dir(Parameters)  and "NUM_FIRST_IMAGE" in dir(Parameters) ):
    numpj = Parameters.NUM_LAST_IMAGE+1-Parameters.NUM_FIRST_IMAGE
 elif("NUM_IMAGE_2" in dir(Parameters)):
    numpj = Parameters.NUM_IMAGE_2 
 else:
	raise Exception, "you should define at least NUM_IMAGE_2 or NUM_LAST_IMAGE+NUM_FIRST_IMAGE in input file" 
 axis_corrections  = Numeric.zeros([numpj ],"f")
 axis_correctionsL = Numeric.zeros([numpj ],"f")



 PIECE_MARGE=1
 Parameters.DO_AXIS_LONGITUDINAL_CORRECTION = 0 
 logger.info("Parameters.DO_AXIS_LONGITUDINAL_CORRECTION %i", Parameters.DO_AXIS_LONGITUDINAL_CORRECTION)
 if( Parameters.DO_AXIS_CORRECTION ):
     s=open(Parameters.AXIS_CORRECTION_FILE,"r").read()
     sl=string.split(s)
     sl2 = string.split(s,"\n")
     if len(  string.split(sl2[0]) )==2:
         Parameters.DO_AXIS_LONGITUDINAL_CORRECTION = 1
     if Parameters.DO_AXIS_LONGITUDINAL_CORRECTION :
         for i in range(numpj) :
             axis_corrections[i], axis_correctionsL[i]   = map( string.atof, string.split(sl2[i]))
             PIECE_MARGE = max(PIECE_MARGE,int( abs( axis_correctionsL[i] )+1) )
     else:
         for i in range(numpj) :
             axis_corrections[i] =string.atof(sl[i])
     logger.info("Parameters.DO_AXIS_LONGITUDINAL_CORRECTION %i", Parameters.DO_AXIS_LONGITUDINAL_CORRECTION)


##########################################################


 if( PARALLEL_MACHINE):
     raise Exception, " PARALLEL_MACHINE not yet implemented for reconstruction from sinograms "


 if Parameters.DO_SINO_FILTER:
      logger.info(" carico filtro ")
      sys.path=[Parameters.SINO_FILTER_PATH] + sys.path
      exec ("from %s import Filter as SINO_FILTER" % Parameters.SINO_FILTER,locals() , globals() )

 sino_pos = 0
 for nsino in range( Parameters.NUM_FIRST_SINOGRAM, Parameters.NUM_LAST_SINOGRAM+1 ):
    if module_mode:
	SINO = sys._getframe(1).f_globals['sinograms'][nsino]
    else:
	name=""
	if(Parameters.NUMBER_LENGTH_VARIES):
    	    name=name+("%d"%nsino)+Parameters.SINOGRAM_POSTFIX
	else:
    	    number = "%d"%nsino
	    if(len(number)< Parameters.LENGTH_OF_NUMERICAL_PART ):
	        number = ("0" *  (Parameters.LENGTH_OF_NUMERICAL_PART-len(number)))+ number
    	    else:
        	number=number[:Parameters.LENGTH_OF_NUMERICAL_PART]

    	    name=name+number+Parameters.SINOGRAM_POSTFIX
	
	name=Parameters.SINOGRAM_PREFIX+name

	# name_out=name[string.rfind(name,"/")+1:] +".vol"

	logger.info(" reading file %s", name)
	start_time = time.time();

	if name.endswith('.edf'):
    	    SINO = EdfFile.EdfFile(name).GetData(0, DataType="FloatValue")
	else:
	    SINO = Numeric.array(Image.open(name))

	edf_timer += time.time() - start_time 
	# SINO = Numeric.transpose(SINO)*Numeric.array([1],"f")

    logger.info( SINO.shape )
    SINO.shape=(1,)+SINO.shape

    if(nsino==Parameters.NUM_FIRST_SINOGRAM):
        name_out=Parameters.OUTPUT_FILE
	SINOWORK=Numeric.empty([Parameters.NUM_LAST_SINOGRAM - Parameters.NUM_FIRST_SINOGRAM + 1,SINO.shape[1], SINO.shape[2]], dtype="f")
	start_time = time.time()
        
        normalise = Numeric.array([1.0],"f")

        pyhst=PyHST_c.PyHST( logger, Parameters.NUM_LAST_SINOGRAM - Parameters.NUM_FIRST_SINOGRAM + 1, name_out ,
                  Parameters.METHOD,
                  Parameters.FFT_OVERSAMPLING_FACTOR,
                  Parameters.DFI_KERNEL_SIZE,
                  Parameters.DFI_KERNEL_POINTS,
                  Parameters.OVERSAMPLING_FACTOR,
                  Parameters.START_VOXEL_1-1 +1,
                  Parameters.START_VOXEL_2-1 +1,
                  Parameters.END_VOXEL_1-Parameters.START_VOXEL_1+1,
                  Parameters.END_VOXEL_2-Parameters.START_VOXEL_2+1,
                  Parameters.ROTATION_AXIS_POSITION,
                  Parameters.ANGLE_OFFSET*DEG2RAD,
                  Parameters.ANGLE_BETWEEN_PROJECTIONS*DEG2RAD,
		  Parameters.NO_SINOGRAM_FILTERING,
                  SINO.shape[2],
                  SINO.shape[1],
                  SINOWORK, axis_corrections ,
                  Parameters.BICUBIC,
                  Parameters.SUMRULE,
                  angles,
                  Parameters.PENTEZONE,
                  Parameters.ZEROOFFMASK,
                  normalise
		  )

        if ( Parameters.FOURIER_FILTER is not None):
              pyhst.setFilterFunct( Parameters() , Parameters.FOURIER_FILTER,Parameters.FOURIER_FILTER_HAS_RAMP )

	init_timer += time.time() - start_time

        WriteInfo( name_out,
               Parameters.END_VOXEL_1-Parameters.START_VOXEL_1+1,
               Parameters.END_VOXEL_2-Parameters.START_VOXEL_2+1,
               -Parameters.NUM_FIRST_SINOGRAM +  Parameters.NUM_LAST_SINOGRAM+1
              )

	if pyhst.astra_scaling:
	    normalise = Numeric.array([ (1.0e6*0.01 / Parameters.IMAGE_PIXEL_SIZE_1)],"f")
	else:
	    normalise = Numeric.array([(math.pi / 2.0) * (1.0 /  SINO.shape[1]) *   (1.0e6*0.01 / Parameters.IMAGE_PIXEL_SIZE_1)],"f")
        logger.debug(" # " * 80)
        logger.debug(" LE FACTEUR DE NORMALISATION EST %s" , normalise)
        logger.debug(" # " * 80)
    if( Parameters.DO_SINO_FILTER ):
       SINO_FILTER(SINO,    Parameters.SINO_FILTER_PARA  )

    SINOWORK[sino_pos,:,:] =SINO[0,:,:] * normalise 
    sino_pos += 1


 start_time = time.time()
 pyhst.calcSlices(sino_pos ,Parameters.OPTIONS)
 c_timer += time.time() - start_time

else:


####################################################################
  # cherche le fichier xml contenant le tag idAc
  idAc = "N_A_"
  file_xml_idac = Parameters.FILE_PREFIX+"_idAc.xml"
  try:
      idAc = simplereadxml(file_xml_idac,"idAc")
  except:
      pass
  Parameters.idAc = idAc



################################################

  numpj = Parameters.NUM_LAST_IMAGE+1-Parameters.NUM_FIRST_IMAGE
  
  axis_corrections=Numeric.zeros([numpj ],"f")
  axis_correctionsL = Numeric.zeros([numpj ],"f")

  Xcorr=0
  Ycorr=0
  PIECE_MARGE=1
  if( Parameters.XYCORRECTIONS):
      
     start_time = time.time();
     Xcorr=EdfFile.EdfFile(Parameters.XYCORRECTIONS+"X.edf")
     Xcorr=(Xcorr.GetData(0,DataType="FloatValue")).astype("f")
     
     Ycorr=EdfFile.EdfFile(Parameters.XYCORRECTIONS+"Y.edf")
     Ycorr=( Ycorr.GetData(0,DataType="FloatValue")).astype("f")
     edf_timer += time.time() - start_time;

     
     PIECE_MARGE = max(PIECE_MARGE, max(max(abs(Ycorr)+1.0)))

  
  Parameters.DO_AXIS_LONGITUDINAL_CORRECTION = 0 
  if( Parameters.DO_AXIS_CORRECTION ):
     s=open(Parameters.AXIS_CORRECTION_FILE,"r").read()
     sl=string.split(s)
     sl2 = string.split(s,"\n")
     
     while( sl2[0][0]=="#" ):
         sl2=sl2[1:]
         
     if len(  string.split(sl2[0]) )==2:
         Parameters.DO_AXIS_LONGITUDINAL_CORRECTION = 1
     if Parameters.DO_AXIS_LONGITUDINAL_CORRECTION :
         for i in range(numpj) :
             axis_corrections[i], axis_correctionsL[i]   = map( string.atof, string.split(sl2[i]))
             PIECE_MARGE = max(PIECE_MARGE, int(abs( axis_correctionsL[i] )+1) )
     else:
         for i in range(numpj) :
             axis_corrections[i]=string.atof(sl[i])
     logger.info("Parameters.DO_AXIS_LONGITUDINAL_CORRECTION %s", Parameters.DO_AXIS_LONGITUDINAL_CORRECTION)
     logger.info(PIECE_MARGE)
     
##########################################################    
  
  if Parameters.DO_CCD_FILTER:
      sys.path=[Parameters.CCD_FILTER_PATH] + sys.path
      logger.info ("loading filter")
      logger.info (sys.path[0])
      exec ("from %s import Filter as CCD_FILTER" % Parameters.CCD_FILTER )
      logger.info(" IMPORTED CCD-FILTER")
      logger.info(CCD_FILTER)
  else:
      CCD_FILTER=0
      
  if Parameters.DO_SINO_FILTER:
      sys.path=[Parameters.SINO_FILTER_PATH] + sys.path
      exec ("from %s import Filter as SINO_FILTER" % Parameters.SINO_FILTER)


  ########################################################################################################  #
  # determine how many slice to read at once
  # --------------------------------------------
  Parameters.nsinos = (Parameters.NUM_LAST_IMAGE - Parameters.NUM_FIRST_IMAGE+1)

  if hasattr(Parameters,'OUTPUT_MEGABYTES'):
    output_limit  = Parameters.OUTPUT_MEGABYTES*1.0e6/(   ( Parameters.END_VOXEL_2 - Parameters.START_VOXEL_2+1 ) * (Parameters.END_VOXEL_1 - Parameters.START_VOXEL_1+1) *4  )
  else:
    output_limit  = 2000*1.0e6/(   ( Parameters.END_VOXEL_2 - Parameters.START_VOXEL_2+1 ) * (Parameters.END_VOXEL_1 - Parameters.START_VOXEL_1+1) *4  )
    
  if(Parameters.ROTATION_VERTICAL):
      input_limit = Parameters.SINOGRAM_MEGABYTES*1.0e6/( Parameters.NUM_IMAGE_1 *Parameters.nsinos  *4  )
  else:
      input_limit = Parameters.SINOGRAM_MEGABYTES*1.0e6/( Parameters.NUM_IMAGE_2 *Parameters.nsinos) 

  Parameters.nslices_atonce = int(min( input_limit, output_limit )) 

  Parameters.nslices_atonce = int(Parameters.nslices_atonce/(4 * PyHST_c.RECONSTRUCTORS)) * (4 * PyHST_c.RECONSTRUCTORS)
  if Parameters.nslices_atonce==0:
    Parameters.nslices_atonce = int(min( input_limit, output_limit )) 

  # =======================================================================================================

  ###############################################################################
  # How many slices we want in total
  # ----------------------------------------------------------------------------

  Parameters.tot_slices = Parameters.END_VOXEL_3 - Parameters.START_VOXEL_3 + 1

  # ///////////////////////////////////////////////////////////////////////////////////////////////////////////////


  #########################################################################################
  # read and store the BACKGROUND in variable background fully dimensioned
  # -------------------------------------------------------------------------------------
  start_time = time.time();
  if Parameters.BINNING is None:
      if(Parameters.SUBTRACT_BACKGROUND):
         logger.info(" reading edf file for background %s", Parameters.BACKGROUND_FILE)
         dark=EdfFile.EdfFile(Parameters.BACKGROUND_FILE)
         background=dark.GetData(0,DataType="FloatValue")
      else:
        background=Numeric.zeros([ Parameters.NUM_IMAGE_2, Parameters.NUM_IMAGE_1 ],"f")

      if(Parameters.ROTATION_VERTICAL==0):
          background=Numeric.transpose(background)
  else:
      background=extract_edf( Parameters.BACKGROUND_FILE, Parameters.ROTATION_VERTICAL , (0,0),
                              (Parameters.NUM_IMAGE_2, Parameters.NUM_IMAGE_1)  , Parameters.BINNING   )
  edf_timer += time.time() - start_time


  logger.info(" backgroundshape = %s" , background.shape)
  logger.info(" Parameters.END_VOXEL_3 %s Parameters.START_VOXEL_3  %s" , Parameters.END_VOXEL_3 , Parameters.START_VOXEL_3)
  logger.info("  Parameters.NUM_IMAGE_2 %s, Parameters.NUM_IMAGE_1 %s",  Parameters.NUM_IMAGE_2, Parameters.NUM_IMAGE_1)



  ###############################################################################################
  # create the big array to contain sinograms
  # --------------------------------------------------------------------------------------------
  tot_passes = int((Parameters.tot_slices-0.001)/ Parameters.nslices_atonce ) + 1
  num_projections =  Parameters.NUM_LAST_IMAGE + 1 - Parameters.NUM_FIRST_IMAGE

  if (Parameters.ROTATION_VERTICAL):
    dim_1 = Parameters.NUM_IMAGE_1
    dim_2 = Parameters.nsinos
    dim_3 = Parameters.nslices_atonce
    if( Parameters.DO_PROJECTION_MEDIAN or Parameters.DO_PROJECTION_MEAN ):
        PROJECTION_MEDIAN = Numeric.zeros( [Parameters.NUM_IMAGE_2,Parameters.NUM_IMAGE_1],"f" ) 
    normalise = Numeric.array([(math.pi / 2.0) * (1.0 / num_projections) *   (1.0e6*0.01 / Parameters.IMAGE_PIXEL_SIZE_1)],"f")
  else:
    dim_1 = Parameters.NUM_IMAGE_2
    dim_2 = Parameters.nsinos
    dim_3 = Parameters.nslices_atonce
    if( Parameters.DO_PROJECTION_MEDIAN or Parameters.DO_PROJECTION_MEAN ):
        PROJECTION_MEDIAN = Numeric.zeros( [Parameters.NUM_IMAGE_1,Parameters.NUM_IMAGE_2],"f" ) 
    normalise = Numeric.array([(math.pi / 2.0) * (1.0 / num_projections) *   (1.0e6*0.01 / Parameters.IMAGE_PIXEL_SIZE_2)],"f")

  if Parameters.BICUBIC:
    normalise=Numeric.array([normalise[0]/ Parameters.BICUBIC],"f")

  if(Parameters.DO_PROJECTION_MEDIAN or Parameters.DO_PROJECTION_MEAN):
    normalise = Numeric.array([1.0],"f")

  if tot_passes > 1:
    multipass = True
  else:
    if Parameters.ROTATION_VERTICAL:
	image_slices =  Parameters.NUM_IMAGE_1
    else:
	image_slices =  Parameters.NUM_IMAGE_2
	
    if 2 * Parameters.tot_slices < image_slices:
	multipass = True
    else:
	multipass = False
    
  suffix3 = Parameters.FILE_POSTFIX[-4:]
  if simulate_data:
    image_reader = SIMULATEReader(Parameters, multipass)
  elif suffix3 == ".edf":
    if parallel_io:
        logger.info("EDFMultiReader")
        image_reader = EDFMultiReader(Parameters, multipass)
    else:
        logger.info("EDFReader")
        image_reader = EDFReader(Parameters, multipass)
  else:
    if use_vips:
        logger.info("VIPSReader")
	image_reader = VIPSReader(Parameters, multipass)
    else:
        logger.info("IMAGEReader")
	image_reader = IMAGEReader(Parameters, multipass)
    
  # >>>>>>>>>>>>>>>> the BIG thing <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  #BIG_SINOS = Numeric.zeros( [dim_3,dim_2,dim_1],"f" )  
  # another reference to access the BIG thing with a different scheme
  #BIG_SINOS_natural_ordering = Numeric.swapaxes(BIG_SINOS ,0,1)

  start_time = time.time()
  pyhst=PyHST_c.PyHST(logger,
                  Parameters.nslices_atonce, 
		  Parameters.OUTPUT_FILE ,
                  Parameters.METHOD,
                  Parameters.FFT_OVERSAMPLING_FACTOR,
                  Parameters.DFI_KERNEL_SIZE,
                  Parameters.DFI_KERNEL_POINTS,
                  Parameters.OVERSAMPLING_FACTOR,
                  Parameters.START_VOXEL_1 -1 +1,
                  Parameters.START_VOXEL_2 -1 +1,
                  Parameters.END_VOXEL_1-Parameters.START_VOXEL_1+1,
                  Parameters.END_VOXEL_2-Parameters.START_VOXEL_2+1,
                  Parameters.ROTATION_AXIS_POSITION,
                  Parameters.ANGLE_OFFSET*DEG2RAD,
                  Parameters.ANGLE_BETWEEN_PROJECTIONS*DEG2RAD,
		  Parameters.NO_SINOGRAM_FILTERING,
                  dim_1,
                  num_projections,
                  [2 if (tot_passes>1) else 1, dim_3, dim_2, dim_1], #BIG_SINOS,
                  axis_corrections ,
                  Parameters.BICUBIC,
                  Parameters.SUMRULE ,
                  (angles),
                  Parameters.PENTEZONE,
                  Parameters.ZEROOFFMASK,
                  normalise)

  if ( Parameters.FOURIER_FILTER is not None):
         pyhst.setFilterFunct( Parameters(), Parameters.FOURIER_FILTER,Parameters.FOURIER_FILTER_HAS_RAMP )
  
  init_timer += time.time() - start_time

  num_files = num_projections + 128
  try:
     resource.setrlimit(resource.RLIMIT_NOFILE, (num_files, num_files))
  except:
     logger.warning("Failed to increase limit of open files")
  limits = resource.getrlimit(resource.RLIMIT_NOFILE)
  
  start_time = time.time()

  if (image_reader.preload) and (limits[0] >= num_files):
     images = []
     for i_pro in range(Parameters.NUM_FIRST_IMAGE, Parameters.NUM_LAST_IMAGE + 1, 1):
        images.insert(i_pro, open_image_file(image_reader, i_pro, Parameters))
  else:
    images = None
    if image_reader.preload:
	logger.warning("The limit on open files prevents pre-parsing of the headers")

  edf_timer += time.time() - start_time
  
  #######################################################################################
  # start the loop over memory bunches
  # --------------------------------------------------------------------------------------
#  while 1:
  for i_pass in range( tot_passes  ):
#  for i_pass in range( tot_passes + 1  ):
#   if i_pass < tot_passes:
    start_time = time.time()
    
    logger.warning(" pass %i of %i", i_pass + 1, tot_passes)

    BIG_SINOS = pyhst.SINOGRAMS
    BIG_SINOS_natural_ordering = Numeric.swapaxes(BIG_SINOS, 0, 1)

    first_slice = i_pass * Parameters.nslices_atonce + Parameters.START_VOXEL_3 - 1
    last_slice  = min(first_slice + Parameters.nslices_atonce - 1, Parameters.END_VOXEL_3 - 1) 

    logger.info(" processing slice %i up to (included) %i", first_slice+1, last_slice+1)

    pre_timer += time.time() - start_time
    start_time = time.time() 
    
    if reconstruct:
        pyhst.startPreprocessing(1 + last_slice - first_slice)
    
    mem_timer += time.time() - start_time
    start_time = time.time() 

    ######################################################################
    # pos_edf and size_edf are passed to the Gobbo's EdfFile routine
    # to extract the interesting part of the edf images
    # -------------------------------------------------------------------
    if(Parameters.ROTATION_VERTICAL==1):
      pos_edf=   [first_slice,0]
      size_edf = [last_slice-first_slice+1, Parameters.NUM_IMAGE_1 ]
  
      if( Parameters.DO_CCD_FILTER or Parameters.DO_AXIS_LONGITUDINAL_CORRECTION or Parameters.XYCORRECTIONS):
        PIECE_MARGE=int(PIECE_MARGE+0.5)
	pos_edf_, size_edf_ = OverlappingLogic (pos_edf,  size_edf, 0, Parameters.NUM_IMAGE_2-1 , PIECE_MARGE )
      else:
        pos_edf_, size_edf_ = pos_edf, size_edf

      logger.info( background.shape )
      logger.info( "%s %s", pos_edf_[0], pos_edf_[0]+size_edf_[0] )
      new_background = background[pos_edf_[0]:pos_edf_[0]+size_edf_[0]]

    if(Parameters.ROTATION_VERTICAL==0):
      pos_edf=   [0,first_slice]
      size_edf = [  Parameters.NUM_IMAGE_2, last_slice-first_slice+1      ]

      if( Parameters.DO_CCD_FILTER   or Parameters.DO_AXIS_LONGITUDINAL_CORRECTION   or Parameters.XYCORRECTIONS  ):
	pos_edf_, size_edf_ = OverlappingLogic (pos_edf,  size_edf, 1,  Parameters.NUM_IMAGE_1-1 ,PIECE_MARGE  )
      else:
        pos_edf_, size_edf_ = pos_edf, size_edf

      new_background = background[pos_edf_[1]:pos_edf_[1]+size_edf_[1]]


    ################################################################################
    # the background that is going to be used for the interesting pat
    # ------------------------------------------------------------------------------
    #
    # new_background = background[first_slice:last_slice+1]
  

    #################################################################################
    # ipro is an index running over the  projection edf-images
    # It can be used directly to access the files by postpending it to the file prefix
    # -------------------------------------------------------------------------------

    if reconstruct:
     if(Parameters.DOUBLEFFCORRECTION):
          ffcorr_  =   EdfFile.EdfFile ( Parameters.DOUBLEFFCORRECTION )
          ffcorr_  =   ffcorr_.GetData(0,DataType="FloatValue")
          
              
          if  ( Parameters.ROTATION_VERTICAL==0):
              ffcorr_  =   Numeric.transpose(ffcorr_)
    
          if ( Parameters.BINNING is not None):
              BINNING= Parameters.BINNING
              s1,s2=ffcorr_.shape
              ffcorr_ = ffcorr_[ 0:  ( s1/BINNING)*BINNING   ,       0:  ( s2/BINNING)*BINNING      ]
              ffcorr_=Numeric.reshape(ffcorr_,[    s1/BINNING, BINNING , s2/BINNING, BINNING ])
              ffcorr_ = Numeric.swapaxes(ffcorr_,1,2)
              ffcorr_=Numeric.sum(Numeric.sum(ffcorr_, axis=-1), axis=-1)/BINNING/BINNING
             
          if( Parameters.TAKE_LOGARITHM):
              ffcorr_[:]  =   Numeric.exp(ffcorr_)

          ffcorr_=ffcorr_.astype("f")

    pre_timer += time.time() - start_time

    for i_pro in range(Parameters.NUM_FIRST_IMAGE, Parameters.NUM_LAST_IMAGE+1,1):
      start_time = time.time();

      if(( i_pro- Parameters.NUM_FIRST_IMAGE)%100 ==0 ):
        logger.info( " READING projection # %i of %i", i_pro- Parameters.NUM_FIRST_IMAGE, Parameters.NUM_LAST_IMAGE+1-Parameters.NUM_FIRST_IMAGE )

      if(Parameters.CORRECT_FLATFIELD):
        if( Parameters.FLATFIELD_CHANGING):
          i_FF, FF_FILE_INTERVAL=ModulesForFF( (i_pro*Parameters.FILE_INTERVAL - Parameters.NUM_FIRST_IMAGE  ),   Parameters.FF_FILE_INTERVAL)
          if( i_FF<Parameters.FILE_INTERVAL ):

            ########################################################
            # convert ipro to the index for the FF
            # ------------------------------------------------------------
            nFF=(i_pro*Parameters.FILE_INTERVAL-i_FF) - Parameters.NUM_FIRST_IMAGE + Parameters.FF_NUM_FIRST_IMAGE
            nFF = min(nFF,  Parameters.FF_NUM_LAST_IMAGE)
            # ////////////////////////////////////////////////////////////////
            ############################################################
            #   read the fisrt  FF of the interpolation range
            # -----------------------------------------------------------
            FF_a = extract_edf_N(Parameters.FF_PREFIX, nFF,
                               Parameters.FF_NUMBER_LENGTH_VARIES,
                               Parameters.FF_LENGTH_OF_NUMERICAL_PART,
                               Parameters.FF_POSTFIX,
                               Parameters.ROTATION_VERTICAL, 
                               pos_edf_, size_edf_, Parameters, rescaleN=0,
                               BINNING=Parameters.BINNING)
	    
#	    if( Parameters.DO_CCD_FILTER):
#		FF_a = Filter_and_Trim( FF_a , pos_edf_, size_edf_, pos_edf, size_edf, CCD_FILTER  , Parameters.CCD_FILTER_PARA )


            # ////////////////////////////////////////////////////////////
            FF_b = FF_a
          if( i_FF>= Parameters.FILE_INTERVAL and i_FF< 2*Parameters.FILE_INTERVAL):
            ########################################################
            # convert ipro to the index for the FF
            # ------------------------------------------------------------
            nFF = (i_pro*Parameters.FILE_INTERVAL-i_FF+1) - 1 - Parameters.NUM_FIRST_IMAGE + Parameters.FF_NUM_FIRST_IMAGE+  FF_FILE_INTERVAL  
            nFF = min(nFF,  Parameters.FF_NUM_LAST_IMAGE)
            # ////////////////////////////////////////////////////////////////
            ############################################################
            # read the second  FF of the interpolation range
            # -----------------------------------------------------------
            FF_b = extract_edf_N(Parameters.FF_PREFIX,
                           nFF ,
                               Parameters.FF_NUMBER_LENGTH_VARIES,
                               Parameters.FF_LENGTH_OF_NUMERICAL_PART,
                               Parameters.FF_POSTFIX,
                               Parameters.ROTATION_VERTICAL, 
                               pos_edf_, size_edf_, Parameters, rescaleN=0,
                               BINNING=Parameters.BINNING)


#	    if( Parameters.DO_CCD_FILTER):
#		FF_b = Filter_and_Trim( FF_b , pos_edf_, size_edf_, pos_edf, size_edf, CCD_FILTER , Parameters.CCD_FILTER_PARA  )

             # ////////////////////////////////////////////////////////////

          ##########################################################################
          # linear interpolation
          # ------------------------------------------------------------------------
          FF = ( FF_a*Numeric.array([FF_FILE_INTERVAL-i_FF]).astype("f") +
                 FF_b *Numeric.array([i_FF]).astype("f")  )/Numeric.array([FF_FILE_INTERVAL]    ).astype("f")
          # /////////////////////////////////////////////////////////////////////////////////////////////////////////////
        else:
            if(ipro==Parameters.NUM_FIRST_IMAGE):
               ###################################################
               # read at the beginning of the scan the only
               # FF available
               # -------------------------------------------------
	       start_time = time.time();
               FF = extract_edf(Parameters.FLATFIELD_FILE,
                                 Parameters.ROTATION_VERTICAL, 
                                 pos_edf_, size_edf_,
                               BINNING=Parameters.BINNING)
	       edf_timer += time.time() - start_time;

#	       if( Parameters.DO_CCD_FILTER):
#		  FF = Filter_and_Trim( FF , pos_edf_, size_edf_, pos_edf, size_edf, CCD_FILTER , Parameters.CCD_FILTER_PARA  )



               # /////////////////////////////////////////////////

        if( Parameters.SUBTRACT_BACKGROUND):
          newFF = FF - new_background
        else:
          newFF = FF

        # clipping newFF 
        newFF=Numeric.maximum(newFF,Numeric.array([0.1],"f") )

      pre_timer += time.time() - start_time;

      start_time = time.time()
      #newitem =  extract_edf_N(Parameters.FILE_PREFIX, i_pro, Parameters.NUMBER_LENGTH_VARIES, Parameters.LENGTH_OF_NUMERICAL_PART, Parameters.FILE_POSTFIX, Parameters.ROTATION_VERTICAL, pos_edf_, size_edf_, Parameters, BINNING=Parameters.BINNING)
      if images is not None:
        newitem = read_image_data(image_reader, i_pro - Parameters.NUM_FIRST_IMAGE, images[i_pro - Parameters.NUM_FIRST_IMAGE], Parameters.ROTATION_VERTICAL, pos_edf_, size_edf_, Parameters.BINNING)
      else:
        newitem = read_image_data(image_reader, i_pro - Parameters.NUM_FIRST_IMAGE, open_image_file(image_reader, i_pro, Parameters), Parameters.ROTATION_VERTICAL, pos_edf_, size_edf_, Parameters.BINNING)
      edf_timer += time.time() - start_time

      start_time = time.time()
      if reconstruct:
       if( Parameters.SUBTRACT_BACKGROUND):
          newitem -= new_background
        
       if( Parameters.CORRECT_FLATFIELD):
        newitem=newitem+ Numeric.array([ Parameters.OFFSETRADIOETFF],"f") 
        newFF=newFF+ Numeric.array([ Parameters.OFFSETRADIOETFF],"f") 
        newitem /=  newitem 

       if( Parameters.TAKE_LOGARITHM):
        # clipping  newitem
        newitem=Numeric.maximum(newitem,Numeric.array([Parameters.ZEROCLIPVALUE],"f") )
        if( Parameters.ONECLIPVALUE is not None):
            newitem=Numeric.minimum(newitem,Numeric.array([Parameters.ONECLIPVALUE] ,"f") )


       if(Parameters.DOUBLEFFCORRECTION):
          
          ffcorr  =   ffcorr_ [pos_edf_[0]:pos_edf_[0]+size_edf_[0]]
          if( Parameters.TAKE_LOGARITHM):
              newitem *= ffcorr
          else:
              newitem -= ffcorr

       if( Parameters.DO_CCD_FILTER or Parameters.DO_AXIS_LONGITUDINAL_CORRECTION or Parameters.XYCORRECTIONS):
          if Parameters.DO_AXIS_LONGITUDINAL_CORRECTION  or Parameters.XYCORRECTIONS:
              from  interpola_filter import Filter as CCD_AXIS_LONGITUDINAL_CORRECTION

              if(axis_correctionsL):
                  correctionL = axis_correctionsL[i_pro]
              else:
                  correctionL = 0
          else:
              CCD_AXIS_LONGITUDINAL_CORRECTION = 0 
              correctionL = 0

          if Ycorr:
              if(Parameters.ROTATION_VERTICAL==1 ):
                  XcorrF   =  Xcorr [pos_edf[0]:pos_edf[0]+size_edf[0]]
                  YcorrF   =  Ycorr [pos_edf[0]:pos_edf[0]+size_edf[0]]
              else:
                  XcorrF   =  Xcorr [:, pos_edf[1]:pos_edf[1]+size_edf[1]]
                  YcorrF   =  Ycorr [:, pos_edf[1]:pos_edf[1]+size_edf[1]]

              XcorrF=Numeric.array(XcorrF)
              YcorrF=Numeric.array(YcorrF)
          else:
              XcorrF=0
              YcorrF=0

          newitem  = Filter_and_Trim( newitem , pos_edf_, size_edf_, pos_edf, size_edf, CCD_FILTER , CCD_AXIS_LONGITUDINAL_CORRECTION , correctionL ,Xcorr, Ycorr, Parameters.CCD_FILTER_PARA)


       if( Parameters.TAKE_LOGARITHM):
        newitem = - Numeric.log(newitem)

      pre_timer += time.time() - start_time

      start_time = time.time()

      if reconstruct: 
        pyhst.transposeSlices(newitem, i_pro - Parameters.NUM_FIRST_IMAGE)

      # throwing new data into the BIG thing
      #throw_start=first_slice -Parameters.START_VOXEL_3 +1
      #throw_end  =last_slice  -Parameters.START_VOXEL_3 +1
#       if( not Parameters.DO_PROJECTION_MEDIAN and (not Parameters.DO_PROJECTION_MEAN)):
#           newitem *= normalise
#       BIG_SINOS_natural_ordering[i_pro-Parameters.NUM_FIRST_IMAGE , 0 :throw_end+1 -throw_start ] = newitem

      newitem = None

      mem_timer += time.time() - start_time
    # calculate slices from first_slice up to last_slice

    # ///////////////////////////////////////////////////////////////////////////////////////  

    start_time = time.time()

    if reconstruct: 
        pyhst.waitPreprocessing()

    mem_timer += time.time() - start_time

#    continue
    start_time = time.time()
    if( Parameters.DO_SINO_FILTER ):
       SINO_FILTER(BIG_SINOS[:last_slice+1-first_slice,:,:],    Parameters.SINO_FILTER_PARA  )

    pre_timer += time.time() - start_time

    if( Parameters.OUTPUT_SINOGRAMS==0 and Parameters.DO_PROJECTION_MEDIAN==0 and Parameters.DO_PROJECTION_MEAN==0):
        logger.debug(" calc slices ")
	start_time = time.time()
        pyhst.calcSlices(last_slice+1-first_slice, Parameters.OPTIONS)
	c_timer += time.time() - start_time
        logger.debug(" calcslices ok ")
    else:
        if( Parameters.OUTPUT_SINOGRAMS):
            for isn in range(first_slice,last_slice+1 ):
                nomesinofile = Parameters.OUTPUT_FILE+"_sino_"+str(isn+1)+".edf"
                logger.debug(nomesinofile)
                f=EdfFile.EdfFile(nomesinofile)
                logger.debug(BIG_SINOS[isn-first_slice,:,:].shape)
                f.WriteImage({}, BIG_SINOS[isn-first_slice,:,:] , Append=0 )
                f=None
                logger.info(" written a sinogram on file %s",nomesinofile)
                
        if( Parameters.DO_PROJECTION_MEDIAN or Parameters.DO_PROJECTION_MEAN ):
	    start_time = time.time()
            PROJECTION_MEDIAN[first_slice:last_slice+1] = pyhst.calcMedian(last_slice+1-first_slice, Parameters.DO_PROJECTION_MEAN ,Parameters.OPTIONS)
	    c_timer += time.time() - start_time


  start_time = time.time()
  pyhst.waitCompletion()
  c_timer += time.time() - start_time

  if(Parameters.OUTPUT_SINOGRAMS==0 and Parameters.DO_PROJECTION_MEDIAN==0 and Parameters.DO_PROJECTION_MEAN==0 ):
        
      # close the object
      WriteInfo( Parameters. OUTPUT_FILE  ,
                 Parameters.END_VOXEL_1-Parameters.START_VOXEL_1+1,
                 Parameters.END_VOXEL_2-Parameters.START_VOXEL_2+1,
                 Parameters.tot_slices
                 )

      WriteInfoXML( Parameters. OUTPUT_FILE  ,Parameters
                    )

  if( Parameters.DO_PROJECTION_MEDIAN==1 or  Parameters.DO_PROJECTION_MEAN==1):

      nomePMfile=Parameters.PROJECTION_MEDIAN_FILENAME

      f=EdfFile.EdfFile(nomePMfile)

      if (   Parameters.ROTATION_VERTICAL ==  0  ):
          PROJECTION_MEDIAN= Numeric.swapaxes( PROJECTION_MEDIAN,0,1)

      
      f.WriteImage({}, PROJECTION_MEDIAN  , Append=0 )
      f=None
      logger.info(" written a sinogram on file %s",nomePMfile)

  start_time = time.time()
  pyhst.close()
  c_timer += time.time() - start_time

  start_time = time.time()
  images = []
  edf_timer += time.time() - start_time

  if((Parameters.DO_HISTOGRAM)or(Parameters.SAVE_JPEG_SLICES)    and ( Parameters.OUTPUT_SINOGRAMS==0 and Parameters.DO_PROJECTION_MEDIAN==0) ):
    logger.info(" Getting min and max value in float volume ... ")
    datamax=0
    datamin=10e6
    size_pixel=4
    sizex,sizey=(Parameters.END_VOXEL_1-Parameters.START_VOXEL_1+1),(Parameters.END_VOXEL_2-Parameters.START_VOXEL_2+1)
    for sliceno in range(0, Parameters.tot_slices):
       Data=read_one_slice(Parameters.OUTPUT_FILE, sizex, sizey, size_pixel, sliceno, Parameters.tot_slices)
       #getting min and max value in image

       data_tmp=Numeric.reshape(Data, [Data.shape[0]*Data.shape[1]])
       max_tmp=max(data_tmp)
       min_tmp=min(data_tmp)
       
       if (max_tmp > datamax): datamax=max_tmp
       if (min_tmp < datamin): datamin=min_tmp
       
    logger.info(" ... done! Min: %s Max: %s", datamin, datamax)
    #to make sure avoiding any rounding errors
    if (datamin<0): datamin=datamin*1.02
    else: datamin=datamin*0.98
    if (datamax>0): datamax=datamax*1.02
    else: datamax=datamax*0.98


  if(Parameters.DO_HISTOGRAM):
    logger.info(" Collecting values in float volume ... ")
    size_pixel=4
    sizex,sizey=(Parameters.END_VOXEL_1-Parameters.START_VOXEL_1+1),(Parameters.END_VOXEL_2-Parameters.START_VOXEL_2+1)
    histogram=Numeric.zeros([1001],"u")
    for sliceno in range(0, Parameters.tot_slices):
       Data=read_one_slice(Parameters.OUTPUT_FILE, sizex, sizey, size_pixel, sliceno, Parameters.tot_slices)

       data_tmp=Numeric.reshape(Data, [Data.shape[0]*Data.shape[1]])
       
       data_tmp=data_tmp - datamin

       data_tmp=data_tmp *(1000.0/ ( (datamax-datamin) *1.00000001 ) )

       data_tmp=data_tmp.astype("i")

       histogram=Numeric.histogram(data_tmp)
#       histogram=arrayfns.histogram(data_tmp)
                           

                
    histofilename=Parameters.OUTPUT_FILE+".histo"
    f=open(histofilename,"w")
    for value in range(0,1001):
     #print ("%f" %((value*(max-min)/100)-abs(min))) + (" %d"%histogram[value])
     if (datamin<0): stringer=("%f" %((value*(datamax-datamin)/1000)-abs(datamin))) + (" %d\n"%histogram[value])
     else: stringer=("%f" %(value*(datamax-datamin)/1000)) + (" %d\n"%histogram[value])
     f.write(stringer)
    f.close()
    logger.info(" ... done! ")

  if (Parameters.SAVE_JPEG_SLICES):
    logger.info(" Converting RAW float-volume file in stack of char-JPEGs ... ")
    size_pixel=4
    sizex,sizey=(Parameters.END_VOXEL_1-Parameters.START_VOXEL_1+1),(Parameters.END_VOXEL_2-Parameters.START_VOXEL_2+1)
    for sliceno in range(0, Parameters.tot_slices):
       Data=read_one_slice(Parameters.OUTPUT_FILE, sizex, sizey, size_pixel, sliceno, Parameters.tot_slices)
       if (datamin<0): Data=Data+abs(datamin)
       Data=Data/Numeric.array([((datamax-datamin)/255.0)], "f")
       DataImage=Data.astype("1")
       im=Image.frombuffer("L", (DataImage.shape[1],DataImage.shape[0]), DataImage.tostring())
       #im.show()
       im = im.transpose(Image.ROTATE_270)
       if ((sliceno+Parameters.START_VOXEL_3)<1000):
           jpgname=Parameters.OUTPUT_FILE+".slice0"+("%d"%(sliceno+Parameters.START_VOXEL_3))+".jpg"
       if ((sliceno+Parameters.START_VOXEL_3)<100):
           jpgname=Parameters.OUTPUT_FILE+".slice00"+("%d"%(sliceno+Parameters.START_VOXEL_3))+".jpg"
       if ((sliceno+Parameters.START_VOXEL_3)<10):
           jpgname=Parameters.OUTPUT_FILE+".slice000"+("%d"%(sliceno+Parameters.START_VOXEL_3))+".jpg"
       if ((sliceno+Parameters.START_VOXEL_3)>=1000):
           jpgname=Parameters.OUTPUT_FILE+".slice"+("%d"%(sliceno+Parameters.START_VOXEL_3))+".jpg"
       im.save(jpgname, "JPEG", quality=Parameters.JPEG_QUALITY)
    logger.info(" ... done! ")
  logger.info(" arrivato alla fin e ")


#print pyhst.recon


max_len = 30
def FormatTimerTitle(title):
    diff = max_len - len(title)
    if diff > 0:
	return title + " " * diff
    else:
	return title

all_timer = time.time() - all_timer

if hasattr(pyhst, "recon_timer"):
    if pyhst.recon != None:
	recons = deepcopy(pyhst.recon)
    
    pyhst_io_timer = pyhst.io_timer
    pyhst_recon_timer = pyhst.recon_timer
    pyhst_comp_timer = pyhst.comp_timer

    start_time = time.time()
    pyhst = None
    #print sys.getrefcount(PyHST_c)
    del sys.modules["PyHST_c_CPU"]
    del PyHST_c
    
    clean_timer = time.time() - start_time
	
    all_timer += clean_timer

    if recons != None:
	for recon in recons:
	  try:
	    if recon.has_key('timer_names'):
		for timer in enumerate(recon['timer_names']):
		    timer_len = len(timer)
		    if timer_len > max_len:
			max_len = timer_len
	  except:
	    continue
    max_len += 2


    recon_timer = 0
    recon_timer_min = all_timer
    recon_timer_max = 0
    recon_timer_n = 0
    if recons != None:
	for recon in recons:
	    try:
		if recon.has_key('slices') and recon.has_key('timers') and (recon["timers"][0] > 0):
	    	    recon_timer += recon["timers"][0]
		    recon_timer_n = recon_timer_n + 1
		    if recon["timers"][0] > recon_timer_max:
			recon_timer_max = recon["timers"][0]
		    if recon["timers"][0] < recon_timer_min:
			recon_timer_min = recon["timers"][0]
	    except:
		continue
	if recon_timer_n > 1:
	    recon_timer = recon_timer / recon_timer_n

    print ""
    print "Timing Information"
    print "------------------"
    print FormatTimerTitle("Overal execution time:             "), all_timer
    print FormatTimerTitle("  Initialization/Cleanup (C-code): "), init_timer + clean_timer
    print FormatTimerTitle("  Loading Data (full):             "), edf_timer
    print FormatTimerTitle("  Preprocessing Data:              "), pre_timer
    print FormatTimerTitle("  Transpose (sched, C-code):       "), mem_timer
    print FormatTimerTitle("  Reconstruction (sched, C-code):  "), c_timer
    print FormatTimerTitle("  Rest (Unknown):                  "), all_timer - init_timer - clean_timer - edf_timer - pre_timer - mem_timer - c_timer
    print ""
    print FormatTimerTitle("Input/Output:                      "), edf_timer + pyhst_io_timer
    print FormatTimerTitle("  Time spent loading EDF files:    "), edf_timer
    print FormatTimerTitle("  Time spent storing results:      "), pyhst_io_timer
    print ""
    print FormatTimerTitle("Time spent within HST library:     "), pyhst_recon_timer
    print FormatTimerTitle("  Reconstruction Schedulling:      "), pyhst_comp_timer
    print FormatTimerTitle("  Actual Reconstruction Time:      "), recon_timer_max
    print ""
    print FormatTimerTitle("FBP Reconstruction (Maximum):      "), recon_timer_max
    print FormatTimerTitle("  Minimal FBP:                     "), recon_timer_min
    print FormatTimerTitle("  Average FBP:                     "), recon_timer
    print ""

    if recon != None:
	for recon in recons:
	  try:
	    if not recon.has_key('timers'):
		continue
	    
    	    if recon.has_key('title'):
		title = recon['title']
    	    else:
		title = "Reconstructor"

	    print title
	    print "-" * len(title)
	    
	    if recon.has_key('slices'):
		slices = recon['slices']
		print FormatTimerTitle("Processed slices: "), slices
	    else:
		slices = 1
	    
	    for i,timing in enumerate(recon['timers']):
		if recon.has_key('timer_names') and recon['timer_names'][i] != None:
		    timer = recon['timer_names'][i].capitalize()
		else:
    		    timer = "Timer " + str(i)
		
		if timer[0] == '*':
		    timer = timer[1:].capitalize()
		    print '%s %f sec      ' % (FormatTimerTitle(timer + ": "), timing)
		else:
		    print '%s %f sec/slice' % (FormatTimerTitle(timer + ": "), timing/slices)
	    print
	  except:
	    continue