#/*##########################################################################
# Copyright (C) 2001-2013 European Synchrotron Radiation Facility
#              PyHST2  
#  European Synchrotron Radiation Facility, Grenoble,France
#
# PyHST2 is  developed at
# the ESRF by the Scientific Software  staff.
# Principal author for PyHST2: Alessandro Mirone.
#
# This program is free software; you can redistribute it and/or modify it 
# under the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 2 of the License, or (at your option) 
# any later version.
#
# PyHST2 is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# PyHST2; if not, write to the Free Software Foundation, Inc., 59 Temple Place,
# Suite 330, Boston, MA 02111-1307, USA.
#
# PyHST2 follows the dual licensing model of Trolltech's Qt and Riverbank's PyQt
# and cannot be used as a free plugin for a non-free program. 
#
# Please contact the ESRF industrial unit (industry@esrf.fr) if this license 
# is a problem for you.
#############################################################################*/

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function


from . import mpistatus


import string
import sys
# raise Exception(sys.argv[0][-12:]=="sphinx-build" )

if(sys.argv[0][-12:]!="sphinx-build"):
    from .  import string_six
    from . import EdfFile
    from . import LTSparsePrep

    from six import string_types

    
    if mpistatus.status :
        from . import setCpuSet

import traceback
import sys
import numpy
Numeric = numpy
import math
import copy
import time
from os.path import basename
import os.path

if mpistatus.status==1 :
    import mpi4py.MPI as MPI
    try:
        import h5py
    except:
        print("H5PY NOT FOUND")
    
    myrank = MPI.COMM_WORLD.Get_rank()
    nprocs = MPI.COMM_WORLD.Get_size()
    mypname = MPI.Get_processor_name()
    comm = MPI.COMM_WORLD
else:
    MPI=None
    myrank=0
    nprocs=1

print(" STATUS OK  " )


import os
import errno


def make_sure_path_exists(path):
    try:
        os.makedirs(path)
    except OSError as exception:
        if exception.errno != errno.EEXIST:
            raise

if( mpistatus.status and   sys.argv[0][-12:]!="sphinx-build"):
    mygpus, MemPerProc, coresperproc, cpusperproc = setCpuSet.setCpuSet()
    ncpus =  coresperproc
    
    def genReduce(token, tipo =  MPI.INT , ope = MPI.MIN ):
        if tipo == MPI.INT:
            dtype = "i"
        if tipo == MPI.FLOAT:
            dtype = "f"

        a    = numpy.array( [token],dtype )
        a_r  = numpy.array( [0]    ,   dtype        )
        comm.Allreduce([a, tipo  ], [a_r, tipo  ], op=ope)
        return a_r[0]

    MemPerProc = genReduce(MemPerProc,MPI.FLOAT,MPI.MIN )

else:
    ncpus=1
    mygpus, MemPerProc, coresperproc, cpusperproc = [], 20, 1,1



def notyetwithquotes(line):
    if "\"" not in line:
        return 1
    if "#" not in line:
        return 0
    
    pos=line.find("\"")
    pos2 = line.find("#")
    if pos<pos2:
        return 0
    else:
        return 1

def treat_par_file(s):
    """
      pythonify the old format (fortran)
      input file given as a string s
    """
    # if myrank==0 :
    #     if Parameters.VERBOSITY>0 :  print( "PYTHONIFYING INPUT FILE")
    s=s.replace("!","#")
    lines=s.split("\n")

    new_lines=""
    for line in lines:
        
        if(  (      line.find("FILE ")>=0    and  (line.find("FILE ") < line.find("="))  )  and notyetwithquotes(line) ):
            if("[" not in line or ( line.find("[")> line.find("#")    )):
                line=add_quotes(line)
        elif (  (line.find("FILE=")>=0    and  (line.find("FILE=") < line.find("="))  )  and notyetwithquotes(line) ):
            if("[" not in line or ( line.find("[")> line.find("#")    )):
                line=add_quotes(line)
        elif( line.find("PREFIX")>=0   and  (line.find("PREFIX") < line.find("="))  and notyetwithquotes(line)   ):
             line=add_quotes(line)
        elif( line.find("POSTFIX")>=0  and (line.find("POSTFIX") < line.find("="))  and notyetwithquotes(line)  ):
             line=add_quotes(line)
        elif( line.find("DS_NAME")>=0  and (line.find("DS_NAME") < line.find("="))  and notyetwithquotes(line)  ):
             line=add_quotes(line)
        elif( line.find("CURRENT_NAME")>=0  and (line.find("CURRENT_NAME") < line.find("="))  and notyetwithquotes(line)  ):
             line=add_quotes(line)
        elif( line[:3] == "IF_" ):
            line=add_quotes(line)
             
        new_lines=new_lines+"\n"+line
    # try to execute the input file
    # if myrank==0 :
    #    if Parameters.VERBOSITY>0 : print( "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 ")
    if myrank==0 : print( "INPUT FILE IS GOOD FOR ME")
    return new_lines


def add_quotes(s):
    pos=s.find("=")
    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


def  StringInfo(Parameters):

    dim1 =  Parameters.END_VOXEL_1 - Parameters.START_VOXEL_1 +1 - 2*Parameters.LT_MARGIN
    dim2 =  Parameters.END_VOXEL_2 - Parameters.START_VOXEL_2 +1 - 2*Parameters.LT_MARGIN
    dim3 =  Parameters.END_VOXEL_3 - Parameters.START_VOXEL_3 +1

    s="! PyHST_SLAVE VOLUME INFO FILE\n"
    s=s+"NUM_X =  %d\n" % dim1
    s=s+"NUM_Y =  %d\n" % dim2
    # s=s+"NUM_Z =  %d\n"
    s=s+"NUM_Z =  %d\n" %dim3

    voxelsize  = Parameters.IMAGE_PIXEL_SIZE_1
    s=s+"voxelSize =  %e\n" % voxelsize

    if sys.byteorder=="big":
        s=s+"BYTEORDER = HIGHBYTEFIRST\n"
    else:
        s=s+"BYTEORDER = LOWBYTEFIRST\n"
        
    s=s+"ValMin =  %14.8e \n" % Parameters.aMin
    s=s+"ValMax =  %14.8e \n" % Parameters.aMax
    s=s+"s1 =  %14.8e \n" % Parameters.saturations[0]
    s=s+"s2 =  %14.8e \n" % Parameters.saturations[1]
    s=s+"S1 =  %14.8e \n" % Parameters.saturations[2]
    s=s+"S2 =  %14.8e \n" % Parameters.saturations[3]
        
    return s

def  StringInfoXML(Parameters):
    nframesmovie=None
    import sys
    if sys.byteorder=="big":
        BYTEORDER = "HIGHBYTEFIRST"
    else:
        BYTEORDER = "LOWBYTEFIRST"
        s=""
    s=s+"<!-- PyHST VOLUME XML FILE -->\n"
    s=s+"<tomodb2>\n"
    s=s+"<reconstruction>\n"
    s=s+"<idAc>%s</idAc>\n" % Parameters.idAc
    s=s+"<listSubVolume>\n"
    s=s+"<subVolume>\n"
    s=s+"<SUBVOLUME_NAME>%s</SUBVOLUME_NAME>\n" % basename(Parameters. OUTPUT_FILE)

    sx =  Parameters.END_VOXEL_1 - Parameters.START_VOXEL_1 +1 - 2*Parameters.LT_MARGIN
    sy =  Parameters.END_VOXEL_2 - Parameters.START_VOXEL_2 +1 - 2*Parameters.LT_MARGIN
    sz =  Parameters.END_VOXEL_3 - Parameters.START_VOXEL_3 +1

    if  nframesmovie is None:
        sz =  Parameters.END_VOXEL_3 - Parameters.START_VOXEL_3 +1
    else:
        sz =  nframesmovie

    x =  Parameters.START_VOXEL_1 
    y =   Parameters.START_VOXEL_2 
    z =   Parameters.START_VOXEL_3 

    s=s+"<SIZEX>%d</SIZEX>\n" % sx 
    s=s+"<SIZEY>%d</SIZEY>\n" % sy 
    # s=s+"<SIZEZ>%d</SIZEZ>\n"  
    s=s+"<SIZEZ>%d</SIZEZ>\n"  % sz 
    
    s=s+"<ORIGINX>%d</ORIGINX>\n" % x 
    s=s+"<ORIGINY>%d</ORIGINY>\n" % y 
    # s=s+"<ORIGINZ>%d</ORIGINZ>\n" 
    s=s+"<ORIGINZ>%d</ORIGINZ>\n" % (z    + Parameters.SOURCE_Y)
    dimrec = sx*sy*sz
    s=s+"<DIM_REC>%ld</DIM_REC>\n"   % dimrec
    voxelsize  = Parameters.IMAGE_PIXEL_SIZE_1
    s=s+"<voxelsize>%e</voxelsize>\n" % voxelsize
    s=s+"<BYTE_ORDER>%s</BYTE_ORDER> \n" % BYTEORDER

    s=s+"<ValMin>%14.8e </ValMin>\n" % Parameters.aMin
    s=s+"<ValMax>%14.8e </ValMax>\n" % Parameters.aMax

    s=s+"<s1>%14.8e</s1> \n" % Parameters.saturations[0]
    s=s+"<s2>%14.8e</s2> \n" % Parameters.saturations[1]
    s=s+"<S1>%14.8e</S1> \n" % Parameters.saturations[2]
    s=s+"<S2>%14.8e</S2> \n" % Parameters.saturations[3]
    
    s=s+"</subVolume>\n"
    s=s+"</listSubVolume>\n"
    s=s+"</reconstruction>\n"
    s=s+"</tomodb2>\n"
    return s


def get_proj_seqnum_list(natonce):
    P=Parameters
    
    if P.DZPERPROJ==0:
        fondo=0
        npjs = P.numpjs
        cima = P.numpjs-1
    else:
        # fondo = -   IMAGE_PIXEL_SIZE_2*1.0/ P.DZPERPROJ 
        fondo = -  int( P.NUM_IMAGE_2*1.0/ P.DZPERPROJ) 
        # cima  = (IMAGE_PIXEL_SIZE_2+natonce)*1.0/ P.DZPERPROJ
        if Parameters.VERBOSITY>2 :
            print( "natonce, P.DZPERPROJ", natonce, P.DZPERPROJ)
        cima  = (natonce)*1.0/ P.DZPERPROJ 
        npjs  = int(cima-fondo)+1

    ppproc = int( ( npjs*1.0/nprocs)+0.999999)
    my_nprocs = min( ppproc,  npjs-myrank*ppproc       )
    return range(fondo+myrank*ppproc , fondo+myrank*ppproc+my_nprocs), npjs,  range(int(fondo), int(cima)+1 )      # npjs teorico

def get_name_edf_proj(n,FILE_PREFIX=None, LENGTH_OF_NUMERICAL_PART=None , NUMBER_LENGTH_VARIES=None,FILE_POSTFIX=None  ):
    P=Parameters
    name=FILE_PREFIX
    nlength = LENGTH_OF_NUMERICAL_PART
    if(NUMBER_LENGTH_VARIES):
        name=name+("%d"%n)+FILE_POSTFIX
    else:
        number = "%d"%n
        if(len(number)< nlength ):
                number = ("0" *  (nlength-len(number)))+ number
        else:
           number=number[:nlength]
        name=name+number+FILE_POSTFIX
    return name

# def get_name_edf_proj(n):
#     P=Parameters
#     name=P.FILE_PREFIX
#     nlength = P.LENGTH_OF_NUMERICAL_PART
#     if(P.NUMBER_LENGTH_VARIES):
#         name=name+("%d"%n)+P.FILE_POSTFIX
#     else:
#         number = "%d"%n
#         if(len(number)< nlength ):
#                 number = ("0" *  (nlength-len(number)))+ number
#         else:
#            number=number[:nlength]
#         name=name+number+P.FILE_POSTFIX
#     return name

def ModulesForFF( P,  ip ,  INTERVALS):
    # if ip<0 or ip> P.numpjs-1:
    if ip<0 or ip> P.NUM_LAST_IMAGE:
        return None

    if type(INTERVALS)==type(1):
        return ip%INTERVALS, INTERVALS, ip//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],i
            elif sum+INTERVALS[i] == ip:
                return 0, INTERVALS[i],i
                
            sum=sum+INTERVALS[i]
        raise Exception( " sum of INTERVALS for FF is smaller than projection index ")

def getsequentialpos_or_add(lista,item ):
    if item not in lista :
        lista.append(item)
        return len(lista)-1
    else:
        pos = lista.index(item)
        return pos



