#/*##########################################################################
# 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.
#############################################################################*/
#include<semaphore.h>


#ifndef FROMCU
#include <fftw3.h>
#endif

#include<complex.h>
#include<cufft.h>
#include"edftools.h"


float  quickselect(float * A, int left, int right, int k);


typedef  float __complex__   fcomplex ;
// typedef  fftwf_complex   fcomplex ;

// fftw_complexf






#define MAXNTOKENS 100
#define MAXALLOCATIONS 100000

#define max(a,b) ( (a)>(b)? (a):(b) )
#define min(a,b) ( (a)<(b)? (a):(b) )



 /* --------------------*/
/*      CCarraylist    */

typedef struct {
  int ntokens ;
  float *datatokens[MAXNTOKENS];
  int *islice_tracker[MAXNTOKENS];
} CCarraylist;
void CCarraylist_initialise(CCarraylist * self) ;
void  CCarraylist_appendarray(CCarraylist* self, float * array) ;



/* --------------------*/
/*   CCreading_infos   */

typedef struct {
  char *proj_reading_type;
  char *ff_reading_type;
  int allfile_list_lenght;
  int proj_file_list_lenght;
  char **proj_file_list;
  char **allfile_list;

  int corr_file_list_lenght;
  char **corr_file_list;


  int ff_file_list_lenght;
  char **ff_file_list;
  float *ff_indexes_coefficients;
  int NFF;
  int * my_proj_num_list;
  int * tot_proj_num_list;
  int * proj_num_offset_list;
  int * proj_mpi_offsets;
  int * proj_mpi_numpjs;
  int * file_proj_num_list;

  long int* headerSizes;
  long int* ff_headerSizes;
  long int* corr_headerSizes;


  
  char *  CURRENT_NAME; 
  float * currents ; 
  int    currents_lenght ; 
  
  int datatype ;
  int byteorder;
  int sizeImage;
  int Dim_1     ;
  int Dim_2     ;

  int corr_restrained;
  int corr_start_2;
  int corr_real_2;


  int ff_datatype ;
  int ff_byteorder;
  int ff_sizeImage;

  int corr_datatype ;
  int corr_byteorder;
  int corr_sizeImage;

  int* pos_s   ;
  int* size_s  ;
  int* pos_s_  ;
  int* size_s_ ;

  int * slices_mpi_offsets ;
  int * slices_mpi_numslices ;

  char * proj_h5_dsname;
  char * ff_h5_dsname;


  int nchunks;
  // da cancellare int nchunks  ;

  int MULTIFRAME;
  int NUM_FIRST_IMAGE;

  int numpjs;
  int numpjs_span;


  int NUM_LAST_IMAGE ; 

} CCreading_infos;

void CCreading_infos_PrereadHeaders(CCreading_infos* self,  int tryconstheader, char ** file_list, int FILE_INTERVAL)  ;




//PP.add
typedef struct {
    cufftReal* d_r_sino_error;
    cufftComplex* d_i_sino_error;
    cufftHandle planRamp_forward;
    cufftHandle planRamp_backward;

    cufftComplex* filter_coeffs;
} PRECOND_PARAMS_DL;



typedef struct lt_infos {
  int Ng;  // numero di gaussiane = max(j)+1
  int Ng2s ;  //  numero di gaussiane che vanno sull'immagine centrale = max(J2s)+1
  // come numero di punti del sino gramma si puo prendere Nsigmas*SLD
  //                                 P.numpjs_span * (num_bins +6*P.LT_MAX_STEP)
  // 
  // Come numero di punti della slice si puo prendere  max(I2s)+1
  int Nslice; 
  
  int N;
  float * C;
  int   * II;
  int   * JJ;
  int   SLD;
  int nprojs ; 
  int Nsigmas;
  float *sigmas;

  int N2s;
  float * C2s;
  int   * I2s;
  int   * J2s;
  int dim2s_s;
  int dim2s_g;
  float *slice_comp ;
  int *corre;
  
  float * gaussian_coeffs ; 

  int lt_reb;

  float *fftwR   ;
  fcomplex *fftwC ;
  double force_per_gaussian;
  int lt_max_step;
  int lt_max_step_reb;
#ifndef FROMCU
  fftwf_plan planr2c, planc2r;
#endif
  
} LT_infos;
#define GPU_LT_MAX_NDIFF4thread  12

typedef struct CCspace_struct CCspace;

void free_LT_infos(LT_infos *p, CCspace *self) ;

typedef struct {
  char * LT_KNOWN_REGIONS_FILE;
  int LT_MARGIN;
  int LT_REBINNING;
  LT_infos*   lt_infos_coarse ;
 } LT_Parameters;



