/tomo/pyhst

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

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


#define _FILE_OFFSET_BITS  64


#include"Python.h"


#ifndef _LARGEFILE_SOURCE
#define _LARGEFILE_SOURCE
#endif

#include <stdlib.h>
#include <stdio.h>
#include<string.h>
#include <math.h>

#include<semaphore.h>

#include <glib.h>

#include"structmember.h"

#include "numpy/arrayobject.h"


#include "hst.h"
#include "debug.h"

/*
#ifdef LINUX
 #define _FILE_OFFSET_BITS  64
#endif
*/



/*
 * The error object to expose 
 */

static PyObject *ErrorObject;
#define onError(message)\
  { PyErr_SetString(ErrorObject, message); return NULL;}



typedef struct {
  PyObject_HEAD
  HSTSetup setup;
  HSTContextPtr hst;
  
  FILE *output;              	/*  Identifier for file outp   */ 
#ifdef HST_USE_FASTWRITER  
  fastwriter_t *fw;
#endif /* HST_USE_FASTWRITER */
  PyObject *logger;		/* Message logger object (logging) */
  void *ibuffer1, *ibuffer2;	/* Pinned storage for the sinograms */
  PyObject *SINO, *SINO1, *SINO2;
  PyObject *SINOGRAMS;          /* Storage for the sinogram. Dimensions:  */
                                /* (num_bins,num_projections,dim_slices)  */
  PyObject *axis_corrections;   /*  Dimensions:  */
                                /* (num_projections)  */
  int astra_scaling;
  int BICUBIC;			/* ignored at the moment */

  int do_custom_angles;
  float * angles_data;

  PyObject *FilterFunct  ; /* object coming from python adn defining the filter */
  int FilterFunctHasRamp;
  PyObject * FilterOwnerInstance   ;
#ifdef HW_USE_PARALLEL_IO
  PyObject ** Projections;
#endif /* HW_USE_PARALLEL_IO */

#ifdef PYHST_MEASURE_TIMINGS
  double recon_timer;
  double comp_timer;
  double io_timer;
  PyObject *recon_info;	 //!< Information about used reconstructors (PyListObject)
#endif /* PYHST_MEASURE_TIMINGS */

} PyHST;

staticforward PyTypeObject PyHSTtype;

#define is_PyHST(v) ((v)->ob_type == &tPyHSType )


static void 
PyHST_log_handler(const gchar *log_domain, GLogLevelFlags log_levels, const gchar *message, PyObject *logger) {
//  PyObject *arglist;
//  arglist = Py_BuildValue("(s)", message);

  if (log_levels&G_LOG_LEVEL_ERROR) PyObject_CallMethod(logger, "critical", "(s)", message);
  else if (log_levels&G_LOG_LEVEL_CRITICAL) PyObject_CallMethod(logger, "error", "(s)", message);
  else if (log_levels&G_LOG_LEVEL_WARNING) PyObject_CallMethod(logger, "warning", "(s)", message);
  else if (log_levels&G_LOG_LEVEL_MESSAGE) PyObject_CallMethod(logger, "info", "(s)", message);
  else if (log_levels&G_LOG_LEVEL_INFO) PyObject_CallMethod(logger, "info", "(s)", message);
  else PyObject_CallMethod(logger, "debug", "(s)", message);

//  Py_DECREF(arglist);
}


/*
 * kind of destructor called by garbage collector
 */
static void
PyHST_dealloc(PyHST *self)
{
#ifdef PYHST_MEASURE_TIMINGS
  if (self->recon_info) Py_DECREF(self->recon_info);
#endif /* PYHST_MEASURE_TIMINGS */
  
  if (self->hst) hst_destroy_context(self->hst);
  

#ifdef HST_USE_FASTWRITER  
  if (self->fw) {
    fastwriter_close(self->fw);
    fastwriter_destroy(self->fw);
  }
#endif /* HST_USE_FASTWRITER */
  if (self->output) fclose(self->output);


  if (self->ibuffer2) {
    if (self->SINO2) Py_XDECREF(self->SINO2);
    hst_pinned_free(self->ibuffer2);
  }
  
  if (self->ibuffer1) {
    if (self->SINO1) Py_XDECREF(self->SINO1);
    hst_pinned_free(self->ibuffer1);
  } 
  
  if (self->SINO) Py_XDECREF(self->SINO);
  if (self->axis_corrections) Py_XDECREF(self->axis_corrections);
  if (self->logger) Py_XDECREF(self->logger);

#ifdef HW_USE_PARALLEL_IO
  if (self->Projections) free(self->Projections);
#endif /* HW_USE_PARALLEL_IO */
 
  if(self && 0 )
    free(self);
}

/*
 * creation and initialization of a new 
 */