def get_proj_reading(preanalisi=0):
    P=Parameters

    ## i seqnum_list sono relativi allo zero che sara dato da 
    #                                              proj_num_offset_list
    seqnum_list, P.numpjs_span,   tot_proj_num_list = get_proj_seqnum_list( P.NSLICESATONCE )
    if Parameters.VERBOSITY>2 :print( " DEBUG ", " seqnum_list, P.numpjs_span ", seqnum_list[0],"..",seqnum_list[-1] , P.numpjs_span)
    if Parameters.VERBOSITY>2 :print( " DEBUG ",  """    ## i seqnum_list sono relativi allo zero che sara dato da proj_num_offset_list""")
    if P.DZPERPROJ!=0:
        # se elicoidale la lista dei file  sara completa e bisognera accompagnare 
        # la seqnum_list di un array che da lo shift in termini di numero d' immagine :  proj_num_offset_list
        # Si leggeranno quelle ( che sono fra P.NUM_FIRST_IMAGE e P.NUM_LAST_IMAGE )

        # proj_num_offset_list = numpy.array(P.first_slices_2r/P.DZPERPROJ, "i") #   calcolato dopo. in preanalisi non si ha

        if P.inverti_ordine:
            file_proj_num_list = numpy.arange( P.numpjs )*( -P.FILE_INTERVAL)+ P.NUM_LAST_IMAGE   ## completa
        else:
            file_proj_num_list = numpy.arange( P.numpjs )* P.FILE_INTERVAL+ P.NUM_FIRST_IMAGE   ## completa

        if Parameters.VERBOSITY>2 :print( " DEBUG ", " file_proj_num_list  " , file_proj_num_list[0]," ... ",file_proj_num_list[-1] )
    else:
        # proj_num_offset_list = numpy.zeros( len(P.first_slices_2r)  ,"i")  #   calcolato dopo. in preanalisi non si ha

        file_proj_num_list = numpy.array(seqnum_list )* P.FILE_INTERVAL+ P.NUM_FIRST_IMAGE

        


        
    if P.FILE_PREFIX[-3:]==".h5" or  P.FILE_PREFIX[-4:]==".nxs":
        # the file type is determined by the postfix
        projtype = "h5"
        file_list = [P.FILE_PREFIX]*len(file_proj_num_list)
        # h5 : proj_num_list will addresses slice contained in the h5 file
        my_proj_num_list = file_proj_num_list
        allfile_list=[]
    else:
        projtype="edf"
        # edf : proj_num_list addresses filenames  contained in file_list
        edfargs = {"FILE_PREFIX":P.FILE_PREFIX, "LENGTH_OF_NUMERICAL_PART":P.LENGTH_OF_NUMERICAL_PART, 
                "NUMBER_LENGTH_VARIES":P.NUMBER_LENGTH_VARIES,"FILE_POSTFIX": P.FILE_POSTFIX }
        file_list = [ get_name_edf_proj(i,**edfargs)  for i in file_proj_num_list]
        
        allfile_list = [ get_name_edf_proj(i,**edfargs)  for i in range(   P.NUM_LAST_IMAGE+1) ]


        
        if P.MULTIFRAME==0:
            if Parameters.VERBOSITY>2 :print( " NO MULTIFRAME !!!!!!!!!!!!!!!!!!!!!!!!!!! " )
            file_list = [ get_name_edf_proj(i,**edfargs)  for i in file_proj_num_list]
        else:
            if Parameters.VERBOSITY>2 :print( " MULTIFRAME !!!!!!!!!!!!!!!!!!!!!!!!!!! " )
            file_list = [ P.FILE_PREFIX  for i in file_proj_num_list]
        

    ########################################################
    #    proj_num_list will addresses projections pointed by a  list
    if P.DZPERPROJ!=0:
        my_proj_num_list = numpy.array(   seqnum_list     )
    else:
        my_proj_num_list = numpy.arange(len(seqnum_list  ))
        # con elicoidale sara accompagnato da proj_num_offset_list

    if Parameters.VERBOSITY>2 :print( " DEBUG ", " my_proj_num_list  " , my_proj_num_list[0], " ... ", my_proj_num_list[-1])

            
    if P.CORR_PROJ_FILE is not None:

        corr_file_list = [ ( P.CORR_PROJ_FILE % i    ) for i in file_proj_num_list ]

    else:
        corr_file_list = []
        

    fftype="notset"
    if(P.CORRECT_FLATFIELD):
        if( P.FLATFIELD_CHANGING):
            ff_file_list = []
            ff_indexes_coefficients =[]

            
            for i_pro in file_proj_num_list:

                # res = ModulesForFF(P, ((i_pro-P.NUM_FIRST_IMAGE)//P.FILE_INTERVAL),   P.FF_FILE_INTERVAL)
                res = ModulesForFF(P, ((i_pro-P.NUM_FIRST_IMAGE)             ),   P.FF_FILE_INTERVAL)

                if res is  None:
                    ff_indexes_coefficients.append([-1,-1, 0.0 , 0.0 ,-1,-1 ]    )
                    continue

                i_FF, FF_FILE_INTERVAL, posa = res

                # nFFa=(i_pro*P.FILE_INTERVAL-i_FF) - P.NUM_FIRST_IMAGE + P.FF_NUM_FIRST_IMAGE
                # nFFb = (i_pro*P.FILE_INTERVAL-i_FF+1) - 1 - P.NUM_FIRST_IMAGE + P.FF_NUM_FIRST_IMAGE+  FF_FILE_INTERVAL  
                nFFa=(i_pro-i_FF) - P.NUM_FIRST_IMAGE + P.FF_NUM_FIRST_IMAGE
                nFFb = (i_pro-i_FF+1) - 1 - P.NUM_FIRST_IMAGE + P.FF_NUM_FIRST_IMAGE+  FF_FILE_INTERVAL  
                nFFa = min(nFFa,  P.FF_NUM_LAST_IMAGE)
                nFFb = min(nFFb,  P.FF_NUM_LAST_IMAGE)
                if(nFFb>nFFa): posb=posa+1
                else:          posb=posa
                coeffb = i_FF*1.0/FF_FILE_INTERVAL
                coeffa = 1.0-coeffb


                if P.FILE_PREFIX[-3:] == ".h5" or  P.FILE_PREFIX[-4:] == ".nxs":
                    fftype="h5"
                    if isinstance( P.FF_PREFIX, tuple) or isinstance( P.FF_PREFIX, list):
                        ff_file_list = list(P.FF_PREFIX)
                    else:
                        ff_file_list = [P.FF_PREFIX]
                    if posb==posa:
                        ff_indexes_coefficients.append([posa,-1, 1.0 , 0.0 , (i_pro-i_FF) , -1]    )
                    else:
                        ff_indexes_coefficients.append([posa,posb, coeffa , coeffb , (i_pro-i_FF) ,  (i_pro-i_FF) +  FF_FILE_INTERVAL  ]    )
                else:
                    fftype="edf"
                    edfargs = {"FILE_PREFIX":P.FF_PREFIX, "LENGTH_OF_NUMERICAL_PART":P.FF_LENGTH_OF_NUMERICAL_PART, 
                               "NUMBER_LENGTH_VARIES":P.FF_NUMBER_LENGTH_VARIES,"FILE_POSTFIX": P.FF_POSTFIX }
                    nomeA = get_name_edf_proj(nFFa,**edfargs)
                    nomeB = get_name_edf_proj(nFFb,**edfargs)

                    posA = getsequentialpos_or_add(ff_file_list,nomeA )
                    posB = getsequentialpos_or_add(ff_file_list,nomeB )
                    if posB==posA:
                        ff_indexes_coefficients.append([posA,-1, 1.0 , 0.0 , (i_pro-i_FF)   , -1      ]    )
                    else:
                        ff_indexes_coefficients.append([posA,posB, coeffa , coeffb ,(i_pro-i_FF) , (i_pro-i_FF) +  FF_FILE_INTERVAL ]    )

        else:
            ff_file_list = [P.FLATFIELD_FILE]
            fftype = "h5" if P.FILE_PREFIX[-3:]==".h5" else "edf"
            ff_indexes_coefficients =[[0 ,-1 ,  1.0, 0.0,0,-1]]*len(file_proj_num_list)  # *len(my_proj_num_list)
    else:
        ff_file_list = []
        ff_indexes_coefficients =[[-1 ,-1 ,  0.0, 0.0,-1,-1]]*len(file_proj_num_list) #  *len(my_proj_num_list)

## each proch irank will processes proj_mpi_numpjs[iproc] projections
#  starting from Nproj =  proj_mpi_offsets[irank]

    ppproc = int( (P.numpjs_span*1.0/nprocs)+0.999999)
    proj_mpi_offsets=[]
    proj_mpi_numpjs =[]
    for iproc in range(nprocs):
        its_nprocs = min( ppproc,  P.numpjs_span-iproc*ppproc       )
        proj_mpi_numpjs.append(its_nprocs )
        if iproc==0:
            proj_mpi_offsets.append(0)
        else:
            proj_mpi_offsets.append( proj_mpi_offsets[-1]+proj_mpi_numpjs[-2])

    res={}
    res["proj_reading_type"]= projtype            # "h5" or "edf"
    res["ff_reading_type"]= fftype            # "h5" or "edf"

    if projtype=="h5":
     res["proj_h5_dsname"]=   P.PROJ_DS_NAME           # "h5" or "edf"

    if fftype=="h5":
     res["ff_h5_dsname"]=     P.FF_DS_NAME           # "h5" or "edf"

    # h5 : my_proj_num_list will addresses slice contained in the h5 file  proj_file_list[0]
    # edf : my_proj_num_list addresses filenames  contained in proj_file_list
    # sequences follows local process numbering of workable slices
    res["proj_file_list"]   = file_list
    res["allfile_list"]   = allfile_list
    if Parameters.VERBOSITY>2 :print( "DEBUGP ", len(file_list))
    res["corr_file_list"]   = corr_file_list
    res["ff_file_list"]   = ff_file_list

    res["my_proj_num_list"]    = numpy.array(my_proj_num_list, "i")
    res["tot_proj_num_list"]    = numpy.array(tot_proj_num_list, "i")
    res["file_proj_num_list"]    = numpy.array(file_proj_num_list, "i")  # e' usato da multiframe, bisognera controllare che effetto fa

    res["ff_indexes_coefficients"]    = numpy.array(ff_indexes_coefficients, "f")  


    ## each proch irank will processes proj_mpi_numpjs[iproc] projections
    #  starting from Nproj =  proj_mpi_offsets[irank]
    res["proj_mpi_offsets"] = numpy.array(proj_mpi_offsets, "i")
    res["proj_mpi_numpjs"]  = numpy.array(proj_mpi_numpjs, "i")

    if preanalisi:
        return res

    if P.DZPERPROJ!=0:
        if P.inverti_ordine:
            proj_num_offset_list = numpy.array(P.first_slices_2r_nonsplitted/P.DZPERPROJ, "i")
        else:
            proj_num_offset_list = numpy.array(P.first_slices_2r_nonsplitted/P.DZPERPROJ, "i")
    else:
        proj_num_offset_list = numpy.zeros( len(P.first_slices_2r_nonsplitted)  ,"i")


    if Parameters.VERBOSITY>2 :print( " DEBUG " , "  proj_num_offset_list  " , proj_num_offset_list[0], " ... " , proj_num_offset_list[-1])
    if Parameters.VERBOSITY>2 :print( " DEBUG " , "  P.first_slices_2r_nonsplitted  " ,P.first_slices_2r_nonsplitted)
        

    res["proj_num_offset_list"]    = numpy.array(proj_num_offset_list, "i")



    #########################################################################
    ##  in caso elicoidale i margini non sono clippati e tutto e' delegato 
    ##  alla fase di ricostruzione , se considerare o no una slice
    ###
    res["pos_s"]  = P.pos_s
    res["size_s"] = P.size_s
    res["pos_s_"] = P.pos_s_
    res["size_s_"]= P.size_s_
    if Parameters.VERBOSITY>2 :print( " ======================================= ")
    if Parameters.VERBOSITY>2 :print( " pos_s  " , res["pos_s"] )
    if Parameters.VERBOSITY>2 :print( " pos_s_ " , res["pos_s_"])
    if Parameters.VERBOSITY>2 :print( " size_s  " , res["size_s"] )
    if Parameters.VERBOSITY>2 :print( " size_s_ " , res["size_s_"])


    P.nchunks = len(res["pos_s"] )
    # la lunghezz dei 4 sopra est P.nchunks

    slices_mpi_offsets   = []
    slices_mpi_numslices = []

    first_slices_2r=[]
    last_slices_2r = []

    Zdims = [ S[0] for S in  P.size_s ] 

    m_u_pp = []
    m_d_pp = []
    m_u_pp_w = []
    m_d_pp_w = []

    for dims0, mup, mdown, mup_orig, mdown_orig, f2r, l2r in zip(Zdims ,  P.CONICITY_MARGIN_UP  , P.CONICITY_MARGIN_DOWN,
                                P.CONICITY_MARGIN_UP_original  , P.CONICITY_MARGIN_DOWN_original,
                                P.first_slices_2r_nonsplitted,  P.last_slices_2r_nonsplitted   ):

        totdim = dims0
        # dims0 = dims0 - mup - mdown

        dim2r = l2r-f2r +1
        slices_2r_pproc = int( ( dim2r*1.0/nprocs)+0.999999)

        for iproc in range(nprocs):

            its_nslices_2r = min(slices_2r_pproc ,  dim2r-iproc* slices_2r_pproc      )

            if iproc==0:
                first_slices_2r.append(f2r ) 
                last_slices_2r.append(f2r+its_nslices_2r -1 )
            else:
                if its_nslices_2r:
                    first_slices_2r.append(  last_slices_2r[-1]+1 ) 
                else:
                    first_slices_2r.append(  first_slices_2r[-1] ) 
                last_slices_2r.append(   first_slices_2r[-1]   +its_nslices_2r -1 )


            base = P.START_VOXEL_3-1
            if its_nslices_2r:
                its_nslices       =  P.corr_slice_slicedect[last_slices_2r[-1]-base]-P.corr_slice_slicedect [first_slices_2r[-1]-base]+1 
            else:
                its_nslices       = 0


            slices_mpi_numslices.append( its_nslices)

            if iproc==0:
                slices_mpi_offsets.append(0  + mdown)
            else:
                slices_mpi_offsets.append( slices_mpi_offsets[-1]+slices_mpi_numslices[-2])


            m_d_pp.append(  min( mdown_orig,  slices_mpi_offsets[-1]  )   ) 
            m_u_pp.append(  min( mup_orig, totdim  -  slices_mpi_offsets[-1] - its_nslices  )   ) 

            m_d_pp_w.append(   mdown_orig )   # ils sont utilises apres pour la repartiions en sous tranches
            m_u_pp_w.append(  mup_orig  ) 

            

    P.slices_mpi_offsets = numpy.array(slices_mpi_offsets,"i")
    P.slices_mpi_numslices =  numpy.array(slices_mpi_numslices,"i")

    P.first_slices_2r = numpy.array(first_slices_2r,"i")
    P.last_slices_2r =  numpy.array(last_slices_2r,"i")

    P.first_slices =  numpy.array( P.first_slices ,"i")
    P.last_slices  =  numpy.array( P.last_slices  ,"i")

    P.CONICITY_MARGIN_UP    = numpy.array(m_u_pp,"i")
    P.CONICITY_MARGIN_DOWN  = numpy.array(m_d_pp,"i")

    P.CONICITY_MARGIN_UP_wished    = numpy.array(m_u_pp_w,"i")
    P.CONICITY_MARGIN_DOWN_wished  = numpy.array(m_d_pp_w,"i")

    res["slices_mpi_offsets"]   =  P.slices_mpi_offsets 
    res["slices_mpi_numslices"] =  P.slices_mpi_numslices

    return  res


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



s=""
if(sys.argv[0][-12:]!="sphinx-build"):
    filename=sys.argv[1]
    try:
        f=open(filename,"r")
    except :
        print( " problems reading file ", filename)
        raise Exception(  " problems reading file %s " % filename  ) 

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

NO=0
YES=1

def default_FBFACTORS_FUNCTION(x):
    return 1