/* ------------*/
/*  Cparameters    */
typedef struct {
  int SUBTRACT_BACKGROUND;
  float REFERENCES_HOTPIXEL_THRESHOLD;
  int CORRECT_FLATFIELD;
  float OFFSETRADIOETFF;
  int TAKE_LOGARITHM;
  float    ZEROCLIPVALUE  ;
  float ONECLIPVALUE ;
  float  ROTATION_AXIS_POSITION ;
  char *  DOUBLEFFCORRECTION;
  int DO_PAGANIN;
  int DO_RING_FILTER;
  int DO_CCD_FILTER;
  int NO_SINOGRAM_FILTERING;
  int NUM_IMAGE_1;  /*  # Number of pixels horizontally  (take into account binning : smaller if binning > 1 ) */
  int NUM_IMAGE_2;  /* # Number of pixels vertically    (take into account binning )  */

  float PAGANIN_Lmicron;
  int PAGANIN_MARGE ;
  float  IMAGE_PIXEL_SIZE_2;
  float  IMAGE_PIXEL_SIZE_1;
  float  VOXEL_SIZE;
  float PUS;
  float PUC;
  int UNSHARP_LoG;
  int CCD_FILTER_KIND;
  void *CCD_FILTER_PARA;

  int DO_OUTPUT_PAGANIN ;
  char *OUTPUT_PAGANIN_FILE ;


  int RING_FILTER_KIND;
  void *RING_FILTER_PARA;

  float *NEURAL_FILTERS;
  float *NEURAL_OFFSETS;
  float *NEURAL_WEIGHTS;
  int NEURAL_NFILTERS;
  int NEURAL_FILTERSIZE;
  float NEURAL_MININ;
  float NEURAL_MAXIN;

  int NEURAL_TRAINING_PIXELS_PER_IMAGE;
  int NEURAL_TRAINING_NLINEAR;
  int NEURAL_TRAINING_USEMASK;
  char *NEURAL_TRAINING_RECONSTRUCTION_FILE;
  char *NEURAL_TRAINING_MASK_FILE;

  int nprojs_span;

  char *OUTPUT_FILE;
  char *PROJ_OUTPUT_FILE ;
  char *OUTPUT_FILE_ORIG;
  char *OUTPUT_FILE_HISTOGRAM;
  // char *XMLOUTPUT;
  // char *INFOOUTPUT;

  float FBFILTER;
  float * FBFACTORS;
  int DFI_NOFVALUES  ;
  int DFI_KERNEL_SIZE ;
  int DFI_OVERSAMPLING_RATE ;
  int DFI_R2C_MODE ;
  int DO_PRECONDITION; //ramp-filter in Fourier domain for precontition
  PRECOND_PARAMS_DL precond_params_dl;

  //DFP
  int DFP_NOFVALUES;
  int DFP_KERNEL_SIZE;
  int DFP_OVERSAMPLING_RATE;

  int USE_DFI;
  int USE_DFP;

  //Wavelets
  int W_LEVELS;
  int W_CYCLE_SPIN;
  float W_FISTA_PARAM;
  char* W_WNAME;
  int W_SWT;
  int W_NORMALIZE;
  int W_THRESHOLD_APPCOEFFS;
  int W_THRESHOLD_COUSINS;

  // Fourier-Wavelet (Munch Et Al) filter
  int FW_LEVELS;
  float FW_SIGMA;
  char* FW_WNAME;
  float FW_SCALING;
  float FW_FILTERPARAM;

  // Fluo
  int FLUO_SINO;
  int FLUO_ITERS;

  // Sirt-Filter
  int SF_ITERATIONS;
  char* SF_SAVEDIR;
  float SF_LAMBDA;
  float *SF_FILTER;

  int DATA_IS_FILTERED;
  int ESTIMATE_BETA;

  // TV+L2
  float BETA_L2;

  // Positivity constraint
  int ITER_POS_CONSTRAINT;

  // For double flat field
  int DOUBLEFFCORRECTION_ONTHEFLY;
  float FF2_SIGMA;



  int start_x ;
  int start_y ;
  int num_x   ;
  int num_y   ;
  float ANGLE_BETWEEN_PROJECTIONS;
  float ANGLE_OFFSET;
  int OVERSAMPLING_FACTOR;
  int do_custom_angles;
  int do_ignore_projections;
  float * custom_angles;
  int* ignore_angles;
  int AXIS_TO_THE_CENTER ;
  char *PADDING ;
  int AVOIDHALFTOMO ;
  int OUTPUT_SINOGRAMS;

  float     PENTEZONE ;
  int zerooffmask ;
  int SUMRULE ;
  int ROTATION_VERTICAL;
  int TRYGPU;
  int TRYGPUCCDFILTER;
  int ALLOWBOTHGPUCPU;
  int RAWDATA_MEMORY_REUSE;
  int BINNING;
  int MYGPU;
  int MULTIFRAME;
  int ITERATIVE_CORRECTIONS;
  float SMOOTH_PENALTY_PARAM;
  float LINE_SEARCH_INIT;
  int PENALTY_TYPE;
  int OPTIM_ALGORITHM;
  int VECTORIALITY ;
  int FISTA;
  int STEPFORPATCHES ;

  int DENOISING_TYPE;
  int N_ITERS_DENOISING;
  float  DUAL_GAP_STOP;
  int CALM_ZONE_LEN;
  float NLM_NOISE_GEOMETRIC_RATIO;
  float NLM_NOISE_INITIAL_FACTOR;

  int STRAIGTHEN_SINOS;
  int FILE_INTERVAL;
  int NUM_FIRST_IMAGE;
  int EDFOUTPUT ;
  float BETA_TV;


  float   ITER_RING_HEIGHT;
  float   ITER_RING_BETA;
  float   ITER_RING_SIZE;
  float LIPSCHITZFACTOR;
  int LIPSCHITZ_ITERATIONS; //PP.add
  //int DO_PRECONDITION;
  float RING_ALPHA; //PP.add
  int NUMBER_OF_RINGS;

  int JOSEPHNOCLIP;

  float    DETECTOR_DUTY_RATIO        ;
  int      DETECTOR_DUTY_OVERSAMPLING ;

  float* patches;
  int    patches_N;
  int    patches_size ;
  int    patches_dim ;
  int    patches_vecto ;
  int    patches_epaisseur ;
  char *nome_directory_distribution;
  float peso_overlap;

  int CONICITY ;
  int CONICITY_FAN ;
  float SOURCE_DISTANCE ;
  float DETECTOR_DISTANCE ;
  float SOURCE_X        ;
  float SOURCE_Y      ;
  float DECT_PSI;
  float DECT_TILT;
  

  float DZPERPROJ;
  float DXPERPROJ;

  int * CONICITY_MARGIN_UP ;
  int * CONICITY_MARGIN_DOWN ;
  int * CONICITY_MARGIN_UP_wished ;
  int * CONICITY_MARGIN_DOWN_wished ;
  int   CONICITY_MARGIN ;

  int STEAM_INVERTER ;
  int PROJ_OUTPUT_RESTRAINED;


  int *first_slices_2r  ;
  int *last_slices_2r  ;
  int *first_slices  ;
  int *last_slices  ;
  int *corr_slice_slicedect  ;
  int patch_ep;

  int do_dump_hrings;

  int *intervals_hrings;
  int N_intervals_hrings;

  int verbosity ;

  LT_Parameters *LT_pars;



  float * dist_h, *dist_v;
  int dist_nrow, dist_ncol;
  
  float * BH_LUT_F, *BH_LUT_U;
  int BH_LUT_N;



} Cparameters ;

#ifndef FROMCU