static char PyHST_new_doc[]=""\
"/*"\
" * creation and initialization of a new PyHST. Lots of arguments:"\
"  if(!PyArg_ParseTuple(args,\"siiiiifffiiiOii:PyHST_new\", &outputname, &(res->oversampling),&(res->start_x),&(res->start_y),"\
"		       &(res->num_x),&(res->num_y),"\
"		       &(res->axis_position),&(res->angle_offset),&(res->angle_increment),"\
"			   &(res->no_filtering),"\
"		       &(res->num_bins),&(res->num_projections),"\
"		       &(res->SINOGRAMS),"\
"                       &(res->axis_corrections),"\
"                       &(res->BICUBIC),"\
"                       &(res->SUMRULE)"\
" ))"\
"    return NULL;"\
"   */";

static PyObject *
PyHST_new( PyObject *self, PyObject *args)
     /* return a new instance of edfobject */
{
  int err = 0;
  PyHST* res;
  char * outputname;
  char * method;

  PyObject *custom_angles  ;
  PyObject *result;
  PyObject *normalise;
  int log_level;

#ifdef PYHST_MEASURE_TIMINGS
  int i, j;
  HSTReconstructorConstContextPtr *reconstructors;
  HSTConstString recon_title;
  HSTConstString *timer_names;
  double *timer_values;
  
  PyObject *recon_dict;
  PyObject *timers_list;
  PyObject *pystr;
#endif /* PYHST_MEASURE_TIMINGS */

  res = (PyHST*) PyObject_NEW(PyHST, &PyHSTtype);

  res->hst = NULL;
  res->output = NULL;
#ifdef HST_USE_FASTWRITER  
  res->fw = NULL;
#endif /* HST_USE_FASTWRITER */
  res->ibuffer1 = NULL;
  res->ibuffer2 = NULL;
  
#ifdef PYHST_MEASURE_TIMINGS
  res->recon_info = NULL;
#endif /* PYHST_MEASURE_TIMINGS */

#ifdef PYHST_ASTRA_SCALING
    res->astra_scaling = 1;
#else
    res->astra_scaling = 0;
#endif

  if(!PyArg_ParseTuple(args,"OissiiiiiiiifffiiiOOiiOfiO:PyHST_new", 
		       &(res->logger), 
		       &(res->setup.max_slices),
		       &outputname, 
		       &method,
		       &(res->setup.fft_oversampling),
		       &(res->setup.dfi_kernel_size),
		       &(res->setup.dfi_kernel_points),
		       &(res->setup.oversampling),
		       &(res->setup.start_x),
		       &(res->setup.start_y),
		       &(res->setup.num_x),
		       &(res->setup.num_y),
		       &(res->setup.axis_position),
		       &(res->setup.angle_offset),
		       &(res->setup.angle_increment),
		       &(res->setup.no_filtering),
		       &(res->setup.num_bins),
		       &(res->setup.num_projections),
		       &(res->SINO),
		       &(res->axis_corrections),
		       &(res->BICUBIC),
		       &(res->setup.sum_rule),
		       &custom_angles,
		       &(res->setup.pente_zone),
		       &(res->setup.zerooffmask),
		       &normalise
  )) {
    res->SINO = NULL;
    res->logger = NULL;
    res->axis_corrections = NULL;

    Py_DECREF(res);
    return NULL;
  }

  float *norm = (float*) ((PyArrayObject *) normalise)->data;
  res->setup.norm = norm[0];

  Py_INCREF(res->logger);
  Py_INCREF(res->axis_corrections);
  Py_INCREF(res->SINO);
  
	// Configuring the logging subsystem 
  result = PyObject_CallMethod(res->logger, "getEffectiveLevel", "()");
  if (result) {
     if (PyInt_Check(result)) log_level = PyInt_AsLong(result);
     else log_level = 0;
//     ||(!PyArg_ParseTuple(result, "i", &log_level))) log_level = 0;
     Py_DECREF(result);
  } else log_level = 0;
  
  pyhst_configure_logger(log_level, (GLogFunc)PyHST_log_handler, res->logger);

#ifdef HW_USE_PARALLEL_IO
 res->Projections = (PyObject**)malloc(res->setup.num_projections * sizeof(PyObject*));
 if (!res->Projections) onError("Memory allocation");
#endif /* HW_USE_PARALLEL_IO */

	// Processing various options
  if(PyArray_Check( custom_angles)) {
    res->do_custom_angles = 1;
    if( ((PyArrayObject *)custom_angles  )  ->descr->type_num != PyArray_FLOAT ) onError("arg custom_angles  is not an array of float " ) ;
    res->angles_data = (float *) ((PyArrayObject *)custom_angles  )->data;
    if(  ((PyArrayObject *)custom_angles  )->dimensions[0] <  res->setup.num_projections         ) {
      onError("arg custom_angles  has not the good lenght " ) ;
    }
  } else {
    res->do_custom_angles = 0;
    res->angles_data =0;
  }
  
  res->FilterFunct=0;
  res->FilterOwnerInstance=0;

  if (res->SINO->ob_type == &PyList_Type) {
    int x, y, z, double_buffer = 1;
    
    PyObject *tuple = PyList_AsTuple(res->SINO);
    
    switch (((PyListObject*)(res->SINO))->ob_size) {
	case 3:
	    if ((tuple)&&(PyArg_ParseTuple(tuple, "iii", &x, &y, &z))) break;
	    else err = 1;
	break;
	case 4:
	    if ((tuple)&&(PyArg_ParseTuple(tuple, "iiii", &double_buffer, &x, &y, &z))) break;
	    else err = 1;
	break;
	default:
	    err = 1;
    }
    

    //pyhst_warning("double buffer: %i\n", double_buffer);

    if (err) {	    
	pyhst_error("ERROR : parsing buffer size");
	Py_DECREF(res);
	return NULL;
    } else {
	npy_intp dims[3] = {x, y, z};
	res->ibuffer1 = hst_pinned_malloc(x * y * z * sizeof(float), HST_PINNED_R);
	if (res->ibuffer1) {
	    res->SINO1 = PyArray_SimpleNewFromData(3, dims, PyArray_FLOAT, res->ibuffer1);
	    if (!res->SINO1) err = 1;
	} else err = 1;
	
	if (double_buffer > 1) {
	    res->ibuffer2 = hst_pinned_malloc(x * y * z * sizeof(float), HST_PINNED_R);
	    if (res->ibuffer2) {
		res->SINO2 = PyArray_SimpleNewFromData(3, dims, PyArray_FLOAT, res->ibuffer2);
		if (!res->SINO2) err = 1;
	    } else err = 1;
	}
	
        if (err) {
	    pyhst_error("ERROR : allocating pinned memory");
	    Py_DECREF(res);
	    return NULL;
	} 
    }
  
    res->SINOGRAMS = res->SINO1;
  } else {
    res->SINOGRAMS = res->SINO;
  }

  if (!strcasecmp(method, "FBP")) {
    res->setup.method = HST_METHOD_FBP;
  } else if (!strcasecmp(method, "DFI")) {
    res->setup.method = HST_METHOD_DFI;
  } else {
      pyhst_error("ERROR : could not open file %s", outputname);
      Py_DECREF(res);
      return NULL;
  }

  if(strlen(outputname)>0){  
#ifdef HST_USE_FASTWRITER  
    res->fw = fastwriter_init(outputname, FASTWRITER_FLAGS_OVERWRITE);
    if (res->fw) {
        fastwriter_set_buffer_size(res->fw, HST_USE_FASTWRITER * 1024ll * 1024ll);
        if (fastwriter_open(res->fw, outputname, FASTWRITER_FLAGS_OVERWRITE)) {
            fastwriter_destroy(res->fw);
            res->fw = NULL;
        }
    }
    if (!res->fw) {
#else /* HST_USE_FASTWRITER */
    res->output=fopen(outputname,"w");
    if(!res->output) {
#endif /* HST_USE_FASTWRITER */
      pyhst_error("ERROR : could not open file %s", outputname);
      Py_DECREF(res);
      return NULL;
    }
  } 
  
  res->hst = hst_create_context();
  if ((!res->hst)||(hst_init_context(res->hst, &res->setup, res->do_custom_angles?res->angles_data:NULL, (float*) ((PyArrayObject *) res->axis_corrections)->data))) {
    Py_DECREF(res);
    return NULL;
  }

#ifdef HST_USE_FASTWRITER  
  hst_set_fastwriter(res->hst, res->fw);
#else /* HST_USE_FASTWRITER */
  hst_set_output_file(res->hst, res->output);
#endif /* HST_USE_FASTWRITER */

#ifdef PYHST_MEASURE_TIMINGS
  res->recon_timer = 0;

  reconstructors = hst_get_configured_reconstructors(res->hst);
  for (i = 0; reconstructors[i]; i++);
  
  res->recon_info = PyList_New(i);//(PyArrayObject*)PyArray_FromDims(1, &i, PyArray_OBJECT);
  if (res->recon_info) {
     for (i = 0; reconstructors[i]; i++) {
        recon_dict = PyDict_New();
	if (recon_dict) {
	    recon_title = hst_reconstructor_get_title(reconstructors[i]);
	    pystr = PyString_FromString(recon_title?recon_title:"Reconstructor");
	    if (pystr) {
		PyDict_SetItemString(recon_dict, "title", pystr);
		Py_DECREF(pystr);
	    }
	    
	    timer_names = hst_reconstructor_get_timers(reconstructors[i], NULL);
	    if (timer_names) {
		for (j = 0; timer_names[j]; j++);
		if (j) {
		    npy_intp dims = j;
		    timers_list = PyArray_SimpleNew(1, &dims, PyArray_DOUBLE);
		    if (timers_list) {
			timer_values = (double *) ((PyArrayObject*)timers_list)->data; 
			for (j = 0; timer_names[j]; j++) {
			    timer_values[j] = 0;
			}
			PyDict_SetItemString(recon_dict, "timers", timers_list);
			Py_DECREF(timers_list);
		    }
		    
		    if (timers_list) timers_list = PyList_New(j);
		    
		    if (timers_list) {
			for (j = 0; timer_names[j]; j++) {
			    pystr = PyString_FromString(timer_names[j]);
			    if (pystr) {
				    // A reference is stollen by this function
				PyList_SetItem(timers_list, j, pystr);
			    } else {
				Py_INCREF(Py_None);
				PyList_SetItem(timers_list, j, Py_None);
			    }
			}
			PyDict_SetItemString(recon_dict, "timer_names", timers_list);
			Py_DECREF(timers_list);
		    }
		    
		    timers_list = PyInt_FromLong(0);
		    if (timers_list) {
			PyDict_SetItemString(recon_dict, "slices", timers_list);
			Py_DECREF(timers_list);
		    }
		}
	    }
	    
	    PyList_SetItem(res->recon_info, i, recon_dict);
	    // SetItem steals a reference, we do not need Py_DECREF(recon_dict);
	} else {
	    Py_INCREF(Py_None);
	    PyList_SetItem(res->recon_info, i, Py_None);
	}
    }
  } else {
    Py_INCREF(Py_None);
    res->recon_info = Py_None;
  }
#endif /* PYHST_MEASURE_TIMINGS */

  return (PyObject*) res;
}


static PyObject *
PyHST_close( PyObject *self, PyObject *args)
     /* return a new instance of edfobject */
{
#ifdef PYHST_MEASURE_TIMINGS
  PyHST *ctx;
  GTimer *timer;
  double time;
   
  ctx = (PyHST *)self;
  timer = g_timer_new();
#endif /* PYHST_MEASURE_TIMINGS */

  if (((PyHST *) self)->output) {
     hst_set_output_file(((PyHST*)self)->hst, NULL);
#ifdef HST_USE_FASTWRITER  
    fastwriter_close(((PyHST*)self)->fw);
    fastwriter_destroy(((PyHST*)self)->fw);
#else /* HST_USE_FASTWRITER */
     fclose( ( (PyHST *) self)->output);
     ((PyHST *) self)->output = NULL;
#endif /* HST_USE_FASTWRITER */
  }

#ifdef PYHST_MEASURE_TIMINGS
  g_timer_stop(timer);
  time = g_timer_elapsed(timer, NULL);
  ctx->io_timer += time;
  ctx->recon_timer += time;
//  ctx->comp_timer += time;
  g_timer_destroy(timer);
#endif /* PYHST_MEASURE_TIMINGS */

  Py_INCREF( Py_None);
  return Py_None;

}


void PyHST_setFilter(PyHST *self, int dim_fft, float *FILTER) {
    int j;

    PyObject *FilterFunct;
    PyObject *FilterOwnerInstance;
    int FilterFunctHasRamp;

    PyObject * pyargs, *pyres;

    FilterFunct = self->FilterFunct;
    FilterOwnerInstance = self->FilterOwnerInstance; 
    FilterFunctHasRamp = self->FilterFunctHasRamp;

    if(FilterFunctHasRamp==0) {
      pyargs = Py_BuildValue("(O,f)", FilterOwnerInstance,0.5);
      pyres  = PyEval_CallObject(FilterFunct, pyargs);
    } else {
      pyargs = Py_BuildValue("(O,l,l)", FilterOwnerInstance,dim_fft/2, dim_fft);
      pyres  = PyEval_CallObject(FilterFunct, pyargs);
    }

    if (pyres==NULL) pyhst_fail("Filter have failed");

    if(FilterFunctHasRamp==0) {
      FILTER[1]= PyFloat_AsDouble(pyres)/dim_fft; ;             /* This is the real part of the highest frequency */
    } else {
      FILTER[1]= PyFloat_AsDouble(pyres)*2.0/dim_fft; ;             /* This is the real part of the highest frequency */
    }

    Py_DECREF(pyargs);
    Py_DECREF(pyres );

    if(FilterFunctHasRamp) {
      pyargs = Py_BuildValue("(O,i,i)", FilterOwnerInstance, 0, dim_fft);
      pyres  = PyEval_CallObject(FilterFunct, pyargs);
      FILTER[0]= PyFloat_AsDouble(pyres)*2.0/dim_fft; ;             /* This is the real part of the highest frequency */ 
    } else {
      FILTER[0] = 0;
    }

    for(j=1;j<dim_fft / 2; j++) {
      if(FilterFunctHasRamp==0) {

	pyargs = Py_BuildValue("(O,f)",FilterOwnerInstance,j * 1.0 / dim_fft );
	pyres  = PyEval_CallObject(FilterFunct, pyargs);

	FILTER[2 * j] =   PyFloat_AsDouble(pyres)*    j * 2.0 / dim_fft/ dim_fft     ;
	FILTER[2 * j+1] = FILTER[2 * j];    
      } else {

	pyargs = Py_BuildValue("(O,i,i)",FilterOwnerInstance,j , dim_fft);
	pyres  = PyEval_CallObject(FilterFunct, pyargs);


	FILTER[2 * j] =   PyFloat_AsDouble(pyres)*     2.0/dim_fft     ;
	FILTER[2 * j+1] = FILTER[2 * j];    
      }

      
      Py_DECREF(pyargs);
      Py_DECREF(pyres );

    }
}


static char PyHST_SetFilterFunct_doc[]=""\
"/*  Set the Filter with a python Function."\
"    needs one argument : a function(x)"\
"    X goes from -0.5 to 0.5 "\
"    the function must return a float"\
" */"\
"";


static PyObject * PyHST_SetFilterFunct(PyObject *self_a, PyObject *args)
{
  PyHST *self;
  self = (PyHST *) self_a;
  if(!PyArg_ParseTuple(args,"OOi:PyHST_SetFilterFunct",&self->FilterOwnerInstance  , &self->FilterFunct, &self->FilterFunctHasRamp) ) {
    return NULL;
  }
  Py_XINCREF(self->FilterFunct);
  Py_XINCREF(self->FilterOwnerInstance);

  hst_set_filter(self->hst, self->FilterFunct?(HSTFilterFunction)PyHST_setFilter:NULL, self);

  Py_INCREF( Py_None);
  return Py_None;

}


static char PyHST_calcSlices_doc[]=""\
"/* calculate a number of slices given as argument"\
"  Takes two arguments : "\
"if(!PyArg_ParseTuple(args,\"ii:PyHST_calcSlices\", &num_slices ))"\
"    return NULL;"\
" */"\
"";

/* l ' ordinamento dimensioni per SINOGRAMS est num_slices self->num_projections  self->num_bins 

  sarebbe saggio verificare en passant che corrisponde alle dimensioni di self->SINOGRAMS
*/
static char PyHST_calcMedian_doc[]="";

#define swap_float( a,b) tmp_swap=(a); (a)=(b); (b)=tmp_swap;

 int compare_floats (const void * A,
                           const void * B)
{

  float *a, *b;
  a=(float *)A;
  b=(float *)B;


  if (*a > *b)
    return 1;
  else if (*a < *b)
    return -1;
  else
    return 0;
}

float getmedian(float *room_for_median, int dim2 ) {

  qsort ( room_for_median , dim2 , sizeof(float), compare_floats );


  return 0.5*( room_for_median[ dim2/2 ]
	       + room_for_median[ (int)(dim2/2.0-0.1 )] 
	       ) ;


}
  
static PyObject *
PyHST_calcMedian(PyObject *self_a, PyObject *args)
{
  int num_slices;
  float *SLICE;
  int first_slice;
  int last_slice;
  PyHST *self;
  PyObject *options , *arr_a ; 


  char *padding;
  char *axis_to_the_center;


  int offset, stride;
  float * data, * toorder;
  int iv,ix,i;
  int nd;
  int dims[2];
  int domean;
  double sum;

  pyhst_debug(" sono in median ");

  first_slice=0;
  options = 0 ; 
  padding = "E"; /* E as extrema, other value is 0 (zero) */
  axis_to_the_center = "N"; /* 'Y' or 'N' */

  self = (PyHST *) self_a;

  if(!PyArg_ParseTuple(args,"ii|O:PyHST_calcSlices", &num_slices, &domean, &options))
    return NULL;
  last_slice=first_slice+num_slices-1;

  SLICE=malloc( num_slices * self->setup.num_bins* sizeof(float) );

  if( ((PyArrayObject *) self->SINOGRAMS)->dimensions[1] !=  self->setup.num_projections) {
    pyhst_fail("  ((PyArrayObject *) self->SINOGRAMS)->dimensions[0] !=  self->num_projections   %d %d  ",  ((PyArrayObject *) self->SINOGRAMS)->dimensions[1] ,  self->setup.num_projections );
  }
  if( ((PyArrayObject *) self->SINOGRAMS)->dimensions[2] !=  self->setup.num_bins) {
    pyhst_fail("  ((PyArrayObject *) self->SINOGRAMS)->dimensions[2] !=  self->num_bins  %d %d  ",  ((PyArrayObject *) self->SINOGRAMS)->dimensions[2] ,  self->setup.num_bins );
  }

  /*   if( ((PyArrayObject *) self->SINOGRAMS)->dimensions[0] !=  num_slices ) { */
  /*     pyhst_debug("  ((PyArrayObject *) self->SINOGRAMS)->dimensions[0] !=  num_slices   %d %d  ",  ((PyArrayObject *) self->SINOGRAMS)->dimensions[0] , num_slices  ); */
  /*     exit(0); */
  /*   } */
  

  data = (float*) ((PyArrayObject *) self->SINOGRAMS)->data;
  toorder = malloc(  self->setup.num_projections*  sizeof(float));
  for(iv=0; iv<num_slices; iv++) {
    for(ix=0; ix<self->setup.num_bins; ix++) {
      
      offset = iv*self->setup.num_bins*self->setup.num_projections + ix;
      stride = self->setup.num_bins;
      for(i=0; i<self->setup.num_projections;i++) {
	toorder[i] = data[offset+i*stride];
      }
      if(domean==0){
	SLICE[iv*self->setup.num_bins  + ix] =getmedian(toorder, self->setup.num_projections);
      } else {
	sum=0;
	for(i=0; i<self->setup.num_projections;i++) {
	  sum += toorder[i];
	}
	SLICE[iv*self->setup.num_bins  + ix] = sum / self->setup.num_projections;
      }
    }
  }
 
  nd=2;
  dims[0] =num_slices ;
  dims[1] = self->setup.num_bins;
  arr_a = PyArray_FromDims(nd,dims,PyArray_FLOAT ) ;
  memcpy( ((PyArrayObject *) arr_a)->data, SLICE  ,   num_slices * self->setup.num_bins*sizeof(float) );

  free(toorder);
  free(SLICE);
  return PyArray_Return((PyArrayObject *) arr_a);

}

static PyObject *
PyHST_RECON(PyHST *self, int num_slices, PyObject *options, PyArrayObject *Oslice) {
  int err;
  
  PyObject * dic_key;
  PyObject * dic_value;

  char *padding;
  char *axis_to_the_center;

#ifdef PYHST_MEASURE_TIMINGS
  GTimer *timer;

  int i;
  HSTReconstructorConstContextPtr *reconstructors;
  PyObject *recon_dict;
  PyObject *timers;
  PyObject *slices;
#endif /* PYHST_MEASURE_TIMINGS */

  padding = "E"; /* E as extrema, other value is 0 (zero) */
  axis_to_the_center = "N"; /* 'Y' or 'N' */

  if(options!=0) {
    dic_key   =  PyString_FromString("padding");
    dic_value =  PyDict_GetItem(options, dic_key);
    if( dic_value  ) {
      padding   =  PyString_AsString(dic_value );  
    }
    Py_DECREF(dic_key);
    dic_key   =  PyString_FromString("axis_to_the_center");
    dic_value =  PyDict_GetItem(options, dic_key);
    if( dic_value  ) {
      axis_to_the_center   =  PyString_AsString(dic_value );  
    }
    Py_DECREF(dic_key);
  }

    // padding[0]=='E' pads with extrema, '0' ( zero ) with zeros , defaults to 'E'
  hst_set_padding(self->hst, *padding == '0');

    // /* 'Y' or 'N' */
  hst_set_axis_mode(self->hst, *axis_to_the_center == 'Y');

#ifdef PYHST_MEASURE_TIMINGS
  timer = g_timer_new();
#endif /* PYHST_MEASURE_TIMINGS */

  if ((Oslice)||(!self->ibuffer2)) {
     err = hst_reconstruct_slices(self->hst,
    	num_slices,
        Oslice?(float*)Oslice->data:NULL,
	(float*) ((PyArrayObject *) self->SINOGRAMS)->data
    );
  } else {
     hst_reconstruct_wait(self->hst);
     
     err = hst_reconstruct_start(self->hst, num_slices, (float*) ((PyArrayObject *) self->SINOGRAMS)->data);
     
     if (self->SINOGRAMS == self->SINO1) self->SINOGRAMS = self->SINO2;
     else self->SINOGRAMS = self->SINO1;
  }

  check_code(err, "hst_reconstruct_slices");

#ifdef PYHST_MEASURE_TIMINGS
  if (timer) {
    self->recon_timer += g_timer_elapsed(timer, NULL);
    g_timer_destroy(timer);
  }

  self->comp_timer = hst_get_recon_timer(self->hst);
  self->io_timer = hst_get_io_timer(self->hst);
  
  if (self->recon_info) {
     reconstructors = hst_get_configured_reconstructors(self->hst);
     for (i = 0; reconstructors[i]; i++) {
        recon_dict = PyList_GetItem(self->recon_info, i);
	if (recon_dict) {
	    timers = PyDict_GetItemString(recon_dict, "timers");
	    if (timers) {
		hst_reconstructor_get_timers(reconstructors[i], (double *) ((PyArrayObject*)timers)->data);
	    }
	    
	    slices  = PyInt_FromLong(reconstructors[i]->processed_slices);
	    if (slices) {
		PyDict_SetItemString(recon_dict, "slices", slices);
		Py_DECREF(slices);
	    }
	}
    }
  }
#endif /* PYHST_MEASURE_TIMINGS */

  Py_INCREF( Py_None);
  return Py_None;
}


static PyObject *
PyHST_calcSlices(PyObject *self_a, PyObject *args)
{
  PyHST *self;

  int num_slices;
  PyObject *options ; 

  self = (PyHST *) self_a;

  if(!PyArg_ParseTuple(args,"i|O:PyHST_calcSlices", &num_slices,  &options))
    return NULL;

  return PyHST_RECON(self, num_slices, options, NULL);
}



static char PyHST_startPreprocessing_doc[]="";
static PyObject *
PyHST_startPreprocessing(PyObject *self_a, PyObject *args) {
  PyHST *self;

  int num_slices;

  self = (PyHST *) self_a;

  if(!PyArg_ParseTuple(args,"i:PyHST_calcSlices", &num_slices))
    return NULL;

  hst_preprocess_start(self->hst, num_slices, (float*) ((PyArrayObject *) self->SINOGRAMS)->data);
  
  Py_INCREF( Py_None);
  return Py_None;
}

static char PyHST_waitPreprocessing_doc[]="";
static PyObject *
PyHST_waitPreprocessing(PyObject *self_a, PyObject *args) {
  int i;
  PyHST *self;

  self = (PyHST *) self_a;

  hst_preprocess_wait(self->hst);
  
#ifdef HW_USE_PARALLEL_IO
  for (i = 0; i < self->setup.num_projections; i++) {
    Py_DECREF(self->Projections[i]);
  }
#endif /* HW_USE_PARALLEL_IO */
  
  Py_INCREF( Py_None);
  return Py_None;
}


static char PyHST_transposeSlices_doc[]=""\
"/* transposes a number of slices given as argument"\
" */"\
"";
static PyObject *
PyHST_transposeSlices(PyObject *self_a, PyObject *args)
{
  PyHST *self;

  int err;
  int proj;
  PyObject *item; 

  self = (PyHST *) self_a;

  if(!PyArg_ParseTuple(args,"Oi:PyHST_calcSlices", &item, &proj))
    return NULL;

  Py_INCREF(item);
  
  float *in = (float*) ((PyArrayObject *) item)->data;
  err =  hst_preprocess_projection(self->hst, proj, in);

#ifdef HW_USE_PARALLEL_IO
  self->Projections[proj] = item;
#else /* HW_USE_PARALLEL_IO */
  Py_DECREF(item);
#endif /* HW_USE_PARALLEL_IO */

  Py_INCREF( Py_None);
  return Py_None;
}



static char PyHST_calcSlicesMemory_doc[]=""\
"/* calculate a number of slices given as argument"\
"  Takes three arguments : "\
"if(!PyArg_ParseTuple(args,\"iiO:PyHST_calcSlices\", &num_slices,  &Oslice))"\
"    return NULL;"\
" */"\
"";
static PyObject *
PyHST_calcSlicesMemory(PyObject *self_a, PyObject *args)
{
  PyHST *self;

  int num_slices;
  PyObject *options ; 
  PyArrayObject * Oslice ; 

  self = (PyHST *) self_a;

  if(!PyArg_ParseTuple(args,"iO|O:PyHST_calcSlicesMemory", &num_slices,  &Oslice, &options )  )  
    return NULL;

  return PyHST_RECON(self, num_slices, options, Oslice);
}



static PyObject *
PyHST_waitCompletion(PyObject *self_a, PyObject *args)
{
  PyHST *self;
  self = (PyHST *) self_a;

#ifdef PYHST_MEASURE_TIMINGS
  int i;
  HSTReconstructorConstContextPtr *reconstructors;
  PyObject *recon_dict;
  PyObject *timers;
  PyObject *slices;

  GTimer *timer;

  timer = g_timer_new();
#endif /* PYHST_MEASURE_TIMINGS */

  hst_reconstruct_wait(self->hst);

#ifdef PYHST_MEASURE_TIMINGS
  if (timer) {
    self->recon_timer += g_timer_elapsed(timer, NULL);
    g_timer_destroy(timer);
  }

  self->comp_timer = hst_get_recon_timer(self->hst);
  self->io_timer = hst_get_io_timer(self->hst);
  
  if (self->recon_info) {
     reconstructors = hst_get_configured_reconstructors(self->hst);
     for (i = 0; reconstructors[i]; i++) {
        recon_dict = PyList_GetItem(self->recon_info, i);
	if (recon_dict) {
	    timers = PyDict_GetItemString(recon_dict, "timers");
	    if (timers) {
		hst_reconstructor_get_timers(reconstructors[i], (double *) ((PyArrayObject*)timers)->data);
	    }
	    
	    slices  = PyInt_FromLong(reconstructors[i]->processed_slices);
	    if (slices) {
		PyDict_SetItemString(recon_dict, "slices", slices);
		Py_DECREF(slices);
	    }
	}
    }
  }

#endif /* PYHST_MEASURE_TIMINGS */

  Py_INCREF( Py_None);
  return Py_None;
}


static PyMethodDef PyHST_methods[]={
  {"startPreprocessing",(PyCFunction) PyHST_startPreprocessing, METH_VARARGS,PyHST_startPreprocessing_doc},
  {"waitPreprocessing",(PyCFunction) PyHST_waitPreprocessing, METH_VARARGS,PyHST_waitPreprocessing_doc},
  {"transposeSlices",(PyCFunction) PyHST_transposeSlices, METH_VARARGS,PyHST_transposeSlices_doc},
  {"calcSlices",(PyCFunction) PyHST_calcSlices, METH_VARARGS,PyHST_calcSlices_doc},
  {"calcMedian",(PyCFunction) PyHST_calcMedian, METH_VARARGS,PyHST_calcMedian_doc},
  {"calcSlicesMemory",(PyCFunction) PyHST_calcSlicesMemory, METH_VARARGS,PyHST_calcSlicesMemory_doc},
  {"waitCompletion",(PyCFunction) PyHST_waitCompletion, METH_VARARGS,NULL},
  {"setFilterFunct",(PyCFunction) PyHST_SetFilterFunct, METH_VARARGS,PyHST_SetFilterFunct_doc},
  {"close", (PyCFunction) PyHST_close,  METH_VARARGS, NULL    },
  { NULL, NULL}
};


static struct memberlist PyHST_memberlist[]={
    {"SINOGRAMS", T_OBJECT, offsetof(PyHST, SINOGRAMS)},
    {"astra_scaling", T_INT, offsetof(PyHST, astra_scaling)},
#ifdef PYHST_MEASURE_TIMINGS
    {"recon_timer", T_DOUBLE, offsetof(PyHST, recon_timer)},
    {"comp_timer", T_DOUBLE, offsetof(PyHST, comp_timer)},
    {"io_timer", T_DOUBLE, offsetof(PyHST, io_timer)},
    {"recon", T_OBJECT, offsetof(PyHST, recon_info)},
#endif /* PYHST_MEASURE_TIMINGS */
    { NULL }
};

static PyMethodDef PyHST_functions[] = {
  {"PyHST", PyHST_new, METH_VARARGS, PyHST_new_doc },
  { NULL, NULL}
};



static PyObject *
PyHST_getattr(PyHST *self, char *attr) 
{
  PyObject *res;

  res= Py_FindMethod(PyHST_methods, (PyObject*) self, attr);
  if(NULL !=res)
    return res;
  else {
    PyErr_Clear();
    return PyMember_Get((char*) self,PyHST_memberlist, attr);
  }
}



static PyObject *logger = NULL;

void finiPyHST_c(void)
{
  hst_free();
  if (logger) {
     pyhst_configure_logger(0, (GLogFunc)NULL, NULL);
     Py_DECREF(logger);
  }
}

void initPyHST_c(void)
{
  PyObject *logger_mod;
  PyObject *m, *d;
  size_t n_rec;
  
  m = Py_InitModule("PyHST_c", PyHST_functions);
  d = PyModule_GetDict(m);
  ErrorObject = Py_BuildValue("s","PyHST_c.error");
  PyDict_SetItemString(d,"error", ErrorObject);
  if(PyErr_Occurred())
    Py_FatalError("can't initialize module PyHST_c");
  else {
    logger_mod = PyImport_ImportModule("logger");
    if (logger_mod) {
	logger = PyObject_GetAttrString(logger_mod, "logger");
	Py_DECREF(logger_mod);
	if (logger) {
	    pyhst_configure_logger(0, (GLogFunc)PyHST_log_handler, logger);
	}
    }

    hst_init();
    n_rec = hst_get_balanced_number_of_reconstructors();
    PyModule_AddIntConstant(m, "RECONSTRUCTORS", n_rec);
    
    Py_AtExit(finiPyHST_c);
  }
#ifdef import_array
  import_array();
#endif
}



static PyTypeObject PyHSTtype = {
  PyObject_HEAD_INIT(&PyType_Type)
  0,
  "PyHST",
  sizeof(PyHST),
  0,
  (destructor) PyHST_dealloc,
  0,
  (getattrfunc) PyHST_getattr,
};