class Parameters:
    """ 
    The input variables are read from the input file. The input file is written with python syntax.
    The Input file  is interpreted as python after a preprocessing
    The preprocessing was introduced at the time of  HST-> PYHST transition to maintaing compatibility with:

      *  NO/YES meaning 0/1 
      *  items containing FILE in their name can be initialized
         without using the "" which are otherwise necessary for strings 
 
    The input file has to set (some of) the  following variables.
    To setup easily an input file you are suggested to start from one of the examples given in the doc.
    """

    ##INPUT

    FILE_PREFIX=None
    """
    this tells where the projection can be found. The projection name is obtained appending XXXX+FILE_POSTFIX to the prefix, where XXXX is the projection number
    whose format is given by the options below. When the MULTIFRAME options is activated, then only one data file for radiography is used and
    FILE_PREFIX is its name.
    
    The only accepted file format so far is EDF/JPeg2000 and HDF5.
    For HDF5 format, FILE_PREFIX must end with ".h5". In this case
    the radiographies are read from one single 3D dataset whose name 
    is set by the PROJ_DS_NAME  variable.

    For Jpeg200  decompression is activated when JP2EDF_DIR is not None
    In this case the corresponding projections files (ending with .jp2) 
    will be searched at paths given by FILE_PREFIX postpending .jp2 instead of .edf.
    The decompressed edf files will be written into  JP2EDF_DIR
    """

    JP2EDF_DIR =  None
    """
    in case of jp2 to edf decompression, edf files will be written in this directory.
    They will be deleted at the end if JP2EDF_REMOVE is 1.
    """

    JP2EDF_REMOVE = 1
    """
    Remove or not edf decompressed (from jp2) files at the end.
    """

    PROJ_DS_NAME="data"
    """
    the data set name for projections inside the FILE_PREFIX in case of hdf5.
    The asked hdf5 datatype is native float.
    It works if your data file is float
    but the conversion routine uint to float or ushort to float
    are not implemented in hdf5 and the documentation
    says that user routines for conversion can be complemented to HDF5 library by the user  
    this is not yet done...
    """

    LT_KNOWN_REGIONS_FILE = ""
    """
    If different from a empty string, then it must be a vol file 
    conaining float32 data which represents the known data. The unknown data are represented by Nans.
    Its format must correspond to the one of the output file 
    """
    LT_DIAMETER = 0
    LT_COARSE   = 40
    # LT_FINE     = 4
    LT_MAX_STEP = 200
    LT_REBINNING = 1
    


    PROJ_OUTPUT_FILE = "v2p_%04d.edf"
    """
    When STEAM INVERTER=1 use this pattern to create projection file names. In this case the variable PROJ_OUTPUT_FILE
    determines where the projection are written.
    """
    PROJ_OUTPUT_RESTRAINED = 0
    """
    When the projected volume is smaller in heigth than the original projection, then the projection will be calculated 
    when this option is active, with a size corresponding to the projected volume, not the original projection.
    As compensation more  informations will be written in the images headers that can be used to  localise the smaller
    projections inside the bigger ones. 

    """


    
    NUM_FIRST_IMAGE=0
    """
    the number of the first projection. The value may refer to the position in the file list or the position in a 3D stack (MULTIFRAME option)
    """
    NUM_LAST_IMAGE =0
    """
    the number of the  last projection ,  inclusive : it means that if first,last==0,1999 there are 2000 projections 
    """

    IMAGE_PIXEL_SIZE_1 = 10.000000 
    """
    Pixel size horizontally (microns)
    """
    IMAGE_PIXEL_SIZE_2 = 10.000000 
    """
    Pixel size vertically
    """


    

    FILE_POSTFIX= ".edf"
    """
    the postfix to be appended to FILE_PREFIX+number
    """
    LENGTH_OF_NUMERICAL_PART=4
    """
    The fixed lenght (when it is fixed) of numerical part. The spaces  are filled with zeros
    """
    NUMBER_LENGTH_VARIES = 0
    """
    If 1, the numerical part lenght is not fixed.
    """
    MULTIFRAME=0
    """
    *  If 0 (default) projection and flat-field data are to be found in sequences of edf files : one sequence for the projection, one sequence for flat-fields.
    *  If 2 : one edf file for all the projections, one edf file for all the flat-fields.  An edf file is expected, in this case, to contain a sequence of images, with one header per image.
    *  If 1 : one edf file for all the projections, one edf file for all the flat-fields. An edf file is expected, in this case, to contain a 3D volume, consisting in a pile of images.
    
    """
    CORRECT_FLATFIELD = 0          # Divide by flat-field image
    """
    If set to 1, division by the flat-field (FF) is activated . Otherwise the  FF files are not necessary.
    """

    FF_MATCHER_SPAN = 0
    """
    
    """

    
    FLATFIELD_CHANGING = 0        
    """
    * if set to 0 or NO, Then only one FF is read from  FLATFIELD_FILE, and used to divide all the projections.
    * if set to 1, instead, several FF files are used, which are separated by a given number of projection. Division is done after linear interpolation.
    """
    FLATFIELD_FILE = ""
    """
    This variable is used when FLATFIELD_CHANGING==0. It is the name of the EDF file containing the FF
    """
    FF_PREFIX = "refHST"
    """
    This is the prefix used to compose the FF filename. if it ends with .h5 it will be considered hdf5.
    In the latter case set also FF_DS_NAME which is the name of the 3D dataset from which FF will be read.
    """
    FF_DS_NAME = "data"
    """
    the data set name for flat field inside the FILE_PREFIX in case of hdf5.
    The asked hdf5 datatype is native float : see discussion above for FILE_PREFIX
    """

    FF_NUM_FIRST_IMAGE = 0 
    """
    the number of the first FF. The value  refers to the position in the file list. The MULTFRAME option does not apply
    yet to FFs.
    """
    FF_NUM_LAST_IMAGE = 1500  
    """
    The number of the last projection for which we have FF data.
    """
    FF_FILE_INTERVAL = 1500 
    """
    The projection interval between FFs. If FF_NUM_FIRST_IMAGE = 0 , FF_NUM_LAST_IMAGE = 1500, and FF_FILE_INTERVAL = 1500
    we have one FF for the first (number 0) projection and another for the projection number 1500. For the projection in between 
    a linear interpolation is done.
    """
 
    FF_LENGTH_OF_NUMERICAL_PART = 4 
    """
    The fixed lenght (when it is fixed) of numerical part. The spaces  are filled with zeros
    """

    FF_NUMBER_LENGTH_VARIES = NO
    """
    If 1, the numerical part lenght is not fixed.
    """

    FF_POSTFIX = ".edf"
    """
    the postfix to be appended to FILE_PREFIX+number
    """


    SUBTRACT_BACKGROUND=0
    """
    if you want to subtract a background from the projections set this to YES and set the following variable
    """
    BACKGROUND_FILE=""
    """
    If the SUBTRACT_BACKGROUND option is active, set it to a file name , for example dark.edf.
    If it ends by .h5 it is considered as hdf5. Set in this case BACKGROUND_DS_NAME.
    """

    BACKGROUND_DS_NAME = "data"
    """
    the data set name for BACKGROUND  inside the FILE_PREFIX in case of hdf5
    The asked hdf5 datatype is native float : see discussion above for FILE_PREFIX
    """


    REFERENCES_HOTPIXEL_THRESHOLD = 0
    """
    if different than zero, then dark and references will be filtered for hotpixel 
    exceding their 3x3 median by more than such value .
    This option is effective for for darks, while for refs it is effective only on hdf5 (or nxs)
    files  when
    they contains multiple slices for the same reference. It is applied after median
    
    """
    
    CURRENT_NAME = ""
    """
    if this variable is set then the current is read in the header of each radiography 
    for the necessary normalisations.
    """
    

    DO_DISTORTION = 0
    DISTORTION_FILE = None
    """ 
    This option defaults to None. No distortion correction will be made for the detector.
    If instead you pass a prefix for this option, let's say  ::

             DISTORTIONS_FILE = mydist

    then the horizontal distortion will be read from the file  *mydist_H.edf*
    and the vertical one from *mydist_V.edf*

    If only one of the file is present the other will be assmed to be zero everywhere.
    """

    DO_BH_LUT = 0
    BH_LUT_FILE = None
    """ 
    This option defaults to None. No beam hardening correction will be made for the detector.
    If instead you pass a filename for this option, let's say  ::

            BH_LUT_FILE  = bhlut.txt

    then two ascii columns will be read :
         * the first contain the integral of mu
         * The second contains a correction factor: after the logarithm, the absorption will be multiplied by such factor for correction

    """





    
    DOUBLEFFCORRECTION=0
    """
    The double-ff-correction corrects the projection for the residual error that remains
    after FF correction. The origin of this residual error can be the effect of 
    several causes, for example a detector non-linearity.

    If let to zero, no correction is done. If set to 1 then 
    one image is read from the double-ff EDF file  name pyhst_dff_1.23.edf residing in the same directory as the output file.
    If the file is not there, it is generated first and then used for the reconstruction 
    The image must contain the logarithm of the correction factor. The data is
    subtracted from the treated projection.
    """

    DOUBLEFFCORRECTION_ONTHEFLY = 0
    """
    This parameter is similar to DOUBLEFFCORRECTION, but the double flat field is computed on the fly.
    Mind that if DOUBLEFFCORRECTION_ONTHEFLY > 0, then we should have DOUBLEFFCORRECTION = 0.
    In the current version, this parameter does not work for helical tomography.
    If DOUBLEFFCORRECTION_ONTHEFLY is not zero, all the projections are averaged.
    A file "pyhst_dff_0.12.edf" is generated.
    This file is used as a double flat-field correction.

    If FF2_SIGMA > 0, the average is high-pass filtered, where the filter is a complementary Gaussian function of STD FF2_SIGMA.
    """

    FF2_SIGMA = 0.0
    """
    This parameter is relevant for DOUBLEFFCORRECTION_ONTHEFLY > 0.
    If FF2_SIGMA > 0, the average of the projections is high-pass filtered, where the filter is a complementary Gaussian function of STD FF2_SIGMA.
    This filtered average is used as a double flat field instead of the plain average.
    """


    ZEROCLIPVALUE=1.0e-9
    """
    After division by the FF, and before the logarithm, the data are clipped to this minimum
    """
    ONECLIPVALUE =1.0e1
    """
    After division by the FF, and before the logarithm, the data are clipped to this maximum
    """

    OFFSETRADIOETFF=0.0  # non la commento per ora


    FILE_INTERVAL =      1 ## riaggiungi
    """
    The reconstruction is done with all the projection from NUM_FIRST_IMAGE to NUM_LAST_IMAGE with a step given by FILE_INTERVAL
    """

    ##GEOMETRY

    NUM_IMAGE_1 = 1
    """
     Number of pixels horizontally  (take into account binning : smaller if binning > 1 )
    """
    NUM_IMAGE_2 = 0
    """
    Number of pixels vertically    (take into account binning )
    """

    ROTATION_VERTICAL=1
    """
    If horizontal data are rotated at reading time
    Then all parameter as if axis were vertical
    """

    ANGLES_FILE=None
    """
    If set to a filename, the projection angles(degree)  are read from that file
    instead of being created as an equispaced array with step=ANGLE_BETWEEN_PROJECTIONS
    """
    IGNORE_FILE=None
    """
    If set to a filename, projections will be ignored.
    The file contains *index* of projections to be ignored, separated by a white space (tab, space, new line).
    Example: 314 315 316 317
    """


    ANGLE_OFFSET=0
    """
    use this if you want to obtain a rotated reconstructed slice.
    """

    ANGLE_BETWEEN_PROJECTIONS=0
    """
    The angular step(degree) between projections.
    """

    ROTATION_AXIS_POSITION=0
    """
    Axis position in pixels

    """

    DZPERPROJ=0.0
    """
    If this parameter is different from zero, the helicoidal scan reconstruction is activated. A positive DZPERPROJ
    means that the rotation axis is going down. In other words the detector is going up respect to the axis.
Units are pixels per projection.
    """
    DXPERPROJ=0.0
    """
    Set this parameter  is the translation is not parallel to the rotation axis. Units are pixels per projection.
    """


    START_VOXEL_1 = 1 
    """
    X-start of reconstruction volume
    This value and the others define the  reconstructed region
    """
    END_VOXEL_1   = 1 
    """
    X-end (inclusive) of reconstruction volume
    """
    START_VOXEL_2 = 1 
    """
    Y-start of reconstruction volume
    This value and the others define the  reconstructed region
    """
    END_VOXEL_2   = 1 
    """
    Y-end (inclusive) of reconstruction volume
    """
    START_VOXEL_3 = 1
    """
    Z-start of reconstruction volume
    This value and the others define the  reconstructed region
    """
    END_VOXEL_3   = 1
    """
    Z-end (inclusive) of reconstruction volume
    """

    OVERSAMPLING_FACTOR=4
    """
    Oversampling factor : the to-be-backprojected data are oversampled.
    Then nearest neighbour approximation is used. When using GPU
    this parameter has no effect : linear interpolation is done at
    the hardware level
    """

    RECONSTRUCT_FROM_SINOGRAMS=0 # non serve


    NO_SINOGRAM_FILTERING=0
    FBFILTER = 0
    """
      * FBFILTER = -1 : no fourier filtering is done (plain back-projection)
      * FBFILTER =  0 : naif ramp filter  ( default)
      * FBFILTER =  1 : Ramp filter for the discretized case    http://www.mathematica-journal.com/issue/v6i2/article/murrell/murrell.pdf
      * FBFILTER =  2 : A filter used for interferometry : an integration is done. Then the discretized ramp filter is used
      * FBFILTER =  3 : A filter which goes like abs(sin(x*PI)) with x=0 for the 0 frequency, 0.5 for the niquyst
      * FBFILTER >  4 and FBFILTER <= 5  : A filter which goes like |x|-alpha*x*x. X being 1 at Niquyst and alpha=FBFILTER-4
      * FBFILTER = 10 : DFI method implemented by Roman Shkarin
    """

    DATA_IS_FILTERED = 0
    """
      Set to 1 if the projection data is already filtered.
      Beware that iterative methods are much less efficient on already filtered data.
    """

    USE_DFI = 0
    '''
      If set to 1, uses the Direct Fourier Inversion method (Fourier-Slice theorem) instead of Filtered Backprojection.
      The process is much faster but the result might be less accurate.
      WARNING : for iterative methods, if USE_DFI is set to 1, then USE_DFP should also be set to 1 (since the projector and backprojector should be adjoint of eachother)
    '''
    USE_DFP = 0
    '''
      If set to 1, uses the Direct Fourier Projection method (Fourier-Slice theorem) instead of projection.
      The process is faster but the result might be less accurate.
      WARNING : for iterative methods, if USE_DFP is set to 1, then USE_DFI should also be set to 1 (since the projector and backprojector should be adjoint of eachother)

    '''

    FBFACTORS_FUNCTION = [default_FBFACTORS_FUNCTION]
    """
    By default f( x)= 1
    """

    
    DFI_KERNEL_SIZE = 9
    """
    The length of kernel which will be used on the step of interpolation (usually from 7 to 9, use odd values).
    """
    DFI_NOFVALUES   = 2047
    """
    Number of presampled values which will be used to calculate 
    DFI_KERNEL_SIZE kernel coefficients (usually equals to 1024 or 2048).
    """
    DFI_OVERSAMPLING_RATE   = 2
    """
    Impacts on the amount of zeros which will be added in the center of the cyclic shifted sinogram 
    (just after the step of truncation of a sinogram) (better to use 1 or 2). 
    """
    DFI_R2C_MODE   = 0
    """
    Turns on or off usage of Real-To-Complex transformation on the step applying 1D FFT to each projection (YES/NO).
    """


    DO_CCD_FILTER=0
    """ Setting this parameter to 1 activates a filtering on each projection.
    The treatement is applied after : -background subtraction  -division by the flat-field -double FF correction ... if these
    options are activated.
    """
    CCD_FILTER=""
    """ 
       This variable is a string containing the name of the method. The available options are:
          *  "CCD_Filter" :
             a median filter is applied on the 3X3 neighbour of every pixel. If the pixel value 
             ( the higher the more photons) exceed the median value more then the "threshold" parameter
             then the pixel value is replaced by the median value.
    """
    CCD_FILTER_PARA ={}
    """
     This parameter is a python dictionary containing all the necessary parameters :
            * for option "CCD_Filter", only threshold is necessary. The dictionary will then be 
              of the form : {"threshold": 0.040000 }
    """

    DO_RING_FILTER=0
    """
    When this option is on (1) a filtering is applied to the sinogramme trying to remove ring artefacts.
    This options was previously named DO_SINO_FILTER in pyhst1.  The old names for this option are still compatible.
    """
    RING_FILTER=""
    """
       This variable is a string and is meaningful when DO_RING_FILTER=1
         *  RING_FILTER = "RING_Filter"
               in this case the sinogram is averaged over all the projections. A high-band pass filter is applied to the average
               and the result is subtracted from all the projections. The details of the frequencies filtering are specified
               through variable RING_FILTER_PARA
         *  RING_FILTER = "RING_Filter_THRESHOLDED"
               This case is similar to the previous one. But the averaging is done not on the sinogram itself, instead
               we reconstruct first the slice, we threshold the most intense part of the slice, and the thresholded
               part of the slice is projected out from the sinogram before averaging.
         *  RING_FILTER = "RING_Filter_SG"
               A Savitsky-Golay  filtering by Dominique Bernard. But the documentation still needs to be written
    """
    RING_FILTER_PARA={}
    """
         *  A typical setup in case of RING_FILTER = "RING_Filter" ::

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

            in this case one suppose that the sinogram has 1300 pixels. The FFT has therefore 1300 frequencie
            and the ar array must provide all the factors to the pass or not these frequencies.
            In the above case one lets all the frequencies  except the fundamental.

            Other forms of filtering can be customized using numpy arrays ::

               ar = numpy.ones(1300,'f')
               ar[0]=0.0
               xx=numpy.arange(600-2)
               ar[2:600]=1.0/(1+numpy.exp((2-xx)/5.0)   )

            the frequency layout is the same as in the original HST code : ar[0] the fundamental,
            ar[1] the highest frequency. Then two by two the sin and cos with increasing frequency.


         * in the  RING_FILTER = "RING_Filter" case the   RING_FILTER_PARA dictionary must be completed
           with the threshold value ::

               RING_FILTER_PARA = {"FILTER": ar, "threshold":0.15 }

    """

    SINO_FILTER_PARA={}

    DO_AXIS_CORRECTION=0
    """
    When set to 1 ( or YES) this option activates the axis translation correction.
    """
    AXIS_CORRECTION_FILE=""
    """
    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.
    The datas are the axis displacements. Positive value means positive displacement. In which directions? The CCD axis directions.
    Data are given in CCD pixel units.
    """
    OPTIONS= { "axis_to_the_center":1,"avoidhalftomo":0 , "padding":"E" }   
    """
    This is a python dictionary which may contain the following keys :
        *  axis_to_the_center: 1 or 0 
                 If 1 is chosen the reconstructed region is centered over the rotation axis.
        *  avoidhalftomo : 1 -1 or 0
                when set to zero, half-tomo is activated when the following condition is satisfied ::
        
                   ANGLE_BETWEEN_PROJECTIONS*num_projections = 2*PI +/-0.001

                where num_projections is the number of projection which is deduced  from  NUM_FIRST_IMAGE, NUM_LAST_IMAGE
                and FILE_INTERVAL

                When it is set to -1, it is always activated

        *  padding  :  "E" or "0"  (mind the quotes)
                this options concerns the Fourier filtering in the Filtered Back Projection(FBP) method.
                In FBP a convolution kernel is applied ( whose kind is selected through
                the FBFILTER parameter ) to the data. At the border of the data region (detector has a limited size )
                the tails of the kernels goes outside the data region, and there are two options to pad it:

                    * "E":
                       the data is padded extending the border pixel value beyond the border limit.
                    * "0"
                       the data is padded with zeros

                The total extent of the padded data after padding is given by the smallest power of 2 
                which is greater or equal to 2*num_bins, where num_bins is the detector width.
                The obtained number is again multiplied by two to host extra symmetric padding in
                the case of half-tomo reconstruction.

    *Padding for Half-Tomo: How it works*

         When "E" padding is selected, in the half-tomo case  the padding is done :
             *  beyond the close-to-the-axis border using the theta+180 projection (or the projection with the nearest angle
                when theta+180 is missing)
             *  beyond the far-from-the-axis border, extending the projection border value.

    """

    PENTEZONE=10.0  #####################################################################
    """
    Doing Half-tomo, each projection has an overlapping region, around the rotation axis,
    whose rotated duplicate can be found at the theta+180 angle. When doing backprojection
    the core of this region is multiplied by 0.5 to compensate the fact that
    this region is counted twice, if one considers also the duplicate. The periferical part
    keeps its wheight=1, while we introduce through the PENTEZONE parametre 
    two transition regions where the weight goes from 1 to 0.5 on one side of the core region,
    and from 0.5 on the other size. Both transition have width=PENTEZONE
    """

    SUMRULE=0
    """
    if SUMRULE=1 then, after reconstruction, a constant value is added to the slice suche that the sum over
    the slice be equal to the sum over the sinogram
    """

    STRAIGTHEN_SINOS      = 0
    """
    This option can be used when your reconstructed slice has a strong cupping that you want to remove.
    In general this problem occurs in local-tomography. When this option is activated each sinogram slice
    is fitted with a sum of Chebychev polynomia up to order 2, the result of the fit is subtracted from the line.
    """

    BINNING = 1
    """
    This option, when >1, activated the binning of your data. (mind that this is a barely tested option)
    """
    ZEROOFFMASK=1
    """
    After reconstruction, sets to zero the region outside the reconstruction mask. The reconstruction mask
    is the set of pixel which is define as:

          * in the half-tomo case, those pixels which  always have  a projection points which falls on the detector for 
            at-least half of the angles
          * in the standard case those pixels whose projection point always falls on the detector  

    """


    # FOURIER_FILTER=None
    # FOURIER_FILTER_HAS_RAMP = NO

    OUTPUT_SINOGRAMS=0
    """
    if set to 1, no reconstruction is done but a sinogram is outputted per each slice
    in an edf file ( the edf output will countain the sinogram instead than the slice )
    """


    # BICUBIC=0

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

    DO_AXIS_LONGITUDINAL_CORRECTION=0  # si setta da sola

    # PROJECTION_MEDIAN_FILENAME = "projectionmedian.edf"


    # XYCORRECTIONS=0

    """
    bla bla
    """

    # SINOGRAM_STEP=1

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


    MULTI_PAGANIN_PARS = None
    CORR_PROJ_FILE     = None


    DO_PAGANIN=0
    """ The Paganin filtering (for homogenuous objects) is activated setting to 1 this  parameter.

         *Paganin, D., Mayo, S. C., Gureyev, T. E., Miller, P. R. and Wilkins, S. W. (2002), Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object. Journal of Microscopy, 206: 33-40. doi: 10.1046/j.1365-2818.2002.01010.x* 
    """
    PAGANIN_Lmicron  =0
    r""" For an infinitely distant point source and a detector at   :math:`R_2`  distance

        .. math::  L^2 = 4 \pi^2  R_2  \frac{\delta}{\mu} =  \pi \lambda R_2 \frac{ \delta}{ \beta }
          
        where :math:`\mu` is the linear absorption coefficient.
        Note that this definition of  :math:`L^2` differes by a factor :math:`4 \pi^2`  from Paganin one. Such factor
        is compensated by the definition of frequencies (without    :math:`2 \pi`). The expression :math:`\pi \lambda*R_2 \frac{ \delta}{ \beta }`
        is what one can literally read in the octave fasttomo scripts.

        The lenght L  enters the factors which are applied  in Fourier space :: 

         1.0/((1.0+FF*PAGANIN_Lmicron*PAGANIN_Lmicron)*N1*N2)

        where  ::


           FF = (f2*f2)+(f1*f1)[:,None] and 
           f1=fft_frequencies(N1)/PIXSZ_1
           f2=fft_frequencies(N2)/PIXSZ_2 
           PIXSZ_1 = Parameters.IMAGE_PIXEL_SIZE_2        
           PIXSZ_2 = Parameters.IMAGE_PIXEL_SIZE_1

        and fft_frequencies  is the following function ::

           def fft_frequencies(Nx):    
              fx = numpy.zeros(Nx,"f")
              fx[0:(Nx)//2+1] = numpy.arange((Nx)//2+1)
              fx[(Nx)//2+1:Nx] = -(Nx)//2+1 + numpy.arange(Nx-(Nx)//2-1)
              return fx /Nx

        Mind that in the expressions above both  PAGANIN_Lmicron and .IMAGE_PIXEL_SIZE_1/2 are in microns

    """
    PAGANIN_MARGE=10
    """  When applying FFT the considered data have a larger extent than the ROI because the Paganin kernel has an effective lenght.
         PAGANIN_MARGE, is expressed in pixels. It should be of the same order of  PAGANIN_Lmicron expressed in pixels ( Mind that PAGANIN_Lmicron
         is in micron )
    """

    DO_OUTPUT_PAGANIN=0
    """ 
    Set this to 1 to write preprocessed projection after Paganin.
    Each projection will be written to an edf file whose prefix is given by OUTPUT_PAGANIN_FILE input variable.
    No further processing is done when this option is on. To use properly this option the best way
    is selecting just the middle slice for reconstruction and then setting a high PAGANIN_MARGE marge
    such to encompass a whole radiography plus  a margin for the paganin kernel.
    The output is the data after Paganin, before logarithm and before interpolation for longitudinal correction.
    """
    OUTPUT_PAGANIN_FILE=""
    """
    Each projection will be written to an edf file whose prefix is this
    """


    PUS=0
    """ 
    Unsharp is part of the paganin alghorithm. it is activsated only inconjunction
    with Paganin. It consists in an additional filter
    that is applied after the Paganin filter
    For unsharp you must specify a non-zero Sigma  ::

         PUS=YourSigma   <<<<<<<   !!!!! SIGMA is in pixel units
         PUC=YourCoefficient    <<<<<<<<<

    The image will be :: 

        (1+Parameters.PUC)*originalPaganinImage -Parameters.PUC*ConvolutedImage 

    """
    PUC=0
    """ See above description for PUS
    """
    UNSHARP_LOG=0
    """   IN THIS CASE THE CONVOLUTION KERNEL IS THE 
          laplacian of the gaussian ( of PUS sigma ) and not the gaussian itself
          the image will be  ::

            originalPaganinImage +Parameters.PUC*ConvolutedImage         

    """

    TRYGPU=1
    TRYGPUCCDFILTER = 1 
    """
    When GPU is available the CCD median filter for removing hot-spots will be done on GPU.
    This is the default option. 
    However it may be interesting to set it to zero to check if (on CPU with many cores )
    The CPU implementation is faster.
    """

    NEEDGPU=0
    
    # KEEPONGOING=0
    # ONLINEVISU=0
    # MPIWITHMASTER=0

    # EDF2H5_FILE=None
    # MAXMEMORY = 0

    OUTPUT_FILE = "dum.vol"
    """  Choice of the output format
        
            * if OUTPUT_FILE has postfix .vol then the reconstructed volume is written in binary form, as an array  of float32s
              to the specified file, with z being the slowest dimension, x the fastest 
            * if OUTPUT_FILE has postfix .edf, or no postfix then each slice is written separately to an edf file, numbered with the 
              z coordinate. 

         besides this, at the end of the calculation, two info files are written :  OUTPUT_FILE.info and  OUTPUT_FILE.xml,
         both containing basically the same informations : NUM_X, NUM_Y, NUM_Z, voxelSize, BYTEORDER, ValMin, ValMax, s1, s2, S1, S2.
         An histogram, on  (10)^6 
         bins is written to histogram_OUTPUT_FILE.edf  ( in edf format )
 

         s1 and s2 are the minimum and maximum values of the range that excludes 0.001% of the voxels on both sides of the histo.
         S1 and S2, instead, exclude 0.2% on both sides
    """
    AXIS_TO_THE_CENTER = 1
    PADDING = "E"
    AVOIDHALFTOMO=0




    ITERATIVE_CORRECTIONS = 0
    """
      The number of iterative correction loops. If zero no iterative correction (default).
    """
    FISTA=1
    """
      the update of the slice, in the iterative procedure,  is done with FISTA
    """
    PENALTY_TYPE=0
    """
      the type of penalty function :
        * 0 : standard l1 norm
        * 1 : modified penalty function as in [Elad, Zibulevski, L1-L2 optimization in signal and image proc, May 2010, p 82]
    """
    
    OPTIM_ALGORITHM=1
    """
      the algorithm that solves the convex optimization problem :
          * 1 : FISTA (with or without smooth penalty depending on PENALTY_TYPE)
          * 2 : Modified conjugate gradient algo. PENALTY_TYPE MUST be set to 1.
          * 3 : Chambolle-Pock algorithm ( Only with total variation )
          * 5 : Conjugate SubGradient Algorithm (only with dictionary learning)
    """
    
    
    SMOOTH_PENALTY_PARAM=0
    """
      The s parameter for modified penalty function as in [Elad, Zibulevski, L1-L2 optimization in signal and image proc, May 2010, p 82]
    """
    LINE_SEARCH_INIT=1.0
    """
      The initial step for line_search() in the modified Conjugate Gradient algo
    """

    DENOISING_TYPE=1
    """
    The kind of denoising used at each step of iterative correction :
      * DENOISING_TYPE=1   -> TV denoising   (translated to gpu from Emmanuelle Gouillart https://github.com/emmanuelle/tomo-tv/)

        * related parameters : ITERATIVE_CORRECTIONS, DO_PRECONDITION=1, BETA_TV, N_ITERS_DENOISING,  DUAL_GAP_STOP, 
 
      * DENOISING_TYPE=4   -> Dictionary patching with Overlapping Patches using FISTA + L1 norm. 
        Following  http://arxiv.org/abs/1305.1256

        * related parameters : ITERATIVE_CORRECTIONS, BETA_TV,  WEIGHT_OVERLAP, PATCHES_FILE, STEPFORPATCHES

   
      * DENOISING_TYPE=5 and ITERATIVE_CORRECTIONS>0  -> NN-FBP: Reconstruct using a trained NN-FBP network
      
        * related parameters : NNFBP_FILTERS_FILE, NNFBP_NLINEAR
        
      * DENOISING_TYPE=6  and ITERATIVE_CORRECTIONS>0  -> NN-FBP: Create a training set to train a NN-FBP network
      
        * related parameters : NNFBP_TRAINING_PIXELS_PER_SLICE, NNFBP_TRAINING_RECONSTRUCTION_FILE, NNFBP_NLINEAR,NNFBP_FILTERS_FILE ,  NNFBP_TRAINING_USEMASK, NNFBP_TRAINING_MASK_FILE , NNFBP_NHIDDEN_NODES

      * DENOISING_TYPE=8   -> Wavelet regularization

        * related parameters : W_WNAME, W_LEVELS, W_FISTA_PARAM, W_CYCLE_SPIN, W_SWT

    """