typedef struct {
  float  * FILTER ;
  float * NEURALFILTER;
  int *minX;
  int * maxX ;
  float *axis_position_corr_s;
  float *cos_s               ;
  float *sin_s               ;
  int fai360;
  float MOVE_X;
  float MOVE_Y;
  fftwf_plan planr2c, planc2r;
  float overlapping , pente_zone, flat_zone, prof_shift, prof_fact;
  float *  angles_per_proj ;
  float prec_gamma;
  int prec_gamma_is_set;
  float Lipschitz_fbdl;
  int dim_fft;
} FBP_prec  ;
#endif


typedef struct {
  float *data;
  int nslices_2r ;
  int nslices_data ;
  int Nfirstslice ;
  int data_start;
  int has_been_set_flag;

} SlicesRequest;

typedef struct {
  int pps_grid_cols; //!< Number of CUDA blocks along X-axis of 1D Fouried sinogram
  int pps_grid_rows; //!< Number of CUDA blocks along Y-axis of 1D Fouried sinogram
  int interp_grid_cols; //!< Number of CUDA blocks along X-axis of 2D Fourier domain
  int interp_grid_rows; //!< Number of CUDA blocks along Y-axis of 2D Fourier domain
  int swap_grid_cols; //!< Number of CUDA blocks along X-axis of 2D Fourier domain
  int swap_grid_rows; //!< Number of CUDA blocks along Y-axis of 2D Fourier domain
  int swap_quad_grid_cols; //!< Number of CUDA blocks along X-axis of 2D Fourier domain
  int swap_quad_grid_rows; //!< Number of CUDA blocks along Y-axis of 2D Fourier domain
  int roi_grid_cols; //!< Number of CUDA blocks along X-axis of ROI
  int roi_grid_rows; //!< Number of CUDA blocks along Y-axis of ROI
  int points_per_thread;  //!< How many points is processed by a single thread (actually square of that)

  int rho_len; //!< Initial length of sinogram  along rho direction
  int rho_ext_len; //!< Extended length of sinogram along rho direction
  int nprojs_span; //!< Number of projections in sinogram (the length along theta direction)
  float L; //!< Length of one side of interpolation kernel
  float L2; //!< A half of length of interpolation kernel
  int ktbl_len; //!< Number of presampled kernel values
  int ktbl_len2; //!< A half of number of presampled kernel values
  int raster_size; //!< Length of side of output (2D Fourier domain)
  int raster_size2; //!< A half of a length of number of presampled kernel values
  int roi_start_x; //!< Starting coordinate of ROI along X-axis
  int roi_start_y; //!< Starting coordinate of ROI along Y-axis
  int roi_x; //!< The length of ROI along X-axis
  int roi_y; //!< The length of ROI along Y-axis
  float table_spacing; //!< An interval between samples of presampled kernel values
  float angle_step_rad; //!< A step between projections (in degrees)
  float theta_max; //!< Max value of theta along Y-axis
  float rho_max; //!< Max value of rho along X-axis
  int oversampling; //!< An oversampling coefficient
  float scale;

  float *ktbl; //!< Temp array of presampled kernel values
  int spectrum_offset_y; //!< Number of dropped values along Y-axis (due to reducing 2D Fourier domain radius)
  float max_radius; //!< Maximum radius of 2D Fourier domain

  int fft_sino_dim;

  cufftHandle fft1d_plan; //!< Complex plan for forward fourier transformations
  cufftHandle ifft2d_plan;  //!< Complex plan for inverse fourier transformations

  void *gpu_ktbl; //!< Array of presampled values at GPU
  float *gpu_truncated_sino; //!< Sinogram after truncation

  float *gpu_zeropadded_sino; //!< Zeropadded sinogram using real values
  cufftComplex *gpu_input; //!< Zeropadded sinogram using complex values
  void *gpu_input_cua; //!< Zeropadded sinogram using complex values

  cufftComplex *gpu_spectrum; //!< Array of reconstructed 2D Fourier domain for GPU
  cufftComplex *gpu_swapped_spectrum; //!< Array of swapped 2D Fourier domain for GPU

  float *gpu_output; //!< Reconstructed slice is going here
  float *gpu_c2r_result; //!< Result of 2D IFFT - Complex-to-Real transform

} DFI_params;

typedef struct {
  int rho_len; //!< Initial length of sinogram along rho direction
  int rho_len2; //!< Initial length of sinogram along rho direction

  int theta_len; //!< Initial length of sinogram along theta direction

  int rho_ext_len; //!< Nearest value power of 2 to the initial known sinogram length

  int slice_size_x; //!< Length along x-axis of the slice
  int slice_size_y; //!< Length along y-axis of the slice

  float L; //!< Length of one side of interpolation kernel
  float L2; //!< A half of length of interpolation kernel

  int ktbl_len; //!< Number of presampled kernel values
  int ktbl_len2; //!< A half of number of presampled kernel values

  float table_spacing; //!< An interval between samples of presampled kernel values
  float angle_step_rad; //!< A step between projections (in degrees)

  float *ktbl; //!< Temp array of presampled kernel values
  void *gpu_ktbl; //!< Array of presampled values at GPU

  int offset_x; //!< Offset from x-axis border on zeropadding
  int offset_y; //!< Offset from y-axis border on zeropadding

  int dfp_grid_cols; //!< Number of CUDA blocks along rho-axis of sinogram
  int dfp_grid_rows; //!< Number of CUDA blocks along theta-axis of sinogram

  int swap_grid_cols; //!< Number of CUDA blocks along x-axis of zeropadded slice
  int swap_grid_rows; //!< Number of CUDA blocks along y-axis of zeropadded slice

  int shift_grid_cols; //!< Number of CUDA blocks along rho-axis of sinogram
  int shift_grid_rows; //!< Number of CUDA blocks along theta-axis of sinogram

  int zeropad_grid_cols; //!< Number of CUDA blocks along x-axis of zeropadded slice
  int zeropad_grid_rows; //!< Number of CUDA blocks along y-axis of zeropadded slice

  int crop_grid_cols; //!< Number of CUDA blocks along rho-axis of sinogram
  int crop_grid_rows; //!< Number of CUDA blocks along theta-axis of sinogram

  cufftComplex *gpu_zeropadded_slice; //!< Result of slice zeropadding
  cufftComplex *gpu_reconstructed_sinogram; //!< Complex numbers of the reconstruction of sinogram obtained with 1D-IFFT
  float *gpu_reconstructed_sinogram_real; //!< Real part of the reconstruction of sinogram
  float *gpu_reconstructed_crop_sinogram_real; //!< Real part of the reconstruction of sinogram

  cufftHandle ifft1d_plan; //!< Complex plan for inverse 1D fourier transformations
  cufftHandle fft2d_plan;  //!< Complex plan for forward 2D fourier transformations

} DFP_params;
void * getLibNameHandle(       const char *dirname,   const char *prefix) ;