#       * DENOISING_TYPE=2   ->  L1 norm of the pixels values.
#
#      * DENOISING_TYPE=3   -> Dictionary patching OMP (WARNING : not really tested )
#        * related parameters : ITERATIVE_CORRECTIONS, ITERATIVE_CORRECTIONS_NOPREC, BETA_TV,  PATCHES_FILE,STEPFORPATCHES
#      * DENOISING_TYPE=30  -> Non local means ( postprocessing )
#
#        * related parameters : ITERATIVE_CORRECTIONS, CALM_ZONE_LEN, NLM_NOISE_GEOMETRIC_RATIO, NLM_NOISE_INITIAL_FACTOR



    VECTORIALITY=1
    """
    Let this always to 1 unless you are doing differential phase tomography.
    In this case the sinogram data must be the horizontal displacement of the beam,
    which is a measure of the phase gradient parallel to the detector horizontal axis.
    PyHST2 reconstructs in this case two images, one obtained by multiplying the data
    by cosinus(angle) , the other by sinus(angle).
    Make sure in this case that ZEROCLIPVALUE is set to a large negative value or your (multiplied) 
    data will be clipped.
    If you use this option with Dictionary learning, provide a vectorial dictionary (key DERIVATIVES: 1
    in the yaml input file for ksvd program which generates the patches )
    """


    FLUO_SINO = 0
    """
    When set to 1, PyHST processes the sinogram as a fluorescence sinogram.
    The statistics is estimated in the iterative reconstruction.
    This only works for the Chambolle-Pock TV solver (DENOISING_TYPE = 1 and OPTIM_ALGORITHM = 3).
    Related parameters: FLUO_ITERS
    """

    FLUO_ITERS = 0
    """
    When FLUO_SINO is set to 1, this determines how many times the TV solver is called.
    The principle is the following. A first "regular" (absorption) reconstruction is performed with the standard TV method.
    The resulting slice is projected, and its square root is taken as the initial statistics.
    Then, FLUO_ITERS iterations of the TV solver are run, each time with an updated statistics.
    """



    WEIGHT_OVERLAP=1.0
    """
    When DENOISING_TYPE=4, the rho of the paper  http://arxiv.org/abs/1305.1256
    """



    ITER_RING_HEIGHT=0.0
    """ Maximum height for rings correction, let it to zero if you want no rings correction.
        or very high value if you want them conpletely free.
    """
    ITER_RING_SIZE=0.0
    """
    (for Dictionary Learning only) estimated width of the rings
    """
    ITER_RING_BETA=0.0
    """ L1 penalisation of rings correction
    """
    LIPSCHITZFACTOR=1.1
    """
       Factor applied to largest eigenvalue(power method)  to estimate Lipschitz factor
    """
    RING_ALPHA = 1.0
    """
       Variable change in rings
    """
    NUMBER_OF_RINGS = 1
    """
       Number of rings to use
    """
    W_LEVELS = 4
    """
      Number of decomposition levels for the Wavelet transform (for DENOISING_TYPE = 8)
    """
    W_CYCLE_SPIN = 1
    """
      Performs a random shift of the image at each iteration. This is a trick to achieve translation invariance, yielding much better visual results.
      When this option is activated, the energy is oscillating (but still globally decreasing).
      Another way of achieving translation invariance is setting W_SWT to 1.
    """
    W_FISTA_PARAM = 4.0
    """
      The revisited FISTA algorithm has the following update equation : y = x + (k-1)/(k+a)*(x-x_old).
      The parameter a > 2  tunes the convergence rate. See https://hal-polytechnique.archives-ouvertes.fr/hal-01060130/document
      Note: this parameter must be greater than 2.
    """
    W_WNAME = "haar"
    """
      Name of the wavelet used for denoising. Availables wavelets are listed at http://wavelets.pybytes.com
    """
    W_SWT = 0
    """
      If activated, use the Stationary Wavelet Transform instead of the (critically sampled) Discrete Wavelet Transform.
      This yields better denoising results, at the expanse of a slower processing.
      When this option is activated, it is not relevant to use W_CYCLE_SPIN.
    """
    W_NORMALIZE = 0
    """
      If set to 1, the wavelet thresholding is normalized accross the scales.
      It usually yields better results, but the regularization parameter has to be increased (approximatively by a factor of 2)
    """
    W_THRESHOLD_APPCOEFFS = 0
    """
      When set to 1, the approximation coefficients of the Wavelet transform are also thresholded.
      Setting this parameter to 1 can be useful at removing even further the noise.
      Default is 0.
    """
    W_THRESHOLD_COUSINS = 0
    """
      When set to 1, the coefficients are thresholded following a "cousinhood" pattern.
      This can be useful to avoid thresholding artefacts when Poisson noise is prominent in the data.
      Default is 0.
    """


    FW_LEVELS = 0
    """
    Wavelets decomposition levels for the Fourier-Wavelet destriping algorithm. 0 means that the de-striper is disabled.
    """

    FW_SIGMA = 1.0
    """
    Gaussian std for the Fourier-Wavelets dampening factor
    """
    FW_WNAME = "db15"
    """
    Wavelet name for the Fourier-Wavelets sinogram de-striping algorithm
    """
    FW_SCALING = 0.0
    """
    If not null, the following process is applied on the sinogram.
    The sum of the sinogram along the columns is computed.
    This "sum line", divided by FW_SCALING, is then subtracted from each sinogram row.
    The Python equivalent code is:
        Rowsum = sino.sum(axis=0)
        sino -= Rowsum/FW_SCALING
    If the Munch filter is enabled (FW_LEVELS > 0), this process is done before the Munch filter.
    """
    FW_FILTERPARAM = 0.0
    """
    If not null, The "row sum" (explained above) is filtered to only keep the high frequencies.
    The high pass filter has the form 1-exp(-nu/FW_FILTERPARAM) where mu = 0, ..., 0.5 is the frequency vector
    normlized to 0.5.
    The Python equivalent code is:
        Rowsum = sino.sum(axis=0)
        Rowsum_filtered = np.fft.irfft(np.fft.rfft(Rowsum)*(1-np.exp(-nu/FW_FILTERPARAM)))
        sino -= Rowsum_filtered/FW_SCALING
    The higher FW_FILTERPARAM, the more low frequencies are kept.
    A good value is 0.01.
    """

    CG_MU = 1e-5
    """
      Parameter tuning the approximation of Total Variation when using the Conjugate Gradient solver.
      The smaller, the closer approximation.
    """

    SF_ITERATIONS = 0
    """
      Number of iterations of the SIRT Filter.
      This option superseeds FB_FILTER one.
      The filter is computed only once with a given number of iterations (and the projection/slice geometry) and can be re-used for any dataset with same geometry.
    """
    SF_FILTER_FILE = "/tmp"
    """
      Directory in which the SIRT filter will be saved.
      The filter is computed only once with a given number of iterations (and the projection/slice geometry) and can be re-used for any dataset with same geometry.
    """
    SF_LAMBDA = 0.0
    """
      Regularization parameter for a Tikhonov (L2 squared) regularization.
    """

    SF_FILTER = numpy.array([1], "f")






    DO_DUMP_HRINGS = 0
    """
       For debugging iterative rings corrections : if 1 dumps at each iteration the corrections on a raw file rings.dat
    """
    INTERVALS_HRINGS = [0,10000000]
    """
       A flat list of intervals where the rings correction are kept. For difficult cases it may be useful to restrain them.
    """


    
    LIPSCHITZ_ITERATIONS=10
    """
     N of iteration in the power method
    """
    DO_PRECONDITION=1
    """
    the preconditioning of the backprojector is turned off by default. Set it
    to zero if you want to perform filtered-backprojections in the iteration loops.
    Mind that BETA_TV will be much different than in the usual case ( much smaller)
    """
    BETA_TV=1.0
    r"""
    * When  DENOISING_TYPE=1, the beta entering in the below formula :

     .. math::
         \frac{1}{2} \left\| Sino - P \cdot Slice \right\|_2^2 + \beta \, TV(Slice)


    whose value is minimized by the iterative procedure

    * When  DENOISING_TYPE=4, the beta of the paper  http://arxiv.org/abs/1305.1256


    """

    # * When  DENOISING_TYPE=2, the beta entering in the below formula :

    #   .. math:: \frac{1}{2} \left\| Sinogramme - P  * Slice  \right\| + \beta ~L_1( Slice) 

    # * When  DENOISING_TYPE=3, the integer part of beta is the number of components used in OMP, 
    #   the fractionary part a in-tolerance to accept components. For example 3.000 will take always 3 components,
    #  3.5 will be selective , 3.9999 being the most selective


    ITER_POS_CONSTRAINT = 0
    """
    Set this value to 1 to enforce a positivity constraint in the reconstructed slice.
    At the moment, this only works for the Chambolle-Pock TV reconstruction (DENOISING_TYPE = 1 and OPTIM_ALGORITHM = 3).
    """


    ESTIMATE_BETA = 0
    """
    When set to 1, PyHST will (coarsely) guess the regularization parameter.
    """

    BETA_L2 = 0.0
    """
    Adds a L2 regularization to the Total Variation, somehow smoothing the result to avoid blocky images.
    Disabled by default.
    """

    N_ITERS_DENOISING = 500
    """
    The number of  iterations when DENOISING_TYPE=1,    for the denoising procedure.
    For DENOISING_TYPE=1 it is the maximum number of iterations, the stop criteria being dual gap.
    """
    DUAL_GAP_STOP = 0.5e-6
    """
    the stop criteria for DENOISING_TYPE=1
    """

    PATCHES_FILE = "patches.h5"
    """
    In the case DENOISING_TYPE=4, the file where square patches are found. It is a  hdf5 file
    containing a 5D dataset named /data.
    The dimesions of this dataset are :

                    * number of patches
                    * number of vectorial components ( 1 for normal thomography, 2 for phase gradient reconstruction)
                    * z dimension of a patch ( 1 if you are doing slice by slice regularisation)
                    * y dimension
                    * x dimension

    The z dimension must be odd ( 1, 3 or more ). if it is n with n>1, n slices will be reconstructed
    and regularised at the same time and the middle one will be taken. 
    The patches can be generated using :
    
          git clone https://gitlab.esrf.Fr/mirone/ksvd
    """

    STEPFORPATCHES=2
    """
      Determines the degree of overlapping between patches. If bigger then 
      patch size then patches are only contiguous. At the other limit,
      when the step is one the overlapping is maximum.
    """
    CALM_ZONE_LEN = 31
    # """" 
    # When  DENOISING_TYPE=30, Non-Local-means denoising is applied. The noise sigma is estimated  by searching 
    # in the image the patch, between all the patches of size CALM_ZONE_LEN, which has the smallest RMSD. 
    # This value, multiplied by NLM_NOISE_INITIAL_FACTOR  is taken as the noise value in the NLM algorithm.
    # Several application of the NLM denoising routine are done, after each application the nois parameter is divided by 
    # NLM_NOISE_GEOMETRIC_RATIO
    # """
    NLM_NOISE_GEOMETRIC_RATIO =2.0
    # """
    # after each application of the filter the noise parameter is divided by this factor, in view of the 
    # next application.
    # """
    NLM_NOISE_INITIAL_FACTOR  =0.5
    # """
    # The RMSD otained from the image is multiplied by this factor
    # """
       
    DETECTOR_DUTY_RATIO        =1.0
    """
    This is the ratio of aquisition time over total time.
    """
    DETECTOR_DUTY_OVERSAMPLING =1
    """
    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.
    """

    CONICITY = 0
    """
    Set to 1 to activate Conicity
    """

    CONICITY_FAN = 0
    """
    Set to 1 to activate Conicity fan geometry ( only horizontal corrections ).
    In this case reconstructions is still slice by slice ( at variance with CONICITY=1
    where the projection of a slices hits sever<al lines ).
    """

    SOURCE_DISTANCE = 0
    """
    Used with CONICITY=1.
    Source-Axis distance in meters
    """
    DETECTOR_DISTANCE = 0
    """
    Used with CONICITY=1.
    DETECTOR-Axis distance in meters
    """
    SOURCE_X        = None
    """
    The X coordinate of the detector pixel where the beam is normal .
    This is a float. The first pixel is 0 ( not 1).
    If you let this variable unset, it will be set equal to  ROTATION_AXIS_POSITION
    """
    SOURCE_Y        = 0
    """
    The Y coordinate of the detector pixel where the beam is normal 
    This is a float. The first pixel is 0 ( not 1).
    Unused if CONICITY_FAN=1
    """


    DECT_PSI  = 0.0
    DECT_TILT = 0.0

    STEAM_INVERTER = 0
    """
    When STEAM INVERTER=1  PyHST goes from volume to projections.
    """



    NSLICESATONCE   =        None
    NPBUNCHES       =        None
    ncpus_expansion =        None
    ALLOWBOTHGPUCPU      = 0
    MAXNTOKENS           = 1
    TRYEDFCONSTANTHEADER = 1
    """
    When enabled (default), PyHST assumes that all the projections have the same header size, if the first and the last do.
    This may lead to a slight gain of performance.
    If the projections do not have the same header size, then set this option to 0.
    If CURRENT_NAME is set, then TRYEDFCONSTANTHEADER is set to zero beacaus all projections need to be read.
    """
    RAWDATA_MEMORY_REUSE = 3#1



    IF_INCLUDE = None
    IF_EXCLUDE = None



    

    JOSEPHNOCLIP=1

    NNFBP_FILTERS_FILE = "neuralfilters.npz"
    """File with trained filter weights, as created during training"""
    NNFBP_NHIDDEN_NODES = 10
    """
    Number of hidden nodes in the neural network
    """
    NNFBP_NLINEAR=4
    """Number of linear steps to take before exponential binning. Higher values will result in higher quality reconstructions, but slower training.
    If DENOISING_TYPE=5, the value for NNFBP_NLINEAR should be equal to the one used when creating the training data with DENOISING_TYPE=6.
    """
    NNFBP_TRAINING_PIXELS_PER_SLICE=1000
    """Number of random pixels to pick per image. Higher values will result in higher quality reconstructions, but slower training. Typical values should result
    in around 10^5 to 10^6 total number of pixels, so if there are 1024 slices, this should be set between 100 and 1000."""
    NNFBP_TRAINING_RECONSTRUCTION_FILE="reconstruction.vol"
    """File with a high quality reconstruction of the object, used during training. The slices of this file should match the slices of the dataset.
    Note that only the .vol extension is supported at the moment, not the .edf extension."""
    NNFBP_TRAINING_USEMASK=0
    """Whether to use an optional mask file (see NNFBP_TRAINING_MASK_FILE)."""
    NNFBP_TRAINING_MASK_FILE="mask.vol"
    """Optional mask file to use during training. The mask file is a binary file of unsigned 8 bit integers, with one value for each pixel of a single slice.
    A value of 0 specifies that this pixel can not be picked during training, and a value>0 specifies that this pixel can be picked during training. The mask can
    be used to train NN-FBP on a specific region of the object, such that it will be able to provide better reconstructions for this region, at the cost of
    reconstruction quality outside the region."""


    ## Parameters for unsharping
    DO_V3D_UNSHARP = 0
    """
    If activated, PyHST2 will read an already reconstructed volume and produce an unsharped one 
    """

    V3D_TIFFREAD = 0
    """
    By default the input volume is a raw data file. If this option is activated, instead, 
    it is a sequence of tiff files
    """


    V3D_TIFFOUTPUT =0
    """
    By default the output volume is written as a raw data file.
    By default this option is set to 0 and no stack of images is written
    to represent the volume.
    When this option is set to 1, the output volume is written
    as a sequence of tiff images.
    When it is set to 2 a sequence of edf images will be written.

    """
    
    
    V3D_CONVZ = 1
    """
    By default unsharping convolution is done on the 3 directions.
    If this option is set to zero, instead, the unsharp is done in the plane only.
    """

    
    V3D_INPUTFILE = ""
    """
    The input filename. When V3D_TIFFREAD=0It can be output.vol as an example.
    When V3D_TIFFREAD=1 or(2) it can be of the like ::

         V3D_INPUTFILE = mytiff%04d.tif

    in this case the numerical index would take 4 digits padded with zeros.

    """
    
    V3D_OUTPUTFILE= ""
    """
    The output filename. when V3D_TIFFOUTPUT =0 it can be output.vol as an example.
    When V3D_TIFFOUTPUT=1 or(2) it can be of the like ::

         V3D_OUTPUTFILE = mytiff%04d.tif

    in this case the numerical index would take 4 digits padded with zeros.

    """
   
    V3D_VOLDIMS = [0,0,0]
    """
    The size of the original volume in the order xyz.
    X being the fastest dimension, Z the slowest one.
    """

    V3D_ROICORNER=[0,0,0]
    """
    The unsharp will be applied to a ROI subvolume whose
    corner the nearest to the origin is given by this parameter.
    """
    
    V3D_ROISIZE=[0,0,0]
    """
    The dimensions of the unsharped ROI in pixel units.
    """
    
    V3D_PUS = 0.0
    """
    
    The convolution SIGMA is in pixel units

    The volume will be :: 

        (1+V3D_PUC)*originalVolume -V3D_PUC*ConvolutedVolume 

    """

    V3D_PUC =0.0
    """

    The volume will be :: 

        (1+V3D_PUC)*originalVolume -V3D_PUC*ConvolutedVolume 


    """


    VERBOSITY = 10
    """
    Level of verbosity
    """

    
    if(sys.argv[0][-12:]!="sphinx-build"):
        print( (" reading " ))
        print( s)
        exec(s)


    DOUBLERUN4DFF=0
    if DOUBLEFFCORRECTION:

        if  not isinstance(DOUBLEFFCORRECTION, string_types):
            dname = os.path.dirname(OUTPUT_FILE)
            dd_path = os.path.join( dname,"pyhst_dff_%.2f.edf"%FF2_SIGMA)
            print( " CONTROLLO ", dd_path)
            if not os.path.exists(dd_path):
                DOUBLEFFCORRECTION = 0
                DOUBLEFFCORRECTION_ONTHEFLY = 1
                DOUBLERUN4DFF = 1
                END_VOXEL_3   = START_VOXEL_3
            else:
                DOUBLEFFCORRECTION = dd_path
        else:
                DOUBLEFFCORRECTION =  DOUBLEFFCORRECTION
                DOUBLEFFCORRECTION_ONTHEFLY = 0
            


    if LT_KNOWN_REGIONS_FILE != "":

        ZEROOFFMASK=0


        
        assert(END_VOXEL_1 - START_VOXEL_1==
               END_VOXEL_2 - START_VOXEL_2  ) 
        assert( OPTIONS["axis_to_the_center"]=='Y' or OPTIONS["axis_to_the_center"]==1 )

        LT_MARGIN = int( (END_VOXEL_1 - START_VOXEL_1+1)*0.414/2 +1.0)

        assert( LT_DIAMETER >= (END_VOXEL_1 - START_VOXEL_1+1)*1.5  ) 
        # assert(LT_COARSE>=10)

        END_VOXEL_1 = END_VOXEL_1 + LT_MARGIN
        END_VOXEL_2 = END_VOXEL_2 + LT_MARGIN
        START_VOXEL_1 = START_VOXEL_1 - LT_MARGIN
        START_VOXEL_2 = START_VOXEL_2 - LT_MARGIN

        LT_MAX_STEP_REBIN = LT_MAX_STEP//LT_REBINNING


        
    else:
        LT_MARGIN = 0




P=Parameters
P.ncpus=ncpus


if len(P.CURRENT_NAME) :
    P.TRYEDFCONSTANTHEADER = 0


if  P.BH_LUT_FILE is not None:
    P.DO_BH_LUT=1
    tmp.numpy.loadtxt( P.BH_LUT_FILE).astpe("f")
    P.BH_LUT_U = numpy.array(tmp[:,0])
    P.BH_LUT_F = numpy.array(tmp[:,1])
else:
    P.DO_BH_LUT=0
    P.BH_LUT_U = numpy.zeros([2],"f")
    P.BH_LUT_F = numpy.zeros([2],"f")
    

if  P.DISTORTION_FILE is not None:

    fnH = P.DISTORTION_FILE+"_H.edf"  
    fnV = P.DISTORTION_FILE+"_V.edf"
    
    shape_distortion = None

    DIST_H = None
    DIST_V = None
    
    if os.path.exists(fnH):
        # DIST_H = numpy.loadtxt(fnH ).astype("f")
        
        tmpf   = EdfFile.EdfFile(  fnH    )
        DIST_H = tmpf.GetData(0,DataType="FloatValue")
        
        shape_distortion = DIST_H.shape
        P.ROTATION_AXIS_POSITION = P.ROTATION_AXIS_POSITION - (  DIST_H[:, int(P.ROTATION_AXIS_POSITION)] ).mean()
        
    if os.path.exists(fnV):
        
        # DIST_V = numpy.loadtxt(fnV ).astype("f")

        tmpf   = EdfFile.EdfFile(   fnV   )
        DIST_H = tmpf.GetData(0,DataType="FloatValue")
        
        
        if shape_distortion is  None:
            shape_distortion = DIST_V.shape
        else:
            assert( shape_distortion == DIST_V.shape   )


    P.DIST_H = numpy.zeros([ P.NUM_IMAGE_2, P.NUM_IMAGE_1 ],"f")
    P.DIST_V = numpy.zeros([ P.NUM_IMAGE_2, P.NUM_IMAGE_1 ],"f")

    if DIST_H is not None:
        P.DIST_H[:] = DIST_H
    if DIST_V is not None:
        P.DIST_V[:] = DIST_V
    
        
    P.DO_DISTORTION = 1

    if P.BINNING >1 :
        P.DIST_H[:]    =  P.DIST_H /      P.BINNING
        P.DIST_V[:]    =  P.DIST_V /      P.BINNING

    P.DIST_H_max = int( abs(P.DIST_H).max() +1 )
    P.DIST_V_max = int( abs(P.DIST_V).max() +1 )

        