typedef struct Gpu_Context_struct Gpu_Context;

int gpu_mainInit(Gpu_Context * self, float *filter );
int gpu_main(Gpu_Context * self, float *WORK , float * SLICE, int do_precondition,
         float    DETECTOR_DUTY_RATIO        ,
         int    DETECTOR_DUTY_OVERSAMPLING,
         float *fidelity,
         int npj_offset, int forcePImultipl);


int gpu_main_2by2(Gpu_Context * self, float *WORK , float * SLICE, int do_precondition,
		  int npj_offset, int forcePImultipl ) ;




int dfi_gpu_main(Gpu_Context * self, float *WORK , float * SLICE, int memisonhost);


float gpu_lt_fit_sino(Gpu_Context * self, int Ng, float *gaussian_coeffs,   int   Ns , float * data, int N,
				       int  *csc_II,  int * csc_JJ  , float * csc_C , int csc_nitems_atonce,
				       int  *csr_II,  int *  csr_JJ , float * csr_C , int csr_nitems_atonce,
				       int SLD,  int Nsigmas,   float *sigmas,  LT_infos* lt, LT_infos *lt2s
				       );


int       gpu_main_conicity(Gpu_Context * self,  float * SLICE, float *WORK_perproje ,
                 int Nfirstslice, int nslices, int data_start, int nslices_data,
                 float source_distance, float detector_distance,
                 float v2x ,   float v2z,float voxel_size,
                float SOURCE_X , float SOURCE_Z );

int       pro_gpu_main_conicity(Gpu_Context * self,  float * SLICE, float *data,
                 int Nfirstslice, int nslices, int data_start, int nslices_data,
                 float source_distance, float detector_distance,
                 float v2x ,   float v2z,float voxel_size,
                float SOURCE_X , float SOURCE_Z );





void fb_dl(Gpu_Context * self, float *WORK , float * SLICE, int precondition,
       float    DETECTOR_DUTY_RATIO        ,
       int    DETECTOR_DUTY_OVERSAMPLING,
       float weight,
       float *patches, int npatches, int dim_patches,  int VECTORIALITY, float *lip,
     void  (*proietta_drings) (void *, float* , int , float)
       ) ;

void search_lipschitz(Gpu_Context * self, float *WORK , float * SLICE, int precondition,
       float    DETECTOR_DUTY_RATIO        ,
       int    DETECTOR_DUTY_OVERSAMPLING,
       float weight,
       float *patches, int npatches, int dim_patches,  int VECTORIALITY, float *lip,
     void  (*proietta_drings) (void *, float* , int , float)
       ) ;


float tv_denoising_fistagpu(Gpu_Context * self, int dim0,int dim1,float  *img,float  *result , float weight,
                float eps , int n_iter_max,   int check_gap_frequency ) ;

float tv_denoising_patches_L1(Gpu_Context * self, int dim0,int dim1,float  *img,float  *result , float weight,
                  float *patches, int npatches, int dim_patches, int N_ITERS_DENOISING,
                  int vectoriality) ;

float tv_denoising_patches_OMP(Gpu_Context * self, int dim0,int dim1,float  *img,float  *result , float weight,
               float *patches, int npatches, int dim_patches) ;


void nonlocalmeans(    Gpu_Context * self,
               float *result,
               float *image,
               int H,
               int W,
               float bruit,
               float mix
               );



void  C_HST_PROJECT_1OVER(
              int  num_bins,        /* Number of bins in each sinogram  */
   int  nprojs_span, /* Number of projections in  each sinogram */
   float*  angles_per_project,
   float   axis_position,
   float* SINOGRAMS ,
   float *SLICE,
   int   dimslice,
   float * axis_corrections,
   float cpu_offset_x,
   float cpu_offset_y
   );


void  C_HST_PROJECT_1OVER_GPU(
                  void *gpuctx,
                  int  num_bins,        /* Number of bins in each sinogram  */
                  int  nprojs_span, /* Number of projections in  each sinogram */
                  float*  angles_per_project,
                  float   axis_position,
                  float* SINOGRAMS ,
                  float *SLICE,
                  int   dimslice,
                  float * axis_corrections,
                  float cpu_offset_x,
                  float cpu_offset_y,
                  int josephnoclip,
                  float    DETECTOR_DUTY_RATIO        ,
                  int    DETECTOR_DUTY_OVERSAMPLING ,
                  int memisonhost,
                  float fan_factor,
                  float source_x
                  );


typedef  int (*gpu_project_Symbol )(
                    void *gpuctx,
                    int  num_bins,        /* Number of bins in each sinogram  */
                    int  nprojs_span, /* Number of projections in  each sinogram */
                    float*  angles_per_project,
                    float   axis_position,
                    float* SINOGRAMS ,
                    float *SLICE,
                    int   dimslice,
                    float * axis_corrections,
                    float cpu_offset_x,
                    float cpu_offset_y,
                    int josephnoclip,
                    float    DETECTOR_DUTY_RATIO        ,
                    int    DETECTOR_DUTY_OVERSAMPLING,
                    int memisonhost,
                    float fan_factor,
                    float source_x
                    );
//PP.add :
typedef  int (*gpu_backproject_Symbol )(
    Gpu_Context * self,
    float *d_S ,
    float * SLICE,
    int do_precondition,
    float DETECTOR_DUTY_RATIO,
    int DETECTOR_DUTY_OVERSAMPLING,
    int doppio,
    int memisonhost
    );

typedef int (*cp_driver_Symbol) (
    Gpu_Context* self,
    float* data,
    float* SLICE,
    float DETECTOR_DUTY_RATIO,
    int DETECTOR_DUTY_OVERSAMPLING,
    float beta_tv,
    float beta_rings,
    float rings_height,
    float alpha_rings
    );

typedef int (*wavelets_driver_Symbol) (
    Gpu_Context* self,
    float* data,
    float* SLICE,
    float DETECTOR_DUTY_RATIO,
    int DETECTOR_DUTY_OVERSAMPLING,
    float beta_tv,
    float beta_rings,
    float rings_height,
    float alpha_rings
    );


typedef int (*sinofilter_driver_Symbol) (
    Gpu_Context* gpuctx,
    float* data
    );




typedef int (*sirtfilter_driver_Symbol) (
    Gpu_Context* self,
    float* data,
    float* SLICE,
    float DETECTOR_DUTY_RATIO,
    int DETECTOR_DUTY_OVERSAMPLING
    );

/*
typedef int (*conjgrad_driver_Symbol) (
    Gpu_Context* self,
    float* data,
    float* SLICE,
    float DETECTOR_DUTY_RATIO,
    int DETECTOR_DUTY_OVERSAMPLING,
    float Lambda
    );
*/


typedef int (*gpu_backproject_dfi_Symbol) (
    Gpu_Context * self,
    float *WORK ,
    float * SLICE
    );



typedef  int (* gpu_mainInit_Symbol )(Gpu_Context * self, float *filter) ;
typedef  int (* gpu_main_Symbol     )(Gpu_Context * self, float *WORK , float * SLICE, int precondition,
                      float    DETECTOR_DUTY_RATIO        ,
                      int    DETECTOR_DUTY_OVERSAMPLING,
                      float *fidelity,
                      int npj_offset, int forcePImultipl) ;


typedef  int (* gpu_main_2by2_Symbol     )(Gpu_Context * self, float *WORK , float * SLICE, int precondition,
                      int npj_offset, int forcePImultipl) ;


typedef float ( *gpu_lt_fit_sino_Symbol) (Gpu_Context * self, int Ng, float *gaussian_coeffs,   int   Ns , float * data, int N,
					  int *csc_II,  int * csc_JJ , float * csc_C ,int csc_nitems_atonce,
					  int *csr_II,  int * csr_JJ , float * csr_C ,int csr_nitems_atonce,
					  int SLD,  int Nsigmas,   float *sigmas,  LT_infos* lt, LT_infos *lt2s
					  );



typedef  int (* dfi_gpu_main_Symbol     )(Gpu_Context * self, float *WORK , float * SLICE) ;




typedef  int (* gpu_main_conicity_Symbol     )(Gpu_Context * self,  float * SLICE, float *WORK_perproje ,
                         int Nfirstslice, int nslices, int data_start, int nslices_data,
                         float source_distance, float detector_distance,
                         float v2x ,   float v2z,float voxel_size,
                         float SOURCE_X , float SOURCE_Z );

typedef  int (* pro_gpu_main_conicity_Symbol     )(Gpu_Context * self,  float * SLICE, float *WORK_perproje ,
                         int Nfirstslice, int nslices, int data_start, int nslices_data,
                         float source_distance, float detector_distance,
                         float v2x ,   float v2z,float voxel_size,
                         float SOURCE_X , float SOURCE_Z );





typedef  void  (* fb_dl_Symbol     )(Gpu_Context * self, float *WORK , float * SLICE, int precondition,
                     float    DETECTOR_DUTY_RATIO        ,
                     int    DETECTOR_DUTY_OVERSAMPLING,
                     float weight,
                     float *patches, int npatches, int dim_patches, int VECTORIALITY, float *lip,
             void  (*proietta_drings) (void *, float* , int , float)
                     ) ;

typedef  void  (* search_lipschitz_Symbol     )(Gpu_Context * self, float *WORK , float * SLICE, int precondition,
                     float    DETECTOR_DUTY_RATIO        ,
                     int    DETECTOR_DUTY_OVERSAMPLING,
                     float weight,
                     float *patches, int npatches, int dim_patches, int VECTORIALITY, float *lip,
             void  (*proietta_drings) (void *, float* , int , float)
                     ) ;

typedef  float (* tv_denoising_fistagpu_Symbol     )(Gpu_Context * self,int dim0,int dim1,float  *img,float  *result , float weight,
                             float eps , int n_iter_max,   int check_gap_frequency ) ;

typedef  float (* tv_denoising_patches_OMP_Symbol     )(Gpu_Context * self,int dim0,int dim1,float  *img,float  *result , float weight,
                            float *patches, int npatches, int dim_patches) ;

typedef  float (* tv_denoising_patches_L1_Symbol     )(Gpu_Context * self,int dim0,int dim1,float  *img,float  *result , float weight,
                               float *patches, int npatches, int dim_patches, int N_ITERS_DENOISING,
                               int vectoriality) ;


typedef  void (*   nonlocalmeans_Symbol   ) (    Gpu_Context * self,
               float *result,
               float *image,
               int H,
               int W,
               float bruit,
               float mix
               );




void *  AllocPinned(Gpu_Context * self,size_t size );
void FreePinned(Gpu_Context * self,void *ptr);


typedef    void *(*AllocPinned_Symbol)(Gpu_Context * self,size_t size );
typedef    void  (* FreePinned_Symbol)(Gpu_Context * self,void * SLICE  );


struct Gpu_Context_struct {
  int inuse ;
  gpu_lt_fit_sino_Symbol  gpu_lt_fit_sino;
  gpu_main_Symbol  gpu_main;
  gpu_main_2by2_Symbol  gpu_main_2by2;
  dfi_gpu_main_Symbol  dfi_gpu_main;
  gpu_main_conicity_Symbol  gpu_main_conicity;
  pro_gpu_main_conicity_Symbol  pro_gpu_main_conicity;
  fb_dl_Symbol  fb_dl;
  gpu_backproject_Symbol gpu_backproj;
  cp_driver_Symbol cp_driver;
//  conjgrad_driver_Symbol conjgrad_driver;
  wavelets_driver_Symbol wavelets_driver;
  sinofilter_driver_Symbol sinofilter_driver;
  sirtfilter_driver_Symbol sirtfilter_driver;
  gpu_backproject_dfi_Symbol gpu_backproj_dfi;
  search_lipschitz_Symbol  search_lipschitz;