else:
    P.DIST_H =numpy.zeros([2,2],"f")
    P.DIST_V =numpy.zeros([2,2],"f")
    P.DO_DISTORTION = 0
    P.DIST_H_max = 0
    P.DIST_V_max = 0

        
if Parameters.ITER_RING_HEIGHT==0:
    Parameters.ALPHA_RINGS=0
 
Parameters.INTERVALS_HRINGS = numpy.array(Parameters.INTERVALS_HRINGS).astype("i")

num_bins = Parameters.NUM_IMAGE_1

two_power = (int)(math.log((2. * num_bins - 1)) / math.log(2.0) + 0.9999)
dim_fft = 1
for i in range(two_power):
    dim_fft*=2
dim_fft*=2

Parameters.FBFACTORS=numpy.zeros(  dim_fft//2  +1 ,"f")
xs = (numpy.arange( dim_fft//2+1 ).astype("f")/(  dim_fft//2 +1  ))
for i in range( dim_fft//2 +1):
    #print( Parameters.FBFACTORS_FUNCTION)
    Parameters.FBFACTORS[i] = Parameters.FBFACTORS_FUNCTION[0](   xs[i] )


# print( " Parameters.FBFACTORS  " ,   Parameters.FBFACTORS)


if Parameters.CONICITY:
    Parameters.SOURCE_Y = int(Parameters.SOURCE_Y)
else:
    Parameters.SOURCE_Y = 0    


Parameters.START_VOXEL_3 =    Parameters.START_VOXEL_3-Parameters.SOURCE_Y
Parameters.END_VOXEL_3 =    Parameters.END_VOXEL_3-Parameters.SOURCE_Y

if Parameters.DO_PRECONDITION is None:
    Parameters.DO_PRECONDITION = 0


if(Parameters.SOURCE_X  is None):
    Parameters.SOURCE_X  = Parameters.ROTATION_AXIS_POSITION


    
if Parameters.NO_SINOGRAM_FILTERING:
    Parameters.FBFILTER=-1

if Parameters.NPBUNCHES is None and mpistatus.status:
    Parameters.NPBUNCHES =    int( genReduce(coresperproc, MPI.INT,MPI.MIN))

if Parameters.TRYGPU==1 and mpistatus.status :
    minNgpus=genReduce(len(mygpus),MPI.INT)
    if minNgpus==0 :
        Parameters.TRYGPU=0
        
if Parameters.ncpus_expansion is None:
    if Parameters.TRYGPU==1 :
        Parameters.ncpus_expansion=2
    else:
        Parameters.ncpus_expansion=1

if Parameters.VERBOSITY>0 : print( " USING  ", Parameters.NPBUNCHES, " bunches " )
## print( " Parameters.TRYGPU ", Parameters.TRYGPU )
## print( " Parameters.ncpus_expansion ", Parameters.ncpus_expansion)
   
P=Parameters

if P.DZPERPROJ<0:
    P.ANGLE_BETWEEN_PROJECTIONS  =  -P.ANGLE_BETWEEN_PROJECTIONS
    P.DZPERPROJ=abs(P.DZPERPROJ)
    P.inverti_ordine = 1
    myprojections = numpy.arange( P.NUM_LAST_IMAGE, P.NUM_FIRST_IMAGE-1, -P.FILE_INTERVAL)
    P.numpjs = len(myprojections)
else:
    myprojections = numpy.arange( P.NUM_FIRST_IMAGE, P.NUM_LAST_IMAGE+1, P.FILE_INTERVAL)
    P.numpjs = len(myprojections)

    P.inverti_ordine = 0

  
P.axis_corrections=numpy.zeros([P.numpjs ],"f")
P.axis_correctionsL = numpy.zeros([P.numpjs ],"f")

P.PIECE_MARGE=1+P.DIST_V_max
check_axis_corrections()


def memory_estimator( P, NSLICESATONCE, n_projs_max, n_ff_max, PIECE_MARGE_touse, MAX_NSLICESATONCE, c_margin, RAWDATA_MEMORY_REUSE=1)    :
    """ This is another test for checking sphinx ability to catch docstrings from functions  into documentation
    """

    sum=0
    sizeV = NSLICESATONCE+2.0*PIECE_MARGE_touse+2*c_margin
    if (P.DOUBLEFFCORRECTION_ONTHEFLY):
        sizeV = P.NUM_IMAGE_2

    # if sizeV>   MAX_NSLICESATONCE : sizeV=MAX_NSLICESATONCE
    if sizeV>  P.NUM_IMAGE_2  : sizeV=P.NUM_IMAGE_2
    ## raw
    sum=sum+ P.NUM_IMAGE_1*n_projs_max*( sizeV  )*1.0/RAWDATA_MEMORY_REUSE
    ## ff
    sum=sum+ P.NUM_IMAGE_1*n_ff_max*( sizeV  )
    ## treated data
    sum=sum+ P.NUM_IMAGE_1*n_projs_max* (NSLICESATONCE + 2*c_margin)
    ## transposed data
    sum=sum+ (P.NUM_IMAGE_1*n_projs_max* (NSLICESATONCE + 2*c_margin))*4
    ## 2 because there are transmit and receive buffers to get it
    
    ## --------------------- C -------------------------------
    ##   request.data

    sum=sum+ P.NUM_IMAGE_1*P.numpjs_span  * (NSLICESATONCE  *1.0/nprocs+2*c_margin )*2

    ## accounting for extra buffers 
    sum = sum + P.NUM_IMAGE_2*P.NUM_IMAGE_1*8* P.NPBUNCHES
    sum = sum + P.numpjs_span*P.NUM_IMAGE_1*8* P.NPBUNCHES


    # sum = sum + P.NUM_IMAGE_2*P.NUM_IMAGE_1*sizeV*4
    # sum = sum + P.NUM_IMAGE_1*( P.numpjs_span  )*2* P.NPBUNCHES

    
    if P.CONICITY:
        sum=sum+ (1+P.END_VOXEL_1-P.START_VOXEL_1)*(1+P.END_VOXEL_1-P.START_VOXEL_1)*(NSLICESATONCE*1.0/nprocs+4)
        sum=sum+ P.NUM_IMAGE_1*P.numpjs_span*(2*c_margin*P.NPBUNCHES+NSLICESATONCE*1.0/nprocs+4)
        # sum=sum+ P.NUM_IMAGE_1*P.numpjs*(NSLICESATONCE*1.0/nprocs+4)

    return sum*4*P.MAXNTOKENS







if not Parameters.DO_V3D_UNSHARP  and Parameters.NSLICESATONCE is None and (sys.argv[0][-12:]!="sphinx-build") :


    if P.CONICITY:
        f_tmp  = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE))    # ??? /P.IMAGE_PIXEL_SIZE_2
        VOXEL_SIZE = P.IMAGE_PIXEL_SIZE_1/f_tmp
        Fact = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE))*VOXEL_SIZE*1.0/P.IMAGE_PIXEL_SIZE_2
        
        DFa  = +(abs(P.DECT_PSI)+abs(P.DECT_TILT) )*3.2/180.0*2
        
        DR1  = (P.END_VOXEL_1 - P.START_VOXEL_1)/2.0
        DR2  = (P.END_VOXEL_2 - P.START_VOXEL_2)/2.0
        RADIUS = math.sqrt( DR1*DR1+DR2*DR2   ) 
        
        alpha = (max (abs(P.START_VOXEL_3),abs(P.END_VOXEL_3) ) *VOXEL_SIZE/1.0e6)/(P.SOURCE_DISTANCE)
        
        C_MARGIN = int(  RADIUS*(
                         alpha*(Fact)+DFa )
                         +0.01)+1
        
        if Parameters.VERBOSITY>1 : print( " C_MARGIN " , C_MARGIN)
    else:
        C_MARGIN=0




    if not Parameters.DO_V3D_UNSHARP  and   P.DZPERPROJ==0 and mpistatus.status  :    
        read_infos   = get_proj_reading(preanalisi=1)
        n_projs_max  =  genReduce(   read_infos["proj_mpi_numpjs"][myrank]     ,MPI.INT, MPI.MAX)
        if(P.CORRECT_FLATFIELD):
            nmax=0
            nmin=1000000
            for  l in read_infos["ff_indexes_coefficients"]:
                
                if(l[0]>-1):
                    nmax=max(nmax,l[0])
                    nmin=min(nmin,l[0])
                if(l[1]>-1):
                    nmax=max(nmax,l[1])
                    nmin=min(nmin,l[1])
            n_ff_max =      genReduce(   nmax+1- nmin    ,MPI.INT, MPI.MAX)
        else:
            n_ff_max=0
    else:
        n_ff_max=0
        
        
    MAX_NSLICESATONCE=P.END_VOXEL_3 +1-P.START_VOXEL_3
    NSLICESATONCE=MAX_NSLICESATONCE

    PIECE_MARGE_touse=  P.PIECE_MARGE
    if P.DO_PAGANIN:
        PIECE_MARGE_touse = PIECE_MARGE_touse + P.PAGANIN_MARGE
    PIECE_MARGE_touse=int(PIECE_MARGE_touse+0.5)

    if (P.DOUBLEFFCORRECTION_ONTHEFLY):
        PIECE_MARGE_touse = P.NUM_IMAGE_2

    
    RAWDATA_MEMORY_REUSE=Parameters.RAWDATA_MEMORY_REUSE
    max_reuse=10
    if Parameters.DO_OUTPUT_PAGANIN:
        assert(P.DZPERPROJ==0)
        max_reuse=n_projs_max
        
    while RAWDATA_MEMORY_REUSE<max_reuse  and  mpistatus.status:
        NSLICESATONCE=MAX_NSLICESATONCE
        while NSLICESATONCE>0:

            if P.DZPERPROJ!=0.0 :
                P.NSLICESATONCE = NSLICESATONCE
                read_infos   = get_proj_reading(preanalisi=1)
                n_projs_max  =  genReduce(   read_infos["proj_mpi_numpjs"][myrank]     ,MPI.INT, MPI.MAX)

            req = memory_estimator(P,NSLICESATONCE,n_projs_max, n_ff_max, PIECE_MARGE_touse, MAX_NSLICESATONCE, C_MARGIN , RAWDATA_MEMORY_REUSE)
            req1 = memory_estimator(P, NSLICESATONCE,n_projs_max, n_ff_max, PIECE_MARGE_touse, MAX_NSLICESATONCE, C_MARGIN , RAWDATA_MEMORY_REUSE+1)
            if Parameters.VERBOSITY>2 : print(  NSLICESATONCE , req, MemPerProc,coresperproc , Parameters.VERBOSITY )# too verbose ?
            if( req <0.75*  MemPerProc*coresperproc     ):
                break
            NSLICESATONCE -=1
        if Parameters.VERBOSITY>0 :  print( " with RAWDATA_MEMORY_REUSE = " , RAWDATA_MEMORY_REUSE , "  NSLICESATONCE  " , NSLICESATONCE)
        if NSLICESATONCE<MAX_NSLICESATONCE and NSLICESATONCE< PIECE_MARGE_touse//2:
            RAWDATA_MEMORY_REUSE +=1
            if Parameters.VERBOSITY>0 :  print( " increased RAWDATA_MEMORY_REUSE" , RAWDATA_MEMORY_REUSE)
        else:
            break
    else:
        if NSLICESATONCE==0:
            raise Exception("PROBLEM: NSLICESATONCE HAS BEEN CALCULATED TO 0 ")
        else:
            print( "WARNING NSLICESATONCE HAS BEEN CALCULATED TO A SMALL VALUE %d COMAPRED TO PIECE_MARGE_touse %d" %( NSLICESATONCE,PIECE_MARGE_touse  ) )


    Parameters.NSLICESATONCE = NSLICESATONCE
    Parameters.RAWDATA_MEMORY_REUSE =  RAWDATA_MEMORY_REUSE


if Parameters.VERBOSITY>0 : print( " Parameters.NSLICESATONCE  " , Parameters.NSLICESATONCE)

if hasattr( Parameters , "DO_SINO_FILTER"   ):
    Parameters.DO_RING_FILTER = Parameters.DO_SINO_FILTER
    if Parameters.DO_RING_FILTER:
        if hasattr(Parameters, "SINO_FILTER"    ):
            Parameters.RING_FILTER="RING_Filter"
            ar = Parameters.SINO_FILTER_PARA["FILTER"]
            Parameters.RING_FILTER_PARA={"FILTER":ar}
         
outputfile = Parameters.OUTPUT_FILE



# if ".vol" in outputfile:
#     outputfile=outputfile[:-4]+"_%04d_%04d"+".vol"
# else:
#     outputfile=outputfile[:-4]+"_%04d_%04d"

    
if "padding" in Parameters.OPTIONS:
    Parameters.PADDING = Parameters.OPTIONS["padding"]
    
if "axis_to_the_center" in Parameters.OPTIONS:
    if  Parameters.OPTIONS["axis_to_the_center"]=='Y':
        Parameters.AXIS_TO_THE_CENTER = 1
    else:
        Parameters.AXIS_TO_THE_CENTER = 0

if "avoidhalftomo" in Parameters.OPTIONS:
    
    if  Parameters.OPTIONS["avoidhalftomo"] in ['Y',1,'1']:
        Parameters.AVOIDHALFTOMO = 1
        if Parameters.VERBOSITY>0 : print( " AVOID  HALF ")
    elif Parameters.OPTIONS["avoidhalftomo"] in ['F',-1,'-1']:
        Parameters.AVOIDHALFTOMO = -1
    else:
        Parameters.AVOIDHALFTOMO = 0
        if Parameters.VERBOSITY>0 : print( " ALLOW HALF " )






if Parameters.PUS!=0:
    Parameters.PUC = Parameters.PUC*max( Parameters.PUS *Parameters.PUS  ,1)  
    Parameters.PUS = Parameters.PUS/2

DEG2RAD=math.pi/180.0

Parameters.ANGLE_BETWEEN_PROJECTIONS = Parameters.ANGLE_BETWEEN_PROJECTIONS *DEG2RAD
Parameters.ANGLE_OFFSET= Parameters.ANGLE_OFFSET*DEG2RAD

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

Parameters.angles=angles




if( Parameters.IGNORE_FILE is not None  ):
    ignore_ = open( Parameters.IGNORE_FILE).read().split( )
    ignore=[]
    n_ignore = 0
    for ign in ignore_:
        try:
            ignore.append(   string.atof(ign)     )
        except:
            break
    n_ignore = len(ignore)
    ignore=(numpy.array(ignore)).astype("i")
else:
    ignore=0
    n_ignore = 0

Parameters.ignore=ignore
Parameters.n_ignore = n_ignore




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 Postprocess(Cspace):

    ####################################################################
    # cherche le fichier xml contenant le tag idAc
    idAc = "N_A_"

    if  not Parameters.RECONSTRUCT_FROM_SINOGRAMS:
        toc= Parameters.FILE_PREFIX
    else:
        toc=  Parameters.SINOGRAM_PREFIX

    file_xml_idac = toc+"_idAc.xml"
    try:
        idAc = simplereadxml(file_xml_idac,"idAc")
    except:
        pass
    Parameters.idAc = idAc
    ################################################

    aMin = Cspace.get_aMin()
    aMax = Cspace.get_aMax()

    val=numpy.array([aMin ])
    res=numpy.array([aMin])
    comm.Allreduce([val, MPI.DOUBLE], [res, MPI.DOUBLE], op=MPI.MIN)
    aMin = res[0]

    val=numpy.array([aMax ])
    res=numpy.array([aMax ])
    comm.Allreduce([val, MPI.DOUBLE], [res, MPI.DOUBLE], op=MPI.MAX)
    aMax = res[0]

    Parameters.aMin = aMin
    Parameters.aMax = aMax
    # print( " getSaturations " )

    sat1,sat2 , Sat1, Sat2= Cspace.getSaturations(aMin,aMax)

    # print( sat1,sat2 , Sat1, Sat2)

    # print( " strings", myrank )


    Parameters.saturations =  (sat1,sat2 , Sat1, Sat2)
    Parameters.INFOOUTPUT = StringInfo   (Parameters)
    Parameters. XMLOUTPUT = StringInfoXML(Parameters)
    
    if myrank==0:
        if Parameters.STEAM_INVERTER==0 :
            print( " scrivo infos " )
            f=open("%s.info"%Parameters.OUTPUT_FILE_ORIG,"w")
            f.write(Parameters.INFOOUTPUT)
            f=open("%s.xml" %Parameters.OUTPUT_FILE_ORIG,"w")
            f.write(Parameters. XMLOUTPUT)
        
        if Parameters.DENOISING_TYPE==6 :
            print( " NEURAL NETWORK TRAINING : ")
            from . import NNFBPTrain
            import glob
            # print( 'Usage:',sys.argv[0],'NFilters ParentPathToTrainingEDFs ParentPathToValidationEDFs outputNPZFile')
            
            edfT = glob.glob('./*_nnfbp_*.edf')
            # edfV = glob.glob(sys.argv[3] + './*_nnfbp_*.edf')
            print( " LERNING FROM files found in current working directory : " , edfT)
            edfV = edfT
            
            n = NNFBPTrain.Network( Parameters.NNFBP_NHIDDEN_NODES ,edfT,edfV)
            n.train()
            n.saveToDisk(   Parameters.NNFBP_FILTERS_FILE  )

   



# if ".vol" in outputfile:
#     outputfile=outputfile[:-4]+"_%04d_%04d"+".vol"
# else:
#     outputfile=outputfile[:-4]+"_%04d_%04d"

if  (sys.argv[0][-12:]!="sphinx-build")   :
    
    Parameters.OUTPUT_FILE_ORIG = outputfile

    histo_outputfile=outputfile

    pos = histo_outputfile.rfind("/")
    
    namepart_histo_outputfile = histo_outputfile[(pos+1):]
    namepart_histo_outputfile=string.replace(namepart_histo_outputfile,".vol",".edf")
    Parameters.OUTPUT_FILE_HISTOGRAM = histo_outputfile[:(pos+1)]+"histogram_"+namepart_histo_outputfile


    Parameters.OUTPUT_FILE_HISTOGRAM = histo_outputfile[:(pos+1)]+"histogram_"+namepart_histo_outputfile

    if Parameters.VECTORIALITY==1:
        vformat = ""
    else:
        vformat =  "_%s_"


    if ".vol" in outputfile  and not (P.OUTPUT_SINOGRAMS or (P.ITERATIVE_CORRECTIONS>0 and P.DENOISING_TYPE==6)):
        outputfile=outputfile + vformat
        P.EDFOUTPUT=0
        if myrank==0  and P.STEAM_INVERTER==0    :
            dimz=P.END_VOXEL_3 +1-P.START_VOXEL_3
            dimy=P.END_VOXEL_2 +1-P.START_VOXEL_2
            dimx=P.END_VOXEL_1 +1-P.START_VOXEL_1
            size = dimz*dimy*dimx
            f=open(outputfile,"w")
            f.seek(size*4-1)
            f.write(" ")
            f.close()
    else:
        P.EDFOUTPUT=1
        if not P.OUTPUT_SINOGRAMS:
            outputfile=outputfile+vformat+"_%04d"+".edf"
        else:
            outputfile=outputfile+"_sinogram_%04d"+".edf"
    Parameters.OUTPUT_FILE=outputfile

    if Parameters.ITERATIVE_CORRECTIONS>0 and Parameters.DENOISING_TYPE==5:
        Parameters.neural_params = {}
        filt_file = numpy.load(Parameters.NNFBP_FILTERS_FILE)
        #filters = numpy.fft.ifftshift(filt_file['filters'],[1])
        #fftfilters = numpy.zeros((filters.shape[0],8192))
        #fftfilters[:,0:1+filters.shape[1]//2]=filters[:,0:1+filters.shape[1]//2]
        #fftfilters[:,8192-filters.shape[1]/2:8192]=filters[:,filters.shape[1]-filters.shape[1]//2:filters.shape[1]]
        #Parameters.neural_params['filters'] = numpy.real(numpy.fft.fft(fftfilters,axis=1)).astype(numpy.float32)
        Parameters.neural_params['filters'] = filt_file['filters'].astype(numpy.float32)
#        import pylab as p
#        p.plot(numpy.fft.ifft(Parameters.neural_params['filters'][0]))
#        p.show()
        Parameters.neural_params['offsets'] = filt_file['offsets'].astype(numpy.float32)
        Parameters.neural_params['weights'] = filt_file['weights'].astype(numpy.float32)
        Parameters.neural_params['minIn'] = float(filt_file['minIn'])
        Parameters.neural_params['maxIn'] = float(filt_file['maxIn'])
    if Parameters.ITERATIVE_CORRECTIONS>0 and Parameters.DENOISING_TYPE==6:
        if os.path.exists(Parameters.NNFBP_TRAINING_RECONSTRUCTION_FILE)==False:
            raise IOError("NNFBP_TRAINING_RECONSTRUCTION_FILE not found: " + Parameters.NNFBP_TRAINING_RECONSTRUCTION_FILE)
        if Parameters.NNFBP_TRAINING_USEMASK==1:
            if os.path.exists(Parameters.NNFBP_TRAINING_MASK_FILE)==False:
                raise IOError("NNFBP_TRAINING_MASK_FILE not found: " + Parameters.NNFBP_TRAINING_MASK_FILE)



class derivedParameters:
    pass


def  OverlappingLogic (pos_edf,  size_edf, axis, last_slice  , marge  ):
        pos_edf_= copy.copy(pos_edf)
        if(pos_edf_[axis]>=marge  or P.DZPERPROJ!=0):
          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)    or   P.DZPERPROJ!=0):
           dsize=dsize+marge
        else:
           dsize=dsize+  last_slice - ((size_edf[axis]-1)+pos_edf[axis])
        if Parameters.VERBOSITY>2 :  print( " qui dsize est " , dsize)
        if Parameters.VERBOSITY>2 : print( " last_slice " , last_slice)
        if Parameters.VERBOSITY>2 : print( " size_edf[axis]" , size_edf[axis])
        if Parameters.VERBOSITY>2 : print( " pos_edf[axis]" , )
        size_edf_=copy.copy(size_edf)
        size_edf_[axis]= size_edf[axis]+dsize
        return pos_edf_, size_edf_


def retrieve_patches(namelist):
    if type(namelist)!=type([]):
        namelist=[namelist]
    res=[]
    for name in namelist:
        h5 = h5py.File(name,"r")
        patches = numpy.array(h5["data"][:],"f")
        res.append(patches)
    patches= numpy.concatenate(res)
    return patches

def readjust_marges(pos, size):
    P.PIECE_MARGE_touse=  P.PIECE_MARGE 
    if P.DO_PAGANIN:
        P.PIECE_MARGE_touse = P.PIECE_MARGE_touse + P.PAGANIN_MARGE
    P.PIECE_MARGE_touse=int(P.PIECE_MARGE_touse+0.5)

    if (P.DOUBLEFFCORRECTION_ONTHEFLY):
        P.PIECE_MARGE_touse = P.NUM_IMAGE_2
        print(("Warning: requested to compute the average of all projections (DOUBLEFFCORRECTION_ONTHEFLY = 1), I am loading all the projections !"))

    pos_, size_ = OverlappingLogic (
        pos,  size, 0, P.NUM_IMAGE_2-1 ,
        P.PIECE_MARGE_touse )
    return pos_, size_