  tv_denoising_fistagpu_Symbol tv_denoising_fistagpu;
  tv_denoising_patches_L1_Symbol tv_denoising_patches_L1;
  tv_denoising_patches_OMP_Symbol tv_denoising_patches_OMP;
  nonlocalmeans_Symbol  nonlocalmeans ;
  gpu_mainInit_Symbol  gpu_mainInit;
  gpu_project_Symbol gpu_project;

  AllocPinned_Symbol AllocPinned;
  FreePinned_Symbol FreePinned;

  // set by user
  int dim_fft;
  int nprojs_span;
  int num_bins ;

  // set by init
  fcomplex  *dev_iWork;
  fcomplex  *dev_iWork_copy;
  float *dev_Filter;
  int planR2C,planC2R, planC2C_2by2;
  void  * a_Proje_voidptr;
  void  * a_cProje_voidptr;
  float *dev_Work_perproje;
  float *d_work ;
  float* dev_rWork;

  void *gpu_streams;

  // set by user
  float overlapping;
  float pente_zone ;
  float flat_zone ;
  float prof_shift;
  float prof_fact ;
  float *axis_position_s  ;
  float  axis_position   ;
  int fai360;

  // set by init
  int NblocchiPerLinea32;
  int NblocchiPerColonna32   ;
  int NblocchiPerLinea;
  int NblocchiPerColonna   ;
  int dimrecx, dimrecy;
  float * d_SLICE ;
  float * d_cos_s, * d_sin_s,  * d_axis_s ;
  
  // set by user
  int num_x ;
  int num_y ;
  float * cos_s            ;
  float * sin_s            ;

  float gpu_offset_x;
  float gpu_offset_y;

  int FBFILTER ;
  int DFI_NOFVALUES  ;
  int DFI_KERNEL_SIZE ;
  int DFI_OVERSAMPLING_RATE ;
  int DFI_R2C_MODE ;

  int DFP_NOFVALUES  ;
  int DFP_KERNEL_SIZE ;
  int DFP_OVERSAMPLING_RATE ;
  int USE_DFI;
  int USE_DFP;


  int DO_PRECONDITION;
  PRECOND_PARAMS_DL precond_params_dl;


  void *gpuctx ;
  int MYGPU;

  int STEPFORPATCHES ;

  float fb_dl_Lipschitz;
  int VECTORIALITY ;
  float * angles_per_proj;
  float *axis_corrections ;


  
  int JOSEPHNOCLIP ;
  int ITERATIVE_CORRECTIONS;
  float SMOOTH_PENALTY_PARAM;
  float LINE_SEARCH_INIT;
  int PENALTY_TYPE;
  int OPTIM_ALGORITHM;
  float peso_overlap;

  int CONICITY ;
  int CONICITY_FAN ;
  float FAN_FACTOR;
  float SOURCE_DISTANCE ;
  float DETECTOR_DISTANCE ;
  float SOURCE_X        ;
  float SOURCE_Y      ;

  float tilt_psi;
  float tilt_tilt;

  int * CONICITY_MARGIN_UP ;
  int * CONICITY_MARGIN_DOWN ;

  float ITER_RING_HEIGHT ;
  float ITER_RING_SIZE ;
  float RING_ALPHA;
  float LIPSCHITZFACTOR;
  int LIPSCHITZ_ITERATIONS;
  float RING_BETA;
  int NUMBER_OF_RINGS;

  //Wavelets
  int W_LEVELS;
  int W_CYCLE_SPIN;
  float W_FISTA_PARAM;
  char* W_WNAME;
  int W_SWT;
  int W_NORMALIZE;
  int W_THRESHOLD_APPCOEFFS;
  int W_THRESHOLD_COUSINS;

  // Fourier-Wavelet (Munch Et Al) filter
  int FW_LEVELS;
  float FW_SIGMA;
  char* FW_WNAME;
  float FW_SCALING;
  float FW_FILTERPARAM;
  cufftHandle* fw_plans1d_r2c;
  cufftHandle* fw_plans1d_c2r;
  cufftHandle fw_plan1d_r2c;
  cufftHandle fw_plan1d_c2r;
  char fw_plans_ok;

  //
  int ESTIMATE_BETA;
  //
  int do_ignore_projections;
  int* ignore_angles;

  // Fluo
  int FLUO_SINO;
  int FLUO_ITERS;
  float* d_Sigma;
  int FLUO_step;

  // Sirt-filter
  int SF_ITERATIONS;
  char* SF_SAVEDIR;
  float SF_LAMBDA;

  int DATA_IS_FILTERED;

  // TV+L2
  float BETA_L2;

  // Positivity constraint
  int ITER_POS_CONSTRAINT;



  DFI_params dfi_params ;
  DFP_params dfp_params ;

  void * void_ccspace_ptr;
  int num_slice;
  int verbosity;


  float DZPERPROJ;
  int *tot_proj_num_list;
  int numpjs;

  
  LT_infos * lt_infos_coarse;
  int lt_planr2c, lt_planc2r; 

  
};


void  proietta_drings(void * void_ccspace_ptr, float* rings_tmp, int doppio , float ) ;


typedef struct Gpu_pag_Context_struct Gpu_pag_Context;


void  gpu_pagCtxDestroy (Gpu_pag_Context * self);
void  gpu_pagCtxCreate (Gpu_pag_Context * self);
void  gpu_pag( Gpu_pag_Context *   self , float * auxbuffer );
void  gpu_pagFree( Gpu_pag_Context *   self );
void  gpu_pagInit( Gpu_pag_Context *   self ) ;






typedef  int (* gpu_pagInit_Symbol )(Gpu_pag_Context * self) ;
typedef  int (* gpu_pag_Symbol     )(Gpu_pag_Context * self, float * mataux ) ;
typedef  int (* gpu_pagCtxCreate_Symbol     )(Gpu_pag_Context * self) ;
typedef  int (* gpu_pagCtxDestroy_Symbol     )(Gpu_pag_Context * self) ;
typedef  int (* gpu_pagFree_Symbol )(Gpu_pag_Context * self) ;


struct Gpu_pag_Context_struct {

  int inuse ;

  // set by user inizio
  void *gpuctx ;
  int MYGPU;
  gpu_pagInit_Symbol  gpu_pagInit;
  gpu_pag_Symbol  gpu_pag ;
  gpu_pagCtxCreate_Symbol  gpu_pagCtxCreate ;
  gpu_pagCtxDestroy_Symbol  gpu_pagCtxDestroy ;
  gpu_pagFree_Symbol  gpu_pagFree;