if  (sys.argv[0][-12:]!="sphinx-build")  and not P.DO_V3D_UNSHARP  :


    P.first_slices_2r_nonsplitted = numpy.arange(  P.START_VOXEL_3 -1,
                                    P.END_VOXEL_3-1 +1,
                                    P.NSLICESATONCE)

    P.last_slices_2r_nonsplitted  = numpy.minimum(P.first_slices_2r_nonsplitted + P.NSLICESATONCE -1,
                                            P.END_VOXEL_3-1)

    P.patch_ep = 0
    if P.STEAM_INVERTER==0 and mpistatus.status:
        if P.ITERATIVE_CORRECTIONS>0 and P.DENOISING_TYPE in [4]:
                

            patches = retrieve_patches(  P.PATCHES_FILE   ) 


            assert(len(patches.shape) in [4,5])

            if(len(patches.shape)==4):
                patches.shape = list(patches.shape[:2])+[1] + list(patches.shape[2:])


            if  patches.shape[2]!=1:
                ep = patches.shape[2]
                assert( ep%2==1)
                ep=ep//2
                P.patch_ep = ep
                P.last_slices_2r_nonsplitted  = numpy.minimum( P.END_VOXEL_3-1   , P.last_slices_2r_nonsplitted+ep)
                P.first_slices_2r_nonsplitted  = numpy.maximum( P.START_VOXEL_3 -1, P.first_slices_2r_nonsplitted-ep)


    P.first_slices = P.first_slices_2r_nonsplitted
    P.last_slices  = P.last_slices_2r_nonsplitted 

    P.CONICITY_MARGIN_DOWN =  numpy.zeros( len( P.first_slices )  ,"i") 
    P.CONICITY_MARGIN_UP   =  numpy.zeros( len( P.first_slices )  ,"i") 

    # P.corr_slice_slicedect = numpy.arange(  P.START_VOXEL_3 -1, P.END_VOXEL_3-1 +1, 1)
    CHANGEME = 20000
    P.corr_slice_slicedect = numpy.arange(  0, P.SOURCE_Y+P.END_VOXEL_3-1 +1   + CHANGEME, 1)
    P.VOXEL_SIZE = P.IMAGE_PIXEL_SIZE_1  

    if P.CONICITY :

        f_tmp  = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE))    # ??? /P.IMAGE_PIXEL_SIZE_2
        P.VOXEL_SIZE = P.IMAGE_PIXEL_SIZE_1/f_tmp

        Fact = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE))*P.VOXEL_SIZE*1.0/P.IMAGE_PIXEL_SIZE_2 
        DFa  = +(abs(P.DECT_PSI)+abs(P.DECT_TILT) )*3.2/180.0*2


        
        DR1  = (P.END_VOXEL_1 - P.START_VOXEL_1)/2.0
        DR2  = (P.END_VOXEL_2 - P.START_VOXEL_2)/2.0
        RADIUS = math.sqrt( DR1*DR1+DR2*DR2   ) 

        for i in range(len(P.first_slices) ):



            if P.DZPERPROJ!=0:
                alpha = ((0-P.SOURCE_Y)*P.VOXEL_SIZE/1.0e6)/(P.SOURCE_DISTANCE)
            else:
                alpha = (abs(P.first_slices[i])*P.VOXEL_SIZE/1.0e6)/(P.SOURCE_DISTANCE)
                
            P.CONICITY_MARGIN_DOWN[i] = int(  RADIUS*(
                                              alpha*(Fact) +DFa)
                                              +0.01)+1

            if P.DZPERPROJ!=0:
                alpha = ((P.NUM_IMAGE_2-P.SOURCE_Y)*P.VOXEL_SIZE/1.0e6)/(P.SOURCE_DISTANCE)
            else:
                alpha = (abs(P.last_slices[i])*P.VOXEL_SIZE/1.0e6)/(P.SOURCE_DISTANCE)


            P.CONICITY_MARGIN_UP[i]   = int(  RADIUS*(
                                              alpha*(Fact) +DFa)
                                              +0.01)+1


        P.first_slices  = P.SOURCE_Y + (P.first_slices)                 * Fact +0.49999
        P.last_slices  = P.SOURCE_Y + (P.last_slices)                   * Fact +0.49999
        P.first_slices  = P.first_slices.astype("i")  
        P.last_slices   = P.last_slices .astype("i")   
    
        P.corr_slice_slicedect  = P.SOURCE_Y + (P.corr_slice_slicedect-P.SOURCE_Y) * Fact +0.49999


    P.corr_slice_slicedect  = P.corr_slice_slicedect.astype("i")
    P.CONICITY_MARGIN_UP_original   =  P.CONICITY_MARGIN_UP  
    P.CONICITY_MARGIN_DOWN_original =  P.CONICITY_MARGIN_DOWN

    # P.last_slices  = numpy.minimum( P.last_slices  + P.CONICITY_MARGIN_UP  ,  P.END_VOXEL_3-1 )
    tmp_array                =   P.last_slices
    if P.DZPERPROJ==0:
        P.last_slices            =   numpy.minimum( P.last_slices  + P.CONICITY_MARGIN_UP  ,  P.NUM_IMAGE_2-1 )
    else:
        P.last_slices            =    P.last_slices  + P.CONICITY_MARGIN_UP  
    P.CONICITY_MARGIN_UP      =   P.last_slices-tmp_array

    tmp_array                =   P.first_slices

    if P.DZPERPROJ==0:
        P.first_slices           =   numpy.maximum( P.first_slices - P.CONICITY_MARGIN_DOWN  ,   0 )
    else:
        P.first_slices           =    P.first_slices - P.CONICITY_MARGIN_DOWN 

    P.CONICITY_MARGIN_DOWN   =  -P.first_slices+tmp_array
    

    if(min(P.CONICITY_MARGIN_UP )<0 or min(P.CONICITY_MARGIN_DOWN )<0) :
        raise Exception( "RECONSTRUCTON REGION PROJECTS OUT OF DETECTOR LIMITS" )



    P.CONICITY_MARGIN =   int(numpy.max(P.CONICITY_MARGIN_UP_original +P.CONICITY_MARGIN_DOWN_original  ) )

    P.pos_s   =  [  [first_slice,0]      
                    for first_slice in P.first_slices ]

    P.size_s  =  [  [   last_slice-first_slice+1     , P.NUM_IMAGE_1 ]
                    for first_slice, last_slice in zip(P.first_slices, P.last_slices) ] 

    P.pos_s_  = [None]*len(P.pos_s )
    P.size_s_ = [None]*len(P.pos_s )
    for pos, size , count in zip( P.pos_s ,P.size_s , range(len(P.pos_s ))      ):
        P.pos_s_[count], P.size_s_[count] = readjust_marges(pos,size)

    # print( " _____________________________________________________________")
    print( P.pos_s_, P.size_s_)

    P.pos_s = numpy.array(P.pos_s, dtype=numpy.int32)
    P.size_s    = numpy.array( P.size_s , dtype=numpy.int32 )
    P.pos_s_    = numpy.array( P.pos_s_ , dtype=numpy.int32 )
    P.size_s_     = numpy.array( P.size_s_ , dtype=numpy.int32  )

    if Parameters.VERBOSITY>2 : print( " DEBUG ",  "   P.pos_s_, P.size_s_               " ,P.pos_s_, P.size_s_)
    if Parameters.VERBOSITY>2 : print( " DEBUG ",  "   P.pos_s, P.size_s               " ,P.pos_s, P.size_s)
    if Parameters.VERBOSITY>2 : print( " DEBUG ", " P.first_slices ", P.first_slices)
    if Parameters.VERBOSITY>2 : print( " DEBUG ", " P.last_slices ", P.last_slices)
    if Parameters.VERBOSITY>2 : print( " DEBUG ", "  P.corr_slice_slicedect ", P.corr_slice_slicedect[0], "  ..  " , P.corr_slice_slicedect[-1])
    



    if Parameters.LT_MARGIN :
        Parameters.LT_INFOS_COARSE = LTSparsePrep.Get_Sparse_LT(Parameters)
        # Parameters.LT_INFOS_FINE   = LTSparsePrep.Get_Sparse_LT(Parameters, finegrid=True)


    
# proj_reading_dict generated by  
def create_arrays_space_as_dictionary(proj_reading_dict):
    P=Parameters

    dimsP = numpy.max(P.size_s_, axis=0)
    num_myprojs = len(proj_reading_dict["my_proj_num_list"])
    num_myprojs_reduced = int( num_myprojs *1.0/ P.RAWDATA_MEMORY_REUSE+0.99999999999)

    if num_myprojs_reduced%P.NPBUNCHES > 0 :
        num_myprojs_reduced =   num_myprojs_reduced  + (P.NPBUNCHES -  num_myprojs_reduced%P.NPBUNCHES     )

    dims = [ num_myprojs_reduced      , dimsP[0], dimsP[1] ]

    nmax=0
    for l in proj_reading_dict["ff_indexes_coefficients"]:
        nmax=max(nmax,l[0])
        nmax=max(nmax,l[1])
        
    
    ff_dims = [ int(1+nmax) , dimsP[0], dimsP[1] ]

    
    ntokensraw=min(P.MAXNTOKENS,len(P.pos_s_))

    if Parameters.VERBOSITY>2 : print( " P.MAXNTOKENS,len(P.pos_s_)  " , P.MAXNTOKENS,len(P.pos_s_))


    if Parameters.VERBOSITY>2 : print( "  rawdatatokens " )
    if Parameters.VERBOSITY>2 : print( dims)

    rawdatatokens =  [ numpy.zeros( dims  ,"f")   for i in range(ntokensraw)  ]  
    ff_rawdatatokens =  [ numpy.zeros(ff_dims  ,"f")   for i in range(ntokensraw)  ]
    
    dims =  numpy.max(P.size_s, axis=0)
    dims =  [  len( proj_reading_dict["my_proj_num_list"] ) , dims[0] , dims[1]   ]

    ntokens   = min( P.MAXNTOKENS,len(P.pos_s_))
    if Parameters.VERBOSITY>2 : print( " P.MAXNTOKENS,len(P.pos_s_)  " , P.MAXNTOKENS,len(P.pos_s_) ,  ntokens, ntokensraw)
    datatokens   =  [ numpy.zeros( dims  ,"f")   for i in range(ntokens)  ]
 
    # slices_mpi_offsets   = []
    # slices_mpi_numslices = []


    dims =  numpy.max(P.size_s, axis=0)
    dims =  [   P.slices_mpi_numslices[myrank] + P.CONICITY_MARGIN  ,   P.numpjs_span     ,dims[1]     ]

    transposeddatatokens   =  [ numpy.zeros( dims  ,"f")   ]
    
    #########################################################################################
    # read and store the BACKGROUND in variable background fully dimensioned
    # -------------------------------------------------------------------------------------
    if(P.SUBTRACT_BACKGROUND):
        if P.BACKGROUND_FILE[-3:]==".h5" or P.BACKGROUND_FILE[-4:]==".nxs":
            h5=h5py.File( P.BACKGROUND_FILE , "r")
            print(h5.keys())
            background = h5[P.BACKGROUND_DS_NAME][:]
#            if len(background.shape) == 3 :
#                print(" DETECTED STACK FOR BACKGROUND, MEDIAN WILL BE APPLIED " )
            background = background.astype("f")
                
        else:
            if Parameters.VERBOSITY>0 :  print( " reading edf file for background ", P.BACKGROUND_FILE)
            dark=EdfFile.EdfFile(P.BACKGROUND_FILE)
            background=dark.GetData(0,DataType="FloatValue")
            
        if P.BINNING >1 :
            if len(background.shape) == 2 :
                
                sh0,sh1 = background.shape
                sh0=(sh0//P.BINNING)*P.BINNING
                sh1=(sh1//P.BINNING)*P.BINNING
                background=numpy.reshape(background[:sh0,:sh1],[sh0 //P.BINNING   , P.BINNING,  sh1//P.BINNING  , P.BINNING ])
                background=numpy.sum(numpy.sum(background , axis=-1), axis=1)//(P.BINNING**2)
                
                
                
            else:
                assert( len(background.shape) == 3   )
                
                shz, sh0,sh1 = background.shape
                
                sh0=(sh0//P.BINNING)*P.BINNING
                sh1=(sh1//P.BINNING)*P.BINNING
                background=numpy.reshape(background[:shz, :sh0,:sh1],[shz, sh0 //P.BINNING   , P.BINNING,  sh1//P.BINNING  , P.BINNING ])
                background=numpy.sum(numpy.sum(background , axis=-1), axis=2)//(P.BINNING**2)


        if len(background.shape) == 3 :       
            background = numpy.swapaxes(background,0,2)
            background = numpy.swapaxes(background,0,1)
            background = numpy.ascontiguousarray(background)
            
        if(Parameters.ROTATION_VERTICAL==0):
            background=numpy.array(numpy.transpose(background), copy=True)
    else:
        background=numpy.zeros([ P.NUM_IMAGE_2, P.NUM_IMAGE_1 ],"f")
    assert(   background.flags['C_CONTIGUOUS'] )
        
    if(P.DOUBLEFFCORRECTION):
        if Parameters.VERBOSITY>0 :  print( " reading edf file for DOUBLEFFCORRECTION  ", P.DOUBLEFFCORRECTION)
        ffcorredf =EdfFile.EdfFile(P.DOUBLEFFCORRECTION)
        ffcorr =ffcorredf.GetData(0,DataType="FloatValue")
        if P.BINNING is not None:
            sh0,sh1 =ffcorr .shape
            sh0=(sh0//P.BINNING)*P.BINNING
            sh1=(sh1//P.BINNING)*P.BINNING
            
            ffcorr =numpy.reshape(ffcorr[:sh0,:sh1],[ sh0//P.BINNING   , P.BINNING,  sh1//P.BINNING  , P.BINNING ])
            ffcorr=numpy.sum(numpy.sum(ffcorr , axis=-1), axis=1)/(P.BINNING**2)
            
        if(Parameters.ROTATION_VERTICAL==0):
            ffcorr =numpy.array(numpy.transpose(ffcorr), copy=True)
    else:
        ffcorr =numpy.zeros([ P.NUM_IMAGE_2, P.NUM_IMAGE_1 ],"f")

    if( Parameters.TAKE_LOGARITHM):
        ffcorr[:]  =   numpy.exp(+ffcorr)
    else:
        ffcorr[:]  =   -ffcorr



    #####################################################################################################
    ## 
    if P.ITERATIVE_CORRECTIONS>0 and P.DENOISING_TYPE in [4]:

        
        patches = retrieve_patches(  P.PATCHES_FILE   ) 



        assert(len(patches.shape)in [4,5])
        if(len(patches.shape)==4):
            patches.shape = list(patches.shape[:2])+[1] + list(patches.shape[2:])
        
        if P.DENOISING_TYPE==4:
            add = numpy.zeros_like( patches[0])
            add[:] = numpy.sqrt(1.0/add.size)
            patches=numpy.concatenate( ([add],patches), axis=0)

    else:
        patches=numpy.array([[[[[0.0]]]]],"f")


    return {"ff_rawdatatokens":ff_rawdatatokens ,"rawdatatokens":rawdatatokens ,
            "datatokens":datatokens, "background":background,"transposeddatatokens":transposeddatatokens,
            "doubleffcorrection":ffcorr ,"patches":patches}



Parameters.ORIGINAL_FILE_PREFIX =  Parameters.FILE_PREFIX
if Parameters.JP2EDF_DIR is not None:
    Parameters.FILE_PREFIX = Parameters.JP2EDF_DIR+"/" + os.path.basename(Parameters.FILE_PREFIX)
    make_sure_path_exists( Parameters.JP2EDF_DIR  )

def get_name_edf_proj(n,FILE_PREFIX=None, LENGTH_OF_NUMERICAL_PART=None , NUMBER_LENGTH_VARIES=None,FILE_POSTFIX=None  ):
    P=Parameters
    name=FILE_PREFIX
    nlength = LENGTH_OF_NUMERICAL_PART
    if(NUMBER_LENGTH_VARIES):
        name=name+("%d"%n)+FILE_POSTFIX
    else:
        number = "%d"%n
        if(len(number)< nlength ):
                number = ("0" *  (nlength-len(number)))+ number
        else:
           number=number[:nlength]
        name = name+number+FILE_POSTFIX
    return name