  // set by user each chunck
  int size_pa0;
  int size_pa1;
  float * kernelbuffer;
  // set by init
  void * d_fftwork, * d_kernelbuffer;
  void * FFTplan_ptr;

  int gpu_pagCtx_initialised ; 
  int NBunches_todo;
  int dones[64000];
  
}   ;


// ---------------------------------------------------------
typedef struct Gpu_med_Context_struct Gpu_med_Context;

void  gpu_medCtxCreate (Gpu_med_Context * self);
void  gpu_medCtxDestroy (Gpu_med_Context * self);
void  gpu_med( Gpu_med_Context * self,int ny, int nx, float *data, float *result, int hwy, int hwx, float threshold);

typedef  void (* gpu_medCtxCreate_Symbol )(Gpu_med_Context * self) ;
typedef  void (* gpu_medCtxDestroy_Symbol )(Gpu_med_Context * self) ;
typedef  void (* gpu_med_Symbol     )(Gpu_med_Context * self,int ny, int nx, float *data, float *result, int hwy, int hwx, float threshold ) ;

struct Gpu_med_Context_struct {

  int inuse ;

  gpu_medCtxCreate_Symbol  gpu_medCtxCreate ;
  gpu_medCtxDestroy_Symbol  gpu_medCtxDestroy ;
  gpu_med_Symbol          gpu_med          ;

  // set by user inizio
  void *gpuctx ;
  int MYGPU;
}   ;




// #ifndef FROMCU
 /* int gpu_main(int num_y, int num_x,  float * SLICE, int num_proj, int num_bins, float *WORK_perproje ,  */
 /*           float axis_position, float * axis_position_s, float * cos_s, float *sin_s , float gpu_offset_x, float gpu_offset_y) ; */

/* ------------*/
/* CCspace     */
struct CCspace_struct{
  CCarraylist* rawdatatokens ;
  CCarraylist* ff_rawdatatokens ;
  CCarraylist* datatokens   ;
  CCarraylist*  transposeddatatokens  ;

  float *background;
  float *axis_corrections;
  float *axis_correctionsL;
  
  float *ffcorr;
  CCreading_infos reading_infos;
  int iproc;
  int nprocs;
  hid_t shared_h5id ; 
  hid_t shared_h5DSid ; 
  hid_t shared_h5DSid_current ; 
  
  int *snXtoken;
  int **ff_read_status;
  sem_t ff_sem;
  sem_t fftw_sem         ;
  sem_t slicesrequest_sem;
  sem_t savgol_sem;
  sem_t fbp_sem;
  sem_t gpustat_sem;
  sem_t gpustat_pag_sem;
  sem_t gpudones_pag_sem;
  sem_t gpustat_med_sem;
  sem_t islicetracker_sem;
  sem_t proupdate_sem;
  sem_t filereading_sem;
  sem_t hdf5sequence_sem;

  
  Gpu_Context * gpu_context;
  Gpu_pag_Context * gpu_pag_context;
  Gpu_med_Context * gpu_med_context;

  // double FF
  float* ff2_localmean_threads;
  float* ff2_localmean_current_thread;
  float* ff2_localmean_process;
  float* ff2_globalmean;
  int* ff2_nels_threads;
  int* ff2_nels_process;
  int ff2_status;
  int ff2_done;
  // ---


  void *sharedHandle;


  SlicesRequest* slicesrequest ;
  int *packet_completion ;
  int gpu_is_apriori_usable;


  Cparameters params;
  float * Coeff_fil;

  sem_t minmax_sem ;
  double aMin, aMax;
  long int *histovalues;
  int invertedsteam_output_created;
  int * packet_has_flown;

  void * allocations[MAXALLOCATIONS];
  long int allocations_sizes[MAXALLOCATIONS];
  
#ifndef FROMCU
  FBP_prec fbp_precalculated;
#endif
} ;
#ifndef FROMCU

void nnfbp_train(CCspace *self,int dim_fft,int num_bins, float *data_orig, int dimslice,float *SLICEcalc,
		 float * SINOGRAMMA, float **WORK, float *WORKbis, float *dumf, fcomplex *dumfC, float *WORK_perproje,
		 float *OVERSAMPLE, int oversampling, int ncpus,float cpu_offset_x,float  cpu_offset_y, int Nfirstslice, int islice,
		 char *nomeout) ;



void nnfbp_reconstruct(CCspace *self,int dim_fft,int num_bins, float *data_orig, int dimslice,float *SLICEcalc,float *SLICE,
		       float * SINOGRAMMA, float **WORK, float *WORKbis, float *dumf, fcomplex *dumfC, float *WORK_perproje, float *OVERSAMPLE,
		       int oversampling, int ncpus,float cpu_offset_x,float  cpu_offset_y);




void   CCspace_getSaturations(CCspace * self, double aMin,double aMax,
                  double *sat1,double *sat2,double *Sat1,double *Sat2  );

void CCspace_initialise     (CCspace * self  ) ;
void CCspace_addRawDataSpace(CCspace * self , float *token );
void CCspace_add2DataSpace(CCspace * self , float *token,const char *key ) ;
void CCspace_read_chunk    (CCspace * self , int sn, int  ntok , int npbunches, int rotation_vertical, int binning,
                int reduced_case, int red_start, int red_end, int red_scope);
void CCspace_tranpose_chunk    (CCspace * self , int sn, int  ntoktreated, int  ntoktransposed , int npbunches, int STEAM_DIRECTION);
void CCspace_dispense_chunk   ( CCspace * self  , int sn,  int  ntoktransposed , int   npbunches  );

void CCspace_preprocess_chunk(CCspace *  self ,int sn,int  ntok, int  ntokt , int npbunches,int ncpus, int doff2 );
void CCspace_InterleavedReadingPreProcessing_chunk(CCspace *  self , int sn, int  ntok, int  ntokt , int npbunches, int ncpus , int do_ff2);

void CCspace_reconstruct(CCspace *  self ,int sn, int npbunches,int ncpus , int STEAM_DIRECTION);
void CCspace_reconstructSHARED(CCspace *  self ,int sn, int npbunches,int ntoktransposed ,int ncpus , int STEAM_DIRECTION);

void CCspace_prepare_concurrent_ff_reading(CCspace * self);
void CCspace_ffstatus_dealloc(CCspace * self);
void   CCspace_end(CCspace * self ) ;

void CCspace_set_nchunks( CCspace * self   ,int   nchunks ) ;




/* for ff_read_status */
#define NOT_ACQUIRED    0
#define IN_ACQUISITION  1
#define ACQUIRED        2

void       Paganin(CCspace *  self,  float * Rawptr,
		   int Pos0, int  Pos1, int  Size0, int  Size1,
		   int pos0, int  pos1, int  size0, int  size1,
		   Cparameters *P , int ncpus,   sem_t* fftw_sem,
		   int pstart, int pend, int poffset,
		   int p_num_offset,
		   int mystart,
		   int npbunches,
		   int ibunch) ;

/* for ccd_filter */
#define CCD_FILTER_NONE_ID       0
#define CCD_Filter_ID            1

typedef struct {
  float threshold;
} CCD_Filter_PARA_struct  ;

/* for ring_filter */
#define RING_FILTER_NONE_ID       0
#define RING_Filter_ID            1
#define RING_Filter_SG_ID            2
#define RING_Filter_THRESHOLDED_ID  3


typedef struct {
  float *FILTER;
  float threshold;
} RING_Filter_PARA_struct  ;

typedef struct {
  int LF ;
  int I_Slope;
  int Lfen;
  float Eps1;
  float Eps2;
  float RProt;
} RING_Filter_SG_PARA_struct  ;



void    CCspace_RING_Filter_implementation(CCspace *   self, float *data,
                       RING_Filter_PARA_struct*  RING_FILTER_PARA ,
                       int nslices,
                       int nprojs_span,
                       int size1 ,
                       int ncpus,int *itrack);


void    CCspace_RING_Filter_SG_implementation(CCspace *   self, float *data,
                          RING_Filter_SG_PARA_struct*  RING_FILTER_PARA ,
                          int nslices,
                          int nprojs_span,
                          int size1 ,
                          int ncpus,int *itrack) ;



void CCspace_Sino_2_Slice( CCspace *    self, float * dataA, int nslices, int nslices_data,int Nfirstslice, int ncpus , int data_start, int steam_direction, int npj_offset, int ibunch, int npbunches);

void CCspace_Sino_2_Slice_conicity( CCspace *    self, float * dataA, int nslices, int nslices_data,
                    int Nfirstslice, int ncpus , int data_start, int steam_direction);



void pro_conic_driver(  CCspace *   self, float *  SINOGRAMMA , float *SLICE, int  nslices,
             int nslices_data,
             int Nfirstslice,
             int data_start,
             float cpu_offset_x,
             float cpu_offset_y
            );

void conic_driver(  CCspace *   self,  float *  data, float *SLICE, int  nslices,
             int nslices_data,
             int Nfirstslice,
             int data_start,
             float *dumf,
             fcomplex *dumfC,
             float *WORK_perproje,
             float **WORK,
             float *WORKbis,
             float cpu_offset_x,
             float cpu_offset_y
            ) ;




#endif

void  CCD_Filter_Implementation(float * Tptr, float * Rawptr,
                int Size0, int Size1 ,
                float threshold,
                       int ncpus );


void SAVGOL(float *Coeff_filtre, int LF  );



float FindNoise(CCspace *  self,float *SLICE_a, int CALM_ZONE_LEN, int size_pa0, int  size_pa1 ) ;

void  nlm_driver(CCspace *self,int dim0,int dim1,float  *img,float  *result , float bruit) ;

float rec_driver(CCspace *self,  float **WORK,float *WORKbis,float *SLICEcalc, int num_bins, int dim_fft,
         float * dumf, fcomplex *dumfC , float *WORK_perproje,float *OVERSAMPLE,int oversampling,
         float * data, int ncpus,float cpu_offset_x,float  cpu_offset_y, int preco, int is_tomolocal,
		 float *fidelity, int npj_offset, int do_2by2);

// Typedef for plugins
typedef float (*rec_driver_func)(CCspace*,  float **,float *,float *, int, int,
         float * , fcomplex * , float *,float *,int ,
         float * , int ,float ,float  , int , int ,
         float *) ;


void pro_driver(
         CCspace * self,int num_bins, float *angles_per_proj,
         float * SINOGRAMMA, float *SLICE, int dimslice, float cpu_offset_x, float cpu_offset_y );

// Typedef for plugins
typedef void (*pro_driver_func)(CCspace *,int, float *, float *, float *, int , float , float );

void fb_dl_driver(CCspace *self, float *data, int num_bins,
      float *SLICEcalc,    int do_precondition, int VECTORIALITY, float w, int num_slice) ;


void correct_slice(CCspace *self,float *SLICE,float *SLICEcalc);

void prepare_correction(CCspace *self, float * data, float * data_origin, float *SINOGRAMMA,int num_bins ) ;

void passeggiaSINO(CCspace *self, float **WORK, float **WORKbis, int  num_bins,  int dim_fft,
           float **dumf, fcomplex **dumfC , float * WORK_perproje, float **OVERSAMPLE,
           int oversampling, int  ncpus   , float cpu_offset_x, float cpu_offset_y,
           float *angles_per_proj,  int dimslice,
           float *SLICE, float *data, float * SINOGRAMMA, int solution, int nsteps);

void passeggiaSLICE(CCspace *self, float **WORK, float **WORKbis, int  num_bins,  int dim_fft,
            float **dumf, fcomplex **dumfC , float * WORK_perproje, float **OVERSAMPLE,
            int oversampling, int  ncpus   , float cpu_offset_x, float cpu_offset_y,
            float *angles_per_proj,  int dimslice,
            float *SLICE, float *data, float * SINOGRAMMA, int nsteps);

float denoising_driver(CCspace *self,int dim0,int dim1,float  *img,float  *result , float weight )  ;
void rotational2zero(CCspace *self,float *SLICE,float *SLICEres) ;

float * medianX(  float *A   , int  dim3 , int  dim2 , int  dim1   , int nproc );
void  CCD_Filter_Implementation(float * Tptr, float * Rawptr,
				int Size0, int Size1 ,
				float threshold, int ncpus) ;